Force constant infinity if variance zero

pycgtool/bondset.py

@@ -26,7 +26,7 @@ class Bond:
     """
     Class holding the properties of a single bonded term.
 
-    Bond lengths, angles and dihedrals are all equivalent, distinguised by the number of atoms present.
+    Bond lengths, angles and dihedrals are all equivalent, distinguished by the number of atoms present.
     """
     __slots__ = ["atoms", "atom_numbers", "values", "eqm", "fconst"]
 
@@ -69,7 +69,8 @@ class Bond:
         try:
             self.fconst = conv[len(self.atoms)]()
         except FloatingPointError:
-            self.fconst = 0
+            # Happens when variance is 0, i.e. infinitely sharp peak
+            self.fconst = float("inf")
 
     def r_squared(self):
         raise NotImplementedError("Bond r-squared is not yet implemented")

pycgtool/mapping.py

@@ -155,6 +155,7 @@ class Mapping:
         aa_residues = (aares for aares in frame if select_predicate(aares))
 
         if cgframe is None:
+            # Frame needs initialising
             aa_residues = list(aa_residues)
             cgframe = self._cg_frame_setup(aa_residues, frame.name)
 
@@ -177,7 +178,7 @@ class Mapping:
                         if self._map_center == "mass":
                             e.args = ("Error with mapping type 'mass', did you provide an itp file?",)
                         else:
-                            e.args = ("Error, unknown mapping type '{0}'".format(e.args[0]),)
+                            e.args = ("Error with mapping type '{0}', unknown mapping type.".format(e.args[0]),)
                         raise
                     bead.coords = calc_coords_weight(ref_coords, coords, cgframe.box, weights)
 

pycgtool/pycgtool.py

@@ -32,8 +32,14 @@ def main(args, config):
         cgframe.output(config.output_name + ".gro", format=config.output)
         logger.info("Mapping will be performed")
     else:
+        cgframe = frame
         logger.info("Mapping will not be performed")
 
+    # Only measure bonds from GRO frame if no XTC is provided
+    # Allows the user to get a topology from a single snapshot
+    if args.bnd and args.xtc is None:
+        bonds.apply(cgframe)
+
     # Main loop - perform mapping and measurement on every frame in XTC
     def main_loop():
         nonlocal cgframe

test/test_bondset.py

@@ -102,4 +102,4 @@ class BondSetTest(unittest.TestCase):
         bondset.boltzmann_invert()
         for bond in bondset.get_bond_lengths("ETH", True):
             self.assertAlmostEqual(1., bond.eqm)
-            self.assertEqual(0., bond.fconst)
+            self.assertEqual(float("inf"), bond.fconst)