Charmm parameters writer

gmso/formats/__init__.py

@@ -10,6 +10,7 @@
 from .lammpsdata import write_lammpsdata
 from .mcf import write_mcf
 from .mol2 import read_mol2
+from .prm_writer import write_prm
 from .top import write_top
 from .xyz import read_xyz, write_xyz
 

gmso/formats/prm_writer.py

@@ -0,0 +1,124 @@
+"""CHARMM Par format."""
+
+from gmso.formats.formats_registry import saves_as
+from gmso.utils.units import LAMMPS_UnitSystems
+
+
+@saves_as(".prm", ".par")
+def write_prm(topology, filename):
+    """Write CHARMM Parameter .prm file given a parametrized structure.
+
+    Notes
+    -----
+    Follows format according to
+    https://www.ks.uiuc.edu/Training/Tutorials/namd/namd-tutorial-unix-html/
+    node25.html
+    Furthermore, ParmEd can support writing CHARMM par, rtf, str files
+    by converting the parmed.Structure into parmed.CharmmParameterSet
+
+    Parmed stores rmin/2 in "rmin"
+    """
+    unit_system = LAMMPS_UnitSystems("real")  # ang, kcal/mol, amu
+    # ATOMS
+    with open(filename, "w") as f:
+        f.write("ATOMS\n")
+        unique_atoms = set()
+        for site in topology.sites:
+            unique_atoms.add(
+                (
+                    site.atom_type.name,
+                    unit_system.convert_parameter(
+                        site.atom_type.mass, n_decimals=6, name="mass"
+                    ),
+                )
+            )
+        for atom in unique_atoms:
+            f.write("MASS -1 {:8s} {:8s}\n".format(atom[0], atom[1]))
+
+        # TODO: Works for harmonic bonds
+        f.write("\nBONDS\n")
+        for bond in topology.bond_types:
+            atom1, atom2 = bond.member_types
+            f.write(
+                "{:8s} {:8s} {:.5f} {:.5f}\n".format(
+                    atom1, atom2, bond.parameters["k"], bond.parameters["r_eq"]
+                )
+            )
+
+        f.write("\nANGLES\n")
+        for angle in topology.angle_types:
+            atom1, atom2, atom3 = angle.member_types
+            if angle.name == "UreyBradleyAnglePotential":
+                f.write(
+                    "{:8s} {:8s} {:8s} {:.5f} {:.5f} {:.5f} {:.5f}\n".format(
+                        atom1,
+                        atom2,
+                        atom3,
+                        angle.parameters["k"],
+                        angle.parameters["theta_eq"],
+                        angle.parameters["kub"],
+                        angle.parameters["r_eq"],
+                    )
+                )
+            else:  # assume harmonic style:
+                f.write(
+                    "{:8s} {:8s} {:8s} {:.5f} {:.5f}\n".format(
+                        atom1,
+                        atom2,
+                        atom3,
+                        angle.parameters["k"],
+                        angle.parameters["theta_eq"],
+                    )
+                )
+
+        # These dihedrals need to be PeriodicTorsion Style (Charmm style)
+        f.write("\nDIHEDRALS\n")
+        for dihedral in topology.dihedral_types:
+            # works for PeriodicTorsion
+            atom1, atom2, atom3, atom4 = dihedral.member_types
+            f.write(
+                "{:8s} {:8s} {:8s} {:8s} {:.5f} {:5f} {:.5f}\n".format(
+                    atom1,
+                    atom2,
+                    atom3,
+                    atom4,
+                    dihedral.parameters["k"][0],
+                    dihedral.parameters["n"][0],
+                    dihedral.parameters["phi_eq"][0],
+                )
+            )
+
+        # TODO: No support for harmonic impropers
+
+        f.write("\nIMPROPER\n")
+        for improper in topology.improper_types:
+            atom1, atom2, atom3, atom4 = improper.member_types
+            f.write(
+                "{:8s} {:8s} {:8s} {:8s} {:.5f} {:5f} {:.5f}\n".format(
+                    atom1,
+                    atom2,
+                    atom3,
+                    atom4,
+                    improper.parameters["phi_k"][0],
+                    improper.parameters["n"][0],
+                    improper.parameters["phi_eq"][0],
+                )
+            )
+
+        f.write("\nNONBONDED\n")
+        nonbonded14 = topology.scaling_factors[0][2]
+
+        for atype in topology.atom_types:
+            # atype, 0.0, epsilon, rmin/2, 0.0, epsilon(1-4), rmin/2 (1-4)
+            f.write(
+                "{:8s} {:8.3f} {:8.3f} {:8.3f} {:8.3f} {:8.3f} {:8.3f}\n".format(
+                    atype.name,
+                    0.0,  # ignored,
+                    atype.parameters["epsilon"],
+                    atype.parameters["sigma"],
+                    0.0,
+                    atype.parameters["epsilon"] * nonbonded14,
+                    atype.parameters["sigma"] * nonbonded14,
+                )
+            )
+        f.write("\nEND")

