Merge pull request #141 from peastman/master

Added option to increase hydrogen mass

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

@@ -127,7 +127,7 @@ class AmberPrmtopFile(object):
 
     def createSystem(self, nonbondedMethod=ff.NoCutoff, nonbondedCutoff=1.0*unit.nanometer,
                      constraints=None, rigidWater=True, implicitSolvent=None, soluteDielectric=1.0, solventDielectric=78.5, removeCMMotion=True,
-                     ewaldErrorTolerance=0.0005):
+                     hydrogenMass=None, ewaldErrorTolerance=0.0005):
         """Construct an OpenMM System representing the topology described by this prmtop file.
 
         Parameters:
@@ -141,6 +141,8 @@ class AmberPrmtopFile(object):
          - soluteDielectric (float=1.0) The solute dielectric constant to use in the implicit solvent model.
          - solventDielectric (float=78.5) The solvent dielectric constant to use in the implicit solvent model.
          - removeCMMotion (boolean=True) If true, a CMMotionRemover will be added to the System
+         - hydrogenMass (mass=None) The mass to use for hydrogen atoms bound to heavy atoms.  Any mass added to a hydrogen is
+           subtracted from the heavy atom to keep their total mass the same.
          - ewaldErrorTolerance (float=0.0005) The error tolerance to use if nonbondedMethod is Ewald or PME.
         Returns: the newly created System
         """
@@ -178,6 +180,14 @@ class AmberPrmtopFile(object):
         sys = amber_file_parser.readAmberSystem(prmtop_loader=self._prmtop, shake=constraintString, nonbondedCutoff=nonbondedCutoff,
                                                  nonbondedMethod=methodMap[nonbondedMethod], flexibleConstraints=False, gbmodel=implicitString,
                                                  soluteDielectric=soluteDielectric, solventDielectric=solventDielectric, rigidWater=rigidWater)
+        if hydrogenMass is not None:
+            for atom1, atom2 in self.topology.bonds():
+                if atom1.element == elem.hydrogen:
+                    (atom1, atom2) = (atom2, atom1)
+                if atom2.element == elem.hydrogen and atom1.element not in (elem.hydrogen, None):
+                    transferMass = hydrogenMass-sys.getParticleMass(atom2.index)
+                    sys.setParticleMass(atom2.index, hydrogenMass)
+                    sys.setParticleMass(atom1.index, sys.getParticleMass(atom1.index)-transferMass)
         for force in sys.getForces():
             if isinstance(force, mm.NonbondedForce):
                 force.setEwaldErrorTolerance(ewaldErrorTolerance)

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

@@ -272,7 +272,7 @@ class ForceField(object):
             self.atoms = [x-self.index for x in self.atoms]
 
     def createSystem(self, topology, nonbondedMethod=NoCutoff, nonbondedCutoff=1.0*unit.nanometer,
-                     constraints=None, rigidWater=True, removeCMMotion=True, **args):
+                     constraints=None, rigidWater=True, removeCMMotion=True, hydrogenMass=None, **args):
         """Construct an OpenMM System representing a Topology with this force field.
 
         Parameters:
@@ -284,6 +284,8 @@ class ForceField(object):
            Allowed values are None, HBonds, AllBonds, or HAngles.
          - rigidWater (boolean=True) If true, water molecules will be fully rigid regardless of the value passed for the constraints argument
          - removeCMMotion (boolean=True) If true, a CMMotionRemover will be added to the System
+         - hydrogenMass (mass=None) The mass to use for hydrogen atoms bound to heavy atoms.  Any mass added to a hydrogen is
+           subtracted from the heavy atom to keep their total mass the same.
          - args Arbitrary additional keyword arguments may also be specified.  This allows extra parameters to be specified that are specific to
            particular force fields.
         Returns: the newly created System
@@ -336,6 +338,17 @@ class ForceField(object):
         for atom in topology.atoms():
             sys.addParticle(self._atomTypes[data.atomType[atom]][1])
         
+        # Adjust masses.
+        
+        if hydrogenMass is not None:
+            for atom1, atom2 in topology.bonds():
+                if atom1.element == elem.hydrogen:
+                    (atom1, atom2) = (atom2, atom1)
+                if atom2.element == elem.hydrogen and atom1.element not in (elem.hydrogen, None):
+                    transferMass = hydrogenMass-sys.getParticleMass(atom2.index)
+                    sys.setParticleMass(atom2.index, hydrogenMass)
+                    sys.setParticleMass(atom1.index, sys.getParticleMass(atom1.index)-transferMass)
+
         # Set periodic boundary conditions.
 
         boxSize = topology.getUnitCellDimensions()

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

@@ -427,7 +427,7 @@ class GromacsTopFile(object):
                     top.addBond(atoms[int(fields[0])-1], atoms[int(fields[1])-1])
 
     def createSystem(self, nonbondedMethod=ff.NoCutoff, nonbondedCutoff=1.0*unit.nanometer,
-                     constraints=None, rigidWater=True, implicitSolvent=None, soluteDielectric=1.0, solventDielectric=78.5, ewaldErrorTolerance=0.0005, removeCMMotion=True):
+                     constraints=None, rigidWater=True, implicitSolvent=None, soluteDielectric=1.0, solventDielectric=78.5, ewaldErrorTolerance=0.0005, removeCMMotion=True, hydrogenMass=None):
         """Construct an OpenMM System representing the topology described by this prmtop file.
 
         Parameters:
