Use GROMACS potential type ids in itp file output

Author: jag1g13

pycgtool/bondset.py

@@ -31,7 +31,7 @@ class Bond:
 
     Bond lengths, angles and dihedrals are all equivalent, distinguished by the number of atoms present.
     """
-    __slots__ = ["atoms", "atom_numbers", "values", "eqm", "fconst", "_func_form"]
+    __slots__ = ["atoms", "atom_numbers", "values", "eqm", "fconst", "gromacs_type_id", "_func_form"]
 
     def __init__(self, atoms, atom_numbers=None, func_form=None):
         """
@@ -48,6 +48,7 @@ def __init__(self, atoms, atom_numbers=None, func_form=None):
         self.fconst = None
 
         self._func_form = func_form
+        self.gromacs_type_id = func_form.gromacs_type_id_by_natoms(len(atoms))
 
     def __len__(self):
         return len(self.atoms)
@@ -273,10 +274,9 @@ def write_bond_angle_dih(bonds, section_header, itp, print_fconst=True, multipli
             if bonds:
                 print("\n[ {0:s} ]".format(section_header), file=itp)
             for bond in bonds:
-                # Factor is usually 1, unless doing correction
                 line = " ".join(["{0:4d}".format(atnum + 1) for atnum in bond.atom_numbers])
                 eqm = math.degrees(bond.eqm) if rad2deg else bond.eqm
-                line += " {0:4d} {1:12.5f}".format(1, eqm)
+                line += " {0:4d} {1:12.5f}".format(bond.gromacs_type_id, eqm)
                 if print_fconst:
                     line += " {0:12.5f}".format(bond.fconst)
                 if multiplicity is not None:

pycgtool/functionalforms.py

@@ -67,7 +67,6 @@ def eqm(values, temp):
         """
         return np.nanmean(values)
 
-    @staticmethod
     @abc.abstractstaticmethod
     def fconst(values, temp):
         """
@@ -78,10 +77,34 @@ def fconst(values, temp):
         :param temp: Temperature of simulation
         :return: Calculated force constant
         """
-        pass
+        raise NotImplementedError
+
+    @abc.abstractproperty
+    def gromacs_type_ids(self):
+        """
+        Return tuple of GROMACS potential type ids when used as length, angle, dihedral.
+        
+        :return tuple[int]: Tuple of GROMACS potential type ids
+        """
+        raise NotImplementedError
+
+    @classmethod
+    def gromacs_type_id_by_natoms(cls, natoms):
+        """
+        Return the GROMACS potential type id for this functional form when used with natoms.
+        
+        :param int natoms: 
+        :return int: GROMACS potential type id 
+        """
+        tipe = cls.gromacs_type_ids[natoms - 2]
+        if tipe is None:
+            raise TypeError("The functional form {0} does not have a defined GROMACS potential type when used with {1} atoms.".format(cls.__name__, natoms))
+        return tipe
 
 
 class Harmonic(FunctionalForm):
+    gromacs_type_ids = (1, 1, 1)  # Consider whether to use improper (type 2) instead, it is actually harmonic
+
     @staticmethod
     def fconst(values, temp):
         rt = 8.314 * temp / 1000.
@@ -90,6 +113,8 @@ def fconst(values, temp):
 
 
 class CosHarmonic(FunctionalForm):
+    gromacs_type_ids = (None, 2, None)
+
     @staticmethod
     def fconst(values, temp):
         rt = 8.314 * temp / 1000.
@@ -99,18 +124,24 @@ def fconst(values, temp):
 
 
 class MartiniDefaultLength(FunctionalForm):
+    gromacs_type_ids = (1, None, None)
+
     @staticmethod
     def fconst(values, temp):
         return 1250.
 
 
 class MartiniDefaultAngle(FunctionalForm):
+    gromacs_type_ids = (None, 2, None)
+
     @staticmethod
     def fconst(values, temp):
         return 25.
 
 
 class MartiniDefaultDihedral(FunctionalForm):
+    gromacs_type_ids = (None, None, 1)
+
     @staticmethod
     def fconst(values, temp):
         return 50.

test/data/sugar_out.itp

@@ -24,12 +24,12 @@ ALLA    1
    5    6    1      0.20597  55319.39933
 
 [ angles ]
-   1    2    3    1     77.88207    503.72544
-   2    3    4    1    116.08126    838.60326
-   3    4    5    1    111.02889    733.45020
-   4    5    6    1     83.28583    946.13165
-   5    6    1    1    143.48577    771.65933
-   6    1    2    1     99.29372    800.62783
+   1    2    3    2     77.88207    503.72544
+   2    3    4    2    116.08126    838.60326
+   3    4    5    2    111.02889    733.45020
+   4    5    6    2     83.28583    946.13165
+   5    6    1    2    143.48577    771.65933
+   6    1    2    2     99.29372    800.62783
 
 [ dihedrals ]
    1    2    3    4    1    -82.85912    253.69305    1

test/test_bondset.py

@@ -24,26 +24,26 @@ class DummyOptions:
 
 
 class BondSetTest(unittest.TestCase):
-    # Columns are: eqm value, standard fc, defaults fc, mixed fc
+    # Columns are: eqm value, standard fc, MARTINI defaults fc
     invert_test_ref_data = [
-        ( 0.220474419132,  72658.18163, 1250, 1520530.416),
-        ( 0.221844516963,  64300.01188, 1250, 1328761.015),
-        ( 0.216313356504,  67934.93368, 1250, 1474281.672),
-        ( 0.253166204438,  19545.27388, 1250,  311446.690),
-        ( 0.205958461836,  55359.06367, 1250, 1322605.992),
-        ( 0.180550961226, 139643.66601, 1250, 4334538.941),
-        ( 1.359314249473,    503.24211,   25,     481.527),
-        ( 2.026002746003,    837.76904,   25,     676.511),
-        ( 1.937848056214,    732.87969,   25,     639.007),
-        ( 1.453592079716,    945.32633,   25,     933.199),
-        ( 2.504189933347,    771.63691,   25,     273.207),
-        ( 1.733002945619,    799.82825,   25,     779.747),
-        (-1.446051810383,    253.75691,   50,    1250),
-        ( 1.067436470329,    125.04591,   50,    1250),
-        (-0.373528903861,    135.50927,   50,    1250),
-        ( 0.927837103158,     51.13975,   50,    1250),
-        (-1.685096988856,     59.38162,   50,    1250),
-        ( 1.315458354592,    279.80889,   50,    1250)
+        ( 0.220474419132,  72658.18163, 1250),
+        ( 0.221844516963,  64300.01188, 1250),
+        ( 0.216313356504,  67934.93368, 1250),
+        ( 0.253166204438,  19545.27388, 1250),
+        ( 0.205958461836,  55359.06367, 1250),
+        ( 0.180550961226, 139643.66601, 1250),
+        ( 1.359314249473,    503.24211,   25),
+        ( 2.026002746003,    837.76904,   25),
+        ( 1.937848056214,    732.87969,   25),
+        ( 1.453592079716,    945.32633,   25),
+        ( 2.504189933347,    771.63691,   25),
+        ( 1.733002945619,    799.82825,   25),
+        (-1.446051810383,    253.75691,   50),
+        ( 1.067436470329,    125.04591,   50),
+        (-0.373528903861,    135.50927,   50),
+        ( 0.927837103158,     51.13975,   50),
+        (-1.685096988856,     59.38162,   50),
+        ( 1.315458354592,    279.80889,   50)
     ]
 
     def test_bondset_create(self):
@@ -119,11 +119,11 @@ class DefaultOptions(DummyOptions):
         measure.boltzmann_invert()
         self.support_check_mean_fc(measure["ALLA"], 2)
 
-    def test_bondset_boltzmann_invert_func_forms(self):
+    def test_bondset_boltzmann_invert_manual_default_fc(self):
         class FuncFormOptions(DummyOptions):
-            length_form = "CosHarmonic"
-            angle_form = "Harmonic"
-            dihedral_form = "MartiniDefaultLength"
+            length_form = "MartiniDefaultLength"
+            angle_form = "MartiniDefaultAngle"
+            dihedral_form = "MartiniDefaultDihedral"
 
         measure = BondSet("test/data/sugar.bnd", FuncFormOptions)
         frame = Frame("test/data/sugar.gro", xtc="test/data/sugar.xtc")
@@ -135,7 +135,7 @@ class FuncFormOptions(DummyOptions):
             measure.apply(cgframe)
 
         measure.boltzmann_invert()
-        self.support_check_mean_fc(measure["ALLA"], 3)
+        self.support_check_mean_fc(measure["ALLA"], 2)
 
     @unittest.skipIf(not mdtraj_present, "MDTRAJ or Scipy not present")
     def test_bondset_boltzmann_invert_mdtraj(self):

test/test_functionalforms.py

@@ -13,6 +13,8 @@ def test_functional_form(self):
 
     def test_functional_form_new(self):
         class TestFunc(FunctionalForm):
+            gromacs_type_ids = (None, None, None)
+
             @staticmethod
             def eqm(values, temp):
                 return "TestResultEqm"