ForceField supports DrudeForce

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

@@ -158,18 +158,20 @@ class ForceField(object):
                 if typeAttrib in attrib:
                     raise ValueError('Tag specifies both a type and a class for the same atom: '+etree.tostring(node))
                 if attrib[classAttrib] not in self._atomClasses:
-                    return None # Unknown atom class
-                types.append(self._atomClasses[attrib[classAttrib]])
+                    types.append(None) # Unknown atom class
+                else:
+                    types.append(self._atomClasses[attrib[classAttrib]])
             else:
                 if typeAttrib not in attrib or attrib[typeAttrib] not in self._atomTypes:
-                    return None # Unknown atom type
-                types.append([attrib[typeAttrib]])
+                    types.append(None) # Unknown atom type
+                else:
+                    types.append([attrib[typeAttrib]])
         return types
 
     def _parseTorsion(self, node):
         """Parse the node defining a torsion."""
         types = self._findAtomTypes(node, 4)
-        if types is None:
+        if None in types:
             return None
         torsion = PeriodicTorsion(types)
         attrib = node.attrib
@@ -416,6 +418,12 @@ class ForceField(object):
         if removeCMMotion:
             sys.addForce(mm.CMMotionRemover())
         
+        # Let generators do postprocessing
+        
+        for force in self._forces:
+            if 'postprocessSystem' in dir(force):
+                force.postprocessSystem(sys, data, args)
+        
         # Execute scripts found in the XML files.
         
         for script in self._scripts:
@@ -483,9 +491,10 @@ def _findAtomMatches(atoms, template, bondedTo, externalBonds, matches, hasMatch
     if position == len(atoms):
         return True
     elem = atoms[position].element
+    name = atoms[position].name
     for i in range(len(atoms)):
         atom = template.atoms[i]
-        if (atom.element == elem or atom.element is None) and not hasMatch[i] and len(atom.bondedTo) == len(bondedTo[position]) and atom.externalBonds == externalBonds[position]:
+        if (atom.element == elem or (atom.element is None and atom.name == name)) and not hasMatch[i] and len(atom.bondedTo) == len(bondedTo[position]) and atom.externalBonds == externalBonds[position]:
             # See if the bonds for this identification are consistent
             
             allBondsMatch = all((bonded > position or matches[bonded] in atom.bondedTo for bonded in bondedTo[position]))
@@ -521,7 +530,7 @@ class HarmonicBondGenerator:
         ff._forces.append(generator)
         for bond in element.findall('Bond'):
             types = ff._findAtomTypes(bond, 2)
-            if types is not None:
+            if None not in types:
                 generator.types1.append(types[0])
                 generator.types2.append(types[1])
                 generator.length.append(float(bond.attrib['length']))
@@ -569,7 +578,7 @@ class HarmonicAngleGenerator:
         ff._forces.append(generator)
         for angle in element.findall('Angle'):
             types = ff._findAtomTypes(angle, 3)
-            if types is not None:
+            if None not in types:
                 generator.types1.append(types[0])
                 generator.types2.append(types[1])
                 generator.types3.append(types[2])
@@ -748,11 +757,11 @@ class RBTorsionGenerator:
         ff._forces.append(generator)
         for torsion in element.findall('Proper'):
             types = ff._findAtomTypes(torsion, 4)
-            if types is not None:
+            if None not in types:
                 generator.proper.append(RBTorsion(types, [float(torsion.attrib['c'+str(i)]) for i in range(6)]))
         for torsion in element.findall('Improper'):
             types = ff._findAtomTypes(torsion, 4)
-            if types is not None:
+            if None not in types:
                 generator.improper.append(RBTorsion(types, [float(torsion.attrib['c'+str(i)]) for i in range(6)]))
     
     def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
@@ -850,7 +859,7 @@ class CMAPTorsionGenerator:
             generator.maps.append(values)
         for torsion in element.findall('Torsion'):
             types = ff._findAtomTypes(torsion, 5)
-            if types is not None:
+            if None not in types:
                 generator.torsions.append(CMAPTorsion(types, int(torsion.attrib['map'])))
     
     def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
@@ -931,7 +940,7 @@ class NonbondedGenerator:
                 raise ValueError('Found multiple NonbondedForce tags with different 1-4 scales') 
         for atom in element.findall('Atom'):
             types = ff._findAtomTypes(atom, 1)
-            if types is not None:
+            if None not in types:
                 values = (float(atom.attrib['charge']), float(atom.attrib['sigma']), float(atom.attrib['epsilon']))
                 for t in types[0]:
                     generator.typeMap[t] = values
@@ -952,7 +961,14 @@ class NonbondedGenerator:
                 force.addParticle(values[0], values[1], values[2])
             else:
                 raise ValueError('No nonbonded parameters defined for atom type '+t)
-        # Create exceptions based on bonds and virtual sites.
+        force.setNonbondedMethod(methodMap[nonbondedMethod])
+        force.setCutoffDistance(nonbondedCutoff)
+        if 'ewaldErrorTolerance' in args:
+            force.setEwaldErrorTolerance(args['ewaldErrorTolerance'])
+        sys.addForce(force)
+    
+    def postprocessSystem(self, sys, data, args):
+        # Create exceptions based on bonds, virtual sites, and Drude particles.
         bondIndices = []
         for bond in data.bonds:
             bondIndices.append((bond.atom1, bond.atom2))
@@ -961,12 +977,26 @@ class NonbondedGenerator:
                 site = sys.getVirtualSite(i)
                 for j in range(site.getNumParticles()):
                     bondIndices.append((i, site.getParticle(j)))
-        force.createExceptionsFromBonds(bondIndices, self.coulomb14scale, self.lj14scale)
-        force.setNonbondedMethod(methodMap[nonbondedMethod])
-        force.setCutoffDistance(nonbondedCutoff)
-        if 'ewaldErrorTolerance' in args:
-            force.setEwaldErrorTolerance(args['ewaldErrorTolerance'])
-        sys.addForce(force)
+        drude = [f for f in sys.getForces() if isinstance(f, mm.DrudeForce)]
+        if len(drude) > 0:
+            drude = drude[0]
+            # For purposes of creating exceptions, a Drude particle is "bonded" to anything
+            # its parent atom is bonded to.
+            drudeMap = {}
+            for i in range(drude.getNumParticles()):
+                params = drude.getParticleParameters(i)
+                drudeMap[params[1]] = params[0]
+            for atom1, atom2 in bondIndices:
+                drude1 = drudeMap[atom1] if atom1 in drudeMap else None
+                drude2 = drudeMap[atom2] if atom2 in drudeMap else None
+                if drude1 is not None:
+                    bondIndices.append((drude1, atom2))
+                    if drude2 is not None:
+                        bondIndices.append((drude1, drude2))
+                if drude2 is not None:
+                    bondIndices.append((atom1, drude2))
+        nonbonded = [f for f in sys.getForces() if isinstance(f, mm.NonbondedForce)][0]
+        nonbonded.createExceptionsFromBonds(bondIndices, self.coulomb14scale, self.lj14scale)        
 
 parsers["NonbondedForce"] = NonbondedGenerator.parseElement
 
@@ -989,7 +1019,7 @@ class GBSAOBCGenerator:
             generator = existing[0]
         for atom in element.findall('Atom'):
             types = ff._findAtomTypes(atom, 1)
-            if types is not None:
+            if None not in types:
                 values = (float(atom.attrib['charge']), float(atom.attrib['radius']), float(atom.attrib['scale']))
                 for t in types[0]:
                     generator.typeMap[t] = values
@@ -1046,7 +1076,7 @@ class GBVIGenerator:
         ff._forces.append(generator)
         for atom in element.findall('Atom'):
             types = ff._findAtomTypes(atom, 1)
-            if types is not None:
+            if None not in types:
                 values = (float(atom.attrib['charge']), float(atom.attrib['radius']), float(atom.attrib['gamma']))
                 for t in types[0]:
                     generator.typeMap[t] = values
@@ -1129,7 +1159,7 @@ class CustomBondGenerator:
             generator.perBondParams.append(param.attrib['name'])
         for bond in element.findall('Bond'):
             types = ff._findAtomTypes(bond, 2)
-            if types is not None:
+            if None not in types:
                 generator.types1.append(types[0])
                 generator.types2.append(types[1])
                 generator.paramValues.append([float(bond.attrib[param]) for param in generator.perBondParams])
@@ -1177,7 +1207,7 @@ class CustomAngleGenerator:
             generator.perAngleParams.append(param.attrib['name'])
         for angle in element.findall('Angle'):
             types = ff._findAtomTypes(angle, 3)
-            if types is not None:
+            if None not in types:
                 generator.types1.append(types[0])
                 generator.types2.append(types[1])
                 generator.types3.append(types[2])
@@ -1238,11 +1268,11 @@ class CustomTorsionGenerator:
             generator.perTorsionParams.append(param.attrib['name'])
         for torsion in element.findall('Proper'):
             types = ff._findAtomTypes(torsion, 4)
-            if types is not None:
+            if None not in types:
                 generator.proper.append(CustomTorsion(types, [float(torsion.attrib[param]) for param in generator.perTorsionParams]))
         for torsion in element.findall('Improper'):
             types = ff._findAtomTypes(torsion, 4)
-            if types is not None:
+            if None not in types:
                 generator.improper.append(CustomTorsion(types, [float(torsion.attrib[param]) for param in generator.perTorsionParams]))
     
     def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
@@ -1329,7 +1359,7 @@ class CustomGBGenerator:
             generator.perParticleParams.append(param.attrib['name'])
         for atom in element.findall('Atom'):
             types = ff._findAtomTypes(atom, 1)
-            if types is not None:
+            if None not in types:
                 values = [float(atom.attrib[param]) for param in generator.perParticleParams]
                 for t in types[0]:
                     generator.typeMap[t] = values
@@ -1431,7 +1461,7 @@ class AmoebaBondGenerator:
         forceField._forces.append(generator)
         for bond in element.findall('Bond'):
             types = forceField._findAtomTypes(bond, 2)
-            if types is not None:
+            if None not in types:
                 generator.types1.append(types[0])
                 generator.types2.append(types[1])
                 generator.length.append(float(bond.attrib['length']))
@@ -1544,7 +1574,7 @@ class AmoebaAngleGenerator:
         forceField._forces.append(generator)
         for angle in element.findall('Angle'):
             types = forceField._findAtomTypes(angle, 3)
-            if types is not None:
+            if None not in types:
 
                 generator.types1.append(types[0])
                 generator.types2.append(types[1])
@@ -2029,7 +2059,7 @@ class AmoebaTorsionGenerator:
 
         for torsion in element.findall('Torsion'):
             types = forceField._findAtomTypes(torsion, 4)
-            if types is not None:
+            if None not in types:
 
                 generator.types1.append(types[0])
                 generator.types2.append(types[1])
@@ -2140,7 +2170,7 @@ class AmoebaPiTorsionGenerator:
 
         for piTorsion in element.findall('PiTorsion'):
             types = forceField._findAtomTypes(piTorsion, 2)
-            if types is not None:
+            if None not in types:
                 generator.types1.append(types[0])
                 generator.types2.append(types[1])
                 generator.k.append(float(piTorsion.attrib['k']))
@@ -2259,7 +2289,7 @@ class AmoebaTorsionTorsionGenerator:
 
         for torsionTorsion in element.findall('TorsionTorsion'):
             types = forceField._findAtomTypes(torsionTorsion, 5)
-            if types is not None:
+            if None not in types:
 
                 generator.types1.append(types[0])
                 generator.types2.append(types[1])
@@ -2496,7 +2526,7 @@ class AmoebaStretchBendGenerator:
 
         for stretchBend in element.findall('StretchBend'):
             types = forceField._findAtomTypes(stretchBend, 3)
-            if types is not None:
+            if None not in types:
 
                 generator.types1.append(types[0])
                 generator.types2.append(types[1])
@@ -2653,7 +2683,7 @@ class AmoebaVdwGenerator:
 
         for atom in element.findall('Vdw'):
             types = forceField._findAtomTypes(atom, 1)
-            if types is not None:
+            if None not in types:
 
                 values = [float(atom.attrib['sigma']), float(atom.attrib['epsilon']), float(atom.attrib['reduction'])]
 
@@ -2926,7 +2956,7 @@ class AmoebaMultipoleGenerator:
 
         for atom in element.findall('Multipole'):
             types = forceField._findAtomTypes(atom, 1)
-            if types is not None:
+            if None not in types:
 
                 # k-indices not provided default to 0
 
@@ -2983,7 +3013,7 @@ class AmoebaMultipoleGenerator:
  
         for atom in element.findall('Polarize'):
             types = forceField._findAtomTypes(atom, 1)
-            if types is not None:
+            if None not in types:
 
                 classIndex = atom.attrib['type']
                 polarizability = float(atom.attrib['polarizability'])
@@ -3511,7 +3541,7 @@ class AmoebaWcaDispersionGenerator:
 
         for atom in element.findall('WcaDispersion'):
             types = forceField._findAtomTypes(atom, 1)
-            if types is not None:
+            if None not in types:
 
                 values = [float(atom.attrib['radius']), float(atom.attrib['epsilon'])]
                 for t in types[0]:
@@ -3855,7 +3885,7 @@ class AmoebaUreyBradleyGenerator:
         forceField._forces.append(generator)
         for bond in element.findall('UreyBradley'):
             types = forceField._findAtomTypes(bond, 3)
-            if types is not None:
+            if None not in types:
 
                 generator.types1.append(types[0])
                 generator.types2.append(types[1])
@@ -3899,3 +3929,87 @@ class AmoebaUreyBradleyGenerator:
 parsers["AmoebaUreyBradleyForce"] = AmoebaUreyBradleyGenerator.parseElement
 
 #=============================================================================================
+
+
+## @private
+class DrudeGenerator:
+    """A DrudeGenerator constructs a DrudeForce."""
+    
+    def __init__(self):
+        self.typeMap = {}
+
+    @staticmethod
+    def parseElement(element, ff):
+        existing = [f for f in ff._forces if isinstance(f, DrudeGenerator)]
+        if len(existing) == 0:
+            generator = DrudeGenerator()
+            ff._forces.append(generator)
+        else:
+            # Multiple <DrudeForce> tags were found, probably in different files.  Simply add more types to the existing one.
+            generator = existing[0]
+        for particle in element.findall('Particle'):
+            types = ff._findAtomTypes(particle, 5)
+            if None not in types[:2]:
+                aniso12 = 0.0
+                aniso34 = 0.0
+                if 'aniso12' in particle.attrib:
+                    aniso12 = float(particle.attrib['aniso12'])
+                if 'aniso34' in particle.attrib:
+                    aniso34 = float(particle.attrib['aniso34'])
+                values = (types[1], types[2], types[3], types[4], float(particle.attrib['charge']), float(particle.attrib['polarizability']), aniso12, aniso34, float(particle.attrib['thole']))
+                for t in types[0]:
+                    generator.typeMap[t] = values
+    
+    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
+        force = mm.DrudeForce()
+        if not any(isinstance(f, mm.NonbondedForce) for f in sys.getForces()):
+            raise ValueError('<DrudeForce> must come after <NonbondedForce> in XML file')
+        
+        # Add Drude particles.
+        
+        drudeMap = {}
+        parentMap = {}
+        for atom in data.atoms:
+            t = data.atomType[atom]
+            if t in self.typeMap:
+                # Find other atoms in the residue that affect the Drude particle.
+                p = [-1, -1, -1, -1]
+                values = self.typeMap[t]
+                for atom2 in atom.residue.atoms():
+                    type2 = data.atomType[atom2]
+                    if type2 in values[0]:
+                        p[0] = atom2.index
+                    elif values[1] is not None and type2 in values[1]:
+                        p[1] = atom2.index
+                    elif values[2] is not None and type2 in values[2]:
+                        p[2] = atom2.index
+                    elif values[3] is not None and type2 in values[3]:
+                        p[3] = atom2.index
+                drudeIndex = force.addParticle(atom.index, p[0], p[1], p[2], p[3], values[4], values[5], values[6], values[7])
+                drudeMap[atom.index] = p[0]
+                parentMap[p[0]] = (atom.index, drudeIndex)
+        sys.addForce(force)
+    
+    def postprocessSystem(self, sys, data, args):
+        # For every nonbonded exclusion between Drude particles, add a screened pair.
+        
+        drude = [f for f in sys.getForces() if isinstance(f, mm.DrudeForce)][0]
+        nonbonded = [f for f in sys.getForces() if isinstance(f, mm.NonbondedForce)][0]
+        particleMap = {}
+        for i in range(drude.getNumParticles()):
+            particleMap[drude.getParticleParameters(i)[0]] = i
+        for i in range(nonbonded.getNumExceptions()):
+            (particle1, particle2, charge, sigma, epsilon) = nonbonded.getExceptionParameters(i)
+            if charge == 0 and epsilon == 0:
+                # This is an exclusion.
+                if particle1 in particleMap and particle2 in particleMap:
+                    # It connects two Drude particles, so add a screened pair.
+                    drude1 = particleMap[particle1]
+                    drude2 = particleMap[particle2]
+                    type1 = data.atomType[data.atoms[drude1]]
+                    type2 = data.atomType[data.atoms[drude2]]
+                    thole1 = self.typeMap[type1][8]
+                    thole2 = self.typeMap[type2][8]
+                    drude.addScreenedPair(drude1, drude2, thole1+thole2)
+
+parsers["DrudeForce"] = DrudeGenerator.parseElement
\ No newline at end of file

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

@@ -34,7 +34,7 @@ __author__ = "Peter Eastman"
 __version__ = "1.0"
 
 from simtk.openmm.app import Topology, PDBFile, ForceField
-from simtk.openmm.app.forcefield import HAngles, _createResidueSignature, _matchResidue
+from simtk.openmm.app.forcefield import HAngles, _createResidueSignature, _matchResidue, DrudeGenerator
 from simtk.openmm.app.topology import Residue
 from simtk.openmm.vec3 import Vec3
 from simtk.openmm import System, Context, NonbondedForce, VerletIntegrator, LocalEnergyMinimizer
@@ -817,6 +817,15 @@ class Modeller(object):
                 bondedToAtomNoEP[atom1.index].add(atom2.index)
                 bondedToAtomNoEP[atom2.index].add(atom1.index)
         
+        # If the force field has a DrudeForce, record the types of Drude particles and their parents since we'll
+        # need them for picking particle positions.
+        
+        drudeTypeMap = {}
+        for force in forcefield._forces:
+            if isinstance(force, DrudeGenerator):
+                for type in force.typeMap:
+                    drudeTypeMap[type] = force.typeMap[type][0]
+        
         # Create the new Topology.
         
         newTopology = Topology()
@@ -894,6 +903,12 @@ class Modeller(object):
                                         position = site.weights[0]*templateAtomPositions[index+site.atoms[0]] + site.weights[1]*templateAtomPositions[index+site.atoms[1]] + site.weights[2]*templateAtomPositions[index+site.atoms[2]]
                                     elif site.type == 'outOfPlane':
                                         position = site.weights[0]*templateAtomPositions[index+site.atoms[0]] + site.weights[1]*templateAtomPositions[index+site.atoms[1]] + site.weights[2]*templateAtomPositions[index+site.atoms[2]]
+                            if position is None and atom.type in drudeTypeMap:
+                                # This is a Drude particle.  Put it on top of its parent atom.
+                                
+                                for atom2, pos in zip(template.atoms, templateAtomPositions):
+                                    if atom2.type in drudeTypeMap[atom.type]:
+                                        position = deepcopy(pos)
                             if position is None:
                                 # We couldn't figure out the correct position.  As a wild guess, just put it at the center of the residue
                                 # and hope that energy minimization will fix it.