gmso/lib/jsons/UreyBradleyAnglePotential.json

@@ -0,0 +1,11 @@
+{
+  "name": "UreyBradleyPotential",
+  "expression": "0.5 * k * (theta-theta_eq)**2 + 0.5 * kub * (r-r_eq)**2",
+  "independent_variables": ["theta", "r"],
+  "expected_parameters_dimensions": {
+    "k":"energy/angle**2",
+    "theta_eq": "angle",
+    "kub": "energy/length**2",
+    "r_eq": "length"
+  }
+}

gmso/tests/test_par.py

@@ -0,0 +1,45 @@
+from gmso import ForceField, Topology
+from gmso.parameterization import apply
+from gmso.tests.base_test import BaseTest
+from gmso.utils.io import get_fn
+
+# TODO: Sorting of atomtypes info in connection types
+# TODO: Test of correct potential forms
+# TODO: Handle iterating over unique types
+# TODO: Handle multiple values in charmm dihedral k as np array
+
+
+class TestPar(BaseTest):
+    def test_save_charmm(self):
+        top = Topology.load(get_fn("charmm_dihedral.mol2"))
+        for site in top.sites:
+            site.name = f"_{site.name}"
+
+        ff = ForceField(get_fn("charmm_truncated.xml"))
+        ptop = apply(
+            top, ff, identify_connections=True, ignore_params=["dihedral", "improper"]
+        )
+        ptop.save("charmm_dihedral.prm")
+
+    def test_par_parameters(self, ethane, oplsaa_forcefield):
+        from gmso.parameterization import apply
+
+        ptop = apply(ethane, oplsaa_forcefield)
+
+        ptop.save(filename="ethane-opls.par")
+        from parmed.charmm import CharmmParameterSet
+
+        pset = CharmmParameterSet.load_set(pfile="ethane-opls.par")
+        assert len(pset.bond_types) == 3
+        assert len(pset.angle_types) == 3
+        assert len(pset.atom_types) == 2
+
+    def test_parameter_forms(self, ethane):
+        # test that the values are correct
+        # write mbuild file and compare to gmso file
+        pass
+
+    def test_raise_errors(self):
+        # only takes harmonic bonds
+        # only takes harmonic of ub angles
+        pass

gmso/utils/files/charmm_dihedral.mol2

@@ -0,0 +1,17 @@
+@<TRIPOS>MOLECULE
+DPPC
+4 3 1 0 1
+SMALL
+USER_CHARGES
+@<TRIPOS>CRYSIN
+   13.9528    17.1622    31.1323    90.0000    90.0000    90.0000  1  1
+@<TRIPOS>ATOM
+      1 CTL1          -0.1558     6.4170     4.4541 CTL1            1 DPPC       1.500000
+      2 OSL         -0.7278     7.6038     3.8783 OSL           1 DPPC      -0.780000
+      3 PL         -0.5108     6.0856     5.8561 PL           1 DPPC      -0.780000
+      4 OSLP        -0.2550     5.1063     3.5611 OSLP          1 DPPC      -0.570000
+
+@<TRIPOS>BOND
+        1 1 2 1
+        2 2 3 1
+        3 1 4 1

