diff --git a/pycgtool/bondset.py b/pycgtool/bondset.py
index 1c436689890b0e887a3ad32cc4db6caaea05ce90..ce30d97efb97d10630960df9f24c3bf54525f8b4 100644
--- a/pycgtool/bondset.py
+++ b/pycgtool/bondset.py
@@ -15,7 +15,7 @@ try:
except ImportError:
pass
-from .util import stat_moments, sliding, dist_with_pbc, transpose_and_sample
+from .util import sliding, dist_with_pbc, transpose_and_sample
from .util import extend_graph_chain, backup_file
from .util import vector_len, vector_cross, vector_angle, vector_angle_signed
from .parsers.cfg import CFG
@@ -64,8 +64,6 @@ class Bond:
raise ValueError("No bonds were measured between atoms {0}".format(self.atoms))
values = np.array(self.values)
- if len(self.atoms) > 2:
- values = np.radians(values)
with np.errstate(divide="raise"):
self.eqm = self._func_form.eqm(values, temp)
@@ -75,9 +73,6 @@ class Bond:
# Happens when variance is 0, i.e. we only have one value
self.fconst = float("inf")
- if len(self.atoms) > 2:
- self.eqm = math.degrees(self.eqm)
-
def __repr__(self):
try:
return "<Bond containing atoms {0} with r_0 {1:.3f} and force constant {2:.3e}>".format(
@@ -269,13 +264,14 @@ class BondSet:
self._populate_atom_numbers(mapping)
backup_file(filename)
- def write_bond_angle_dih(bonds, section_header, itp, print_fconst=True, multiplicity=None):
+ def write_bond_angle_dih(bonds, section_header, itp, print_fconst=True, multiplicity=None, rad2deg=False):
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])
- line += " {0:4d} {1:12.5f}".format(1, bond.eqm)
+ eqm = math.degrees(bond.eqm) if rad2deg else bond.eqm
+ line += " {0:4d} {1:12.5f}".format(1, eqm)
if print_fconst:
line += " {0:12.5f}".format(bond.fconst)
if multiplicity is not None:
@@ -312,8 +308,8 @@ class BondSet:
), file=itp)
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_angles(mol), "angles", itp, rad2deg=True)
+ write_bond_angle_dih(self.get_bond_dihedrals(mol), "dihedrals", itp, multiplicity=1, rad2deg=True)
write_bond_angle_dih(self.get_bond_length_constraints(mol), "constraints", itp, print_fconst=False)
def apply(self, frame):
@@ -329,7 +325,7 @@ class BondSet:
def calc_angle(atoms):
veca = dist_with_pbc(atoms[0].coords, atoms[1].coords, frame.box)
vecb = dist_with_pbc(atoms[1].coords, atoms[2].coords, frame.box)
- return math.degrees(math.pi - vector_angle(veca, vecb))
+ return math.pi - vector_angle(veca, vecb)
def calc_dihedral(atoms):
veca = dist_with_pbc(atoms[0].coords, atoms[1].coords, frame.box)
@@ -339,7 +335,7 @@ class BondSet:
c1 = vector_cross(veca, vecb)
c2 = vector_cross(vecb, vecc)
- return math.degrees(vector_angle_signed(c1, c2, vecb))
+ return vector_angle_signed(c1, c2, vecb)
calc = {2: calc_length,
3: calc_angle,
@@ -389,6 +385,7 @@ class BondSet:
except ValueError:
pass
+ # TODO add test
def dump_values(self, target_number=10000):
"""
Output measured bond values to files for length, angles and dihedrals.
@@ -396,9 +393,11 @@ class BondSet:
:param target_number: Approx number of sample measurements to output. If None, all samples will be output
"""
- def write_bonds_to_file(bonds, filename):
+ def write_bonds_to_file(bonds, filename, rad2deg=False):
with open(filename, "w") as f:
for row in transpose_and_sample((bond.values for bond in bonds), n=target_number):
+ if rad2deg:
+ row = [math.degrees(val) for val in row]
print((len(row) * "{:12.5f}").format(*row), file=f)
for mol in self._molecules:
@@ -410,11 +409,11 @@ class BondSet:
bonds = self.get_bond_angles(mol)
if bonds:
- write_bonds_to_file(bonds, "{0}_angle.dat".format(mol))
+ write_bonds_to_file(bonds, "{0}_angle.dat".format(mol), rad2deg=True)
bonds = self.get_bond_dihedrals(mol)
if bonds:
- write_bonds_to_file(bonds, "{0}_dihedral.dat".format(mol))
+ write_bonds_to_file(bonds, "{0}_dihedral.dat".format(mol), rad2deg=True)
def __len__(self):
return len(self._molecules)
diff --git a/test/test_bondset.py b/test/test_bondset.py
index 62a7ca433b0193e1be80d8348f95bd9a84714bd4..c9a1e18d41396a798e80b0ace05853461e16954a 100644
--- a/test/test_bondset.py
+++ b/test/test_bondset.py
@@ -1,6 +1,7 @@
import unittest
import logging
+import math
from pycgtool.bondset import BondSet
from pycgtool.frame import Frame
@@ -25,24 +26,24 @@ class DummyOptions:
class BondSetTest(unittest.TestCase):
# Columns are: eqm value, standard fc, defaults fc, mixed 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),
- ( 77.882969526805, 503.24211, 25, 481.527),
- (116.081406627900, 837.76904, 25, 676.511),
- (111.030514958715, 732.87969, 25, 639.007),
- ( 83.284691301386, 945.32633, 25, 933.199),
- (143.479514279933, 771.63691, 25, 273.207),
- ( 99.293754667718, 799.82825, 25, 779.747),
- (-82.852665692244, 253.75691, 50, 1250),
- ( 61.159604648237, 125.04591, 50, 1250),
- (-21.401629717440, 135.50927, 50, 1250),
- ( 53.161150086611, 51.13975, 50, 1250),
- (-96.548945531698, 59.38162, 50, 1250),
- ( 75.370211843364, 279.80889, 50, 1250)
+ ( 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)
]
def test_bondset_create(self):
@@ -61,12 +62,14 @@ class BondSetTest(unittest.TestCase):
delta=0.2225376 / 500)
# Second six are angles
self.assertEqual(1, len(measure["ALLA"][6].values))
- self.assertAlmostEqual(77.22779289, measure["ALLA"][6].values[0],
- delta=77.22779289 / 500)
+ expected = math.radians(77.22779289)
+ self.assertAlmostEqual(expected, measure["ALLA"][6].values[0],
+ delta=expected / 500)
# Final six are dihedrals
self.assertEqual(1, len(measure["ALLA"][12].values))
- self.assertAlmostEqual(-89.5552903, measure["ALLA"][12].values[0],
- delta=89.552903 / 500)
+ expected = math.radians(-89.5552903)
+ self.assertAlmostEqual(expected, measure["ALLA"][12].values[0],
+ delta=abs(expected) / 500)
def test_bondset_remove_triangles(self):
bondset = BondSet("test/data/triangle.bnd", DummyOptions)