@@ -442,6 +442,8 @@ class GromacsTopFile(object):
          - solventDielectric (float=78.5) The solvent dielectric constant to use in the implicit solvent model.
          - ewaldErrorTolerance (float=0.0005) The error tolerance to use if nonbondedMethod is Ewald or PME.
          - removeCMMotion (boolean=True) If true, a CMMotionRemover will be added to the System
+         - hydrogenMass (mass=None) The mass to use for hydrogen atoms bound to heavy atoms.  Any mass added to a hydrogen is
+           subtracted from the heavy atom to keep their total mass the same.
         Returns: the newly created System
         """
         # Create the System.
@@ -734,6 +736,17 @@ class GromacsTopFile(object):
         nb.setCutoffDistance(nonbondedCutoff)
         nb.setEwaldErrorTolerance(ewaldErrorTolerance)
         
+        # Adjust masses.
+        
+        if hydrogenMass is not None:
+            for atom1, atom2 in self.topology.bonds():
+                if atom1.element == elem.hydrogen:
+                    (atom1, atom2) = (atom2, atom1)
+                if atom2.element == elem.hydrogen and atom1.element not in (elem.hydrogen, None):
+                    transferMass = hydrogenMass-sys.getParticleMass(atom2.index)
+                    sys.setParticleMass(atom2.index, hydrogenMass)
+                    sys.setParticleMass(atom1.index, sys.getParticleMass(atom1.index)-transferMass)
+
         # Add a CMMotionRemover.
 
         if removeCMMotion:

wrappers/python/tests/TestAmberPrmtopFile.py

@@ -3,6 +3,7 @@ from validateConstraints import *
 from simtk.openmm.app import *
 from simtk.openmm import *
 from simtk.unit import *
+import simtk.openmm.app.element as elem
 
 class TestAmberPrmtopFile(unittest.TestCase):
 
@@ -132,5 +133,20 @@ class TestAmberPrmtopFile(unittest.TestCase):
                     self.assertTrue(found_matching_solvent_dielectric and 
                                     found_matching_solute_dielectric)
 
+    def test_HydrogenMass(self):
+        """Test that altering the mass of hydrogens works correctly."""
+        
+        topology = self.prmtop1.topology
+        hydrogenMass = 4*amu
+        system1 = self.prmtop1.createSystem()
+        system2 = self.prmtop1.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)
+
 if __name__ == '__main__':
     unittest.main()

wrappers/python/tests/TestForceField.py

@@ -3,6 +3,7 @@ from validateConstraints import *
 from simtk.openmm.app import *
 from simtk.openmm import *
 from simtk.unit import *
+import simtk.openmm.app.element as elem
 
 class TestForceField(unittest.TestCase):
     """Test the ForceField.createSystem() method."""
@@ -104,6 +105,20 @@ class TestForceField(unittest.TestCase):
         self.assertTrue(found_matching_solvent_dielectric and 
                         found_matching_solute_dielectric)
 
+    def test_HydrogenMass(self):
+        """Test that altering the mass of hydrogens works correctly."""
+        
+        topology = self.pdb1.topology
+        hydrogenMass = 4*amu
+        system1 = self.forcefield1.createSystem(topology)
+        system2 = self.forcefield1.createSystem(topology, 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)
 
 
 if __name__ == '__main__':

wrappers/python/tests/TestGromacsTopFile.py

@@ -3,6 +3,7 @@ from validateConstraints import *
 from simtk.openmm.app import *
 from simtk.openmm import *
 from simtk.unit import *
+import simtk.openmm.app.element as elem
 
 class TestGromacsTopFile(unittest.TestCase):
 
@@ -103,6 +104,21 @@ class TestGromacsTopFile(unittest.TestCase):
         self.assertTrue(found_matching_solvent_dielectric and 
                         found_matching_solute_dielectric)
 
+    def test_HydrogenMass(self):
+        """Test that altering the mass of hydrogens works correctly."""
+        
+        topology = self.top1.topology
+        hydrogenMass = 4*amu
+        system1 = self.top1.createSystem()
+        system2 = self.top1.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)
+
 if __name__ == '__main__':
     unittest.main()