Merge pull request #2625 from peastman/drudemass

Added drudeMass argument for setting mass of Drude particles

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

@@ -804,7 +804,8 @@ class CharmmPsfFile(object):
                      ewaldErrorTolerance=0.0005,
                      flexibleConstraints=True,
                      verbose=False,
-                     gbsaModel=None):
+                     gbsaModel=None,
+                     drudeMass=0.4*u.amu):
         """Construct an OpenMM System representing the topology described by the
         prmtop file. You MUST have loaded a parameter set into this PSF before
         calling createSystem. If not, AttributeError will be raised. ValueError
@@ -862,6 +863,9 @@ class CharmmPsfFile(object):
         gbsaModel : str=None
             Can be ACE (to use the ACE solvation model) or None. Other values
             raise a ValueError
+        drudeMass : mass=0.4*amu
+            The mass to use for Drude particles.  Any mass added to a Drude particle is
+            subtracted from its parent atom to keep their total mass the same.
         """
         # Load the parameter set
         self.loadParameters(params)
@@ -1460,6 +1464,17 @@ class CharmmPsfFile(object):
                     drude2 = ia2 + 1
                     drudeforce.addScreenedPair(particleMap[drude1], particleMap[drude2], thole1+thole2)
 
+            # Set the masses of Drude particles.
+            if not u.is_quantity(drudeMass):
+                drudeMass *= u.dalton
+            for i in range(drudeforce.getNumParticles()):
+                params = drudeforce.getParticleParameters(i)
+                particle = params[0]
+                parent = params[1]
+                transferMass = drudeMass-system.getParticleMass(particle)
+                system.setParticleMass(particle, drudeMass)
+                system.setParticleMass(parent, system.getParticleMass(parent)-transferMass)
+
         # If we needed a CustomNonbondedForce, map all of the exceptions from
         # the NonbondedForce to the CustomNonbondedForce
         if has_nbfix_terms:

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

@@ -1082,7 +1082,7 @@ class ForceField(object):
 
     def createSystem(self, topology, nonbondedMethod=NoCutoff, nonbondedCutoff=1.0*unit.nanometer,
                      constraints=None, rigidWater=None, removeCMMotion=True, hydrogenMass=None, residueTemplates=dict(),
-                     ignoreExternalBonds=False, switchDistance=None, flexibleConstraints=False, **args):
+                     ignoreExternalBonds=False, switchDistance=None, flexibleConstraints=False, drudeMass=0.4*unit.amu, **args):
         """Construct an OpenMM System representing a Topology with this force field.
 
         Parameters
@@ -1124,6 +1124,9 @@ class ForceField(object):
             Lennard-Jones interactions. If this is None, no switching function will be used.
         flexibleConstraints : boolean=False
             If True, parameters for constrained degrees of freedom will be added to the System
+        drudeMass : mass=0.4*amu
+            The mass to use for Drude particles.  Any mass added to a Drude particle is
+            subtracted from its parent 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
@@ -1136,6 +1139,7 @@ class ForceField(object):
         """
         args['switchDistance'] = switchDistance
         args['flexibleConstraints'] = flexibleConstraints
+        args['drudeMass'] = drudeMass
         data = ForceField._SystemData(topology)
         rigidResidue = [False]*topology.getNumResidues()
 
@@ -5803,6 +5807,19 @@ class DrudeGenerator(object):
                     thole2 = self.typeMap[type2][8]
                     drude.addScreenedPair(drude1, drude2, thole1+thole2)
 
+        # Set the masses of Drude particles.
+
+        drudeMass = args['drudeMass']
+        if not unit.is_quantity(drudeMass):
+            drudeMass *= unit.dalton
+        for i in range(drude.getNumParticles()):
+            params = drude.getParticleParameters(i)
+            particle = params[0]
+            parent = params[1]
+            transferMass = drudeMass-sys.getParticleMass(particle)
+            sys.setParticleMass(particle, drudeMass)
+            sys.setParticleMass(parent, sys.getParticleMass(parent)-transferMass)
+
 parsers["DrudeForce"] = DrudeGenerator.parseElement
 
 #=============================================================================================

