diff --git a/pycgtool/bondset.py b/pycgtool/bondset.py index 0ce9e97..e8184db 100644 --- a/pycgtool/bondset.py +++ b/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: diff --git a/pycgtool/functionalforms.py b/pycgtool/functionalforms.py index 987d96d..7c35f78 100644 --- a/pycgtool/functionalforms.py +++ b/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. diff --git a/test/data/sugar_out.itp b/test/data/sugar_out.itp index 2a9a5b7..60a9f10 100644 --- a/test/data/sugar_out.itp +++ b/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 diff --git a/test/test_bondset.py b/test/test_bondset.py index 7394917..66eb451 100644 --- a/test/test_bondset.py +++ b/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): diff --git a/test/test_functionalforms.py b/test/test_functionalforms.py index 78a4464..b68da28 100644 --- a/test/test_functionalforms.py +++ b/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"