Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fixes in LAMMPS export #238

Merged
merged 3 commits into from
Jun 30, 2021
Merged

Fixes in LAMMPS export #238

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
c7a564e
BUG: Possibly fix improper torsions in LAMMPS export
mattwthompson Jun 23, 2021
cf25601
LINT: Remove old checks for cvff compatability
mattwthompson Jun 23, 2021
f240eeb
Merge remote-tracking branch 'upstream/master' into lammps-impropers
mattwthompson Jun 30, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion openff/interchange/drivers/lammps.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"]
Expand Down
47 changes: 38 additions & 9 deletions openff/interchange/interop/internal/lammps.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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"]
Expand All @@ -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"]
Expand All @@ -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"]
Expand All @@ -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"]
Expand All @@ -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"]
Expand All @@ -215,23 +221,42 @@ def _write_improper_coeffs(lmp_file: IO, openff_sys: Interchange):
idivf = int(params["idivf"])
k = k / idivf

if (phase != 180) or (n != 2):
# 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))

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 "
"at the moment, see PR #126"
"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"
)

# See https://lammps.sandia.gov/doc/improper_fourier.html
# cos(n * x - pi) == - cos(n * x)
# 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")
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()})
Expand Down Expand Up @@ -265,6 +290,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"]
Expand Down Expand Up @@ -298,6 +324,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,
Expand Down Expand Up @@ -331,6 +358,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,
Expand Down Expand Up @@ -366,6 +394,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,
Expand Down
15 changes: 4 additions & 11 deletions openff/interchange/tests/interop/internal/test_lammps.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,21 +21,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
Expand Down Expand Up @@ -84,6 +77,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,
},
)