diff --git a/mbuild/formats/lammpsdata.py b/mbuild/formats/lammpsdata.py
index 47e3c174d..382933dbc 100644
--- a/mbuild/formats/lammpsdata.py
+++ b/mbuild/formats/lammpsdata.py
@@ -41,6 +41,7 @@ def write_lammpsdata(
use_urey_bradleys=False,
use_rb_torsions=True,
use_dihedrals=False,
+ zero_dihedral_weighting_factor=False,
):
"""Output a LAMMPS data file.
@@ -81,6 +82,9 @@ def write_lammpsdata(
use_dihedrals:
If True, will treat dihedrals as CHARMM-style dihedrals while looking
for `structure.dihedrals`
+ zero_dihedral_weighting_factor:
+ If True, will set weighting parameter to zero in CHARMM-style dihedrals.
+ This should be True if the CHARMM dihedral style is used in non-CHARMM forcefields.
Notes
-----
@@ -240,13 +244,6 @@ def write_lammpsdata(
"Forcefield XML and structure"
)
- # Check impropers
- for dihedral in structure.dihedrals:
- if dihedral.improper:
- raise ValueError(
- "Amber-style impropers are currently not supported"
- )
-
bonds = [[b.atom1.idx + 1, b.atom2.idx + 1] for b in structure.bonds]
angles = [
[angle.atom1.idx + 1, angle.atom2.idx + 1, angle.atom3.idx + 1]
@@ -261,6 +258,7 @@ def write_lammpsdata(
dihedrals = [
[d.atom1.idx + 1, d.atom2.idx + 1, d.atom3.idx + 1, d.atom4.idx + 1]
for d in structure.dihedrals
+ if not d.improper
]
else:
dihedrals = []
@@ -268,6 +266,14 @@ def write_lammpsdata(
[i.atom1.idx + 1, i.atom2.idx + 1, i.atom3.idx + 1, i.atom4.idx + 1]
for i in structure.impropers
]
+ imp_dihedrals = [
+ [d.atom1.idx + 1, d.atom2.idx + 1, d.atom3.idx + 1, d.atom4.idx + 1]
+ for d in structure.dihedrals
+ if d.improper
+ ]
+
+ if impropers and imp_dihedrals:
+ raise ValueError("Use of multiple improper styles is not supported")
if bonds:
if len(structure.bond_types) == 0:
@@ -288,9 +294,18 @@ def write_lammpsdata(
epsilon_conversion_factor,
)
+ if imp_dihedrals:
+ (
+ imp_dihedral_types,
+ unique_imp_dihedral_types,
+ ) = _get_improper_dihedral_types(structure, epsilon_conversion_factor)
if dihedrals:
dihedral_types, unique_dihedral_types = _get_dihedral_types(
- structure, use_rb_torsions, use_dihedrals, epsilon_conversion_factor
+ structure,
+ use_rb_torsions,
+ use_dihedrals,
+ epsilon_conversion_factor,
+ zero_dihedral_weighting_factor,
)
if impropers:
@@ -305,7 +320,9 @@ def write_lammpsdata(
data.write("{:d} bonds\n".format(len(bonds)))
data.write("{:d} angles\n".format(len(angles)))
data.write("{:d} dihedrals\n".format(len(dihedrals)))
- data.write("{:d} impropers\n\n".format(len(impropers)))
+ data.write(
+ "{:d} impropers\n\n".format(len(impropers) + len(imp_dihedrals))
+ )
data.write("{:d} atom types\n".format(len(set(types))))
if atom_style in ["full", "molecular"]:
@@ -321,6 +338,10 @@ def write_lammpsdata(
data.write(
"{:d} improper types\n".format(len(set(improper_types)))
)
+ elif imp_dihedrals:
+ data.write(
+ "{:d} improper types\n".format(len(set(imp_dihedral_types)))
+ )
data.write("\n")
# Box data
@@ -669,6 +690,30 @@ def write_lammpsdata(
params[5],
)
)
+ elif imp_dihedrals:
+ # Improper dihedral coefficients
+ sorted_imp_dihedral_types = {
+ k: v
+ for k, v in sorted(
+ unique_imp_dihedral_types.items(),
+ key=lambda item: item[1],
+ )
+ }
+ data.write("\nImproper Coeffs # cvff\n")
+ data.write("#K, d, n\n")
+ for params, idx in sorted_imp_dihedral_types.items():
+ data.write(
+ "{}\t{:.5f}\t{:d}\t{:d}\t# {}\t{}\t{}\t{}\n".format(
+ idx,
+ params[0],
+ params[1],
+ params[2],
+ params[5],
+ params[6],
+ params[7],
+ params[8],
+ )
+ )
# Atom data
data.write("\nAtoms\n\n")
@@ -749,6 +794,22 @@ def write_lammpsdata(
improper[3],
)
)
+ elif imp_dihedrals:
+ data.write("\nImpropers\n\n")
+ for i, improper in enumerate(imp_dihedrals):
+ # The atoms are written central-atom third in LAMMPS data file.
+ # This is correct for AMBER impropers even though
+ # LAMMPS documentation implies central-atom-first.
+ data.write(
+ "{:d}\t{:d}\t{:d}\t{:d}\t{:d}\t{:d}\n".format(
+ i + 1,
+ imp_dihedral_types[i],
+ improper[0],
+ improper[1],
+ improper[2],
+ improper[3],
+ )
+ )
def _get_bond_types(
@@ -886,7 +947,11 @@ def _get_angle_types(
def _get_dihedral_types(
- structure, use_rb_torsions, use_dihedrals, epsilon_conversion_factor
+ structure,
+ use_rb_torsions,
+ use_dihedrals,
+ epsilon_conversion_factor,
+ zero_dihedral_weighting_factor,
):
lj_unit = 1 / epsilon_conversion_factor
if use_rb_torsions:
@@ -940,7 +1005,10 @@ def _get_dihedral_types(
structure.join_dihedrals()
for dihedral in structure.dihedrals:
if not dihedral.improper:
- weight = 1 / len(dihedral.type)
+ if zero_dihedral_weighting_factor:
+ weight = 0.0
+ else:
+ weight = 1.0 / len(dihedral.type)
for dih_type in dihedral.type:
charmm_dihedrals.append(
(
@@ -969,6 +1037,44 @@ def _get_dihedral_types(
return dihedral_types, unique_dihedral_types
+def _get_improper_dihedral_types(structure, epsilon_conversion_factor):
+ lj_unit = 1 / epsilon_conversion_factor
+ improper_dihedrals = []
+ for dihedral in structure.dihedrals:
+ if dihedral.improper:
+ dih_type = dihedral.type
+ phase = abs(int(round(dih_type.phase, 0)))
+ if not (phase == 0 or phase == 180):
+ raise ValueError("Improper dihedral phase must be 0 or 180")
+ if phase:
+ d = -1
+ else:
+ d = 1
+ improper_dihedrals.append(
+ (
+ round(dih_type.phi_k * lj_unit, 3),
+ d,
+ int(round(dih_type.per, 0)),
+ round(dih_type.scee, 1),
+ round(dih_type.scnb, 1),
+ dihedral.atom1.type,
+ dihedral.atom2.type,
+ dihedral.atom3.type,
+ dihedral.atom4.type,
+ )
+ )
+ unique_imp_dihedral_types = dict(enumerate(set(improper_dihedrals)))
+ unique_imp_dihedral_types = OrderedDict(
+ [(y, x + 1) for x, y in unique_imp_dihedral_types.items()]
+ )
+ imp_dihedral_types = [
+ unique_imp_dihedral_types[dihedral_info]
+ for dihedral_info in improper_dihedrals
+ ]
+
+ return imp_dihedral_types, unique_imp_dihedral_types
+
+
def _get_impropers(structure, epsilon_conversion_factor):
lj_unit = 1 / epsilon_conversion_factor
unique_improper_types = dict(
diff --git a/mbuild/tests/test_lammpsdata.py b/mbuild/tests/test_lammpsdata.py
index 0b57083e8..1bdf482a0 100755
--- a/mbuild/tests/test_lammpsdata.py
+++ b/mbuild/tests/test_lammpsdata.py
@@ -22,6 +22,10 @@ def test_save_forcefield(self, ethane, unit_style):
@pytest.mark.skipif(not has_foyer, reason="Foyer package not installed")
def test_save_charmm(self):
+ from foyer import Forcefield
+
+ from mbuild.formats.lammpsdata import write_lammpsdata
+
cmpd = mb.load(get_fn("charmm_dihedral.mol2"))
for i in cmpd.particles():
i.name = "_{}".format(i.name)
@@ -29,16 +33,12 @@ def test_save_charmm(self):
box=cmpd.get_boundingbox(),
residues=set([p.parent.name for p in cmpd.particles()]),
)
-
- from foyer import Forcefield
-
ff = Forcefield(forcefield_files=[get_fn("charmm_truncated.xml")])
structure = ff.apply(structure, assert_dihedral_params=False)
-
- from mbuild.formats.lammpsdata import write_lammpsdata
-
write_lammpsdata(structure, "charmm_dihedral.lammps")
out_lammps = open("charmm_dihedral.lammps", "r").readlines()
+ found_angles = False
+ found_dihedrals = False
for i, line in enumerate(out_lammps):
if "Angle Coeffs" in line:
assert "# charmm" in line
@@ -47,12 +47,129 @@ def test_save_charmm(self):
in out_lammps[i + 1]
)
assert len(out_lammps[i + 2].split("#")[0].split()) == 5
+ found_angles = True
+ elif "Dihedral Coeffs" in line:
+ assert "# charmm" in line
+ assert "#k, n, phi, weight" in out_lammps[i + 1]
+ assert len(out_lammps[i + 2].split("#")[0].split()) == 5
+ found_dihedrals = True
+ else:
+ pass
+ assert found_angles
+ assert found_dihedrals
+
+ @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed")
+ def test_singleterm_charmm(self):
+ from foyer import Forcefield
+
+ from mbuild.formats.lammpsdata import write_lammpsdata
+
+ cmpd = mb.load(get_fn("charmm_dihedral.mol2"))
+ for i in cmpd.particles():
+ i.name = "_{}".format(i.name)
+ structure = cmpd.to_parmed(
+ box=cmpd.get_boundingbox(),
+ residues=set([p.parent.name for p in cmpd.particles()]),
+ )
+ ff = Forcefield(
+ forcefield_files=[get_fn("charmm_truncated_singleterm.xml")]
+ )
+ structure = ff.apply(structure, assert_dihedral_params=False)
+ write_lammpsdata(structure, "charmm_dihedral_singleterm.lammps")
+ out_lammps = open("charmm_dihedral_singleterm.lammps", "r").readlines()
+ found_dihedrals = False
+ for i, line in enumerate(out_lammps):
+ if "Dihedral Coeffs" in line:
+ assert "# charmm" in line
+ assert "#k, n, phi, weight" in out_lammps[i + 1]
+ assert len(out_lammps[i + 2].split("#")[0].split()) == 5
+ assert float(
+ out_lammps[i + 2].split("#")[0].split()[4]
+ ) == float("1.0")
+ found_dihedrals = True
+ else:
+ pass
+ assert found_dihedrals
+
+ @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed")
+ def test_charmm_improper(self):
+ from foyer import Forcefield
+
+ import mbuild as mb
+ from mbuild.formats.lammpsdata import write_lammpsdata
+
+ system = mb.Compound()
+ first = mb.Particle(name="_CTL2", pos=[-1, 0, 0])
+ second = mb.Particle(name="_CL", pos=[0, 0, 0])
+ third = mb.Particle(name="_OBL", pos=[1, 0, 0])
+ fourth = mb.Particle(name="_OHL", pos=[0, 1, 0])
+
+ system.add([first, second, third, fourth])
+
+ system.add_bond((first, second))
+ system.add_bond((second, third))
+ system.add_bond((second, fourth))
+
+ ff = Forcefield(forcefield_files=[get_fn("charmm36_cooh.xml")])
+ struc = ff.apply(
+ system,
+ assert_angle_params=False,
+ assert_dihedral_params=False,
+ assert_improper_params=False,
+ )
+
+ write_lammpsdata(struc, "charmm_improper.lammps")
+ out_lammps = open("charmm_improper.lammps", "r").readlines()
+ found_impropers = False
+ for i, line in enumerate(out_lammps):
+ if "Improper Coeffs" in line:
+ assert "# harmonic" in line
+ assert "k, phi" in out_lammps[i + 1]
+ assert len(out_lammps[i + 2].split("#")[0].split()) == 3
+ found_impropers = True
+ assert found_impropers
+
+ @pytest.mark.skipif(not has_foyer, reason="Foyer package not installed")
+ def test_amber(self):
+ from foyer import Forcefield
+
+ from mbuild.formats.lammpsdata import write_lammpsdata
+
+ cmpd = mb.load("C1(=CC=CC=C1)F", smiles=True)
+
+ ff = Forcefield(forcefield_files=[get_fn("gaff_test.xml")])
+ structure = ff.apply(cmpd)
+
+ write_lammpsdata(
+ structure, "amber.lammps", zero_dihedral_weighting_factor=True
+ )
+ out_lammps = open("amber.lammps", "r").readlines()
+ found_angles = False
+ found_dihedrals = False
+ found_impropers = False
+ for i, line in enumerate(out_lammps):
+ if "Angle Coeffs" in line:
+ assert "# harmonic" in line
+ assert "#\treduced_k\t\ttheteq(deg)" in out_lammps[i + 1]
+ assert len(out_lammps[i + 2].split("#")[0].split()) == 3
+ found_angles = True
elif "Dihedral Coeffs" in line:
assert "# charmm" in line
assert "#k, n, phi, weight" in out_lammps[i + 1]
assert len(out_lammps[i + 2].split("#")[0].split()) == 5
+ assert float(out_lammps[i + 2].split("#")[0].split()[4]) == 0.0
+ found_dihedrals = True
+ elif "Improper Coeffs" in line:
+ assert "# cvff" in line
+ assert "#K, d, n" in out_lammps[i + 1]
+ assert len(out_lammps[i + 2].split("#")[0].split()) == 4
+ assert out_lammps[i + 2].split("#")[0].split()[2] == "-1"
+ found_impropers = True
else:
pass
+ assert found_angles
+ assert found_dihedrals
+ assert found_impropers
@pytest.mark.skipif(not has_foyer, reason="Foyer package not installed")
@pytest.mark.parametrize("unit_style", ["real", "lj"])
diff --git a/mbuild/utils/reference/charmm36_cooh.xml b/mbuild/utils/reference/charmm36_cooh.xml
new file mode 100644
index 000000000..322f9bc42
--- /dev/null
+++ b/mbuild/utils/reference/charmm36_cooh.xml
@@ -0,0 +1,51 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/mbuild/utils/reference/charmm_truncated_singleterm.xml b/mbuild/utils/reference/charmm_truncated_singleterm.xml
new file mode 100644
index 000000000..1869566a7
--- /dev/null
+++ b/mbuild/utils/reference/charmm_truncated_singleterm.xml
@@ -0,0 +1,43 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/mbuild/utils/reference/gaff_test.xml b/mbuild/utils/reference/gaff_test.xml
new file mode 100644
index 000000000..f6d6885cb
--- /dev/null
+++ b/mbuild/utils/reference/gaff_test.xml
@@ -0,0 +1,33 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+