Merge pull request #481 from peastman/exclusions

Force field XML files can specify CustomNonbondedForces

docs/usersguide/application.rst

@@ -2186,6 +2186,49 @@ second atom has class OS and the third has class P:
 
     <Proper class1="" class2="OS" class3="P" class4="" per="3" phase="0.0" k="0.66944"/>
 
+<CustomNonbondedForce>
+===============
+
+To add a CustomNonbondedForce to the System, include a tag that looks like this:
+
+.. code-block:: xml
+
+    <CustomNonbondedForce energy="scale*epsilon1*epsilon2*((sigma1+sigma2)/r)^12" bondCutoff="3">
+     <GlobalParameter name="scale" defaultValue="1"/>
+     <PerParticleParameter name="sigma"/>
+     <PerParticleParameter name="epsilon"/>
+     <Atom type="0" sigma="0.3249" epsilon="0.7112"/>
+     <Atom type="1" sigma="0.1069" epsilon="0.0656"/>
+     <Atom type="2" sigma="0.3399" epsilon="0.4577"/>
+     ...
+    </CustomNonbondedForce>
+
+The energy expression for the CustomNonbondedForce is specified by the
+:code:`energy` attribute.  This is a mathematical expression that gives the
+energy of each pairwise interaction as a function of the distance *r*\ .  It
+also may depend on an arbitrary list of global or per-particle parameters.  Use
+a :code:`<GlobalParameter>` tag to define a global parameter, and a
+:code:`<PerParticleParameter>` tag to define a per-particle parameter.
+
+Exclusions are created automatically based on the :code:`bondCutoff` attribute.
+After setting the nonbonded parameters for all atoms, the force field calls
+:code:`createExclusionsFromBonds()` on the CustomNonbondedForce, passing in this
+value as its argument.  To avoid creating exclusions, set :code:`bondCutoff` to 0.
+
+Each :code:`<Atom>` tag specifies the parameters for one atom type
+(specified with the :code:`type` attribute) or atom class (specified with
+the :code:`class` attribute).  It is fine to mix these two methods, having
+some tags specify a type and others specify a class.  However you do it, you
+must make sure that a unique set of parameters is defined for every atom type.
+The remaining attributes are the values to use for the per-atom parameters. All
+per-atom parameters must be specified for every :code:`<Atom>` tag, and the
+attribute name must match the name of the parameter.  For instance, if there is
+a per-atom parameter with the name “radius”, then every :code:`<Atom>` tag
+must include an attribute called :code:`radius`\ .
+
+CustomNonbondedForce also allows you to define tabulated functions.  See section
+:ref:`tabulated-functions` for details.
+
 <CustomGBForce>
 ===============
 
@@ -2243,30 +2286,8 @@ attribute name must match the name of the parameter.  For instance, if there is
 a per-atom parameter with the name “radius”, then every :code:`<Atom>` tag
 must include an attribute called :code:`radius`\ .
 
-CustomGBForce also allows you to define tabulated functions.  To define a
-function, include a :code:`<Function>` tag inside the
-:code:`<CustomGBForce>` tag:
-
-.. code-block:: xml
-
-    <Function name="myfn" min="-5" max="5">
-    0.983674857694 -0.980096396266 -0.975743130031 -0.970451936613 -0.964027580076
-    -0.956237458128 -0.946806012846 -0.935409070603 -0.921668554406 -0.905148253645
-    -0.885351648202 -0.861723159313 -0.833654607012 -0.800499021761 -0.761594155956
-    -0.716297870199 -0.664036770268 -0.604367777117 -0.537049566998 -0.46211715726
-    -0.379948962255 -0.291312612452 -0.197375320225 -0.099667994625 0.0
-    0.099667994625 0.197375320225 0.291312612452 0.379948962255 0.46211715726
-    0.537049566998 0.604367777117 0.664036770268 0.716297870199 0.761594155956
-    0.800499021761 0.833654607012 0.861723159313 0.885351648202 0.905148253645
-    0.921668554406 0.935409070603 0.946806012846 0.956237458128 0.964027580076
-    0.970451936613 0.975743130031 0.980096396266 0.983674857694 0.986614298151
-    0.989027402201
-    </Function>
-
-The tag’s attributes define the name of the function and the range of values for
-which it is defined.  The tabulated values are listed inside the body of the
-tag, with successive values separated by white space.  Again, see the API
-documentation for more details.
+CustomGBForce also allows you to define tabulated functions.  See section
+:ref:`tabulated-functions` for details.
 
 Writing Custom Expressions
 ==========================
