Correctly set OpenMM constraint tolerance, fix its tests

--testall

moldesign/_tests/helpers.py

@@ -146,9 +146,11 @@ def assert_something_resembling_minimization_happened(p0, e0, traj, mol):
 def assert_something_resembling_dynamics_happened(traj, mol, p0, t0, duration):
     """ Checks that dynamics routines have produced a consistent trajectory
     """
-    frame_interval = mol.integrator.frame_interval
+    frame_interval = mol.integrator.params.frame_interval
+    timestep = mol.integrator.params.timestep
+    epsilon_t = 1e-5 * timestep  # numerical noise in time
     if isinstance(frame_interval, int):
-        frame_interval *= mol.integrator.timestep
+        frame_interval *= mol.integrator.params.timestep
 
     assert mol is traj.mol
     assert traj.time[0] == t0
@@ -170,28 +172,44 @@ def assert_something_resembling_dynamics_happened(traj, mol, p0, t0, duration):
 
         if istep == len(traj) - 1:
             # last frame - it just needs
-            assert step.time >= expected_end - 1e-5 * u.timestep
+            assert step.time >= expected_end - epsilon_t
             break
 
         elif istep != 0:
-            assert (step.time - lasttime - frame_interval) < 1e-5 * mol.integrator.timestep
+            assert (step.time - lasttime - frame_interval) < epsilon_t
         lasttime = step.time
 
     # If frame_interval doesn't divide duration exactly, it's permitted to go beyond duration
-    assert duration <= mol.time - traj.time[0] < duration + frame_interval
+    assert duration <= mol.time - traj.time[0] + epsilon_t
+    assert mol.time - traj.time[0] < duration + frame_interval + epsilon_t
 
     if mol.constraints or mol.energy_model.params.get('constrain_hbonds', False):
-        assert_constraints_satisfied(traj, p0, mol)
+        assert_constraints_satisfied(traj, p0, mol, rigorous_constraints=True)
 
 
-def assert_constraints_satisfied(traj, p0, mol):
+def assert_constraints_satisfied(traj, p0, mol, rigorous_constraints=False):
     """ Checks that constraints were satisfied during molecular motion
     """
+    # TODO: check water rigidity, if called for
+    if mol.integrator is not None and mol.integrator.params.get('constrain_hbonds', False):
+        check_hbonds = True
+        assert mdt.geom.HBondsConstraint(mol).satisfied()  # not sure what tolerance should be here
+    else:
+        check_hbonds = False
+
     for constraint in mol.constraints:
         assert constraint.satisfied()
 
-    # TODO: check whole trajectory (not necessarily conserved during minimization however)
-    # TODO: check hbonds, if energy model calls for it
+    if rigorous_constraints:
+        testmol = mol.copy()
+        if check_hbonds:
+            hbonds = mdt.HBondsConstraint(testmol)
+        for frame in traj.frames[1:]:  # ok if starting positions don't obey constraints
+            testmol.positions = frame.positions
+            for constraint in testmol.constraints:
+                assert constraint.satisfied()
+            if check_hbonds:
+                assert hbonds.satisfied()
 
 
 def assert_almost_equal(actual, desired, **kwargs):

moldesign/_tests/test_alignments.py

@@ -8,8 +8,7 @@
 from moldesign.mathutils import normalized
 from moldesign import units as u
 
-from .molecule_fixtures import (benzene, h2, ligand3aid,
-                                pdb3aid, small_molecule, pdb1yu8, ligand_residue_3aid)
+from .molecule_fixtures import *
 from .helpers import assert_almost_equal
 
 

moldesign/models/openmm.py

@@ -66,6 +66,7 @@ def _reset(self):
         self.mm_integrator = None
         self._prepped_integrator = 'uninitialized'
         self._constraints_set = False
+        self._required_tolerance = None
 
     def get_openmm_simulation(self):
         if opm.force_remote:
@@ -133,6 +134,8 @@ def prep(self):
 
         platform, platform_properties = self._get_platform()
 
+        if self._required_tolerance:
+            self.mm_integrator.setConstraintTolerance(float(self._required_tolerance))
         self.sim = app.Simulation(self.mol.ff.parmed_obj.topology,
                                   self.mm_system,
                                   self.mm_integrator,
@@ -221,7 +224,7 @@ def _set_constraints(self):
 
         # openmm uses a global tolerance, calculated as ``constraint_violation/constraint_dist``
         # (since only distance constraints are supported). Here we calculate the necessary value
-        required_tolerance = 1e-5
+        required_tolerance = None
 
         # Constrain atom positions
         for constraint in self.mol.constraints:
@@ -237,6 +240,9 @@ def _set_constraints(self):
                 system.addConstraint(constraint.a1.index,
                                      constraint.a2.index,
                                      opm.pint2simtk(constraint.value))
+
+                if required_tolerance is None:
+                    required_tolerance = 1e-5
                 required_tolerance = min(required_tolerance, constraint.tolerance/constraint.value)
 
             elif constraint.desc == 'hbonds':
@@ -265,10 +271,8 @@ def _set_constraints(self):
                 else:
                     ic += 1
 
-        if self.mm_integrator is not None:
-            self.mm_integrator.setConstraintTolerance(float(required_tolerance))
-
         self._constraints_set = True
+        self._required_tolerance = required_tolerance
 
     @staticmethod
     def _make_dummy_integrator():