Began working on CPU implementation of CustomManyParticleForce

openmmapi/include/openmm/CustomManyParticleForce.h

@@ -36,6 +36,7 @@
 #include "TabulatedFunction.h"
 #include "internal/windowsExport.h"
 #include <set>
+#include <string>
 #include <vector>
 
 namespace OpenMM {

platforms/cpu/include/CompiledExpressionSet.h

+#ifndef OPENMM_COMPILEDEXPRESSIONSET_H_
+#define OPENMM_COMPILEDEXPRESSIONSET_H_
+
+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * 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) 2014 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "lepton/CompiledExpression.h"
+#include "windowsExportCpu.h"
+#include <string>
+#include <vector>
+
+namespace OpenMM {
+
+/**
+ * This class simplifies the management of a set of related CompiledExpressions that share variables.
+ */
+class OPENMM_EXPORT_CPU CompiledExpressionSet {
+public:
+    CompiledExpressionSet();
+    /**
+     * Add a CompiledExpression to the set.
+     */
+    void registerExpression(Lepton::CompiledExpression& expression);
+    /**
+     * Get the index of a particular variable.
+     */
+    int getVariableIndex(const std::string& name);
+    /**
+     * Set the value of a variable on every CompiledExpression.
+     * 
+     * @param index    the index of the variable, as returned by getVariableIndex()
+     * @param value    the value to set it to
+     */
+    void setVariable(int index, double value);
+private:
+    std::vector<Lepton::CompiledExpression*> expressions;
+    std::vector<std::string> variables;
+    std::vector<std::vector<double*> > variableReferences;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_COMPILEDEXPRESSIONSET_H_*/

platforms/cpu/include/CpuCustomManyParticleForce.h

+
+/* Portions copyright (c) 2009-2014 Stanford University and Simbios.
+ * Contributors: Peter Eastman
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+#ifndef OPENMM_CPU_CUSTOM_MANY_PARTICLE_FORCE_H__
+#define OPENMM_CPU_CUSTOM_MANY_PARTICLE_FORCE_H__
+
+#include "ReferenceForce.h"
+#include "ReferenceBondIxn.h"
+#include "CompiledExpressionSet.h"
+#include "openmm/CustomManyParticleForce.h"
+#include "openmm/internal/vectorize.h"
+#include "lepton/CompiledExpression.h"
+#include "lepton/ParsedExpression.h"
+#include <map>
+#include <set>
+#include <vector>
+
+namespace OpenMM {
+
+class CpuCustomManyParticleForce {
+private:
+
+    class ParticleTermInfo;
+    class DistanceTermInfo;
+    class AngleTermInfo;
+    class DihedralTermInfo;
+    int numParticlesPerSet, numPerParticleParameters, numTypes;
+    bool useCutoff, usePeriodic;
+    RealOpenMM cutoffDistance;
+    RealOpenMM periodicBoxSize[3];
+    CompiledExpressionSet expressionSet;
+    Lepton::CompiledExpression energyExpression;
+    std::vector<std::vector<int> > particleParamIndices;
+    std::vector<std::set<int> > exclusions;
+    std::vector<int> particleTypes;
+    std::vector<int> orderIndex;
+    std::vector<std::vector<int> > particleOrder;
+    std::vector<ParticleTermInfo> particleTerms;
+    std::vector<DistanceTermInfo> distanceTerms;
+    std::vector<AngleTermInfo> angleTerms;
+    std::vector<DihedralTermInfo> dihedralTerms;
+
+    void loopOverInteractions(std::vector<int>& particles, int loopIndex, float* posq, std::vector<OpenMM::RealVec>& atomCoordinates,
+                              RealOpenMM** particleParameters, std::vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy, const fvec4& boxSize, const fvec4& invBoxSize);
+
+    /**---------------------------------------------------------------------------------------
+
+       Calculate custom interaction for one set of particles
+
+       @param particles        the indices of the particles
+       @param posq             atom coordinates in float format
+       @param atomCoordinates  atom coordinates
+       @param forces           force array (forces added)
+       @param totalEnergy      total energy
+
+       --------------------------------------------------------------------------------------- */
+
+    void calculateOneIxn(const std::vector<int>& particles, float* posq, std::vector<OpenMM::RealVec>& atomCoordinates,
+                         std::vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy, const fvec4& boxSize, const fvec4& invBoxSize);
+
+    /**
+     * Compute the displacement and squared distance between two points, optionally using
+     * periodic boundary conditions.
+     */
+    void getDeltaR(const fvec4& posI, const fvec4& posJ, fvec4& deltaR, float& r2, const fvec4& boxSize, const fvec4& invBoxSize) const;
+
+    void computeDelta(int atom1, int atom2, RealOpenMM* delta, std::vector<OpenMM::RealVec>& atomCoordinates) const;
+
+    static RealOpenMM computeAngle(RealOpenMM* vec1, RealOpenMM* vec2);
+
+
+public:
+
+    /**---------------------------------------------------------------------------------------
+
+       Constructor
+
+       --------------------------------------------------------------------------------------- */
+
+    CpuCustomManyParticleForce(const OpenMM::CustomManyParticleForce& force);
+
+    /**---------------------------------------------------------------------------------------
+
+       Destructor
+
+       --------------------------------------------------------------------------------------- */
+
+    ~CpuCustomManyParticleForce();
+
+    /**---------------------------------------------------------------------------------------
+
+       Set the force to use a cutoff.
+
+       @param distance            the cutoff distance
+
+       --------------------------------------------------------------------------------------- */
+
+    void setUseCutoff(RealOpenMM distance);
+
+    /**---------------------------------------------------------------------------------------
+
+       Set the force to use periodic boundary conditions.  This requires that a cutoff has
+       already been set, and the smallest side of the periodic box is at least twice the cutoff
+       distance.
+
+       @param boxSize             the X, Y, and Z widths of the periodic box
+
+       --------------------------------------------------------------------------------------- */
+
+    void setPeriodic(OpenMM::RealVec& boxSize);
+
+    /**---------------------------------------------------------------------------------------
+
+       Calculate the interaction
+
+       @param posq               atom coordinates in float format
+       @param atomCoordinates    atom coordinates
+       @param particleParameters particle parameter values (particleParameters[particleIndex][parameterIndex])
+       @param globalParameters   the values of global parameters
+       @param forces             force array (forces added)
+       @param totalEnergy        total energy
+
+       --------------------------------------------------------------------------------------- */
+
+    void calculateIxn(float* posq, std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** particleParameters,
+                      const std::map<std::string, double>& globalParameters,
+                      std::vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy);
+
+// ---------------------------------------------------------------------------------------
+
+};
+
+class CpuCustomManyParticleForce::ParticleTermInfo {
+public:
+    std::string name;
+    int atom, component, variableIndex;
+    Lepton::CompiledExpression forceExpression;
+    ParticleTermInfo(const std::string& name, int atom, int component, const Lepton::CompiledExpression& forceExpression, CompiledExpressionSet& set) :
+            name(name), atom(atom), component(component), forceExpression(forceExpression) {
+        variableIndex = set.getVariableIndex(name);
+    }
+};
+
+class CpuCustomManyParticleForce::DistanceTermInfo {
+public:
+    std::string name;
+    int p1, p2, variableIndex;
+    Lepton::CompiledExpression forceExpression;
+    mutable RealOpenMM delta[ReferenceForce::LastDeltaRIndex];
+    DistanceTermInfo(const std::string& name, const std::vector<int>& atoms, const Lepton::CompiledExpression& forceExpression, CompiledExpressionSet& set) :
+            name(name), p1(atoms[0]), p2(atoms[1]), forceExpression(forceExpression) {
+        variableIndex = set.getVariableIndex(name);
+    }
+};
+
+class CpuCustomManyParticleForce::AngleTermInfo {
+public:
+    std::string name;
+    int p1, p2, p3, variableIndex;
+    Lepton::CompiledExpression forceExpression;
+    mutable RealOpenMM delta1[ReferenceForce::LastDeltaRIndex];
+    mutable RealOpenMM delta2[ReferenceForce::LastDeltaRIndex];
+    AngleTermInfo(const std::string& name, const std::vector<int>& atoms, const Lepton::CompiledExpression& forceExpression, CompiledExpressionSet& set) :
+            name(name), p1(atoms[0]), p2(atoms[1]), p3(atoms[2]), forceExpression(forceExpression) {
+        variableIndex = set.getVariableIndex(name);
+    }
+};
+
+class CpuCustomManyParticleForce::DihedralTermInfo {
+public:
+    std::string name;
+    int p1, p2, p3, p4, variableIndex;
+    Lepton::CompiledExpression forceExpression;
+    mutable RealOpenMM delta1[ReferenceForce::LastDeltaRIndex];
+    mutable RealOpenMM delta2[ReferenceForce::LastDeltaRIndex];
+    mutable RealOpenMM delta3[ReferenceForce::LastDeltaRIndex];
+    mutable RealOpenMM cross1[3];
+    mutable RealOpenMM cross2[3];
+    DihedralTermInfo(const std::string& name, const std::vector<int>& atoms, const Lepton::CompiledExpression& forceExpression, CompiledExpressionSet& set) :
+            name(name), p1(atoms[0]), p2(atoms[1]), p3(atoms[2]), p4(atoms[3]), forceExpression(forceExpression) {
+        variableIndex = set.getVariableIndex(name);
+    }
+};
+
+} // namespace OpenMM
+
+#endif // OPENMM_CPU_CUSTOM_MANY_PARTICLE_FORCE_H__