@@ -2303,6 +2324,55 @@ is exactly equivalent to
 The definition of an intermediate value may itself involve other intermediate
 values.  All uses of a value must appear *before* that value’s definition.
 
+.. _tabulated-functions:
+
+TabulatedFunctions
+==================
+
+Some forces, such as CustomNonbondedForce and CustomGBForce, allow you to define
+tabulated functions.  To define a function, include a :code:`<Function>` tag inside the
+:code:`<CustomNonbondedForce>` or :code:`<CustomGBForce>` tag:
+
+.. code-block:: xml
+
+    <Function name="myfn" type="Continuous1D" min="-5" max="5">
+    0.983674857694 -0.980096396266 -0.975743130031 -0.970451936613 -0.964027580076
+    -0.956237458128 -0.946806012846 -0.935409070603 -0.921668554406 -0.905148253645
+    -0.885351648202 -0.861723159313 -0.833654607012 -0.800499021761 -0.761594155956
+    -0.716297870199 -0.664036770268 -0.604367777117 -0.537049566998 -0.46211715726
+    -0.379948962255 -0.291312612452 -0.197375320225 -0.099667994625 0.0
+    0.099667994625 0.197375320225 0.291312612452 0.379948962255 0.46211715726
+    0.537049566998 0.604367777117 0.664036770268 0.716297870199 0.761594155956
+    0.800499021761 0.833654607012 0.861723159313 0.885351648202 0.905148253645
+    0.921668554406 0.935409070603 0.946806012846 0.956237458128 0.964027580076
+    0.970451936613 0.975743130031 0.980096396266 0.983674857694 0.986614298151
+    0.989027402201
+    </Function>
+
+The tag’s attributes define the name of the function, the type of function, and
+the range of values for which it is defined.  The required set of attributed
+depends on the function type:
+
+.. tabularcolumns:: |l|L|
+
+============  =======================================================
+Type          Required Attributes
+============  =======================================================
+Continuous1D  min, max
+Continuous2D  xmin, ymin, xmax, ymax, xsize, ysize
+Continuous3D  xmin, ymin, zmin, xmax, ymax, zmax, xsize, ysize, zsize
+Discrete1D
+Discrete2D    xsize, ysize
+Discrete3D    xsize, ysize, zsize
+============  =======================================================
+
+
+The "min" and "max" attributes define the range of the independent variables for
+a continuous function.  The "size" attributes define the size of the table along
+each axis.  The tabulated values are listed inside the body of the tag, with
+successive values separated by white space.  See the API documentation for more
+details.
+
 
 Using Multiple Files
 ********************

openmmapi/include/openmm/CustomNonbondedForce.h

@@ -343,6 +343,8 @@ public:
     void setParticleParameters(int index, const std::vector<double>& parameters);
     /**
      * Add a particle pair to the list of interactions that should be excluded.
+     * 
+     * In many cases, you can use createExclusionsFromBonds() rather than adding each exclusion explicitly.
      *
      * @param particle1  the index of the first particle in the pair
      * @param particle2  the index of the second particle in the pair
@@ -365,6 +367,15 @@ public:
      * @param particle2  the index of the second particle in the pair
      */
     void setExclusionParticles(int index, int particle1, int particle2);
+    /**
+     * Identify exclusions based on the molecular topology.  Particles which are separated by up to a specified number of
+     * bonds are added as exclusions.
+     *
+     * @param bonds       the set of bonds based on which to construct exclusions.  Each element specifies the indices of
+     *                    two particles that are bonded to each other.
+     * @param bondCutoff  pairs of particles that are separated by this many bonds or fewer are added to the list of exclusions
+     */
+    void createExclusionsFromBonds(const std::vector<std::pair<int, int> >& bonds, int bondCutoff);
     /**
      * Add a tabulated function that may appear in the energy expression.
      *

openmmapi/src/CustomNonbondedForce.cpp

@@ -174,6 +174,33 @@ void CustomNonbondedForce::setExclusionParticles(int index, int particle1, int p
     exclusions[index].particle1 = particle1;
     exclusions[index].particle2 = particle2;
 }
+
+void CustomNonbondedForce::createExclusionsFromBonds(const vector<pair<int, int> >& bonds, int bondCutoff) {
+    if (bondCutoff < 1)
+        return;
+    vector<set<int> > exclusions(particles.size());
+    vector<set<int> > bonded12(exclusions.size());
+    for (int i = 0; i < (int) bonds.size(); ++i) {
+        int p1 = bonds[i].first;
+        int p2 = bonds[i].second;
+        exclusions[p1].insert(p2);
+        exclusions[p2].insert(p1);
+        bonded12[p1].insert(p2);
+        bonded12[p2].insert(p1);
+    }
+    for (int level = 0; level < bondCutoff-1; level++) {
+        vector<set<int> > currentExclusions = exclusions;
+        for (int i = 0; i < (int) particles.size(); i++) {
+            for (set<int>::const_iterator iter = currentExclusions[i].begin(); iter != currentExclusions[i].end(); ++iter)
+                exclusions[*iter].insert(bonded12[i].begin(), bonded12[i].end());
+        }
+    }
+    for (int i = 0; i < (int) exclusions.size(); ++i)
+        for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter)
+            if (*iter < i)
+                addExclusion(*iter, i);
+}
+
 int CustomNonbondedForce::addTabulatedFunction(const std::string& name, TabulatedFunction* function) {
     functions.push_back(FunctionInfo(name, function));
     return functions.size()-1;

tests/TestFindExclusions.cpp

@@ -30,6 +30,7 @@
  * -------------------------------------------------------------------------- */
 
 #include "openmm/internal/AssertionUtilities.h"
