Do not alter mass of hydrogens in rigid water (#2975)

* Do not alter mass of hydrogens in rigid water

* Do not alter mass of hydrogens in rigid water

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

@@ -200,7 +200,8 @@ class AmberPrmtopFile(object):
         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.
+            total mass the same.  If rigidWater is used to make water molecules
+            rigid, then water hydrogens are not altered.
         ewaldErrorTolerance : float=0.0005
             The error tolerance to use if nonbondedMethod is Ewald, PME, or LJPME.
         switchDistance : float=0*nanometers
@@ -284,6 +285,8 @@ class AmberPrmtopFile(object):
             for atom1, atom2 in self.topology.bonds():
                 if atom1.element == elem.hydrogen:
                     (atom1, atom2) = (atom2, atom1)
+                if rigidWater and atom2.residue.name == 'HOH':
+                    continue
                 if atom2.element == elem.hydrogen and atom1.element not in (elem.hydrogen, None):
                     transferMass = hydrogenMass-sys.getParticleMass(atom2.index)
                     sys.setParticleMass(atom2.index, hydrogenMass)

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

@@ -855,7 +855,8 @@ class CharmmPsfFile(object):
         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.
+            total mass the same.  If rigidWater is used to make water molecules
+            rigid, then water hydrogens are not altered.
         ewaldErrorTolerance : float=0.0005
             The error tolerance to use if the nonbonded method is Ewald, PME, or LJPME.
         flexibleConstraints : bool=True
@@ -1581,6 +1582,8 @@ class CharmmPsfFile(object):
                 if (bond.atom1.type.atomic_number != 1 and
                     bond.atom2.type.atomic_number != 1):
                     continue
+                if rigidWater and _is_bond_in_water(bond):
+                    continue
                 atom1, atom2 = bond.atom1, bond.atom2
                 if atom1.type.atomic_number == 1:
                     atom1, atom2 = atom2, atom1 # now atom2 is hydrogen for sure

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

@@ -1135,7 +1135,8 @@ class ForceField(object):
         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.
+            their total mass the same.  If rigidWater is used to make water molecules
+            rigid, then water hydrogens are not altered.
         residueTemplates : dict=dict()
             Key: Topology Residue object
             Value: string, name of _TemplateData residue template object to use for (Key) residue.
@@ -1211,7 +1212,7 @@ class ForceField(object):
             for atom1, atom2 in topology.bonds():
                 if atom1.element is elem.hydrogen:
                     (atom1, atom2) = (atom2, atom1)
-                if atom2.element is elem.hydrogen and atom1.element not in (elem.hydrogen, None):
+                if atom2.element is elem.hydrogen and atom1.element not in (elem.hydrogen, None) and not rigidResidue[atom2.residue.index]:
                     transferMass = hydrogenMass-sys.getParticleMass(atom2.index)
                     sys.setParticleMass(atom2.index, hydrogenMass)
                     sys.setParticleMass(atom1.index, sys.getParticleMass(atom1.index)-transferMass)

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

@@ -627,7 +627,8 @@ class GromacsTopFile(object):
         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.
+            total mass the same.  If rigidWater is used to make water molecules
+            rigid, then water hydrogens are not altered.
         switchDistance : float=None
             The distance at which the potential energy switching function is turned on for
             Lennard-Jones interactions. If this is None, no switching function will be used.
@@ -1202,6 +1203,8 @@ class GromacsTopFile(object):
             for atom1, atom2 in self.topology.bonds():
                 if atom1.element == elem.hydrogen:
                     (atom1, atom2) = (atom2, atom1)
+                if rigidWater and atom2.residue.name == 'HOH':
+                    continue
                 if atom2.element == elem.hydrogen and atom1.element not in (elem.hydrogen, None):
                     transferMass = hydrogenMass-sys.getParticleMass(atom2.index)
                     sys.setParticleMass(atom2.index, hydrogenMass)

wrappers/python/tests/TestAmberPrmtopFile.py

@@ -145,6 +145,9 @@ class TestAmberPrmtopFile(unittest.TestCase):
         for atom in topology.atoms():
             if atom.element == elem.hydrogen:
                 self.assertNotEqual(hydrogenMass, system1.getParticleMass(atom.index))
+                if atom.residue.name == 'HOH':
+                    self.assertEqual(system1.getParticleMass(atom.index), system2.getParticleMass(atom.index))
+                else:
                     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)

wrappers/python/tests/TestCharmmFiles.py

@@ -91,6 +91,9 @@ class TestCharmmFiles(unittest.TestCase):
         for atom in topology.atoms():
             if atom.element == elem.hydrogen:
                 self.assertNotEqual(hydrogenMass, system1.getParticleMass(atom.index))
+                if atom.residue.name == 'HOH':
+                    self.assertEqual(system1.getParticleMass(atom.index), system2.getParticleMass(atom.index))
+                else:
                     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)

wrappers/python/tests/TestForceField.py

@@ -246,6 +246,9 @@ class TestForceField(unittest.TestCase):
         for atom in topology.atoms():
             if atom.element == elem.hydrogen:
                 self.assertNotEqual(hydrogenMass, system1.getParticleMass(atom.index))
+                if atom.residue.name == 'HOH':
+                    self.assertEqual(system1.getParticleMass(atom.index), system2.getParticleMass(atom.index))
+                else:
                     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)

wrappers/python/tests/TestGromacsTopFile.py

@@ -176,6 +176,9 @@ class TestGromacsTopFile(unittest.TestCase):
         for atom in topology.atoms():
             if atom.element == elem.hydrogen:
                 self.assertNotEqual(hydrogenMass, system1.getParticleMass(atom.index))
+                if atom.residue.name == 'HOH':
+                    self.assertEqual(system1.getParticleMass(atom.index), system2.getParticleMass(atom.index))
+                else:
                     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)