Match equation in paper for BI - factor 2 and sin

Use default fc on bond lengths too if chosen

pycgtool.py

@@ -22,21 +22,22 @@ if __name__ == "__main__":
     input_files.add_argument('--end', type=int, default=-1, help="Frame number to end")
 
     args = parser.parse_args()
-    config = Options([("output_name", "out"),
+    config = Options([
+        ("output_name", "out"),
         ("output", "gro"),
         ("output_xtc", args.outputxtc),
         ("map_only", not bool(args.bnd)),
         ("map_center", "geom"),
         ("constr_threshold", 100000),
         ("dump_measurements", bool(args.bnd) and not bool(args.map)),
-                      ("dump_n_values", 100000),
+        ("dump_n_values", 10000),
         ("output_forcefield", False),
         ("temperature", 310),
-                      ("angle_default_fc", True),
+        ("default_fc", False),
         ("generate_angles", True),
         ("generate_dihedrals", False),
-                      ("empirical_corr", False)],
-                     args)
+        ("empirical_corr", False)
+    ], args)
 
     if not args.map and not args.bnd:
         parser.error("One or both of -m and -b is required.")

pycgtool/bondset.py

@@ -47,23 +47,20 @@ class Bond:
     def __len__(self):
         return len(self.atoms)
 
-    def boltzmann_invert(self, temp=310, angle_default_fc=True):
+    def boltzmann_invert(self, temp=310, default_fc=True):
         """
-        Perform Boltzmann Inversion using measured values of this bond to calculate equilibrium value and force constant.
+        Perform Boltzmann Inversion using measured values of bond to calculate equilibrium value and force constant.
 
         :param temp: Temperature at which the simulation was performed
-        :param angle_default_fc: Use default value of 25 kJ mol-1 rad-2 for angle force constants
+        :param default_fc: Use default value of 1250 kJ mol-1 nm-2 for length fc and 25 kJ mol-1 for angle fc
         """
         mean, var = stat_moments(self.values)
 
         rt = 8.314 * temp / 1000.
         rad2 = np.pi * np.pi / (180. * 180.)
-        conv = {2: lambda: rt / var,
-                3: lambda: 25. if angle_default_fc else rt * np.sin(np.radians(mean))**2 / (var * rad2),
-                4: lambda: rt / (var * rad2)}
-        # conv = {2: lambda: rt / (2 * var),
-        #         3: lambda: 25. if angle_default_fc else rt / (2 * np.sin(np.radians(mean))**2 * var * rad2),
-        #         4: lambda: rt / (var * rad2)}
+        conv = {2: lambda: 1250. if default_fc else rt / (2 * var),
+                3: lambda: 25. if default_fc else rt / (2 * np.sin(np.radians(mean))**2 * var * rad2),
+                4: lambda: rt / (2 * var * rad2)}
 
         self.eqm = mean
         try:
@@ -115,9 +112,9 @@ class BondSet:
             self._temperature = 310.
 
         try:
-            self._angle_default_fc = options.angle_default_fc
+            self._default_fc = options.default_fc
         except AttributeError:
-            self._angle_default_fc = True
+            self._default_fc = False
 
         with CFG(filename) as cfg:
             for mol in cfg:
@@ -237,16 +234,17 @@ class BondSet:
         self._populate_atom_numbers(mapping)
         backup_file(filename, verbose=True)
 
-        def write_bond_angle_dih(bonds, section_header, itp, fconst=True, multiplicity=None):
+        def write_bond_angle_dih(bonds, section_header, itp, print_fconst=True, multiplicity=None):
             if bonds:
                 print("\n[ {0:s} ]".format(section_header), file=itp)
             for bond in bonds:
                 # Factor is usually 1, unless doing correction
                 eqm = bond.eqm * self._empirical_correction_factor
+                fconst = bond.fconst * self._empirical_correction_factor
                 line = " ".join(["{0:4d}".format(atnum + 1) for atnum in bond.atom_numbers])
                 line += " {0:4d} {1:12.5f}".format(1, eqm)
-                if fconst:
-                    line += " {0:12.5f}".format(bond.fconst)
+                if print_fconst:
+                    line += " {0:12.5f}".format(fconst)
                 if multiplicity is not None:
                     line += " {0:4d}".format(multiplicity)
                 print(line, file=itp)
@@ -275,7 +273,7 @@ class BondSet:
                 write_bond_angle_dih(self.get_bond_lengths(mol), "bonds", itp)
                 write_bond_angle_dih(self.get_bond_angles(mol), "angles", itp)
                 write_bond_angle_dih(self.get_bond_dihedrals(mol), "dihedrals", itp, multiplicity=1)
-                write_bond_angle_dih(self.get_bond_length_constraints(mol), "constraints", itp, fconst=False)
+                write_bond_angle_dih(self.get_bond_length_constraints(mol), "constraints", itp, print_fconst=False)
 
     def apply(self, frame):
         """
@@ -343,7 +341,7 @@ class BondSet:
 
         for bond in bond_iter_wrap:
             bond.boltzmann_invert(temp=self._temperature,
-                                  angle_default_fc=self._angle_default_fc)
+                                  default_fc=self._default_fc)
 
     def dump_values(self, target_number=10000):
         """

pycgtool/pycgtool.py