platforms/cpu/include/CpuKernels.h

@@ -33,6 +33,7 @@
  * -------------------------------------------------------------------------- */
 
 #include "CpuBondForce.h"
+#include "CpuCustomManyParticleForce.h"
 #include "CpuCustomNonbondedForce.h"
 #include "CpuGBSAOBCForce.h"
 #include "CpuLangevinDynamics.h"
@@ -302,6 +303,48 @@ private:
     CpuGBSAOBCForce obc;
 };
 
+/**
+ * This kernel is invoked by CustomManyParticleForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CpuCalcCustomManyParticleForceKernel : public CalcCustomManyParticleForceKernel {
+public:
+    CpuCalcCustomManyParticleForceKernel(std::string name, const Platform& platform, CpuPlatform::PlatformData& data) : CalcCustomManyParticleForceKernel(name, platform),
+            data(data), ixn(NULL) {
+    }
+    ~CpuCalcCustomManyParticleForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CustomManyParticleForce this kernel will be used for
+     */
+    void initialize(const System& system, const CustomManyParticleForce& force);
+    /**
+     * Execute the kernel to calculate the forces and/or energy.
+     *
+     * @param context        the context in which to execute this kernel
+     * @param includeForces  true if forces should be calculated
+     * @param includeEnergy  true if the energy should be calculated
+     * @return the potential energy due to the force
+     */
+    double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CustomManyParticleForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomManyParticleForce& force);
+private:
+    CpuPlatform::PlatformData& data;
+    int numParticles;
+    RealOpenMM cutoffDistance;
+    RealOpenMM **particleParamArray;
+    CpuCustomManyParticleForce* ixn;
+    std::vector<std::string> globalParameterNames;
+    NonbondedMethod nonbondedMethod;
+};
+
 /**
  * This kernel is invoked by LangevinIntegrator to take one time step.
  */

