Use radians internally for angles

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)

test/test_bondset.py

 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)