Make CharmmPsfFile ignore H-H water bond that is required for rigid water with...

Make CharmmPsfFile ignore H-H water bond that is required for rigid water with original CHARMM implementation of constrains.

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

@@ -272,6 +272,9 @@ class CharmmPsfFile(object):
                 drudepair_list.append([min(id1,id2), max(id1,id2)])
             elif (atom_list[id1].name[0:2]=='LP' or atom_list[id2].name[0:2]=='LP' or atom_list[id1].name=='OM' or atom_list[id2].name=='OM'):
                 pass
+            # Ignore H-H bond in water if present
+            elif atom_list[id1].name[0]=='H' and atom_list[id2].name[0]=='H' and (atom_list[id1].residue.resname in WATNAMES):
+                pass
             else:
                 bond_list.append(Bond(atom_list[id1], atom_list[id2]))
         bond_list.changed = False

wrappers/python/tests/TestCharmmFiles.py

@@ -395,6 +395,93 @@ class TestCharmmFiles(unittest.TestCase):
         self.assertEqual(sorted(system.getConstraintParameters(i)[:2] for i in range(system.getNumConstraints())),
                          sorted(hBonds_water + hAngles_water + allBonds_methanol + hAngles_methanol))
 
+    def test_Constraints_charmm(self):
+        """Tests that CHARMM and OpenMM implementation of CHARMM force field produce the same constraints and energy"""
+        warnings.filterwarnings('ignore', category=CharmmPSFWarning)
+        psf = CharmmPsfFile('systems/ala3_solv.psf',
+                            unitCellDimensions=Vec3(32.7119500, 32.9959600, 33.0071500) * angstroms)
+        crd = CharmmCrdFile('systems/ala3_solv.crd')
+        params = CharmmParameterSet('systems/par_all36_prot.prm',
+                                    'systems/toppar_water_ions.str')
+        plat = Platform.getPlatformByName('Reference')
+        system_charmm = psf.createSystem(params, nonbondedMethod=PME,
+                                  nonbondedCutoff=8 * angstroms)
+        topology = psf.topology
+        forcefield = ForceField('charmm36.xml', 'charmm36/water.xml')
+        system_openmm = forcefield.createSystem(topology, nonbondedMethod=PME,
+                                  nonbondedCutoff=8 * angstroms)
+        # Test different combinations of constrains/rigidWater parameters
+        system_charmm = psf.createSystem(params, constraints=None, rigidWater=False, nonbondedMethod=PME,
+                                         nonbondedCutoff=8 * angstroms)
+        system_openmm = forcefield.createSystem(topology, constraints=None, rigidWater=False, nonbondedMethod=PME,
+                                                nonbondedCutoff=8 * angstroms)
+        self.assertEqual(system_charmm.getNumConstraints(), 0)
+        self.assertEqual(system_openmm.getNumConstraints(), 0)
+        con_charmm = Context(system_charmm, VerletIntegrator(2 * femtoseconds), plat)
+        con_charmm.setPositions(crd.positions)
+        con_openmm = Context(system_openmm, VerletIntegrator(2 * femtoseconds), plat)
+        con_openmm.setPositions(crd.positions)
+        state_charmm = con_charmm.getState(getEnergy=True, enforcePeriodicBox=True)
+        ene_charmm = state_charmm.getPotentialEnergy().value_in_unit(kilocalories_per_mole)
+        state_openmm = con_openmm.getState(getEnergy=True, enforcePeriodicBox=True)
+        ene_openmm = state_openmm.getPotentialEnergy().value_in_unit(kilocalories_per_mole)
+        self.assertAlmostEqual(ene_charmm, ene_openmm, delta=0.05)
+
+        system_charmm = psf.createSystem(params, constraints=None, rigidWater=True, nonbondedMethod=PME,
+                                         nonbondedCutoff=8 * angstroms)
+        system_openmm = forcefield.createSystem(topology, constraints=None, rigidWater=True, nonbondedMethod=PME,
+                                                nonbondedCutoff=8 * angstroms)
+        self.assertEqual(
+            sorted(system_charmm.getConstraintParameters(i) for i in range(system_charmm.getNumConstraints())),
+            sorted(system_charmm.getConstraintParameters(j) for j in range(system_charmm.getNumConstraints())))
+
+        system_charmm = psf.createSystem(params, constraints=HBonds, rigidWater=False, nonbondedMethod=PME,
+                                         nonbondedCutoff=8 * angstroms)
+        system_openmm = forcefield.createSystem(topology, constraints=HBonds, rigidWater=False, nonbondedMethod=PME,
+                                                nonbondedCutoff=8 * angstroms)
+        self.assertEqual(
+            sorted(system_charmm.getConstraintParameters(i) for i in range(system_charmm.getNumConstraints())),
+            sorted(system_charmm.getConstraintParameters(j) for j in range(system_charmm.getNumConstraints())))
+
+        system_charmm = psf.createSystem(params, constraints=HBonds, rigidWater=True, nonbondedMethod=PME,
+                                         nonbondedCutoff=8 * angstroms)
+        system_openmm = forcefield.createSystem(topology, constraints=HBonds, rigidWater=True, nonbondedMethod=PME,
+                                                nonbondedCutoff=8 * angstroms)
+        self.assertEqual(
+            sorted(system_charmm.getConstraintParameters(i) for i in range(system_charmm.getNumConstraints())),
+            sorted(system_charmm.getConstraintParameters(j) for j in range(system_charmm.getNumConstraints())))
+
+        system_charmm = psf.createSystem(params, constraints=AllBonds, rigidWater=False, nonbondedMethod=PME,
+                                         nonbondedCutoff=8 * angstroms)
+        system_openmm = forcefield.createSystem(topology, constraints=AllBonds, rigidWater=False, nonbondedMethod=PME,
+                                                nonbondedCutoff=8 * angstroms)
+        self.assertEqual(
+            sorted(system_charmm.getConstraintParameters(i) for i in range(system_charmm.getNumConstraints())),
+            sorted(system_charmm.getConstraintParameters(j) for j in range(system_charmm.getNumConstraints())))
+
+        system_charmm = psf.createSystem(params, constraints=AllBonds, rigidWater=True, nonbondedMethod=PME,
+                                         nonbondedCutoff=8 * angstroms)
+        system_openmm = forcefield.createSystem(topology, constraints=AllBonds, rigidWater=True, nonbondedMethod=PME,
+                                                nonbondedCutoff=8 * angstroms)
+        self.assertEqual(
+            sorted(system_charmm.getConstraintParameters(i) for i in range(system_charmm.getNumConstraints())),
+            sorted(system_charmm.getConstraintParameters(j) for j in range(system_charmm.getNumConstraints())))
+
+        system_charmm = psf.createSystem(params, constraints=HAngles, rigidWater=False, nonbondedMethod=PME,
+                                         nonbondedCutoff=8 * angstroms)
+        system_openmm = forcefield.createSystem(topology, constraints=HAngles, rigidWater=False, nonbondedMethod=PME,
+                                                nonbondedCutoff=8 * angstroms)
+        self.assertEqual(
+            sorted(system_charmm.getConstraintParameters(i) for i in range(system_charmm.getNumConstraints())),
+            sorted(system_charmm.getConstraintParameters(j) for j in range(system_charmm.getNumConstraints())))
+
+        system_charmm = psf.createSystem(params, constraints=HAngles, rigidWater=True, nonbondedMethod=PME,
+                                         nonbondedCutoff=8 * angstroms)
+        system_openmm = forcefield.createSystem(topology, constraints=HAngles, rigidWater=True, nonbondedMethod=PME,
+                                                nonbondedCutoff=8 * angstroms)
+        self.assertEqual(
+            sorted(system_charmm.getConstraintParameters(i) for i in range(system_charmm.getNumConstraints())),
+            sorted(system_charmm.getConstraintParameters(j) for j in range(system_charmm.getNumConstraints())))
 
 if __name__ == '__main__':
     unittest.main()