From e43a88a7afda86ddcd90a500d5009b1f8a4e7bcf Mon Sep 17 00:00:00 2001
From: James Graham <J.A.Graham@soton.ac.uk>
Date: Fri, 5 Aug 2016 17:07:53 +0100
Subject: [PATCH] Force constant infinity if variance zero

---
 pycgtool/bondset.py  | 5 +++--
 pycgtool/mapping.py  | 3 ++-
 pycgtool/pycgtool.py | 6 ++++++
 test/test_bondset.py | 2 +-
 4 files changed, 12 insertions(+), 4 deletions(-)

diff --git a/pycgtool/bondset.py b/pycgtool/bondset.py
index 73bbfa8..3a4b568 100644
--- a/pycgtool/bondset.py
+++ b/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")
diff --git a/pycgtool/mapping.py b/pycgtool/mapping.py
index 10a6b75..6235300 100644
--- a/pycgtool/mapping.py
+++ b/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)
 
diff --git a/pycgtool/pycgtool.py b/pycgtool/pycgtool.py
index deae69e..efef173 100755
--- a/pycgtool/pycgtool.py
+++ b/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
diff --git a/test/test_bondset.py b/test/test_bondset.py
index 5ba6348..faf529e 100644
--- a/test/test_bondset.py
+++ b/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)
-- 
GitLab