gmso/utils/files/charmm_truncated.xml

@@ -0,0 +1,44 @@
+<ForceField><AtomTypes>
+    <Type class="CTL1" element="_CTL1" mass="12.011" name="CTL1" def="[_CTL1]"/>
+    <Type class="OSL" element="_OSL" mass="16.008" name="OSL" def="[_OSL]"/>
+    <Type class="PL" element="_PL" mass="15" name="PL" def="[_PL]"/>
+    <Type class="OSLP" element="_OSLP" mass="16.008" name="OSLP" def="[_OSLP]"/>
+  </AtomTypes>
+  <HarmonicBondForce>
+    <Bond k="1.00000" length="1" type1="CTL1" type2="CTL1"/>
+    <Bond k="2" length="2." type1="CTL1" type2="OSL"/>
+    <Bond k="3.00000" length="3" type1="CTL1" type2="PL"/>
+    <Bond k="4.00000" length="4" type1="CTL1" type2="OSLP"/>
+    <Bond k="5.00000" length="5" type1="OSL" type2="PL"/>
+    <Bond k="6.00000" length="6" type1="OSL" type2="OSL"/>
+    <Bond k="7.00000" length="7" type1="OSL" type2="OSLP"/>
+    <Bond k="8.00000" length="8" type1="PL" type2="OSLP"/>
+    <Bond k="9.00000" length="9" type1="PL" type2="PL"/>
+    <Bond k="10.00000" length="10." type1="OSLP" type2="OSLP"/>
+  </HarmonicBondForce>
+  <HarmonicAngleForce>
+    <Angle angle="1.93732" k="418.40000" type1="CTL1" type2="OSL" type3="PL"/>
+    <Angle angle="2.93732" k="518.40000" type1="CTL1" type2="OSL" type3="OSLP"/>
+    <Angle angle="3.93732" k="618.40000" type1="OSLP" type2="OSL" type3="PL"/>
+    <Angle angle="4.93732" k="718.40000" type1="OSLP" type2="CTL1" type3="OSL"/>
+  </HarmonicAngleForce>
+  <AmoebaUreyBradleyForce>
+    <UreyBradley d="0.21400" k="20920.00000" class1="CTL1" class2="OSL" class3="PL"/>
+     </AmoebaUreyBradleyForce>
+  <PeriodicTorsionForce>
+      <!--<Proper k1="69.00000" periodicity1="1" phase1="0.00000" type1="PL" type2="OSL" type3="CTL1" type4="OSLP"/>-->
+      <Proper k1="69.00000" periodicity1="1" phase1="0.00000" type1="PL" type2="OSL" type3="CTL1" type4="OSLP" k2="138.00000" periodicity2="2" phase2="45.00000" k3="138.00000" periodicity3="3" phase3="90.00000"/>
+  </PeriodicTorsionForce>
+  <CustomTorsionForce energy="k*(theta-theta0)^2">
+    <PerTorsionParameter name="k"/>
+    <PerTorsionParameter name="theta0"/>
+    <Improper k="200.83200" theta0="0.00000" type1="PL" type2="OSL" type3="CTL1" type4="OSLP"/>
+    <Improper k="300.83200" theta0="1.00000" type1="OSL" type2="PL" type3="CTL1" type4="OSLP"/>
+ </CustomTorsionForce>
+  <NonbondedForce coulomb14scale="1.0" lj14scale="1.0">
+    <Atom epsilon="1.192464" sigma="1.04000135244450124" type="CTL1" charge="0.0"/>
+    <Atom epsilon="2.192464" sigma="2.04000135244450124" type="OSL" charge="0.0"/>
+    <Atom epsilon="3.192464" sigma="3.04000135244450124" type="PL" charge="0.0"/>
+    <Atom epsilon="4.192464" sigma="4.04000135244450124" type="OSLP" charge="0.0"/>
+    </NonbondedForce>
+  </ForceField>