platforms/cpu/src/CompiledExpressionSet.cpp

+/* Portions copyright (c) 2014 Stanford University and Simbios.
+ * Contributors: Peter Eastman
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+#include "CompiledExpressionSet.h"
+
+using namespace OpenMM;
+using namespace Lepton;
+using namespace std;
+
+CompiledExpressionSet::CompiledExpressionSet() {
+}
+
+void CompiledExpressionSet::registerExpression(Lepton::CompiledExpression& expression) {
+    expressions.push_back(&expression);
+    for (int i = 0; i < (int) variables.size(); i++)
+        if (expression.getVariables().find(variables[i]) != expression.getVariables().end())
+            variableReferences[i].push_back(&expression.getVariableReference(variables[i]));
+}
+
+int CompiledExpressionSet::getVariableIndex(const std::string& name) {
+    for (int i = 0; i < (int) variables.size(); i++)
+        if (variables[i] == name)
+            return i;
+    int index = variables.size();
+    variables.push_back(name);
+    variableReferences.push_back(vector<double*>());
+    for (int i = 0; i < (int) expressions.size(); i++)
+        if (expressions[i]->getVariables().find(name) != expressions[i]->getVariables().end())
+            variableReferences[index].push_back(&expressions[i]->getVariableReference(name));
+    return index;
+}
+
+void CompiledExpressionSet::setVariable(int index, double value) {
+    for (int i = 0; i < (int) variableReferences[index].size(); i++)
+        *variableReferences[index][i] = value;
+}

platforms/cpu/src/CpuCustomManyParticleForce.cpp

+
+/* Portions copyright (c) 2009-2014 Stanford University and Simbios.
+ * Contributors: Peter Eastman
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+#include <string.h>
+#include <sstream>
+#include <utility>
+
+#include "SimTKOpenMMCommon.h"
+#include "SimTKOpenMMLog.h"
+#include "SimTKOpenMMUtilities.h"
+#include "ReferenceForce.h"
+#include "CpuCustomManyParticleForce.h"
+#include "ReferenceTabulatedFunction.h"
+#include "openmm/internal/CustomManyParticleForceImpl.h"
+#include "lepton/CustomFunction.h"
+
+using namespace OpenMM;
+using namespace std;
+
+CpuCustomManyParticleForce::CpuCustomManyParticleForce(const CustomManyParticleForce& force) : useCutoff(false), usePeriodic(false) {
+    numParticlesPerSet = force.getNumParticlesPerSet();
+    numPerParticleParameters = force.getNumPerParticleParameters();
+    
+    // Create custom functions for the tabulated functions.
+
+    map<string, Lepton::CustomFunction*> functions;
+    for (int i = 0; i < (int) force.getNumTabulatedFunctions(); i++)
+        functions[force.getTabulatedFunctionName(i)] = createReferenceTabulatedFunction(force.getTabulatedFunction(i));
+
+    // Parse the expression and create the object used to calculate the interaction.
+
+    map<string, vector<int> > distances;
+    map<string, vector<int> > angles;
+    map<string, vector<int> > dihedrals;
+    Lepton::ParsedExpression energyExpr = CustomManyParticleForceImpl::prepareExpression(force, functions, distances, angles, dihedrals);
+    energyExpression = energyExpr.createCompiledExpression();
+    expressionSet.registerExpression(energyExpression);
+    vector<string> particleParameterNames;
+    if (force.getNonbondedMethod() != CustomManyParticleForce::NoCutoff)
+        setUseCutoff(force.getCutoffDistance());
+
+    // Delete the custom functions.
+
+    for (map<string, Lepton::CustomFunction*>::iterator iter = functions.begin(); iter != functions.end(); iter++)
+        delete iter->second;
+    
+    // Differentiate the energy to get expressions for the force.
+
+    particleParamIndices.resize(numParticlesPerSet);
+    for (int i = 0; i < numParticlesPerSet; i++) {
+        stringstream xname, yname, zname;
+        xname << 'x' << (i+1);
+        yname << 'y' << (i+1);
+        zname << 'z' << (i+1);
+        particleTerms.push_back(CpuCustomManyParticleForce::ParticleTermInfo(xname.str(), i, 0, energyExpr.differentiate(xname.str()).optimize().createCompiledExpression(), expressionSet));
+        particleTerms.push_back(CpuCustomManyParticleForce::ParticleTermInfo(yname.str(), i, 1, energyExpr.differentiate(yname.str()).optimize().createCompiledExpression(), expressionSet));
+        particleTerms.push_back(CpuCustomManyParticleForce::ParticleTermInfo(zname.str(), i, 2, energyExpr.differentiate(zname.str()).optimize().createCompiledExpression(), expressionSet));
+        for (int j = 0; j < numPerParticleParameters; j++) {
+            stringstream paramname;
+            paramname << force.getPerParticleParameterName(j) << (i+1);
+            particleParamIndices[i].push_back(expressionSet.getVariableIndex(paramname.str()));
+        }
+    }
+    for (map<string, vector<int> >::const_iterator iter = distances.begin(); iter != distances.end(); ++iter)
+        distanceTerms.push_back(CpuCustomManyParticleForce::DistanceTermInfo(iter->first, iter->second, energyExpr.differentiate(iter->first).optimize().createCompiledExpression(), expressionSet));
+    for (map<string, vector<int> >::const_iterator iter = angles.begin(); iter != angles.end(); ++iter)
+        angleTerms.push_back(CpuCustomManyParticleForce::AngleTermInfo(iter->first, iter->second, energyExpr.differentiate(iter->first).optimize().createCompiledExpression(), expressionSet));
+    for (map<string, vector<int> >::const_iterator iter = dihedrals.begin(); iter != dihedrals.end(); ++iter)
+        dihedralTerms.push_back(CpuCustomManyParticleForce::DihedralTermInfo(iter->first, iter->second, energyExpr.differentiate(iter->first).optimize().createCompiledExpression(), expressionSet));
+    for (int i = 0; i < particleTerms.size(); i++)
+        expressionSet.registerExpression(particleTerms[i].forceExpression);
+    for (int i = 0; i < distanceTerms.size(); i++)
+        expressionSet.registerExpression(distanceTerms[i].forceExpression);
+    for (int i = 0; i < angleTerms.size(); i++)
+        expressionSet.registerExpression(angleTerms[i].forceExpression);
+    for (int i = 0; i < dihedralTerms.size(); i++)
+        expressionSet.registerExpression(dihedralTerms[i].forceExpression);
+    
+    // Record exclusions.
+    
+    exclusions.resize(force.getNumParticles());
+    for (int i = 0; i < (int) force.getNumExclusions(); i++) {
+        int p1, p2;
+        force.getExclusionParticles(i, p1, p2);
+        exclusions[p1].insert(p2);
+        exclusions[p2].insert(p1);
+    }
+    
+    // Record information about type filters.
+    
+    CustomManyParticleForceImpl::buildFilterArrays(force, numTypes, particleTypes, orderIndex, particleOrder);
+}
+
+CpuCustomManyParticleForce::~CpuCustomManyParticleForce( ){
+}
+
+void CpuCustomManyParticleForce::calculateIxn(float* posq, vector<RealVec>& atomCoordinates, RealOpenMM** particleParameters,
+                                                  const map<string, double>& globalParameters, vector<RealVec>& forces,
+                                                  RealOpenMM* totalEnergy) {
+    for (map<string, double>::const_iterator iter = globalParameters.begin(); iter != globalParameters.end(); ++iter)
+        expressionSet.setVariable(expressionSet.getVariableIndex(iter->first), iter->second);
+    vector<int> particles(numParticlesPerSet);
+    fvec4 boxSize(periodicBoxSize[0], periodicBoxSize[1], periodicBoxSize[2], 0);
+    fvec4 invBoxSize((1/periodicBoxSize[0]), (1/periodicBoxSize[1]), (1/periodicBoxSize[2]), 0);
+    loopOverInteractions(particles, 0, posq, atomCoordinates, particleParameters, forces, totalEnergy, boxSize, invBoxSize);
+}
+
+void CpuCustomManyParticleForce::setUseCutoff(RealOpenMM distance) {
+    useCutoff = true;
+    cutoffDistance = distance;
+}
+
+void CpuCustomManyParticleForce::setPeriodic(RealVec& boxSize) {
+    assert(useCutoff);
+    assert(boxSize[0] >= 2.0*cutoffDistance);
+    assert(boxSize[1] >= 2.0*cutoffDistance);
+    assert(boxSize[2] >= 2.0*cutoffDistance);
+    usePeriodic = true;
+    periodicBoxSize[0] = boxSize[0];
+    periodicBoxSize[1] = boxSize[1];
+    periodicBoxSize[2] = boxSize[2];
+}
+
+void CpuCustomManyParticleForce::loopOverInteractions(vector<int>& particles, int loopIndex, float* posq, vector<OpenMM::RealVec>& atomCoordinates,
+                                                          RealOpenMM** particleParameters, vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy, const fvec4& boxSize, const fvec4& invBoxSize) {
+    int numParticles = atomCoordinates.size();
+    int start = (loopIndex == 0 ? 0 : particles[loopIndex-1]+1);
+    for (int i = start; i < numParticles; i++) {
+        particles[loopIndex] = i;
+        for (int j = 0; j < numPerParticleParameters; j++)
+            expressionSet.setVariable(particleParamIndices[loopIndex][j], particleParameters[i][j]);
+        if (loopIndex == numParticlesPerSet-1)
+            calculateOneIxn(particles, posq, atomCoordinates, forces, totalEnergy, boxSize, invBoxSize);
+        else
+            loopOverInteractions(particles, loopIndex+1, posq, atomCoordinates, particleParameters, forces, totalEnergy, boxSize, invBoxSize);
+    }
+}
+
+void CpuCustomManyParticleForce::calculateOneIxn(const vector<int>& particles, float* posq, vector<RealVec>& atomCoordinates, vector<RealVec>& forces, RealOpenMM* totalEnergy, const fvec4& boxSize, const fvec4& invBoxSize) {
+    // Select the ordering to use for the particles.
+    
+    vector<int> permutedParticles(numParticlesPerSet);
+    if (particleOrder.size() == 1) {
+        // There are no filters, so we don't need to worry about ordering.
+        
+        permutedParticles = particles;
+    }
+    else {
+        int index = 0;
+        for (int i = numParticlesPerSet-1; i >= 0; i--)
+            index = particleTypes[particles[i]]+numTypes*index;
+        int order = orderIndex[index];
+        if (order == -1)
+            return;
+        for (int i = 0; i < numParticlesPerSet; i++)
+            permutedParticles[i] = particles[particleOrder[order][i]];
+    }
+    
+    // Decide whether to include this interaction.
+    
+    for (int i = 0; i < (int) permutedParticles.size(); i++) {
+        int p1 = permutedParticles[i];
+        for (int j = i+1; j < (int) permutedParticles.size(); j++) {
+            int p2 = permutedParticles[j];
+            if (exclusions[p1].find(p2) != exclusions[p1].end())
+                return;
+            if (useCutoff) {
+//                fvec4 deltaR;
+//                fvec4 posI(posq+4*p1);
+//                fvec4 posJ(posq+4*p2);
+//                float r2;
+//                getDeltaR(posI, posJ, deltaR, r2, boxSize, invBoxSize);
+//                if (r2 >= cutoffDistance*cutoffDistance)
+//                    return;
+//
+                RealOpenMM delta[ReferenceForce::LastDeltaRIndex];
+                computeDelta(p1, p2, delta, atomCoordinates);
+                if (delta[ReferenceForce::RIndex] >= cutoffDistance)
+                    return;
+            }
+        }
+    }
+    
+    // Compute all of the variables the energy can depend on.
+
+    for (int i = 0; i < (int) particleTerms.size(); i++) {
+        const ParticleTermInfo& term = particleTerms[i];
+        expressionSet.setVariable(term.variableIndex, atomCoordinates[term.atom][term.component]);
+    }
+    for (int i = 0; i < (int) distanceTerms.size(); i++) {
+        const DistanceTermInfo& term = distanceTerms[i];
+        computeDelta(permutedParticles[term.p1], permutedParticles[term.p2], term.delta, atomCoordinates);
+        expressionSet.setVariable(term.variableIndex, term.delta[ReferenceForce::RIndex]);
+    }
+    for (int i = 0; i < (int) angleTerms.size(); i++) {
+        const AngleTermInfo& term = angleTerms[i];
+        computeDelta(permutedParticles[term.p1], permutedParticles[term.p2], term.delta1, atomCoordinates);
+        computeDelta(permutedParticles[term.p3], permutedParticles[term.p2], term.delta2, atomCoordinates);
+        expressionSet.setVariable(term.variableIndex, computeAngle(term.delta1, term.delta2));
+    }
+    for (int i = 0; i < (int) dihedralTerms.size(); i++) {
+        const DihedralTermInfo& term = dihedralTerms[i];
+        computeDelta(permutedParticles[term.p2], permutedParticles[term.p1], term.delta1, atomCoordinates);
+        computeDelta(permutedParticles[term.p2], permutedParticles[term.p3], term.delta2, atomCoordinates);
+        computeDelta(permutedParticles[term.p4], permutedParticles[term.p3], term.delta3, atomCoordinates);
+        RealOpenMM dotDihedral, signOfDihedral;
+        RealOpenMM* crossProduct[] = {term.cross1, term.cross2};
+        expressionSet.setVariable(term.variableIndex, ReferenceBondIxn::getDihedralAngleBetweenThreeVectors(term.delta1, term.delta2, term.delta3, crossProduct, &dotDihedral, term.delta1, &signOfDihedral, 1));
+    }
+    
+    // Apply forces based on individual particle coordinates.
+    
+    for (int i = 0; i < (int) particleTerms.size(); i++) {
+        const ParticleTermInfo& term = particleTerms[i];
+        forces[permutedParticles[term.atom]][term.component] -= term.forceExpression.evaluate();
+    }
+
+    // Apply forces based on distances.
+
+    for (int i = 0; i < (int) distanceTerms.size(); i++) {
+        const DistanceTermInfo& term = distanceTerms[i];
+        RealOpenMM dEdR = (RealOpenMM) (term.forceExpression.evaluate()/(term.delta[ReferenceForce::RIndex]));
+        for (int i = 0; i < 3; i++) {
+           RealOpenMM force  = -dEdR*term.delta[i];
+           forces[permutedParticles[term.p1]][i] -= force;
+           forces[permutedParticles[term.p2]][i] += force;
+        }
+    }
+
+    // Apply forces based on angles.
+
+    for (int i = 0; i < (int) angleTerms.size(); i++) {
+        const AngleTermInfo& term = angleTerms[i];
+        RealOpenMM dEdTheta = (RealOpenMM) term.forceExpression.evaluate();
+        RealOpenMM thetaCross[ReferenceForce::LastDeltaRIndex];
+        SimTKOpenMMUtilities::crossProductVector3(term.delta1, term.delta2, thetaCross);
+        RealOpenMM lengthThetaCross = SQRT(DOT3(thetaCross, thetaCross));
+        if (lengthThetaCross < 1.0e-06)
+            lengthThetaCross = (RealOpenMM) 1.0e-06;
+        RealOpenMM termA = dEdTheta/(term.delta1[ReferenceForce::R2Index]*lengthThetaCross);
+        RealOpenMM termC = -dEdTheta/(term.delta2[ReferenceForce::R2Index]*lengthThetaCross);
+        RealOpenMM deltaCrossP[3][3];
+        SimTKOpenMMUtilities::crossProductVector3(term.delta1, thetaCross, deltaCrossP[0]);
+        SimTKOpenMMUtilities::crossProductVector3(term.delta2, thetaCross, deltaCrossP[2]);
+        for (int i = 0; i < 3; i++) {
+            deltaCrossP[0][i] *= termA;
+            deltaCrossP[2][i] *= termC;
+            deltaCrossP[1][i] = -(deltaCrossP[0][i]+deltaCrossP[2][i]);
+        }
+        for (int i = 0; i < 3; i++) {
+            forces[permutedParticles[term.p1]][i] += deltaCrossP[0][i];
+            forces[permutedParticles[term.p2]][i] += deltaCrossP[1][i];
+            forces[permutedParticles[term.p3]][i] += deltaCrossP[2][i];
+        }
+    }
+
+    // Apply forces based on dihedrals.
+
+    for (int i = 0; i < (int) dihedralTerms.size(); i++) {
+        const DihedralTermInfo& term = dihedralTerms[i];
+        RealOpenMM dEdTheta = (RealOpenMM) term.forceExpression.evaluate();
+        RealOpenMM internalF[4][3];
+        RealOpenMM forceFactors[4];
+        RealOpenMM normCross1 = DOT3(term.cross1, term.cross1);
+        RealOpenMM normBC = term.delta2[ReferenceForce::RIndex];
+        forceFactors[0] = (-dEdTheta*normBC)/normCross1;
+        RealOpenMM normCross2 = DOT3(term.cross2, term.cross2);
+                   forceFactors[3] = (dEdTheta*normBC)/normCross2;
+                   forceFactors[1] = DOT3(term.delta1, term.delta2);
+                   forceFactors[1] /= term.delta2[ReferenceForce::R2Index];
+                   forceFactors[2] = DOT3(term.delta3, term.delta2);
+                   forceFactors[2] /= term.delta2[ReferenceForce::R2Index];
+        for (int i = 0; i < 3; i++) {
+            internalF[0][i] = forceFactors[0]*term.cross1[i];
+            internalF[3][i] = forceFactors[3]*term.cross2[i];
+            RealOpenMM s = forceFactors[1]*internalF[0][i] - forceFactors[2]*internalF[3][i];
+            internalF[1][i] = internalF[0][i] - s;
+            internalF[2][i] = internalF[3][i] + s;
+        }
+        for (int i = 0; i < 3; i++) {
+            forces[permutedParticles[term.p1]][i] += internalF[0][i];
+            forces[permutedParticles[term.p2]][i] -= internalF[1][i];
+            forces[permutedParticles[term.p3]][i] -= internalF[2][i];
+            forces[permutedParticles[term.p4]][i] += internalF[3][i];
+        }
+    }
+
+    // Add the energy
+
+    if (totalEnergy)
+        *totalEnergy += (RealOpenMM) energyExpression.evaluate();
+}
+
+void CpuCustomManyParticleForce::computeDelta(int atom1, int atom2, RealOpenMM* delta, vector<RealVec>& atomCoordinates) const {
+    if (usePeriodic)
+        ReferenceForce::getDeltaRPeriodic(atomCoordinates[atom1], atomCoordinates[atom2], periodicBoxSize, delta);
+    else
+        ReferenceForce::getDeltaR(atomCoordinates[atom1], atomCoordinates[atom2], delta);
+}
+
+RealOpenMM CpuCustomManyParticleForce::computeAngle(RealOpenMM* vec1, RealOpenMM* vec2) {
+    RealOpenMM dot = DOT3(vec1, vec2);
+    RealOpenMM cosine = dot/SQRT((vec1[ReferenceForce::R2Index]*vec2[ReferenceForce::R2Index]));
+    RealOpenMM angle;
+    if (cosine >= 1)
+        angle = 0;
+    else if (cosine <= -1)
+        angle = PI_M;
+    else
+        angle = ACOS(cosine);
+    return angle;
+}
+
+void CpuCustomManyParticleForce::getDeltaR(const fvec4& posI, const fvec4& posJ, fvec4& deltaR, float& r2, const fvec4& boxSize, const fvec4& invBoxSize) const {
+    deltaR = posJ-posI;
+    if (usePeriodic) {
+        fvec4 base = round(deltaR*invBoxSize)*boxSize;
+        deltaR = deltaR-base;
+    }
+    r2 = dot3(deltaR, deltaR);
+}

platforms/cpu/src/CpuKernelFactory.cpp

@@ -49,6 +49,8 @@ KernelImpl* CpuKernelFactory::createKernelImpl(std::string name, const Platform&
         return new CpuCalcNonbondedForceKernel(name, platform, data);
     if (name == CalcCustomNonbondedForceKernel::Name())
         return new CpuCalcCustomNonbondedForceKernel(name, platform, data);
+    if (name == CalcCustomManyParticleForceKernel::Name())
+        return new CpuCalcCustomManyParticleForceKernel(name, platform, data);
     if (name == CalcGBSAOBCForceKernel::Name())
         return new CpuCalcGBSAOBCForceKernel(name, platform, data);
     if (name == IntegrateLangevinStepKernel::Name())

platforms/cpu/src/CpuKernels.cpp

@@ -835,6 +835,74 @@ void CpuCalcGBSAOBCForceKernel::copyParametersToContext(ContextImpl& context, co
     obc.setParticleParameters(particleParams);
 }
 
+CpuCalcCustomManyParticleForceKernel::~CpuCalcCustomManyParticleForceKernel() {
+    if (particleParamArray != NULL) {
+        for (int i = 0; i < numParticles; i++)
+            delete[] particleParamArray[i];
+        delete[] particleParamArray;
+    }
+    if (ixn != NULL)
+        delete ixn;
+}
+
+void CpuCalcCustomManyParticleForceKernel::initialize(const System& system, const CustomManyParticleForce& force) {
+
+    // Build the arrays.
+
+    numParticles = system.getNumParticles();
+    int numParticleParameters = force.getNumPerParticleParameters();
+    particleParamArray = new double*[numParticles];
+    for (int i = 0; i < numParticles; i++)
+        particleParamArray[i] = new double[numParticleParameters];
+    for (int i = 0; i < numParticles; ++i) {
+        vector<double> parameters;
+        int type;
+        force.getParticleParameters(i, parameters, type);
+        for (int j = 0; j < numParticleParameters; j++)
+            particleParamArray[i][j] = parameters[j];
+    }
+    for (int i = 0; i < force.getNumGlobalParameters(); i++)
+        globalParameterNames.push_back(force.getGlobalParameterName(i));
+    ixn = new CpuCustomManyParticleForce(force);
+    nonbondedMethod = CalcCustomManyParticleForceKernel::NonbondedMethod(force.getNonbondedMethod());
+    cutoffDistance = force.getCutoffDistance();
+}
+
+double CpuCalcCustomManyParticleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
+    vector<RealVec>& posData = extractPositions(context);
+    vector<RealVec>& forceData = extractForces(context);
+    RealOpenMM energy = 0;
+    map<string, double> globalParameters;
+    for (int i = 0; i < (int) globalParameterNames.size(); i++)
+        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
+    if (nonbondedMethod == CutoffPeriodic) {
+        RealVec& box = extractBoxSize(context);
+        double minAllowedSize = 2*cutoffDistance;
+        if (box[0] < minAllowedSize || box[1] < minAllowedSize || box[2] < minAllowedSize)
+            throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
+        ixn->setPeriodic(box);
+    }
+    ixn->calculateIxn(&data.posq[0], posData, particleParamArray, globalParameters, forceData, includeEnergy ? &energy : NULL);
+    return energy;
+}
+
+void CpuCalcCustomManyParticleForceKernel::copyParametersToContext(ContextImpl& context, const CustomManyParticleForce& force) {
+    if (numParticles != force.getNumParticles())
+        throw OpenMMException("updateParametersInContext: The number of particles has changed");
+
+    // Record the values.
+
+    int numParameters = force.getNumPerParticleParameters();
+    vector<double> params;
+    for (int i = 0; i < numParticles; ++i) {
+        vector<double> parameters;
+        int type;
+        force.getParticleParameters(i, parameters, type);
+        for (int j = 0; j < numParameters; j++)
+            particleParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
+    }
+}
+
 CpuIntegrateLangevinStepKernel::~CpuIntegrateLangevinStepKernel() {
     if (dynamics)
         delete dynamics;

platforms/cpu/src/CpuPlatform.cpp

@@ -65,6 +65,7 @@ CpuPlatform::CpuPlatform() {
     registerKernelFactory(CalcRBTorsionForceKernel::Name(), factory);
     registerKernelFactory(CalcNonbondedForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomNonbondedForceKernel::Name(), factory);
+    registerKernelFactory(CalcCustomManyParticleForceKernel::Name(), factory);
     registerKernelFactory(CalcGBSAOBCForceKernel::Name(), factory);
     registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);
     platformProperties.push_back(CpuThreads());