From 9cb0384118cb8606c7299b2342fd428b4d9f3e6a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tobias=20H=C3=83=C2=BCfner?= Date: Wed, 6 Nov 2024 11:02:00 +0100 Subject: [PATCH 1/5] Added functionality to convert simulation box vectors to restricted triclinic box format in lammps --- .../interop/lammps/export/export.py | 44 +++++++++++++++---- 1 file changed, 35 insertions(+), 9 deletions(-) diff --git a/openff/interchange/interop/lammps/export/export.py b/openff/interchange/interop/lammps/export/export.py index 9f25a50ac..f6814082a 100644 --- a/openff/interchange/interop/lammps/export/export.py +++ b/openff/interchange/interop/lammps/export/export.py @@ -68,21 +68,47 @@ 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 = [ax,0,0] + b_new = [bx,by,0] + c_new = [cx,cy,cz] + + 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") From 0fbde383754fc35f0f90b43a1f309e78e8978ab6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tobias=20H=C3=83=C2=BCfner?= Date: Wed, 6 Nov 2024 11:52:14 +0100 Subject: [PATCH 2/5] Added lammps triclinic box test --- .../interop/lammps/export/test_lammps.py | 63 +++++++++++++++++++ 1 file changed, 63 insertions(+) diff --git a/openff/interchange/_tests/unit_tests/interop/lammps/export/test_lammps.py b/openff/interchange/_tests/unit_tests/interop/lammps/export/test_lammps.py index 06eacf209..77a084628 100644 --- a/openff/interchange/_tests/unit_tests/interop/lammps/export/test_lammps.py +++ b/openff/interchange/_tests/unit_tests/interop/lammps/export/test_lammps.py @@ -15,6 +15,69 @@ @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. + """ + mol = MoleculeWithConformer.from_smiles(mol) + mol.conformers[0] -= numpy.min(mol.conformers[0], axis=0) + top = Topology.from_molecules(n_mols * [mol]) + + box = numpy.zeros((3,3), dtype=float) * unit.angstrom + + box[0] = [51.34903463831951, 0, 0] * unit.angstrom + box[1] = [-0.03849979989403723, 50.9134404144338, 0] * unit.angstrom + box[2] = [-2.5907782992729538, 0.3720740833800747, 49.80705567557188] * unit.angstrom + + top.box_vectors = box + + 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 = top.box_vectors + + reference = get_openmm_energies( + interchange=interchange, + round_positions=3, + ) + + lmp_energies = get_lammps_energies( + interchange=interchange, + round_positions=3, + ) + + lmp_energies.compare( + reference, + tolerances={ + "Nonbonded": 1 * unit.kilojoule_per_mole, + "Torsion": 0.02 * unit.kilojoule_per_mole, + }, + ) + @pytest.mark.parametrize("n_mols", [1, 2]) @pytest.mark.parametrize( "mol", From 8aa89210bb8f8767369beb9621904d8cbc1c954f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tobias=20H=C3=83=C2=BCfner?= Date: Wed, 6 Nov 2024 12:08:11 +0100 Subject: [PATCH 3/5] Added test to compare box vectors before and after the transformation. --- .../interchange/interop/lammps/export/export.py | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/openff/interchange/interop/lammps/export/export.py b/openff/interchange/interop/lammps/export/export.py index f6814082a..3e86dfb3c 100644 --- a/openff/interchange/interop/lammps/export/export.py +++ b/openff/interchange/interop/lammps/export/export.py @@ -94,12 +94,17 @@ def to_lammps(interchange: Interchange, file_path: Path | str): xy, xz, yz = bx, cx, cy - a_new = [ax,0,0] - b_new = [bx,by,0] - c_new = [cx,cy,cz] - - import warnings - warnings.warn( + 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( From 93384ed44aad23798252ae06fc348d2d75b4857b Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 6 Nov 2024 11:10:43 +0000 Subject: [PATCH 4/5] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../interop/lammps/export/test_lammps.py | 2 +- .../interop/lammps/export/export.py | 32 ++++++++++--------- 2 files changed, 18 insertions(+), 16 deletions(-) diff --git a/openff/interchange/_tests/unit_tests/interop/lammps/export/test_lammps.py b/openff/interchange/_tests/unit_tests/interop/lammps/export/test_lammps.py index 77a084628..da533c8ba 100644 --- a/openff/interchange/_tests/unit_tests/interop/lammps/export/test_lammps.py +++ b/openff/interchange/_tests/unit_tests/interop/lammps/export/test_lammps.py @@ -41,7 +41,7 @@ def test_to_lammps_single_mols_triclinic( mol.conformers[0] -= numpy.min(mol.conformers[0], axis=0) top = Topology.from_molecules(n_mols * [mol]) - box = numpy.zeros((3,3), dtype=float) * unit.angstrom + box = numpy.zeros((3, 3), dtype=float) * unit.angstrom box[0] = [51.34903463831951, 0, 0] * unit.angstrom box[1] = [-0.03849979989403723, 50.9134404144338, 0] * unit.angstrom diff --git a/openff/interchange/interop/lammps/export/export.py b/openff/interchange/interop/lammps/export/export.py index 3e86dfb3c..75ffdfe7d 100644 --- a/openff/interchange/interop/lammps/export/export.py +++ b/openff/interchange/interop/lammps/export/export.py @@ -81,31 +81,33 @@ def to_lammps(interchange: Interchange, file_path: Path | str): 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)) + 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) + 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]) + 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(a, a_new) check *= numpy.allclose(b, b_new) check *= numpy.allclose(c, c_new) - - if not check: + + 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}") + 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 + ax:.10g} xlo xhi\n" From d23f95344af75f592e0d107a123205e90ac94877 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tobias=20H=C3=83=C2=BCfner?= Date: Wed, 6 Nov 2024 14:22:58 +0100 Subject: [PATCH 5/5] Modified lammps triclinic test. --- .../interop/lammps/export/test_lammps.py | 46 ++++++++++++------- 1 file changed, 30 insertions(+), 16 deletions(-) diff --git a/openff/interchange/_tests/unit_tests/interop/lammps/export/test_lammps.py b/openff/interchange/_tests/unit_tests/interop/lammps/export/test_lammps.py index 77a084628..07eab8740 100644 --- a/openff/interchange/_tests/unit_tests/interop/lammps/export/test_lammps.py +++ b/openff/interchange/_tests/unit_tests/interop/lammps/export/test_lammps.py @@ -12,7 +12,6 @@ rng = numpy.random.default_rng(821) - @needs_lmp class TestLammps: @pytest.mark.parametrize("n_mols", [1, 2]) @@ -36,18 +35,24 @@ def test_to_lammps_single_mols_triclinic( """ 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]) - box = numpy.zeros((3,3), dtype=float) * unit.angstrom + box_t = numpy.zeros((3,3), dtype=float) * unit.angstrom + + 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[0] = [51.34903463831951, 0, 0] * unit.angstrom - box[1] = [-0.03849979989403723, 50.9134404144338, 0] * unit.angstrom - box[2] = [-2.5907782992729538, 0.3720740833800747, 49.80705567557188] * unit.angstrom + box_o = 50. * numpy.eye(3) * unit.angstrom - top.box_vectors = box + top.box_vectors = box_o if n_mols == 1: positions = mol.conformers[0] @@ -58,26 +63,35 @@ def test_to_lammps_single_mols_triclinic( interchange = Interchange.from_smirnoff(sage_unconstrained, top) interchange.positions = positions - interchange.box = top.box_vectors + interchange.box = box_o - reference = get_openmm_energies( + reference_o = get_openmm_energies( interchange=interchange, round_positions=3, ) - lmp_energies = get_lammps_energies( + lmp_energies_o = get_lammps_energies( interchange=interchange, round_positions=3, ) - lmp_energies.compare( - reference, - tolerances={ - "Nonbonded": 1 * unit.kilojoule_per_mole, - "Torsion": 0.02 * unit.kilojoule_per_mole, - }, + 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", @@ -132,7 +146,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, }, )