Merge pull request #1731 from kyleabeauchamp/fixmts

Fix substeps division in MTSIntegrator

wrappers/python/simtk/openmm/mtsintegrator.py

@@ -28,7 +28,7 @@ DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
 OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
 USE OR OTHER DEALINGS IN THE SOFTWARE.
 """
-from __future__ import absolute_import
+from __future__ import absolute_import, division
 __author__ = "Peter Eastman"
 __version__ = "1.0"
 
@@ -90,9 +90,10 @@ class MTSIntegrator(CustomIntegrator):
 
     def _createSubsteps(self, parentSubsteps, groups):
         group, substeps = groups[0]
-        stepsPerParentStep = substeps//parentSubsteps
+        stepsPerParentStep = substeps / parentSubsteps
         if stepsPerParentStep < 1 or stepsPerParentStep != int(stepsPerParentStep):
             raise ValueError("The number for substeps for each group must be a multiple of the number for the previous group")
+        stepsPerParentStep = int(stepsPerParentStep)
         if group < 0 or group > 31:
             raise ValueError("Force group must be between 0 and 31")
         for i in range(stepsPerParentStep):

wrappers/python/tests/TestIntegrators.py

@@ -96,5 +96,29 @@ class TestIntegrators(unittest.TestCase):
             # Take a step
             integrator.step(1)
 
+    def testBadGroups(self):
+        """Test the MTS integrator with bad force group substeps."""
+        # Create a periodic solvated system with PME
+        pdb = PDBFile('systems/alanine-dipeptide-explicit.pdb')
+        ff = ForceField('amber99sbildn.xml', 'tip3p.xml')
+        system = ff.createSystem(pdb.topology, cutoffMethod=PME)
+
+        # Split forces into groups
+        for force in system.getForces():
+            if force.__class__.__name__ == 'NonbondedForce':
+                force.setForceGroup(1)
+                force.setReciprocalSpaceForceGroup(2)
+            else:
+                force.setForceGroup(0)
+
+        with self.assertRaises(ValueError):
+            # Create an integrator
+            integrator = MTSIntegrator(4*femtoseconds, [(2,1), (1,3), (0,8)])
+
+            # Run a few steps of dynamics
+            context = Context(system, integrator)
+            context.setPositions(pdb.positions)
+            integrator.step(10)
+
 if __name__ == '__main__':
     unittest.main()