+#include "openmm/CustomNonbondedForce.h"
 #include "openmm/NonbondedForce.h"
 #include "openmm/OpenMMException.h"
 #include <iostream>
@@ -166,10 +167,62 @@ void testReplaceExceptions() {
     ASSERT(charge == 5.0);
 }
 
+/**
+ * This is the same as testFindExceptions(), except it tests adding exclusions to a CustomNonbondedForce.
+ */
+
+void testFindCustomExclusions() {
+    CustomNonbondedForce nonbonded("r");
+    vector<pair<int, int> > bonds;
+    vector<double> params;
+    for (int i = 0; i < NUM_ATOMS; i++)
+        nonbonded.addParticle(params);
+    // loop over all main-chain atoms (even numbered atoms)
+    for (int i = 0; i < NUM_ATOMS-1; i += 2)
+    {
+        // side-chain bonds
+        bonds.push_back(pair<int, int>(i, i+1));
+        // main-chain bonds
+        if (i < NUM_ATOMS-2) // penultimate atom (NUM_ATOMS-2) has no subsequent main-chain atom
+            bonds.push_back(pair<int, int>(i, i+2));
+    }
+    nonbonded.createExclusionsFromBonds(bonds, 3);
+
+    // Build lists of the expected exclusions.
+
+    vector<set<int> > expectedExclusions(NUM_ATOMS);
+    int totalExclusions = 0;
+    for (int i = 0; i < NUM_ATOMS; i += 2) {
+        addAtomsToExclusions(i, i+1, expectedExclusions, totalExclusions);
+        addAtomsToExclusions(i, i+2, expectedExclusions, totalExclusions);
+        addAtomsToExclusions(i, i+3, expectedExclusions, totalExclusions);
+        addAtomsToExclusions(i, i+4, expectedExclusions, totalExclusions);
+        addAtomsToExclusions(i+1, i+2, expectedExclusions, totalExclusions);
+        addAtomsToExclusions(i, i+5, expectedExclusions, totalExclusions);
+        addAtomsToExclusions(i, i+6, expectedExclusions, totalExclusions);
+        addAtomsToExclusions(i+1, i+3, expectedExclusions, totalExclusions);
+        addAtomsToExclusions(i+1, i+4, expectedExclusions, totalExclusions);
+    }
+    for (int i = 0; i < nonbonded.getNumExclusions(); i++) {
+        int particle1, particle2;
+        nonbonded.getExclusionParticles(i, particle1, particle2);
+    }
+
+    // Compare them to the exceptions that were generated.
+
+    ASSERT_EQUAL(totalExclusions, nonbonded.getNumExclusions());
+    for (int i = 0; i < nonbonded.getNumExclusions(); i++) {
+        int particle1, particle2;
+        nonbonded.getExclusionParticles(i, particle1, particle2);
+        ASSERT(expectedExclusions[particle1].find(particle2) != expectedExclusions[particle1].end());
+    }
+}
+
 int main() {
     try {
         testFindExceptions();
         testReplaceExceptions();
+        testFindCustomExclusions();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;

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

@@ -6,7 +6,7 @@ Simbios, the NIH National Center for Physics-Based Simulation of
 Biological Structures at Stanford, funded under the NIH Roadmap for
 Medical Research, grant U54 GM072970. See https://simtk.org.
 
-Portions copyright (c) 2012-2013 Stanford University and the Authors.
+Portions copyright (c) 2012-2014 Stanford University and the Authors.
 Authors: Peter Eastman, Mark Friedrichs
 Contributors:
 
@@ -1485,6 +1485,104 @@ class CustomTorsionGenerator:
 parsers["CustomTorsionForce"] = CustomTorsionGenerator.parseElement
 
 
+## @private
+class CustomNonbondedGenerator:
+    """A CustomNonbondedGenerator constructs a CustomNonbondedForce."""
+
+    def __init__(self, energy, bondCutoff):
+        self.energy = energy
+        self.bondCutoff = bondCutoff
+        self.typeMap = {}
+        self.globalParams = {}
+        self.perParticleParams = []
+        self.functions = []
+
+    @staticmethod
+    def parseElement(element, ff):
+        generator = CustomNonbondedGenerator(element.attrib['energy'], int(element.attrib['bondCutoff']))
+        ff._forces.append(generator)
+        for param in element.findall('GlobalParameter'):
+            generator.globalParams[param.attrib['name']] = float(param.attrib['defaultValue'])
+        for param in element.findall('PerParticleParameter'):
+            generator.perParticleParams.append(param.attrib['name'])
+        for atom in element.findall('Atom'):
+            types = ff._findAtomTypes(atom, 1)
+            if None not in types:
+                values = [float(atom.attrib[param]) for param in generator.perParticleParams]
+                for t in types[0]:
+                    generator.typeMap[t] = values
+
+    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
+        methodMap = {NoCutoff:mm.CustomNonbondedForce.NoCutoff,
+                     CutoffNonPeriodic:mm.CustomNonbondedForce.CutoffNonPeriodic,
+                     CutoffPeriodic:mm.CustomNonbondedForce.CutoffPeriodic}
+        if nonbondedMethod not in methodMap:
+            raise ValueError('Illegal nonbonded method for CustomNonbondedForce')
+        force = mm.CustomNonbondedForce(self.energy)
+        for param in self.globalParams:
+            force.addGlobalParameter(param, self.globalParams[param])
+        for param in self.perParticleParams:
+            force.addPerParticleParameter(param)
+        for (name, type, values, params) in self.functions:
+            if type == 'Continuous1D':
+                force.addTabulatedFunction(name, mm.Continuous1DFunction(values, params['min'], params['max']))
+            elif type == 'Continuous2D':
+                force.addTabulatedFunction(name, mm.Continuous2DFunction(params['xsize'], params['ysize'], values, params['xmin'], params['xmax'], params['ymin'], params['ymax']))
+            elif type == 'Continuous3D':
+                force.addTabulatedFunction(name, mm.Continuous2DFunction(params['xsize'], params['ysize'], params['zsize'], values, params['xmin'], params['xmax'], params['ymin'], params['ymax'], params['zmin'], params['zmax']))
+            elif type == 'Discrete1D':
+                force.addTabulatedFunction(name, mm.Discrete1DFunction(values))
+            elif type == 'Discrete2D':
+                force.addTabulatedFunction(name, mm.Discrete2DFunction(params['xsize'], params['ysize'], values))
+            elif type == 'Discrete3D':
+                force.addTabulatedFunction(name, mm.Discrete2DFunction(params['xsize'], params['ysize'], params['zsize'], values))
+        for atom in data.atoms:
+            t = data.atomType[atom]
+            if t in self.typeMap:
+                values = self.typeMap[t]
+                force.addParticle(self.typeMap[t])
+            else:
+                raise ValueError('No CustomNonbonded parameters defined for atom type '+t)
+        force.setNonbondedMethod(methodMap[nonbondedMethod])
+        force.setCutoffDistance(nonbondedCutoff)
+        sys.addForce(force)
+
+    def postprocessSystem(self, sys, data, args):
+        # Create exceptions based on bonds, virtual sites, and Drude particles.
+        if self.bondCutoff == 0:
+            return
+        bondIndices = []
+        for bond in data.bonds:
+            bondIndices.append((bond.atom1, bond.atom2))
+        for i in range(sys.getNumParticles()):
+            if sys.isVirtualSite(i):
+                site = sys.getVirtualSite(i)
+                for j in range(site.getNumParticles()):
+                    bondIndices.append((i, site.getParticle(j)))
+        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.CustomNonbondedForce)][0]
+        nonbonded.createExclusionsFromBonds(bondIndices, self.bondCutoff)
+
+parsers["CustomNonbondedForce"] = CustomNonbondedGenerator.parseElement
+
+
 ## @private
 class CustomGBGenerator:
     """A CustomGBGenerator constructs a CustomGBForce."""
@@ -1493,7 +1591,6 @@ class CustomGBGenerator:
         self.typeMap = {}
         self.globalParams = {}
         self.perParticleParams = []
-        self.paramValues = []
         self.computedValues = []
         self.energyTerms = []
         self.functions = []