diff --git a/pycgtool.py b/pycgtool.py
index 2df400c1dff4b7f11d3d29bc2c254d02e67656b8..4fe8dc588c34218e8de7b687f9e5591aa11355ed 100755
--- a/pycgtool.py
+++ b/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"),
- ("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),
- ("output_forcefield", False),
- ("temperature", 310),
- ("angle_default_fc", True),
- ("generate_angles", True),
- ("generate_dihedrals", False),
- ("empirical_corr", False)],
- args)
+ 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", 10000),
+ ("output_forcefield", False),
+ ("temperature", 310),
+ ("default_fc", False),
+ ("generate_angles", True),
+ ("generate_dihedrals", False),
+ ("empirical_corr", False)
+ ], args)
if not args.map and not args.bnd:
parser.error("One or both of -m and -b is required.")
diff --git a/pycgtool/bondset.py b/pycgtool/bondset.py
index a6f6f7b3e207472c75f166108e6b9f097506a597..79928ab322d085f2754896f01ba6fe177875e1e4 100644
--- a/pycgtool/bondset.py
+++ b/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):
"""
diff --git a/pycgtool/pycgtool.py b/pycgtool/pycgtool.py
index efef1731bcf5878d8ffc4c41402bdf046d9ed6d2..f0f447392160270e2a62ce1a604b0fd91eb9a4e1 100755
--- a/pycgtool/pycgtool.py
+++ b/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")
diff --git a/test/data/sugar_out.itp b/test/data/sugar_out.itp
index 1fc051bbd25ecb080da05cba82d823b88d49b704..8e111740a646fe4404daffd95adaba00c3f971d6 100644
--- a/test/data/sugar_out.itp
+++ b/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
diff --git a/test/test_bondset.py b/test/test_bondset.py
index 14f7ce22bdd2383a026df25a7f61983e0dc5867e..d167245aeefe254ba2bc80256dd8322cbe989396 100644
--- a/test/test_bondset.py
+++ b/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))
diff --git a/test/test_pycgtool.py b/test/test_pycgtool.py
index f71f8aca045cddd641614f5a5df0768cfbb72e72..cfd6f6d84b3b60ade34eecaefff9667852158dcb 100644
--- a/test/test_pycgtool.py
+++ b/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)])