wrappers/python/tests/TestCharmmFiles.py

@@ -96,6 +96,33 @@ class TestCharmmFiles(unittest.TestCase):
         totalMass2 = sum([system2.getParticleMass(i) for i in range(system2.getNumParticles())]).value_in_unit(amu)
         self.assertAlmostEqual(totalMass1, totalMass2)
 
+    def test_DrudeMass(self):
+        """Test that setting the mass of Drude particles works correctly."""
+
+        psf = CharmmPsfFile('systems/ala3_solv_drude.psf')
+        crd = CharmmCrdFile('systems/ala3_solv_drude.crd')
+        params = CharmmParameterSet('systems/toppar_drude_master_protein_2013e.str')
+        system = psf.createSystem(params, drudeMass=0)
+        trueMass = [system.getParticleMass(i) for i in range(system.getNumParticles())]
+        drudeMass = 0.3*amu
+        system = psf.createSystem(params, drudeMass=drudeMass)
+        adjustedMass = [system.getParticleMass(i) for i in range(system.getNumParticles())]
+        drudeForce = [f for f in system.getForces() if isinstance(f, DrudeForce)][0]
+        drudeParticles = set()
+        parentParticles = set()
+        for i in range(drudeForce.getNumParticles()):
+            params = drudeForce.getParticleParameters(i)
+            drudeParticles.add(params[0])
+            parentParticles.add(params[1])
+        for i in range(system.getNumParticles()):
+            if i in drudeParticles:
+                self.assertEqual(0*amu, trueMass[i])
+                self.assertEqual(drudeMass, adjustedMass[i])
+            elif i in parentParticles:
+                self.assertEqual(trueMass[i]-drudeMass, adjustedMass[i])
+            else:
+                self.assertEqual(trueMass[i], adjustedMass[i])
+
     def test_NBFIX(self):
         """Tests CHARMM systems with NBFIX Lennard-Jones modifications"""
         warnings.filterwarnings('ignore', category=CharmmPSFWarning)

wrappers/python/tests/TestForceField.py

@@ -251,6 +251,34 @@ class TestForceField(unittest.TestCase):
         totalMass2 = sum([system2.getParticleMass(i) for i in range(system2.getNumParticles())]).value_in_unit(amu)
         self.assertAlmostEqual(totalMass1, totalMass2)
 
+    def test_DrudeMass(self):
+        """Test that setting the mass of Drude particles works correctly."""
+
+        forcefield = ForceField('charmm_polar_2013.xml')
+        pdb = PDBFile('systems/ala_ala_ala.pdb')
+        modeller = Modeller(pdb.topology, pdb.positions)
+        modeller.addExtraParticles(forcefield)
+        system = forcefield.createSystem(modeller.topology, drudeMass=0)
+        trueMass = [system.getParticleMass(i) for i in range(system.getNumParticles())]
+        drudeMass = 0.3*amu
+        system = forcefield.createSystem(modeller.topology, drudeMass=drudeMass)
+        adjustedMass = [system.getParticleMass(i) for i in range(system.getNumParticles())]
+        drudeForce = [f for f in system.getForces() if isinstance(f, DrudeForce)][0]
+        drudeParticles = set()
+        parentParticles = set()
+        for i in range(drudeForce.getNumParticles()):
+            params = drudeForce.getParticleParameters(i)
+            drudeParticles.add(params[0])
+            parentParticles.add(params[1])
+        for i in range(system.getNumParticles()):
+            if i in drudeParticles:
+                self.assertEqual(0*amu, trueMass[i])
+                self.assertEqual(drudeMass, adjustedMass[i])
+            elif i in parentParticles:
+                self.assertEqual(trueMass[i]-drudeMass, adjustedMass[i])
+            else:
+                self.assertEqual(trueMass[i], adjustedMass[i])
+
     def test_Forces(self):
         """Compute forces and compare them to ones generated with a previous version of OpenMM to ensure they haven't changed."""