Merge pull request #156 from rmcgibbo/dmstests

dmstests

wrappers/python/simtk/openmm/app/__init__.py

@@ -20,11 +20,11 @@ from amberinpcrdfile import AmberInpcrdFile
 from dcdfile import DCDFile
 from gromacsgrofile import GromacsGroFile
 from gromacstopfile import GromacsTopFile
-from desmonddmsfile import DesmondDMSFile
 from dcdreporter import DCDReporter
 from modeller import Modeller
 from statedatareporter import StateDataReporter
 from element import Element
+from desmonddmsfile import DesmondDMSFile
 
 # Enumerated values
 

wrappers/python/simtk/openmm/app/desmonddmsfile.py

@@ -111,9 +111,9 @@ class DesmondDMSFile(object):
 
         boxVectors = []
         for x, y, z in self._conn.execute('SELECT x, y, z FROM global_cell'):
-            boxVectors.append(mm.Vec3(x, y, z)*angstrom)
+            boxVectors.append(mm.Vec3(x, y, z))
         unitCellDimensions = [boxVectors[0][0], boxVectors[1][1], boxVectors[2][2]]
-        top.setUnitCellDimensions(unitCellDimensions)
+        top.setUnitCellDimensions(unitCellDimensions*angstrom)
 
         atoms = {}
         lastChain = None
@@ -122,10 +122,12 @@ class DesmondDMSFile(object):
         q = '''SELECT id, name, anum, resname, resid, chain, x, y, z
         FROM particle'''
         for (atomId, atomName, atomNumber, resName, resId, chain, x, y, z) in self._conn.execute(q):
+            newChain = False
             if chain != lastChain:
                 lastChain = chain
                 c = top.addChain()
-            if resId != lastResId:
+                newChain = True
+            if resId != lastResId or newChain:
                 lastResId = resId
                 if resName in PDBFile._residueNameReplacements:
                     resName = PDBFile._residueNameReplacements[resName]
@@ -144,11 +146,12 @@ class DesmondDMSFile(object):
                 atomName = atomReplacements[atomName]
 
             atoms[atomId] = top.addAtom(atomName, elem, r)
-            positions.append(mm.Vec3(x, y, z)*angstrom)
+            positions.append(mm.Vec3(x, y, z))
 
         for p0, p1 in self._conn.execute('SELECT p0, p1 FROM bond'):
             top.addBond(atoms[p0], atoms[p1])
 
+        positions = positions*angstrom
         return top, positions
 
     def createSystem(self, nonbondedMethod=ff.NoCutoff, nonbondedCutoff=1.0*nanometer,

wrappers/python/tests/TestDesmondDMSFile.py

+import os
+import unittest
+from simtk.openmm.app import *
+from simtk.openmm import *
+from simtk.unit import *
+import simtk.openmm.app.element as elem
+
+class TestDesmondDMSFile(unittest.TestCase):
+    def setUp(self):
+        """Set up the tests by loading the input files."""
+
+        # alanine dipeptide with explicit water
+        path = os.path.join(os.path.dirname(__file__), 'systems/alanine-dipeptide-explicit-amber99SBILDN-tip3p.dms')
+        self.dms = DesmondDMSFile(path)
+    
+    def test_NonbondedMethod(self):
+        """Test all five options for the nonbondedMethod parameter."""
+
+        methodMap = {NoCutoff:NonbondedForce.NoCutoff, 
+                     CutoffNonPeriodic:NonbondedForce.CutoffNonPeriodic, 
+                     CutoffPeriodic:NonbondedForce.CutoffPeriodic, 
+                     Ewald:NonbondedForce.Ewald, PME: NonbondedForce.PME}
+        for method in methodMap:
+            system = self.dms.createSystem(nonbondedMethod=method)
+            forces = system.getForces()
+            self.assertTrue(any(isinstance(f, NonbondedForce) and 
+                                f.getNonbondedMethod()==methodMap[method] 
+                                for f in forces))
+                                
+    def test_Cutoff(self):
+        """Test to make sure the nonbondedCutoff parameter is passed correctly."""
+
+        for method in [CutoffNonPeriodic, CutoffPeriodic, Ewald, PME]:
+            system = self.dms.createSystem(nonbondedMethod=method, 
+                                            nonbondedCutoff=2*nanometer)
+            cutoff_distance = 0.0*nanometer
+            cutoff_check = 2.0*nanometer
+            for force in system.getForces():
+                if isinstance(force, NonbondedForce):
+                    cutoff_distance = force.getCutoffDistance()
+            self.assertEqual(cutoff_distance, cutoff_check)
+
+    def test_EwaldErrorTolerance(self):
+        """Test to make sure the ewaldErrorTolerance parameter is passed correctly."""
+
+        for method in [Ewald, PME]:
+            system = self.dms.createSystem(nonbondedMethod=method,
+                                            ewaldErrorTolerance=1e-6)
+            tolerance = 0
+            tolerance_check = 1e-6
+            for force in system.getForces():
+                if isinstance(force, NonbondedForce):
+                    tolerance = force.getEwaldErrorTolerance()
+            self.assertEqual(tolerance, tolerance_check)
+
+    def test_RemoveCMMotion(self):
+        """Test both options (True and False) for the removeCMMotion parameter."""
+
+        for b in [True, False]:
+            system = self.dms.createSystem(removeCMMotion=b)
+            self.assertEqual(any(isinstance(f, CMMotionRemover) for f in system.getForces()), b)
+
+    def test_HydrogenMass(self):
+        """Test that altering the mass of hydrogens works correctly."""
+        
+        topology = self.dms.topology
+        hydrogenMass = 4*amu
+        system1 = self.dms.createSystem()
+        system2 = self.dms.createSystem(hydrogenMass=hydrogenMass)
+        for atom in topology.atoms():
+            if atom.element == elem.hydrogen:
+                self.assertNotEqual(hydrogenMass, system1.getParticleMass(atom.index))
+                self.assertEqual(hydrogenMass, system2.getParticleMass(atom.index))
+        totalMass1 = sum([system1.getParticleMass(i) for i in range(system1.getNumParticles())]).value_in_unit(amu)
+        totalMass2 = sum([system2.getParticleMass(i) for i in range(system2.getNumParticles())]).value_in_unit(amu)
+        self.assertAlmostEqual(totalMass1, totalMass2)