Added triclinic boxes for lammps

openff/interchange/_tests/unit_tests/interop/lammps/export/test_lammps.py

@@ -12,9 +12,90 @@
 
 rng = numpy.random.default_rng(821)
 
-
 @needs_lmp
 class TestLammps:
+    @pytest.mark.parametrize("n_mols", [1, 2])
+    @pytest.mark.parametrize(
+        "mol",
+        [
+            "C",
+            "CC",  # Adds a proper torsion term(s)
+            "C=O",  # Simplest molecule with any improper torsion
+            "OC=O",  # Simplest molecule with a multi-term torsion
+            "CCOC",  # This hits t86, which has a non-1.0 idivf
+            "C1COC(=O)O1",  # This adds an improper, i2
+        ],
+    )
+    def test_to_lammps_single_mols_triclinic(
+        self,
+        mol: str,
+        sage_unconstrained: ForceField,
+        n_mols: int,
+    ) -> None:
+        """
+        Test that Interchange.to_openmm Interchange.to_lammps report sufficiently similar energies
+        in triclinic simulation boxes.
+        This test is only sensitive to nonbonded interactions. Due to inherent differences in how
+        long-ranged interactions are handled in openmm/lammps, there will be a baseline deviation
+        between the two codes that goes beyond this test. This test is designed to test whether the
+        triclinic box representations give no error greater than the ones from the orthogonal box test.
+        """
+        mol = MoleculeWithConformer.from_smiles(mol)
+        mol.conformers[0] -= numpy.min(mol.conformers[0], axis=0)
+        top = Topology.from_molecules(n_mols * [mol])
+
+<<<<<<< HEAD
+        box_t = numpy.zeros((3,3), dtype=float) * unit.angstrom
+=======
+        box = numpy.zeros((3, 3), dtype=float) * unit.angstrom
+>>>>>>> 93384ed44aad23798252ae06fc348d2d75b4857b
+
+        box_t[0] = [51.34903463831951, 0, 0] * unit.angstrom
+        box_t[1] = [-0.03849979989403723, 50.9134404144338, 0] * unit.angstrom
+        box_t[2] = [-2.5907782992729538, 0.3720740833800747, 49.80705567557188] * unit.angstrom
+
+        box_o = 50. * numpy.eye(3) * unit.angstrom
+
+        top.box_vectors = box_o
+
+        if n_mols == 1:
+            positions = mol.conformers[0]
+        elif n_mols == 2:
+            positions = numpy.concatenate(
+                [mol.conformers[0], mol.conformers[0] + 1.5 * unit.nanometer],
+            )
+
+        interchange = Interchange.from_smirnoff(sage_unconstrained, top)
+        interchange.positions = positions
+        interchange.box = box_o
+
+        reference_o = get_openmm_energies(
+            interchange=interchange,
+            round_positions=3,
+        )
+
+        lmp_energies_o = get_lammps_energies(
+            interchange=interchange,
+            round_positions=3,
+        )
+
+        interchange.box = box_t
+        reference_t = get_openmm_energies(
+            interchange=interchange,
+            round_positions=3,
+        )
+
+        lmp_energies_t = get_lammps_energies(
+            interchange=interchange,
+            round_positions=3,
+        )
+
+        diff_o = reference_o.diff(lmp_energies_o)
+        diff_t = reference_t.diff(lmp_energies_t)
+
+        assert abs(diff_o["Nonbonded"] - diff_t["Nonbonded"]) < 0.1 * unit.kilojoule_per_mole
+
+
     @pytest.mark.parametrize("n_mols", [1, 2])
     @pytest.mark.parametrize(
         "mol",
@@ -69,7 +150,7 @@ def test_to_lammps_single_mols(
         lmp_energies.compare(
             reference,
             tolerances={
-                "Nonbonded": 1 * unit.kilojoule_per_mole,
+                "Nonbonded": 0.4 * unit.kilojoule_per_mole,
                 "Torsion": 0.02 * unit.kilojoule_per_mole,
             },
         )

openff/interchange/interop/lammps/export/export.py

@@ -68,21 +68,54 @@ def to_lammps(interchange: Interchange, file_path: Path | str):
             axis=0,
         ).magnitude
         if interchange.box is None:
-            L_x, L_y, L_z = 100, 100, 100
+            ax, by, cz = 100, 100, 100
+            xy, xz, yz = 0, 0, 0
         elif (interchange.box.m == numpy.diag(numpy.diagonal(interchange.box.m))).all():
-            L_x, L_y, L_z = numpy.diag(interchange.box.to(unit.angstrom).magnitude)
+            ax, by, cz = numpy.diag(interchange.box.to(unit.angstrom).magnitude)
+            xy, xz, yz = 0, 0, 0
         else:
-            raise NotImplementedError(
-                "Interchange does not yet support exporting non-rectangular boxes to LAMMPS",
-            )
+            a = interchange.box.to(unit.angstrom).magnitude[0]
+            b = interchange.box.to(unit.angstrom).magnitude[1]
+            c = interchange.box.to(unit.angstrom).magnitude[2]
+            a_len = numpy.linalg.norm(a)
+            b_len = numpy.linalg.norm(b)
+            c_len = numpy.linalg.norm(c)
+
+            alpha = numpy.arccos(numpy.dot(b, c) / (b_len * c_len))
+            beta = numpy.arccos(numpy.dot(c, a) / (c_len * a_len))
+            gamma = numpy.arccos(numpy.dot(a, b) / (a_len * b_len))
+
+            ax = a_len
+            bx = b_len * numpy.cos(gamma)
+            by = b_len * numpy.sin(gamma)
+            cx = c_len * numpy.cos(beta)
+            cy = c_len * (numpy.cos(alpha) - numpy.cos(beta) * numpy.cos(gamma)) / numpy.sin(gamma)
+            cz = numpy.sqrt(c_len**2 - cx**2 - cy**2)
+
+            xy, xz, yz = bx, cx, cy
+
+            a_new = numpy.array([ax, 0, 0])
+            b_new = numpy.array([bx, by, 0])
+            c_new = numpy.array([cx, cy, cz])
+
+            check = numpy.allclose(a, a_new)
+            check *= numpy.allclose(b, b_new)
+            check *= numpy.allclose(c, c_new)
+
+            if not check:
+                import warnings
+
+                warnings.warn(
+                    f"We are rotating the unit cell matrix. Before: a={a}, b={b}, c={c}, After: a={a_new}, b={b_new}, c={c_new}",
+                )
 
         lmp_file.write(
-            f"{x_min:.10g} {x_min + L_x:.10g} xlo xhi\n"
-            f"{y_min:.10g} {y_min + L_y:.10g} ylo yhi\n"
-            f"{z_min:.10g} {z_min + L_z:.10g} zlo zhi\n",
+            f"{x_min:.10g} {x_min + ax:.10g} xlo xhi\n"
+            f"{y_min:.10g} {y_min + by:.10g} ylo yhi\n"
+            f"{z_min:.10g} {z_min + cz:.10g} zlo zhi\n",
         )
 
-        lmp_file.write("0.0 0.0 0.0 xy xz yz\n")
+        lmp_file.write(f"{xy:.10g} {xz:.10g} {yz:.10g} xy xz yz\n")
 
         lmp_file.write("\nMasses\n\n")