@@ -63,6 +63,8 @@ def main(args, config):
         if args.map:
             logger.info("Beginning Boltzmann inversion")
             bonds.boltzmann_invert(progress=True)
+            print("Note: charges are not currently handled")
+            print("If your CG molecule should contain charges the itp will need to be edited")
             if config.output_forcefield:
                 logger.info("Creating GROMACS forcefield directory")
                 ff = ForceField("ff" + config.output_name + ".ff")

test/data/sugar_out.itp

@@ -17,27 +17,27 @@ ALLA    1
    6 P4      1 ALLA O5      6    0.000
 
 [ bonds ]
-   2    3    1      0.21998  35162.10022
-   3    4    1      0.21648  65545.10313
-   4    5    1      0.26454  19022.17121
-   5    6    1      0.20452  53457.20659
+   1    2    1      0.22209  88271.49992
+   2    3    1      0.21998  17581.05011
+   3    4    1      0.21648  32772.55156
+   4    5    1      0.26454   9511.08560
+   5    6    1      0.20452  26728.60329
 
 [ angles ]
-   1    2    3    1     77.69557    383.16060
-   2    3    4    1    115.59599    262.66648
-   3    4    5    1    108.51038    698.56915
-   4    5    6    1     81.23827   1098.41706
-   5    6    1    1    146.49898     48.08599
-   6    1    2    1     99.84045   1571.73513
+   1    2    3    1     77.69557    210.24268
+   2    3    4    1    115.59599    198.52391
+   3    4    5    1    108.51038    431.97481
+   4    5    6    1     81.23827    575.61056
+   5    6    1    1    146.49898    259.05013
+   6    1    2    1     99.84045    833.86832
 
 [ dihedrals ]
-   1    2    3    4    1    -85.45004    234.10405    1
-   2    3    4    5    1     65.21378    117.46132    1
-   3    4    5    6    1    -22.73796     94.20175    1
-   4    5    6    1    1     55.18087     37.92107    1
-   5    6    1    2    1    -96.52349     40.73340    1
-   6    1    2    3    1     73.82477    259.20126    1
+   1    2    3    4    1    -85.45004    117.05203    1
+   2    3    4    5    1     65.21378     58.73066    1
+   3    4    5    6    1    -22.73796     47.10087    1
+   4    5    6    1    1     55.18087     18.96053    1
+   5    6    1    2    1    -96.52349     20.36670    1
+   6    1    2    3    1     73.82477    129.60063    1
 
 [ constraints ]
-   1    2    1      0.22209
    6    1    1      0.17931

test/test_bondset.py

@@ -48,24 +48,24 @@ class BondSetTest(unittest.TestCase):
         self.assertEqual(0, len(angles))
 
     def test_bondset_boltzmann_invert(self):
-        ref = [(0.222817161647, 19116816.3363),
-               (0.216766804211, 55835643.1619),
-               (0.223815134299, 6199518.81029),
-               (0.239439104442, 642593936.491),
-               (0.215497208685, 7366326.46847),
-               (0.179544192953, 20521533.3397),
-               (76.9617550012, 23447621.6268),
-               (116.009786726, 1210005.81229),
-               (109.247878273, 311284.586971),
-               (84.4955293511, 1713127.90333),
-               (149.557001531, 62116.1326728),
-               (99.0524708671, 192028.755165),
-               (-89.5321545897, 182481.962675),
-               (70.040930903, 337940.374373),
-               (-21.1194544981, 1990817.75509),
-               (48.3727747848, 71424.537429),
-               (-85.923133935, 124765.874439),
-               (70.3195444564, 1555761.24713)]
+        ref = [(  0.222817161647,   9560719),
+               (  0.216766804211,  27930081),
+               (  0.223815134299,   3099614),
+               (  0.239439104442, 322294637),
+               (  0.215497208685,   3683331),
+               (  0.179544192953,  10260492),
+               ( 76.9617550012,    13048242),
+               (116.009786726,       927619),
+               (109.247878273,       195899),
+               ( 84.4955293511,      871991),
+               (149.557001531,       471274),
+               ( 99.0524708671,      100943),
+               (-89.5321545897,       91248),
+               ( 70.040930903,       169002),
+               (-21.1194544981,      995309),
+               ( 48.3727747848,       35712),
+               (-85.923133935,        62394),
+               ( 70.3195444564,      778389)]
 
         measure = BondSet("test/data/sugar.bnd", DummyOptions)
         frame = Frame("test/data/sugar.gro", xtc="test/data/sugar.xtc")
@@ -78,7 +78,7 @@ class BondSetTest(unittest.TestCase):
         measure.boltzmann_invert()
         for i, bond in enumerate(measure["ALLA"]):
             # Require accuracy to 0.5%
-            # Allows for slight modifications to algorithm
+            # Allows for slight modifications to code
             self.assertAlmostEqual(ref[i][0], bond.eqm, delta=abs(ref[i][0] / 200))
             self.assertAlmostEqual(ref[i][1], bond.fconst, delta=abs(ref[i][1] / 200))
 

test/test_pycgtool.py

@@ -38,7 +38,7 @@ class PycgtoolTest(unittest.TestCase):
                       ("dump_n_values", 100000),
                       ("output_forcefield", False),
                       ("temperature", 310),
-                      ("angle_default_fc", True),
+                      ("default_fc", False),
                       ("generate_angles", True),
                       ("generate_dihedrals", False)])