From c7a564e67938c112d0deb7666995b9c42816cc14 Mon Sep 17 00:00:00 2001 From: "Matthew W. Thompson" Date: Wed, 23 Jun 2021 17:00:27 -0500 Subject: [PATCH 1/2] BUG: Possibly fix improper torsions in LAMMPS export --- .pre-commit-config.yaml | 2 +- openff/interchange/drivers/lammps.py | 2 +- openff/interchange/interop/internal/lammps.py | 43 +++++++++++++++++-- .../tests/interop/internal/test_lammps.py | 15 ++----- 4 files changed, 45 insertions(+), 17 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index ee1129da6..04556a6ba 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -39,7 +39,7 @@ repos: rev: 1.4.0 hooks: - id: interrogate - args: [--quiet, --fail-under=45, openff/interchange] + args: [--fail-under=45, openff/interchange/] - repo: https://github.com/nbQA-dev/nbQA rev: 0.13.0 hooks: diff --git a/openff/interchange/drivers/lammps.py b/openff/interchange/drivers/lammps.py index be98de66e..cf6e97e57 100644 --- a/openff/interchange/drivers/lammps.py +++ b/openff/interchange/drivers/lammps.py @@ -115,7 +115,7 @@ def _write_lammps_input( fo.write("dihedral_style hybrid fourier\n") if "ImproperTorsions" in off_sys.handlers: if len(off_sys["ImproperTorsions"].potentials) > 0: - fo.write("improper_style hybrid cvff \n") + fo.write("improper_style cvff\n") vdw_hander = off_sys.handlers["vdW"] electrostatics_handler = off_sys.handlers["Electrostatics"] diff --git a/openff/interchange/interop/internal/lammps.py b/openff/interchange/interop/internal/lammps.py index 2f1afbda8..291a5ad7f 100644 --- a/openff/interchange/interop/internal/lammps.py +++ b/openff/interchange/interop/internal/lammps.py @@ -10,6 +10,7 @@ def to_lammps(openff_sys: Interchange, file_path: Union[Path, str]): + """Write an Interchange object to a LAMMPS data file""" if isinstance(file_path, str): path = Path(file_path) @@ -127,6 +128,7 @@ def to_lammps(openff_sys: Interchange, file_path: Union[Path, str]): def _write_pair_coeffs(lmp_file: IO, openff_sys: Interchange, atom_type_map: Dict): + """Write the Pair Coeffs section of a LAMMPS data file""" lmp_file.write("Pair Coeffs\n\n") vdw_handler = openff_sys["vdW"] @@ -143,6 +145,7 @@ def _write_pair_coeffs(lmp_file: IO, openff_sys: Interchange, atom_type_map: Dic def _write_bond_coeffs(lmp_file: IO, openff_sys: Interchange): + """Write the Bond Coeffs section of a LAMMPS data file""" lmp_file.write("Bond Coeffs\n\n") bond_handler = openff_sys.handlers["Bonds"] @@ -161,6 +164,7 @@ def _write_bond_coeffs(lmp_file: IO, openff_sys: Interchange): def _write_angle_coeffs(lmp_file: IO, openff_sys: Interchange): + """Write the Angle Coeffs section of a LAMMPS data file""" lmp_file.write("\nAngle Coeffs\n\n") angle_handler = openff_sys.handlers["Angles"] @@ -179,6 +183,7 @@ def _write_angle_coeffs(lmp_file: IO, openff_sys: Interchange): def _write_proper_coeffs(lmp_file: IO, openff_sys: Interchange): + """Write the Dihedral Coeffs section of a LAMMPS data file""" lmp_file.write("\nDihedral Coeffs\n\n") proper_handler = openff_sys.handlers["ProperTorsions"] @@ -201,6 +206,7 @@ def _write_proper_coeffs(lmp_file: IO, openff_sys: Interchange): def _write_improper_coeffs(lmp_file: IO, openff_sys: Interchange): + """Write the Improper Coeffs section of a LAMMPS data file""" lmp_file.write("\nImproper Coeffs\n\n") improper_handler = openff_sys.handlers["ImproperTorsions"] @@ -221,17 +227,42 @@ def _write_improper_coeffs(lmp_file: IO, openff_sys: Interchange): "at the moment, see PR #126" ) - # See https://lammps.sandia.gov/doc/improper_fourier.html - # cos(n * x - pi) == - cos(n * x) + # See https://lammps.sandia.gov/doc/improper_cvff.html + # E_periodic = k * (1 + cos(n * theta - phase)) + # E_cvff = k * (1 + d * cos(n * theta)) + # k_periodic = k_cvff + # if phase = 0, + # d_cvff = 1 + # n_periodic = n_cvff + # if phase = 180, + # cos(n * x - pi) == - cos(n * x) + # d_cvff = -1 + # n_periodic = n_cvff # k * (1 + cos(n * phi - pi / 2)) == k * (1 - cos(n * phi)) - d = -1 - lmp_file.write(f"{improper_type_idx+1:d} cvff {k:.16g}\t{d:d}\t{n:.16g}\n") + if phase == 0: + k_cvff = k + d_cvff = 1 + n_cvff = n + elif phase == 180: + k_cvff = k + d_cvff = -1 + n_cvff = n + else: + raise UnsupportedExportError( + "Improper exports to LAMMPS are funky and not well-supported, the only compatibility" + "found between periodidic impropers is with improper_style cvff when phase = 0 or 180 degrees" + ) + + lmp_file.write( + f"{improper_type_idx+1:d} {k_cvff:.16g}\t{d_cvff:d}\t{n_cvff:.16g}\n" + ) lmp_file.write("\n") def _write_atoms(lmp_file: IO, openff_sys: Interchange, atom_type_map: Dict): + """Write the Atoms section of a LAMMPS data file""" lmp_file.write("\nAtoms\n\n") atom_type_map_inv = dict({v: k for k, v in atom_type_map.items()}) @@ -265,6 +296,7 @@ def _write_atoms(lmp_file: IO, openff_sys: Interchange, atom_type_map: Dict): def _write_bonds(lmp_file: IO, openff_sys: Interchange): + """Write the Bonds section of a LAMMPS data file""" lmp_file.write("\nBonds\n\n") bond_handler = openff_sys["Bonds"] @@ -298,6 +330,7 @@ def _write_bonds(lmp_file: IO, openff_sys: Interchange): def _write_angles(lmp_file: IO, openff_sys: Interchange): + """Write the Angles section of a LAMMPS data file""" from openff.interchange.components.mdtraj import ( _iterate_angles, _store_bond_partners, @@ -331,6 +364,7 @@ def _write_angles(lmp_file: IO, openff_sys: Interchange): def _write_propers(lmp_file: IO, openff_sys: Interchange): + """Write the Dihedrals section of a LAMMPS data file""" from openff.interchange.components.mdtraj import ( _iterate_propers, _store_bond_partners, @@ -366,6 +400,7 @@ def _write_propers(lmp_file: IO, openff_sys: Interchange): def _write_impropers(lmp_file: IO, openff_sys: Interchange): + """Write the Impropers section of a LAMMPS data file""" from openff.interchange.components.mdtraj import ( _iterate_impropers, _store_bond_partners, diff --git a/openff/interchange/tests/interop/internal/test_lammps.py b/openff/interchange/tests/interop/internal/test_lammps.py index 5582263fd..097d8588b 100644 --- a/openff/interchange/tests/interop/internal/test_lammps.py +++ b/openff/interchange/tests/interop/internal/test_lammps.py @@ -20,21 +20,14 @@ "C", "CC", # Adds a proper torsion term(s) "C=O", # Simplest molecule with any improper torsion - pytest.param( - "OC=O", - marks=pytest.mark.xfail(reason="degenerate impropers"), - ), # Simplest molecule with a multi-term torsion + "OC=O", # Simplest molecule with a multi-term torsion "CCOC", # This hits t86, which has a non-1.0 idivf - pytest.param( - "C1COC(=O)O1", - marks=pytest.mark.xfail(reason="degenerate impropers"), - ), # This adds an improper, i2 + "C1COC(=O)O1", # This adds an improper, i2 ], ) def test_to_lammps_single_mols(mol, n_mols): """ - Test that ForceField.create_openmm_system and Interchange.to_openmm produce - objects with similar energies + Test that Interchange.to_openmm Interchange.to_lammps report sufficiently similar energies. TODO: Tighten tolerances TODO: Test periodic and non-periodic @@ -83,6 +76,6 @@ def test_to_lammps_single_mols(mol, n_mols): "Nonbonded": 100 * omm_unit.kilojoule_per_mole, "Electrostatics": 100 * omm_unit.kilojoule_per_mole, "vdW": 100 * omm_unit.kilojoule_per_mole, - "Torsion": 0.005 * omm_unit.kilojoule_per_mole, + "Torsion": 3e-5 * omm_unit.kilojoule_per_mole, }, ) From cf25601feed7fa09c81418c42ee5d1ad865b4e59 Mon Sep 17 00:00:00 2001 From: "Matthew W. Thompson" Date: Wed, 23 Jun 2021 17:34:56 -0500 Subject: [PATCH 2/2] LINT: Remove old checks for cvff compatability --- openff/interchange/interop/internal/lammps.py | 6 ------ 1 file changed, 6 deletions(-) diff --git a/openff/interchange/interop/internal/lammps.py b/openff/interchange/interop/internal/lammps.py index 291a5ad7f..0f16dcb17 100644 --- a/openff/interchange/interop/internal/lammps.py +++ b/openff/interchange/interop/internal/lammps.py @@ -221,12 +221,6 @@ def _write_improper_coeffs(lmp_file: IO, openff_sys: Interchange): idivf = int(params["idivf"]) k = k / idivf - if (phase != 180) or (n != 2): - raise UnsupportedExportError( - "Improper exports to LAMMPS are funky and not well-supported " - "at the moment, see PR #126" - ) - # See https://lammps.sandia.gov/doc/improper_cvff.html # E_periodic = k * (1 + cos(n * theta - phase)) # E_cvff = k * (1 + d * cos(n * theta))