Began implementing CustomManyParticleForce

olla/include/openmm/kernels.h

@@ -44,6 +44,7 @@
 #include "openmm/CustomHbondForce.h"
 #include "openmm/CustomIntegrator.h"
 #include "openmm/CustomNonbondedForce.h"
+#include "openmm/CustomManyParticleForce.h"
 #include "openmm/CustomTorsionForce.h"
 #include "openmm/GBSAOBCForce.h"
 #include "openmm/GBVIForce.h"
@@ -826,6 +827,46 @@ public:
     virtual void copyParametersToContext(ContextImpl& context, const CustomCompoundBondForce& force) = 0;
 };
 
+/**
+ * This kernel is invoked by CustomManyParticleForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CalcCustomManyParticleForceKernel : public KernelImpl {
+public:
+    enum NonbondedMethod {
+        NoCutoff = 0,
+        CutoffNonPeriodic = 1,
+        CutoffPeriodic = 2
+    };
+    static std::string Name() {
+        return "CalcCustomManyParticleForce";
+    }
+    CalcCustomManyParticleForceKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
+    }
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CustomManyParticleForce this kernel will be used for
+     */
+    virtual void initialize(const System& system, const CustomManyParticleForce& force) = 0;
+    /**
+     * 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
+     */
+    virtual double execute(ContextImpl& context, bool includeForces, bool includeEnergy) = 0;
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CustomManyParticleForce to copy the parameters from
+     */
+    virtual void copyParametersToContext(ContextImpl& context, const CustomManyParticleForce& force) = 0;
+};
+
 /**
  * This kernel is invoked by VerletIntegrator to take one time step.
  */

openmmapi/include/openmm/CustomManyParticleForce.h

+#ifndef OPENMM_CUSTOMTHREEBODYFORCE_H_
+#define OPENMM_CUSTOMTHREEBODYFORCE_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) 2008-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 "Force.h"
+#include "TabulatedFunction.h"
+#include "internal/windowsExport.h"
+#include <set>
+#include <vector>
+
+namespace OpenMM {
+
+/**
+ * This class supports a wide variety of nonbonded N-particle interactions, where N is user specified.  The
+ * interaction energy is determined by an arbitrary, user specified algebraic expression that is evaluated for
+ * every possible set of N particles in the system.  It may depend on the positions of the individual particles,
+ * the distances between pairs of particles, the angles formed by sets of three particles, and the dihedral
+ * angles formed by sets of four particles.
+ * 
+ * Be aware that the cost of evaluating an N-particle interaction increases very rapidly with N.  Values larger
+ * than N=3 are rarely used.
+ * 
+ * We refer to a set of particles for which the energy is being evaluated  as p1, p2, p3, etc.  The energy expression
+ * may depend on the following variables and functions:
+ * 
+ * <ul>
+ * <li>x1, y1, z1, x2, y2, z2, etc.: The x, y, and z coordinates of the particle positions.  For example, x1
+ * is the x coordinate of particle p1, and y3 is the y coordinate of particle p3.</li>
+ * <li>distance(p1, p2): the distance between particles p1 and p2 (where "p1" and "p2" may be replaced by the names
+ * of whichever particles you want to calculate the distance between).</li>
+ * <li>angle(p1, p2, p3): the angle formed by the three specified particles.</li>
+ * <li>dihedral(p1, p2, p3, p4): the dihedral angle formed by the four specified particles.</li>
+ * <li>arbitrary global and per-particle parameters that you define.</li>
+ * </ul>
+ *
+ * To use this class, create a CustomManyParticleForce object, passing an algebraic expression to the constructor
+ * that defines the interaction energy of each set of particles.  Then call addPerParticleParameter() to define per-particle
+ * parameters, and addGlobalParameter() to define global parameters.  The values of per-particle parameters are specified as
+ * part of the system definition, while values of global parameters may be modified during a simulation by calling Context::setParameter().
+ * 
+ * Next, call addParticle() once for each particle in the System to set the values of its per-particle parameters.
+ * The number of particles for which you set parameters must be exactly equal to the number of particles in the
+ * System, or else an exception will be thrown when you try to create a Context.  After a particle has been added,
+ * you can modify its parameters by calling setParticleParameters().  This will have no effect on Contexts that already exist
+ * unless you call updateParametersInContext().
+ * 
+ * Multi-particle interactions can be very expensive to evaluate, so they are usually used with a cutoff distance.  If two particles
+ * are further apart than the cutoff, <i>all</i> sets that include those two particles will be omitted.
+ * 
+ * CustomManyParticleForce also lets you specify "exclusions", particular pairs of particles whose interactions should be
+ * omitted from force and energy calculations.  This is most often used for particles that are bonded to each other.
+ * If you specify a pair of particles as an exclusion, <i>all</i> sets that include those two particles will be omitted.
+ *
+ * As an example, the following code creates a CustomManyParticleForce that implements an Axilrod-Teller potential.  This
+ * is an interaction between three particles that depends on all three distances and angles formed by the particles.
+ *
+ * <tt>CustomManyParticleForce* force = new CustomManyParticleForce(3,
+ *     "C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
+ *     "theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
+ *     "r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
+ * </tt>
+ *
+ * This force depends on one parameter, C.  The following code defines it as a global parameter:
+ *
+ * <tt><pre>
+ * force->addGlobalParameter("C", 1.0);
+ * </pre></tt>
+ * 
+ * The expression <i>must</i> be symmetric with respect to the particles.  It typically will only be evaluated once for
+ * each set of particles, and no guarantee is made about which particle will be identified as "particle 1".  In the above
+ * example, the energy only depends on the products cos(theta1)*cos(theta2)*cos(theta3) and r12*r13*r23, both of which are
+ * unchanged if the labels p1, p2, and p3 are permuted.  If that were not true, the results would be undefined, because
+ * permuting the labels would change the energy.
+ * 
+ * In some cases this requirement is overly restrictive.  When some particles are fundamentally different from others,
+ * the expression may be inherently non-symmetric.  An example would be a water model that involves three particles,
+ * two of which <i>must</i> be hydrogen and one of which <i>must</i> be oxygen.  Cases like this can be implemented
+ * using particle types.
+ * 
+ * A particle type is an integer that you specify when you call addParticle().  (If you omit the argument, it defaults
+ * to 0.)  For the water model, you could specify 0 for all oxygen atoms and 1 for all hydrogen atoms.  You can then
+ * call setTypeFilter() to specify the list of allowed types for each of the N particles involved in an interaction:
+ * 
+ * <tt><pre>
+ * set&lt;int&gt; oxygenTypes, hydrogenTypes;
+ * oxygenTypes.insert(0);
+ * hydrogenTypes.insert(1);
+ * force->setTypeFilter(0, oxygenTypes);
+ * force->setTypeFilter(1, hydrogenTypes);
+ * force->setTypeFilter(2, hydrogenTypes);
+ * </pre></tt>
+ *
+ * This specifies that of the three particles in an interaction, p1 must be oxygen while p2 and p3 must be hydrogen.
+ * The energy expression will only be evaluated for triplets of particles that satisfy those requirements.  It will
+ * still only be evaluated once for each triplet, so it must still be symmetric with respect to p2 and p3.
+ * 
+ * Expressions may involve the operators + (add), - (subtract), * (multiply), / (divide), and ^ (power), and the following
+ * functions: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, erf, erfc, min, max, abs, step, delta.  All trigonometric functions
+ * are defined in radians, and log is the natural logarithm.  step(x) = 0 if x is less than 0, 1 otherwise.  delta(x) = 1 if x is 0, 0 otherwise.
+ * The names of per-particle parameters have the suffix "1", "2", etc. appended to them to indicate the values for the multiple interacting particles.
+ * For example, if you define a per-particle parameter called "charge", then the variable "charge2" is the charge of particle p2.
+ * As seen above, the expression may also involve intermediate quantities that are defined following the main expression, using ";" as a separator.
+ *
+ * In addition, you can call addTabulatedFunction() to define a new function based on tabulated values.  You specify the function by
+ * creating a TabulatedFunction object.  That function can then appear in the expression.
+ */
+
+class OPENMM_EXPORT CustomManyParticleForce : public Force {
+public:
+    /**
+     * This is an enumeration of the different methods that may be used for handling long range nonbonded forces.
+     */
+    enum NonbondedMethod {
+        /**
+         * No cutoff is applied to nonbonded interactions.  The full set of interactions is computed exactly.
+         * This necessarily means that periodic boundary conditions cannot be used.  This is the default.
+         */
+        NoCutoff = 0,
+        /**
+         * Interactions are ignored if any two particles are further apart than the cutoff distance.
+         */
+        CutoffNonPeriodic = 1,
+        /**
+         * Periodic boundary conditions are used, so that each particle interacts only with the nearest periodic copy of
+         * each other particle.  Interactions are ignored if any two particles are further apart than the cutoff distance.
+         */
+        CutoffPeriodic = 2,
+    };
+    /**
+     * Create a CustomManyParticleForce.
+     *
+     * @param particlesPerSet  the number of particles in each set for which the energy is evaluated
+     * @param energy           an algebraic expression giving the interaction energy of each triplet as a function
+     *                         of particle positions, inter-particle distances, angles, and any global and per-particle parameters
+     */
+    explicit CustomManyParticleForce(int particlesPerSet, const std::string& energy);
+    ~CustomManyParticleForce();
+    /**
+     * Get the number of particles in each set for which the energy is evaluated
+     */
+    int getNumParticlesPerSet() const {
+        return particlesPerSet;
+    }
+    /**
+     * Get the number of particles for which force field parameters have been defined.
+     */
+    int getNumParticles() const {
+        return particles.size();
+    }
+    /**
+     * Get the number of particle pairs whose interactions should be excluded.
+     */
+    int getNumExclusions() const {
+        return exclusions.size();
+    }
+    /**
+     * Get the number of per-particle parameters that the interaction depends on.
+     */
+    int getNumPerParticleParameters() const {
+        return particleParameters.size();
+    }
+    /**
+     * Get the number of global parameters that the interaction depends on.
+     */
+    int getNumGlobalParameters() const {
+        return globalParameters.size();
+    }
+    /**
+     * Get the number of tabulated functions that have been defined.
+     */
+    int getNumTabulatedFunctions() const {
+        return functions.size();
+    }
+    /**
+     * Get the algebraic expression that gives the interaction energy of each bond
+     */
+    const std::string& getEnergyFunction() const;
+    /**
+     * Set the algebraic expression that gives the interaction energy of each bond
+     */
+    void setEnergyFunction(const std::string& energy);
+    /**
+     * Get the method used for handling long range nonbonded interactions.
+     */
+    NonbondedMethod getNonbondedMethod() const;
+    /**
+     * Set the method used for handling long range nonbonded interactions.
+     */
+    void setNonbondedMethod(NonbondedMethod method);
+    /**
+     * Get the cutoff distance (in nm) being used for nonbonded interactions.  If the NonbondedMethod in use
+     * is NoCutoff, this value will have no effect.
+     *
+     * @return the cutoff distance, measured in nm
+     */
+    double getCutoffDistance() const;
+    /**
+     * Set the cutoff distance (in nm) being used for nonbonded interactions.  If the NonbondedMethod in use
+     * is NoCutoff, this value will have no effect.
+     *
+     * @param distance    the cutoff distance, measured in nm
+     */
+    void setCutoffDistance(double distance);
+    /**
+     * Add a new per-particle parameter that the interaction may depend on.
+     *
+     * @param name     the name of the parameter
+     * @return the index of the parameter that was added
+     */
+    int addPerParticleParameter(const std::string& name);
+    /**
+     * Get the name of a per-particle parameter.
+     *
+     * @param index     the index of the parameter for which to get the name
+     * @return the parameter name
+     */
+    const std::string& getPerParticleParameterName(int index) const;
+    /**
+     * Set the name of a per-particle parameter.
+     *
+     * @param index          the index of the parameter for which to set the name
+     * @param name           the name of the parameter
+     */
+    void setPerParticleParameterName(int index, const std::string& name);
+    /**
+     * Add a new global parameter that the interaction may depend on.
+     *
+     * @param name             the name of the parameter
+     * @param defaultValue     the default value of the parameter
+     * @return the index of the parameter that was added
+     */
+    int addGlobalParameter(const std::string& name, double defaultValue);
+    /**
+     * Get the name of a global parameter.
+     *
+     * @param index     the index of the parameter for which to get the name
+     * @return the parameter name
+     */
+    const std::string& getGlobalParameterName(int index) const;
+    /**
+     * Set the name of a global parameter.
+     *
+     * @param index          the index of the parameter for which to set the name
+     * @param name           the name of the parameter
+     */
+    void setGlobalParameterName(int index, const std::string& name);
+    /**
+     * Get the default value of a global parameter.
+     *
+     * @param index     the index of the parameter for which to get the default value
+     * @return the parameter default value
+     */
+    double getGlobalParameterDefaultValue(int index) const;
+    /**
+     * Set the default value of a global parameter.
+     *
+     * @param index          the index of the parameter for which to set the default value
+     * @param name           the default value of the parameter
+     */
+    void setGlobalParameterDefaultValue(int index, double defaultValue);
+    /**
+     * Add the nonbonded force parameters for a particle.  This should be called once for each particle
+     * in the System.  When it is called for the i'th time, it specifies the parameters for the i'th particle.
+     *
+     * @param parameters    the list of parameters for the new particle
+     * @param type          the type of the new particle
+     * @return the index of the particle that was added
+     */
+    int addParticle(const std::vector<double>& parameters, int type=0);
+    /**
+     * Get the nonbonded force parameters for a particle.
+     *
+     * @param index       the index of the particle for which to get parameters
+     * @param parameters  the list of parameters for the specified particle
+     * @param type        the type of the specified particle
+     */
+    void getParticleParameters(int index, std::vector<double>& parameters, int& type) const;
+    /**
+     * Set the nonbonded force parameters for a particle.
+     *
+     * @param index       the index of the particle for which to set parameters
+     * @param parameters  the list of parameters for the specified particle
+     * @param type        the type of the specified particle
+     */
+    void setParticleParameters(int index, const std::vector<double>& parameters, int type);
+    /**
+     * 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
+     * @return the index of the exclusion that was added
+     */
+    int addExclusion(int particle1, int particle2);
+    /**
+     * Get the particles in a pair whose interaction should be excluded.
+     *
+     * @param index      the index of the exclusion for which to get particle indices
+     * @param particle1  the index of the first particle in the pair
+     * @param particle2  the index of the second particle in the pair
+     */
+    void getExclusionParticles(int index, int& particle1, int& particle2) const;
+    /**
+     * Set the particles in a pair whose interaction should be excluded.
+     *
+     * @param index      the index of the exclusion for which to set particle indices
+     * @param particle1  the index of the first particle in the pair
+     * @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);
+    /**
+     * Get the allowed particle types for one of the particles involved in the interaction.
+     * If this an empty set (the default), no filter is applied and all interactions are evaluated
+     * regardless of the type of the specified particle.
+     * 
+     * @param index    the index of the particle within the interaction (between 0 and getNumParticlesPerSet())
+     * @param types    the allowed types for the specified particle
+     */
+    void getTypeFilter(int index, std::set<int>& types) const;
+    /**
+     * Set the allowed particle types for one of the particles involved in the interaction.
+     * If this an empty set (the default), no filter is applied and all interactions are evaluated
+     * regardless of the type of the specified particle.
+     * 
+     * @param index    the index of the particle within the interaction (between 0 and getNumParticlesPerSet())
+     * @param types    the allowed types for the specified particle
+     */
+    void setTypeFilter(int index, const std::set<int>& types);
+    /**
+     * Add a tabulated function that may appear in the energy expression.
+     *
+     * @param name           the name of the function as it appears in expressions
+     * @param function       a TabulatedFunction object defining the function.  The TabulatedFunction
+     *                       should have been created on the heap with the "new" operator.  The
+     *                       Force takes over ownership of it, and deletes it when the Force itself is deleted.
+     * @return the index of the function that was added
+     */
+    int addTabulatedFunction(const std::string& name, TabulatedFunction* function);
+    /**
+     * Get a const reference to a tabulated function that may appear in the energy expression.
+     *
+     * @param index     the index of the function to get
+     * @return the TabulatedFunction object defining the function
+     */
+    const TabulatedFunction& getTabulatedFunction(int index) const;
+    /**
+     * Get a reference to a tabulated function that may appear in the energy expression.
+     *
+     * @param index     the index of the function to get
+     * @return the TabulatedFunction object defining the function
+     */
+    TabulatedFunction& getTabulatedFunction(int index);
+    /**
+     * Get the name of a tabulated function that may appear in the energy expression.
+     *
+     * @param index     the index of the function to get
+     * @return the name of the function as it appears in expressions
+     */
+    const std::string& getTabulatedFunctionName(int index) const;
+    /**
+     * Update the per-particle parameters in a Context to match those stored in this Force object.  This method provides
+     * an efficient method to update certain parameters in an existing Context without needing to reinitialize it.
+     * Simply call setParticleParameters() to modify this object's parameters, then call updateParametersInContext()
+     * to copy them over to the Context.
+     * 
+     * This method has several limitations.  The only information it updates is the values of per-particle parameters.
+     * All other aspects of the Force (the energy function, nonbonded method, cutoff distance, etc.) are unaffected and can
+     * only be changed by reinitializing the Context.  Also, this method cannot be used to add new particles, only to change
+     * the parameters of existing ones.
+     */
+    void updateParametersInContext(Context& context);
+protected:
+    ForceImpl* createImpl() const;
+private:
+    class ParticleInfo;
+    class ParticleParameterInfo;
+    class GlobalParameterInfo;
+    class ExclusionInfo;
+    class FunctionInfo;
+    int particlesPerSet;
+    NonbondedMethod nonbondedMethod;
+    double cutoffDistance;
+    std::string energyExpression;
+    std::vector<ParticleParameterInfo> particleParameters;
+    std::vector<GlobalParameterInfo> globalParameters;
+    std::vector<ParticleInfo> particles;
+    std::vector<ExclusionInfo> exclusions;
+    std::vector<FunctionInfo> functions;
+    std::vector<std::set<int> > typeFilters;
+};
+
+
+/**
+ * This is an internal class used to record information about a particle.
+ * @private
+ */
+class CustomManyParticleForce::ParticleInfo {
+public:
+    std::vector<double> parameters;
+    int type;
+    ParticleInfo() {
+    }
+    ParticleInfo(const std::vector<double>& parameters, int type) : parameters(parameters), type(type) {
+    }
+};
+
+/**
+ * This is an internal class used to record information about a per-particle parameter.
+ * @private
+ */
+class CustomManyParticleForce::ParticleParameterInfo {
+public:
+    std::string name;
+    ParticleParameterInfo() {
+    }
+    ParticleParameterInfo(const std::string& name) : name(name) {
+    }
+};
+
+/**
+ * This is an internal class used to record information about a global parameter.
+ * @private
+ */
+class CustomManyParticleForce::GlobalParameterInfo {
+public:
+    std::string name;
+    double defaultValue;
+    GlobalParameterInfo() {
+    }
+    GlobalParameterInfo(const std::string& name, double defaultValue) : name(name), defaultValue(defaultValue) {
+    }
+};
+
+/**
+ * This is an internal class used to record information about an exclusion.
+ * @private
+ */
+class CustomManyParticleForce::ExclusionInfo {
+public:
+    int particle1, particle2;
+    ExclusionInfo() {
+        particle1 = particle2 = -1;
+    }
+    ExclusionInfo(int particle1, int particle2) :
+        particle1(particle1), particle2(particle2) {
+    }
+};
+
+/**
+ * This is an internal class used to record information about a tabulated function.
+ * @private
+ */
+class CustomManyParticleForce::FunctionInfo {
+public:
+    std::string name;
+    TabulatedFunction* function;
+    FunctionInfo() {
+    }
+    FunctionInfo(const std::string& name, TabulatedFunction* function) : name(name), function(function) {
+    }
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_CUSTOMTHREEBODYFORCE_H_*/

openmmapi/include/openmm/internal/CustomManyParticleForceImpl.h

+#ifndef OPENMM_CUSTOMMANYPARTICLEFORCEIMPL_H_
+#define OPENMM_CUSTOMMANYPARTICLEFORCEIMPL_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) 2008-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 "ForceImpl.h"
+#include "openmm/CustomManyParticleForce.h"
+#include "openmm/Kernel.h"
+#include "lepton/CustomFunction.h"
+#include "lepton/ExpressionTreeNode.h"
+#include "lepton/ParsedExpression.h"
+#include <utility>
+#include <map>
+#include <string>
+
+namespace OpenMM {
+
+/**
+ * This is the internal implementation of CustomManyParticleForce.
+ */
+
+class OPENMM_EXPORT CustomManyParticleForceImpl : public ForceImpl {
+public:
+    CustomManyParticleForceImpl(const CustomManyParticleForce& owner);
+    ~CustomManyParticleForceImpl();
+    void initialize(ContextImpl& context);
+    const CustomManyParticleForce& getOwner() const {
+        return owner;
+    }
+    void updateContextState(ContextImpl& context) {
+        // This force field doesn't update the state directly.
+    }
+    double calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups);
+    std::map<std::string, double> getDefaultParameters();
+    std::vector<std::string> getKernelNames();
+    void updateParametersInContext(ContextImpl& context);
+    /**
+     * This is a utility routine that parses the energy expression, identifies the angles and dihedrals
+     * in it, and replaces them with variables.
+     *
+     * @param force     the CustomManyParticleForce to process
+     * @param functions definitions of custom function that may appear in the expression
+     * @param distances on exit, this will contain an entry for each distance used in the expression.  The key is the name
+     *                  of the corresponding variable, and the value is the list of particle indices.
+     * @param angles    on exit, this will contain an entry for each angle used in the expression.  The key is the name
+     *                  of the corresponding variable, and the value is the list of particle indices.
+     * @param dihedrals on exit, this will contain an entry for each dihedral used in the expression.  The key is the name
+     *                  of the corresponding variable, and the value is the list of particle indices.
+     * @return a Parsed expression for the energy
+     */
+    static Lepton::ParsedExpression prepareExpression(const CustomManyParticleForce& force, const std::map<std::string, Lepton::CustomFunction*>& functions, std::map<std::string, std::vector<int> >& distances,
+            std::map<std::string, std::vector<int> >& angles, std::map<std::string, std::vector<int> >& dihedrals);
+private:
+    class FunctionPlaceholder;
+    static Lepton::ExpressionTreeNode replaceFunctions(const Lepton::ExpressionTreeNode& node, std::map<std::string, int> atoms,
+            std::map<std::string, std::vector<int> >& distances, std::map<std::string, std::vector<int> >& angles,
+            std::map<std::string, std::vector<int> >& dihedrals);
+    const CustomManyParticleForce& owner;
+    Kernel kernel;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_CUSTOMMANYPARTICLEFORCEIMPL_H_*/

openmmapi/src/CustomManyParticleForce.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   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) 2008-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 "openmm/Force.h"
+#include "openmm/OpenMMException.h"
+#include "openmm/CustomManyParticleForce.h"
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/internal/CustomManyParticleForceImpl.h"
+#include <cmath>
+#include <map>
+#include <set>
+#include <sstream>
+#include <utility>
+
+using namespace OpenMM;
+using namespace std;
+
+CustomManyParticleForce::CustomManyParticleForce(int particlesPerSet, const string& energy) :
+        particlesPerSet(particlesPerSet), energyExpression(energy), typeFilters(particlesPerSet) {
+}
+
+CustomManyParticleForce::~CustomManyParticleForce() {
+    for (int i = 0; i < (int) functions.size(); i++)
+        delete functions[i].function;
+}
+
+const string& CustomManyParticleForce::getEnergyFunction() const {
+    return energyExpression;
+}
+
+void CustomManyParticleForce::setEnergyFunction(const string& energy) {
+    energyExpression = energy;
+}
+
+CustomManyParticleForce::NonbondedMethod CustomManyParticleForce::getNonbondedMethod() const {
+    return nonbondedMethod;
+}
+
+void CustomManyParticleForce::setNonbondedMethod(NonbondedMethod method) {
+    nonbondedMethod = method;
+}
+
+double CustomManyParticleForce::getCutoffDistance() const {
+    return cutoffDistance;
+}
+
+void CustomManyParticleForce::setCutoffDistance(double distance) {
+    cutoffDistance = distance;
+}
+
+int CustomManyParticleForce::addPerParticleParameter(const string& name) {
+    particleParameters.push_back(ParticleParameterInfo(name));
+    return particleParameters.size()-1;
+}
+
+const string& CustomManyParticleForce::getPerParticleParameterName(int index) const {
+    ASSERT_VALID_INDEX(index, particleParameters);
+    return particleParameters[index].name;
+}
+
+void CustomManyParticleForce::setPerParticleParameterName(int index, const string& name) {
+    ASSERT_VALID_INDEX(index, particleParameters);
+    particleParameters[index].name = name;
+}
+
+int CustomManyParticleForce::addGlobalParameter(const string& name, double defaultValue) {
+    globalParameters.push_back(GlobalParameterInfo(name, defaultValue));
+    return globalParameters.size()-1;
+}
+
+const string& CustomManyParticleForce::getGlobalParameterName(int index) const {
+    ASSERT_VALID_INDEX(index, globalParameters);
+    return globalParameters[index].name;
+}
+
+void CustomManyParticleForce::setGlobalParameterName(int index, const string& name) {
+    ASSERT_VALID_INDEX(index, globalParameters);
+    globalParameters[index].name = name;
+}
+
+double CustomManyParticleForce::getGlobalParameterDefaultValue(int index) const {
+    ASSERT_VALID_INDEX(index, globalParameters);
+    return globalParameters[index].defaultValue;
+}
+
+void CustomManyParticleForce::setGlobalParameterDefaultValue(int index, double defaultValue) {
+    ASSERT_VALID_INDEX(index, globalParameters);
+    globalParameters[index].defaultValue = defaultValue;
+}
+
+int CustomManyParticleForce::addParticle(const vector<double>& parameters, int type) {
+    particles.push_back(ParticleInfo(parameters, type));
+    return particles.size()-1;
+}
+
+void CustomManyParticleForce::getParticleParameters(int index, vector<double>& parameters, int& type) const {
+    ASSERT_VALID_INDEX(index, particles);
+    parameters = particles[index].parameters;
+    type = particles[index].type;
+}
+
+void CustomManyParticleForce::setParticleParameters(int index, const vector<double>& parameters, int type) {
+    ASSERT_VALID_INDEX(index, particles);
+    particles[index].parameters = parameters;
+    particles[index].type = type;
+}
+
+int CustomManyParticleForce::addExclusion(int particle1, int particle2) {
+    exclusions.push_back(ExclusionInfo(particle1, particle2));
+    return exclusions.size()-1;
+}
+void CustomManyParticleForce::getExclusionParticles(int index, int& particle1, int& particle2) const {
+    ASSERT_VALID_INDEX(index, exclusions);
+    particle1 = exclusions[index].particle1;
+    particle2 = exclusions[index].particle2;
+}
+
+void CustomManyParticleForce::setExclusionParticles(int index, int particle1, int particle2) {
+    ASSERT_VALID_INDEX(index, exclusions);
+    exclusions[index].particle1 = particle1;
+    exclusions[index].particle2 = particle2;
+}
+
+void CustomManyParticleForce::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);
+}
+
+void CustomManyParticleForce::getTypeFilter(int index, set<int>& types) const {
+    if (index < 0 || index >= particlesPerSet)
+        throw OpenMMException("CustomManyParticleForce: index to getTypeFilter out of range");
+    types = typeFilters[index];
+}
+
+void CustomManyParticleForce::setTypeFilter(int index, const set<int>& types) {
+    if (index < 0 || index >= particlesPerSet)
+        throw OpenMMException("CustomManyParticleForce: index to setTypeFilter out of range");
+    typeFilters[index] = types;
+}
+
+int CustomManyParticleForce::addTabulatedFunction(const string& name, TabulatedFunction* function) {
+    functions.push_back(FunctionInfo(name, function));
+    return functions.size()-1;
+}
+
+const TabulatedFunction& CustomManyParticleForce::getTabulatedFunction(int index) const {
+    ASSERT_VALID_INDEX(index, functions);
+    return *functions[index].function;
+}
+
+TabulatedFunction& CustomManyParticleForce::getTabulatedFunction(int index) {
+    ASSERT_VALID_INDEX(index, functions);
+    return *functions[index].function;
+}
+
+const string& CustomManyParticleForce::getTabulatedFunctionName(int index) const {
+    ASSERT_VALID_INDEX(index, functions);
+    return functions[index].name;
+}
+
+ForceImpl* CustomManyParticleForce::createImpl() const {
+    return new CustomManyParticleForceImpl(*this);
+}
+
+void CustomManyParticleForce::updateParametersInContext(Context& context) {
+    dynamic_cast<CustomManyParticleForceImpl&>(getImplInContext(context)).updateParametersInContext(getContextImpl(context));
+}

openmmapi/src/CustomManyParticleForceImpl.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   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) 2008-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 "openmm/OpenMMException.h"
+#include "openmm/internal/ContextImpl.h"
+#include "openmm/internal/CustomManyParticleForceImpl.h"
+#include "openmm/kernels.h"
+#include "lepton/Operation.h"
+#include "lepton/Parser.h"
+#include <sstream>
+
+using namespace OpenMM;
+using Lepton::CustomFunction;
+using Lepton::ExpressionTreeNode;
+using Lepton::Operation;
+using Lepton::ParsedExpression;
+using std::map;
+using std::pair;
+using std::vector;
+using std::set;
+using std::string;
+using std::stringstream;
+
+/**
+ * This class serves as a placeholder for angles and dihedrals in expressions.
+ */
+class CustomManyParticleForceImpl::FunctionPlaceholder : public CustomFunction {
+public:
+    int numArguments;
+    FunctionPlaceholder(int numArguments) : numArguments(numArguments) {
+    }
+    int getNumArguments() const {
+        return numArguments;
+    }
+    double evaluate(const double* arguments) const {
+        return 0.0;
+    }
+    double evaluateDerivative(const double* arguments, const int* derivOrder) const {
+        return 0.0;
+    }
+    CustomFunction* clone() const {
+        return new FunctionPlaceholder(numArguments);
+    }
+};
+
+CustomManyParticleForceImpl::CustomManyParticleForceImpl(const CustomManyParticleForce& owner) : owner(owner) {
+}
+
+CustomManyParticleForceImpl::~CustomManyParticleForceImpl() {
+}
+
+void CustomManyParticleForceImpl::initialize(ContextImpl& context) {
+    kernel = context.getPlatform().createKernel(CalcCustomManyParticleForceKernel::Name(), context);
+
+    // Check for errors in the specification of parameters and exclusions.
+
+    const System& system = context.getSystem();
+    if (owner.getNumParticles() != system.getNumParticles())
+        throw OpenMMException("CustomManyParticleForce must have exactly as many particles as the System it belongs to.");
+    vector<set<int> > exclusions(owner.getNumParticles());
+    vector<double> parameters;
+    int type;
+    int numParameters = owner.getNumPerParticleParameters();
+    for (int i = 0; i < owner.getNumParticles(); i++) {
+        owner.getParticleParameters(i, parameters, type);
+        if (parameters.size() != numParameters) {
+            stringstream msg;
+            msg << "CustomManyParticleForce: Wrong number of parameters for particle ";
+            msg << i;
+            throw OpenMMException(msg.str());
+        }
+    }
+    for (int i = 0; i < owner.getNumExclusions(); i++) {
+        int particle1, particle2;
+        owner.getExclusionParticles(i, particle1, particle2);
+        if (particle1 < 0 || particle1 >= owner.getNumParticles()) {
+            stringstream msg;
+            msg << "CustomManyParticleForce: Illegal particle index for an exclusion: ";
+            msg << particle1;
+            throw OpenMMException(msg.str());
+        }
+        if (particle2 < 0 || particle2 >= owner.getNumParticles()) {
+            stringstream msg;
+            msg << "CustomManyParticleForce: Illegal particle index for an exclusion: ";
+            msg << particle2;
+            throw OpenMMException(msg.str());
+        }
+        if (exclusions[particle1].count(particle2) > 0 || exclusions[particle2].count(particle1) > 0) {
+            stringstream msg;
+            msg << "CustomManyParticleForce: Multiple exclusions are specified for particles ";
+            msg << particle1;
+            msg << " and ";
+            msg << particle2;
+            throw OpenMMException(msg.str());
+        }
+        exclusions[particle1].insert(particle2);
+        exclusions[particle2].insert(particle1);
+    }
+    if (owner.getNonbondedMethod() == CustomManyParticleForce::CutoffPeriodic) {
+        Vec3 boxVectors[3];
+        system.getDefaultPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
+        double cutoff = owner.getCutoffDistance();
+        if (cutoff > 0.5*boxVectors[0][0] || cutoff > 0.5*boxVectors[1][1] || cutoff > 0.5*boxVectors[2][2])
+            throw OpenMMException("CustomManyParticleForce: The cutoff distance cannot be greater than half the periodic box size.");
+    }
+    kernel.getAs<CalcCustomManyParticleForceKernel>().initialize(context.getSystem(), owner);
+}
+
+double CustomManyParticleForceImpl::calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups) {
+    if ((groups&(1<<owner.getForceGroup())) != 0)
+        return kernel.getAs<CalcCustomManyParticleForceKernel>().execute(context, includeForces, includeEnergy);
+    return 0.0;
+}
+
+vector<string> CustomManyParticleForceImpl::getKernelNames() {
+    vector<string> names;
+    names.push_back(CalcCustomManyParticleForceKernel::Name());
+    return names;
+}
+
+map<string, double> CustomManyParticleForceImpl::getDefaultParameters() {
+    map<string, double> parameters;
+    for (int i = 0; i < owner.getNumGlobalParameters(); i++)
+        parameters[owner.getGlobalParameterName(i)] = owner.getGlobalParameterDefaultValue(i);
+    return parameters;
+}
+
+ParsedExpression CustomManyParticleForceImpl::prepareExpression(const CustomManyParticleForce& force, const map<string, CustomFunction*>& customFunctions, map<string, vector<int> >& distances,
+        map<string, vector<int> >& angles, map<string, vector<int> >& dihedrals) {
+    CustomManyParticleForceImpl::FunctionPlaceholder custom(1);
+    CustomManyParticleForceImpl::FunctionPlaceholder distance(2);
+    CustomManyParticleForceImpl::FunctionPlaceholder angle(3);
+    CustomManyParticleForceImpl::FunctionPlaceholder dihedral(4);
+    map<string, CustomFunction*> functions = customFunctions;
+    functions["distance"] = &distance;
+    functions["angle"] = &angle;
+    functions["dihedral"] = &dihedral;
+    ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction(), functions);
+    map<string, int> atoms;
+    for (int i = 0; i < force.getNumParticlesPerSet(); i++) {
+        stringstream name;
+        name << 'p' << (i+1);
+        atoms[name.str()] = i;
+    }
+    return ParsedExpression(replaceFunctions(expression.getRootNode(), atoms, distances, angles, dihedrals)).optimize();
+}
+
+ExpressionTreeNode CustomManyParticleForceImpl::replaceFunctions(const ExpressionTreeNode& node, map<string, int> atoms,
+        map<string, vector<int> >& distances, map<string, vector<int> >& angles, map<string, vector<int> >& dihedrals) {
+    const Operation& op = node.getOperation();
+    if (op.getId() != Operation::CUSTOM || (op.getName() != "distance" && op.getName() != "angle" && op.getName() != "dihedral"))
+    {
+        // This is not an angle or dihedral, so process its children.
+
+        vector<ExpressionTreeNode> children;
+        for (int i = 0; i < (int) node.getChildren().size(); i++)
+            children.push_back(replaceFunctions(node.getChildren()[i], atoms, distances, angles, dihedrals));
+        return ExpressionTreeNode(op.clone(), children);
+    }
+    const Operation::Custom& custom = static_cast<const Operation::Custom&>(op);
+
+    // Identify the atoms this term is based on.
+
+    int numArgs = custom.getNumArguments();
+    vector<int> indices(numArgs);
+    for (int i = 0; i < numArgs; i++) {
+        map<string, int>::const_iterator iter = atoms.find(node.getChildren()[i].getOperation().getName());
+        if (iter == atoms.end())
+            throw OpenMMException("CustomManyParticleForce: Unknown particle '"+node.getChildren()[i].getOperation().getName()+"'");
+        indices[i] = iter->second;
+    }
+    
+    // Select a name for the variable and add it to the appropriate map.
+    
+    stringstream variable;
+    if (numArgs == 2)
+        variable << "distance";
+    else if (numArgs == 3)
+        variable << "angle";
+    else
+        variable << "dihedral";
+    for (int i = 0; i < numArgs; i++)
+        variable << indices[i];
+    string name = variable.str();
+    if (numArgs == 2)
+        distances[name] = indices;
+    else if (numArgs == 3)
+        angles[name] = indices;
+    else
+        dihedrals[name] = indices;
+    
+    // Return a new node that represents it as a simple variable.
+    
+    return ExpressionTreeNode(new Operation::Variable(name));
+}
+
+void CustomManyParticleForceImpl::updateParametersInContext(ContextImpl& context) {
+    kernel.getAs<CalcCustomManyParticleForceKernel>().copyParametersToContext(context, owner);
+}

platforms/reference/include/ReferenceBondIxn.h

@@ -82,7 +82,7 @@ class OPENMM_EXPORT ReferenceBondIxn {
       
          --------------------------------------------------------------------------------------- */
       
-      RealOpenMM getNormedDotProduct( RealOpenMM* vector1, RealOpenMM* vector2, int hasREntry ) const;
+      static RealOpenMM getNormedDotProduct( RealOpenMM* vector1, RealOpenMM* vector2, int hasREntry );
       
       /**---------------------------------------------------------------------------------------
       
@@ -98,8 +98,8 @@ class OPENMM_EXPORT ReferenceBondIxn {
       
          --------------------------------------------------------------------------------------- */
       
-      RealOpenMM getAngleBetweenTwoVectors( RealOpenMM* vector1, RealOpenMM* vector2, 
-                                            RealOpenMM* outputDotProduct, int hasREntry ) const;
+      static RealOpenMM getAngleBetweenTwoVectors( RealOpenMM* vector1, RealOpenMM* vector2, 
+                                                   RealOpenMM* outputDotProduct, int hasREntry );
       
       /**---------------------------------------------------------------------------------------
       
@@ -119,10 +119,10 @@ class OPENMM_EXPORT ReferenceBondIxn {
       
          --------------------------------------------------------------------------------------- */
       
-      RealOpenMM getDihedralAngleBetweenThreeVectors( RealOpenMM* vector1, RealOpenMM* vector2, 
-                                                      RealOpenMM* vector3, RealOpenMM** outputCrossProduct, 
-                                                      RealOpenMM* cosineOfAngle, RealOpenMM* signVector, 
-                                                      RealOpenMM* signOfAngle, int hasREntry ) const;
+      static RealOpenMM getDihedralAngleBetweenThreeVectors( RealOpenMM* vector1, RealOpenMM* vector2, 
+                                                             RealOpenMM* vector3, RealOpenMM** outputCrossProduct, 
+                                                             RealOpenMM* cosineOfAngle, RealOpenMM* signVector, 
+                                                             RealOpenMM* signOfAngle, int hasREntry );
       
 };
 

platforms/reference/include/ReferenceCustomManyParticleIxn.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 __ReferenceCustomManyParticleIxn_H__
+#define __ReferenceCustomManyParticleIxn_H__
+
+#include "ReferenceBondIxn.h"
+#include "lepton/ExpressionProgram.h"
+#include "lepton/ParsedExpression.h"
+#include <map>
+#include <vector>
+
+// ---------------------------------------------------------------------------------------
+
+class ReferenceCustomManyParticleIxn {
+
+   private:
+
+      class ParticleTermInfo;
+      class DistanceTermInfo;
+      class AngleTermInfo;
+      class DihedralTermInfo;
+      int numParticlesPerSet, numPerParticleParameters;
+      bool useCutoff, usePeriodic;
+      RealOpenMM cutoffDistance;
+      RealOpenMM periodicBoxSize[3];
+      Lepton::ExpressionProgram energyExpression;
+      std::vector<std::vector<std::string> > particleParamNames;
+      std::vector<ParticleTermInfo> particleTerms;
+      std::vector<DistanceTermInfo> distanceTerms;
+      std::vector<AngleTermInfo> angleTerms;
+      std::vector<DihedralTermInfo> dihedralTerms;
+
+      void loopOverInteractions(std::vector<int>& particles, int loopIndex, std::vector<OpenMM::RealVec>& atomCoordinates,
+                                RealOpenMM** particleParameters, std::map<std::string, double>& variables, std::vector<OpenMM::RealVec>& forces,
+                                RealOpenMM* totalEnergy) const;
+
+      /**---------------------------------------------------------------------------------------
+
+         Calculate custom interaction for one set of particles
+
+         @param particles        the indices of the particles
+         @param atomCoordinates  atom coordinates
+         @param variables        the values of variables that may appear in expressions
+         @param forces           force array (forces added)
+         @param totalEnergy      total energy
+
+         --------------------------------------------------------------------------------------- */
+
+      void calculateOneIxn(const std::vector<int>& particles, std::vector<OpenMM::RealVec>& atomCoordinates,
+                           std::map<std::string, double>& variables, std::vector<OpenMM::RealVec>& forces,
+                           RealOpenMM* totalEnergy) const;
+
+      void computeDelta(int atom1, int atom2, RealOpenMM* delta, std::vector<OpenMM::RealVec>& atomCoordinates) const;
+
+      static RealOpenMM computeAngle(RealOpenMM* vec1, RealOpenMM* vec2);
+
+
+   public:
+
+      /**---------------------------------------------------------------------------------------
+
+         Constructor
+
+         --------------------------------------------------------------------------------------- */
+
+       ReferenceCustomManyParticleIxn(int numParticlesPerSet, const Lepton::ParsedExpression& energyExpression,
+                               const std::vector<std::string>& particleParameterNames, const std::map<std::string, std::vector<int> >& distances,
+                               const std::map<std::string, std::vector<int> >& angles, const std::map<std::string, std::vector<int> >& dihedrals);
+
+      /**---------------------------------------------------------------------------------------
+
+         Destructor
+
+         --------------------------------------------------------------------------------------- */
+
+       ~ReferenceCustomManyParticleIxn();
+
+      /**---------------------------------------------------------------------------------------
+
+         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 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(std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** particleParameters,
+                        const std::map<std::string, double>& globalParameters,
+                        std::vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy) const;
+
+// ---------------------------------------------------------------------------------------
+
+};
+
+class ReferenceCustomManyParticleIxn::ParticleTermInfo {
+public:
+    std::string name;
+    int atom, component;
+    Lepton::ExpressionProgram forceExpression;
+    ParticleTermInfo(const std::string& name, int atom, int component, const Lepton::ExpressionProgram& forceExpression) :
+            name(name), atom(atom), component(component), forceExpression(forceExpression) {
+    }
+};
+
+class ReferenceCustomManyParticleIxn::DistanceTermInfo {
+public:
+    std::string name;
+    int p1, p2;
+    Lepton::ExpressionProgram forceExpression;
+    mutable RealOpenMM delta[ReferenceForce::LastDeltaRIndex];
+    DistanceTermInfo(const std::string& name, const std::vector<int>& atoms, const Lepton::ExpressionProgram& forceExpression) :
+            name(name), p1(atoms[0]), p2(atoms[1]), forceExpression(forceExpression) {
+    }
+};
+
+class ReferenceCustomManyParticleIxn::AngleTermInfo {
+public:
+    std::string name;
+    int p1, p2, p3;
+    Lepton::ExpressionProgram forceExpression;
+    mutable RealOpenMM delta1[ReferenceForce::LastDeltaRIndex];
+    mutable RealOpenMM delta2[ReferenceForce::LastDeltaRIndex];
+    AngleTermInfo(const std::string& name, const std::vector<int>& atoms, const Lepton::ExpressionProgram& forceExpression) :
+            name(name), p1(atoms[0]), p2(atoms[1]), p3(atoms[2]), forceExpression(forceExpression) {
+    }
+};
+
+class ReferenceCustomManyParticleIxn::DihedralTermInfo {
+public:
+    std::string name;
+    int p1, p2, p3, p4;
+    Lepton::ExpressionProgram 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::ExpressionProgram& forceExpression) :
+            name(name), p1(atoms[0]), p2(atoms[1]), p3(atoms[2]), p4(atoms[3]), forceExpression(forceExpression) {
+    }
+};
+
+#endif // __ReferenceCustomManyParticleIxn_H__

platforms/reference/include/ReferenceKernels.h

@@ -44,6 +44,7 @@ class CpuGBVI;
 class ReferenceAndersenThermostat;
 class ReferenceCustomCompoundBondIxn;
 class ReferenceCustomHbondIxn;
+class ReferenceCustomManyParticleIxn;
 class ReferenceBrownianDynamics;
 class ReferenceStochasticDynamics;
 class ReferenceConstraintAlgorithm;
@@ -861,6 +862,46 @@ private:
     std::vector<std::string> globalParameterNames;
 };
 
+/**
+ * This kernel is invoked by CustomManyParticleForce to calculate the forces acting on the system and the energy of the system.
+ */
+class ReferenceCalcCustomManyParticleForceKernel : public CalcCustomManyParticleForceKernel {
+public:
+    ReferenceCalcCustomManyParticleForceKernel(std::string name, const Platform& platform) : CalcCustomManyParticleForceKernel(name, platform), ixn(NULL) {
+    }
+    ~ReferenceCalcCustomManyParticleForceKernel();
+    /**
+     * 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:
+    int numParticles;
+    RealOpenMM cutoffDistance;
+    RealOpenMM **particleParamArray;
+    ReferenceCustomManyParticleIxn* ixn;
+    std::vector<std::string> globalParameterNames;
+    NonbondedMethod nonbondedMethod;
+};
+
 /**
  * This kernel is invoked by VerletIntegrator to take one time step.
  */

platforms/reference/src/ReferenceKernelFactory.cpp

@@ -80,6 +80,8 @@ KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Pla
         return new ReferenceCalcCustomHbondForceKernel(name, platform);
     if (name == CalcCustomCompoundBondForceKernel::Name())
         return new ReferenceCalcCustomCompoundBondForceKernel(name, platform);
+    if (name == CalcCustomManyParticleForceKernel::Name())
+        return new ReferenceCalcCustomManyParticleForceKernel(name, platform);
     if (name == IntegrateVerletStepKernel::Name())
         return new ReferenceIntegrateVerletStepKernel(name, platform, data);
     if (name == IntegrateLangevinStepKernel::Name())

platforms/reference/src/ReferenceKernels.cpp

@@ -47,6 +47,7 @@
 #include "ReferenceCustomGBIxn.h"
 #include "ReferenceCustomHbondIxn.h"
 #include "ReferenceCustomNonbondedIxn.h"
+#include "ReferenceCustomManyParticleIxn.h"
 #include "ReferenceCustomTorsionIxn.h"
 #include "ReferenceHarmonicBondIxn.h"
 #include "ReferenceLJCoulomb14.h"
@@ -67,6 +68,7 @@
 #include "openmm/internal/ContextImpl.h"
 #include "openmm/internal/CustomCompoundBondForceImpl.h"
 #include "openmm/internal/CustomHbondForceImpl.h"
+#include "openmm/internal/CustomManyParticleForceImpl.h"
 #include "openmm/internal/CustomNonbondedForceImpl.h"
 #include "openmm/internal/CMAPTorsionForceImpl.h"
 #include "openmm/internal/NonbondedForceImpl.h"
@@ -1615,6 +1617,91 @@ void ReferenceCalcCustomCompoundBondForceKernel::copyParametersToContext(Context
     }
 }
 
+ReferenceCalcCustomManyParticleForceKernel::~ReferenceCalcCustomManyParticleForceKernel() {
+    disposeRealArray(particleParamArray, numParticles);
+    if (ixn != NULL)
+        delete ixn;
+}
+
+void ReferenceCalcCustomManyParticleForceKernel::initialize(const System& system, const CustomManyParticleForce& force) {
+
+    // Build the arrays.
+
+    numParticles = system.getNumParticles();
+    int numParticleParameters = force.getNumPerParticleParameters();
+    particleParamArray = allocateRealArray(numParticles, 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];
+    }
+
+    // Create custom functions for the tabulated functions.
+
+    map<string, Lepton::CustomFunction*> functions;
+    for (int i = 0; i < 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 energyExpression = CustomManyParticleForceImpl::prepareExpression(force, functions, distances, angles, dihedrals);
+    vector<string> particleParameterNames;
+    for (int i = 0; i < numParticleParameters; i++)
+        particleParameterNames.push_back(force.getPerParticleParameterName(i));
+    for (int i = 0; i < force.getNumGlobalParameters(); i++)
+        globalParameterNames.push_back(force.getGlobalParameterName(i));
+    ixn = new ReferenceCustomManyParticleIxn(force.getNumParticlesPerSet(), energyExpression, particleParameterNames, distances, angles, dihedrals);
+    nonbondedMethod = CalcCustomManyParticleForceKernel::NonbondedMethod(force.getNonbondedMethod());
+    cutoffDistance = force.getCutoffDistance();
+    if (nonbondedMethod != NoCutoff)
+        ixn->setUseCutoff(cutoffDistance);
+
+    // Delete the custom functions.
+
+    for (map<string, Lepton::CustomFunction*>::iterator iter = functions.begin(); iter != functions.end(); iter++)
+        delete iter->second;
+}
+
+double ReferenceCalcCustomManyParticleForceKernel::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(posData, particleParamArray, globalParameters, forceData, includeEnergy ? &energy : NULL);
+    return energy;
+}
+
+void ReferenceCalcCustomManyParticleForceKernel::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]);
+    }
+}
+
 ReferenceIntegrateVerletStepKernel::~ReferenceIntegrateVerletStepKernel() {
     if (dynamics)
         delete dynamics;

platforms/reference/src/ReferencePlatform.cpp

@@ -63,6 +63,7 @@ ReferencePlatform::ReferencePlatform() {
     registerKernelFactory(CalcCustomExternalForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomHbondForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomCompoundBondForceKernel::Name(), factory);
+    registerKernelFactory(CalcCustomManyParticleForceKernel::Name(), factory);
     registerKernelFactory(IntegrateVerletStepKernel::Name(), factory);
     registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);
     registerKernelFactory(IntegrateBrownianStepKernel::Name(), factory);

platforms/reference/src/SimTKReference/ReferenceBondIxn.cpp

@@ -105,7 +105,7 @@ ReferenceBondIxn::~ReferenceBondIxn( ){
    --------------------------------------------------------------------------------------- */
 
 RealOpenMM ReferenceBondIxn::getNormedDotProduct( RealOpenMM* vector1, RealOpenMM* vector2,
-                                                  int hasREntry = 0 ) const {
+                                                  int hasREntry = 0 ) {
 
    // ---------------------------------------------------------------------------------------
 
@@ -187,7 +187,7 @@ RealOpenMM ReferenceBondIxn::getNormedDotProduct( RealOpenMM* vector1, RealOpenM
 
 RealOpenMM ReferenceBondIxn::getAngleBetweenTwoVectors( RealOpenMM* vector1, RealOpenMM* vector2, 
                                                         RealOpenMM* outputDotProduct = NULL,
-                                                        int hasREntry = 0 ) const {
+                                                        int hasREntry = 0 ) {
 
    // ---------------------------------------------------------------------------------------
 
@@ -249,7 +249,7 @@ RealOpenMM ReferenceBondIxn::getDihedralAngleBetweenThreeVectors( RealOpenMM*  v
                                                                   RealOpenMM*  cosineOfAngle       = NULL, 
                                                                   RealOpenMM*  signVector          = NULL, 
                                                                   RealOpenMM*  signOfAngle         = NULL, 
-                                                                   int          hasREntry = 0 ) const {
+                                                                   int          hasREntry = 0 ) {
 
    // ---------------------------------------------------------------------------------------
 

platforms/reference/src/SimTKReference/ReferenceCustomManyParticleIxn.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 "ReferenceCustomManyParticleIxn.h"
+
+using std::map;
+using std::pair;
+using std::string;
+using std::stringstream;
+using std::vector;
+using OpenMM::RealVec;
+
+ReferenceCustomManyParticleIxn::ReferenceCustomManyParticleIxn(int numParticlesPerSet,
+            const Lepton::ParsedExpression& energyExpression, const vector<string>& particleParameterNames,
+            const map<string, vector<int> >& distances, const map<string, vector<int> >& angles, const map<string, vector<int> >& dihedrals) :
+            numParticlesPerSet(numParticlesPerSet), energyExpression(energyExpression.createProgram()), useCutoff(false), usePeriodic(false) {
+    particleParamNames.resize(numParticlesPerSet);
+    numPerParticleParameters = particleParameterNames.size();
+    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(ReferenceCustomManyParticleIxn::ParticleTermInfo(xname.str(), i, 0, energyExpression.differentiate(xname.str()).optimize().createProgram()));
+        particleTerms.push_back(ReferenceCustomManyParticleIxn::ParticleTermInfo(yname.str(), i, 1, energyExpression.differentiate(yname.str()).optimize().createProgram()));
+        particleTerms.push_back(ReferenceCustomManyParticleIxn::ParticleTermInfo(zname.str(), i, 2, energyExpression.differentiate(zname.str()).optimize().createProgram()));
+        for (int j = 0; j < numPerParticleParameters; j++) {
+            stringstream paramname;
+            paramname << particleParameterNames[j] << (i+1);
+            particleParamNames[i].push_back(paramname.str());
+        }
+    }
+    for (map<string, vector<int> >::const_iterator iter = distances.begin(); iter != distances.end(); ++iter)
+        distanceTerms.push_back(ReferenceCustomManyParticleIxn::DistanceTermInfo(iter->first, iter->second, energyExpression.differentiate(iter->first).optimize().createProgram()));
+    for (map<string, vector<int> >::const_iterator iter = angles.begin(); iter != angles.end(); ++iter)
+        angleTerms.push_back(ReferenceCustomManyParticleIxn::AngleTermInfo(iter->first, iter->second, energyExpression.differentiate(iter->first).optimize().createProgram()));
+    for (map<string, vector<int> >::const_iterator iter = dihedrals.begin(); iter != dihedrals.end(); ++iter)
+        dihedralTerms.push_back(ReferenceCustomManyParticleIxn::DihedralTermInfo(iter->first, iter->second, energyExpression.differentiate(iter->first).optimize().createProgram()));
+}
+
+ReferenceCustomManyParticleIxn::~ReferenceCustomManyParticleIxn( ){
+}
+
+void ReferenceCustomManyParticleIxn::calculateIxn(vector<RealVec>& atomCoordinates, RealOpenMM** particleParameters,
+                                                  const map<string, double>& globalParameters, vector<RealVec>& forces,
+                                                  RealOpenMM* totalEnergy) const {
+    map<string, double> variables = globalParameters;
+    vector<int> particles(numParticlesPerSet);
+    loopOverInteractions(particles, 0, atomCoordinates, particleParameters, variables, forces, totalEnergy);
+}
+
+void ReferenceCustomManyParticleIxn::setUseCutoff(RealOpenMM distance) {
+    useCutoff = true;
+    cutoffDistance = distance;
+}
+
+void ReferenceCustomManyParticleIxn::setPeriodic(RealVec& boxSize) {
+    assert(cutoff);
+    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 ReferenceCustomManyParticleIxn::loopOverInteractions(vector<int>& particles, int loopIndex, vector<OpenMM::RealVec>& atomCoordinates,
+                                                          RealOpenMM** particleParameters, map<string, double>& variables, vector<OpenMM::RealVec>& forces,
+                                                          RealOpenMM* totalEnergy) const {
+    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++)
+            variables[particleParamNames[loopIndex][j]] = particleParameters[i][j];
+        if (loopIndex == numParticlesPerSet-1)
+            calculateOneIxn(particles, atomCoordinates, variables, forces, totalEnergy);
+        else
+            loopOverInteractions(particles, loopIndex+1, atomCoordinates, particleParameters, variables, forces, totalEnergy);
+    }
+}
+
+void ReferenceCustomManyParticleIxn::calculateOneIxn(const vector<int>& particles, vector<RealVec>& atomCoordinates,
+                        map<string, double>& variables, vector<RealVec>& forces, RealOpenMM* totalEnergy) const {
+    // Compute all of the variables the energy can depend on.
+
+    for (int i = 0; i < (int) particleTerms.size(); i++) {
+        const ParticleTermInfo& term = particleTerms[i];
+        variables[term.name] = atomCoordinates[term.atom][term.component];
+    }
+    for (int i = 0; i < (int) distanceTerms.size(); i++) {
+        const DistanceTermInfo& term = distanceTerms[i];
+        computeDelta(particles[term.p1], particles[term.p2], term.delta, atomCoordinates);
+        variables[term.name] = term.delta[ReferenceForce::RIndex];
+        if (useCutoff && term.delta[ReferenceForce::RIndex] > cutoffDistance)
+            return;
+    }
+    for (int i = 0; i < (int) angleTerms.size(); i++) {
+        const AngleTermInfo& term = angleTerms[i];
+        computeDelta(particles[term.p1], particles[term.p2], term.delta1, atomCoordinates);
+        computeDelta(particles[term.p3], particles[term.p2], term.delta2, atomCoordinates);
+        variables[term.name] = computeAngle(term.delta1, term.delta2);
+    }
+    for (int i = 0; i < (int) dihedralTerms.size(); i++) {
+        const DihedralTermInfo& term = dihedralTerms[i];
+        computeDelta(particles[term.p2], particles[term.p1], term.delta1, atomCoordinates);
+        computeDelta(particles[term.p2], particles[term.p3], term.delta2, atomCoordinates);
+        computeDelta(particles[term.p4], particles[term.p3], term.delta3, atomCoordinates);
+        RealOpenMM dotDihedral, signOfDihedral;
+        RealOpenMM* crossProduct[] = {term.cross1, term.cross2};
+        variables[term.name] = 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[particles[term.atom]][term.component] -= term.forceExpression.evaluate(variables);
+    }
+
+    // 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(variables)/(term.delta[ReferenceForce::RIndex]));
+        for (int i = 0; i < 3; i++) {
+           RealOpenMM force  = -dEdR*term.delta[i];
+           forces[particles[term.p1]][i] -= force;
+           forces[particles[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(variables);
+        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[particles[term.p1]][i] += deltaCrossP[0][i];
+            forces[particles[term.p2]][i] += deltaCrossP[1][i];
+            forces[particles[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(variables);
+        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[particles[term.p1]][i] += internalF[0][i];
+            forces[particles[term.p2]][i] -= internalF[1][i];
+            forces[particles[term.p3]][i] -= internalF[2][i];
+            forces[particles[term.p4]][i] += internalF[3][i];
+        }
+    }
+
+    // Add the energy
+
+    if (totalEnergy)
+        *totalEnergy += (RealOpenMM) energyExpression.evaluate(variables);
+}
+
+void ReferenceCustomManyParticleIxn::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 ReferenceCustomManyParticleIxn::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;
+}

platforms/reference/tests/TestReferenceCustomManyParticleForce.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   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.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the reference implementation of CustomManyParticleForce.
+ */
+
+#ifdef WIN32
+  #define _USE_MATH_DEFINES // Needed to get M_PI
+#endif
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "ReferencePlatform.h"
+#include "openmm/CustomManyParticleForce.h"
+#include "openmm/System.h"
+#include "openmm/VerletIntegrator.h"
+#include "sfmt/SFMT.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+const double TOL = 1e-5;
+
+void validateAxilrodTeller(CustomManyParticleForce* force, const vector<Vec3>& positions, const vector<const int*>& expectedSets, double boxSize) {
+    // Create a System and Context.
+    
+    int numParticles = force->getNumParticles();
+    CustomManyParticleForce::NonbondedMethod nonbondedMethod = force->getNonbondedMethod();
+    System system;
+    for (int i = 0; i < numParticles; i++)
+        system.addParticle(1.0);
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    system.addForce(force);
+    VerletIntegrator integrator(0.001);
+    ReferencePlatform platform;
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    State state1 = context.getState(State::Forces | State::Energy);
+    double c = context.getParameter("C");
+    
+    // See if the energy matches the expected value.
+    
+    double expectedEnergy = 0;
+    for (int i = 0; i < (int) expectedSets.size(); i++) {
+        int p1 = expectedSets[i][0];
+        int p2 = expectedSets[i][1];
+        int p3 = expectedSets[i][2];
+        Vec3 d12 = positions[p2]-positions[p1];
+        Vec3 d13 = positions[p3]-positions[p1];
+        Vec3 d23 = positions[p3]-positions[p2];
+        if (nonbondedMethod == CustomManyParticleForce::CutoffPeriodic) {
+            for (int j = 0; j < 3; j++) {
+                d12[j] -= floor(d12[j]/boxSize+0.5f)*boxSize;
+                d13[j] -= floor(d13[j]/boxSize+0.5f)*boxSize;
+                d23[j] -= floor(d23[j]/boxSize+0.5f)*boxSize;
+            }
+        }
+        double r12 = sqrt(d12.dot(d12));
+        double r13 = sqrt(d13.dot(d13));
+        double r23 = sqrt(d23.dot(d23));
+        double ctheta1 = d12.dot(d13)/(r12*r13);
+        double ctheta2 = -d12.dot(d23)/(r12*r23);
+        double ctheta3 = d13.dot(d23)/(r13*r23);
+        double rprod = r12*r13*r23;
+        expectedEnergy += c*(1+3*ctheta1*ctheta2*ctheta3)/(rprod*rprod*rprod);
+    }
+    ASSERT_EQUAL_TOL(expectedEnergy, state1.getPotentialEnergy(), 1e-5);
+
+    // Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.
+
+    const vector<Vec3>& forces = state1.getForces();
+    double norm = 0.0;
+    for (int i = 0; i < (int) forces.size(); ++i)
+        norm += forces[i].dot(forces[i]);
+    norm = std::sqrt(norm);
+    const double stepSize = 1e-3;
+    double step = 0.5*stepSize/norm;
+    vector<Vec3> positions2(numParticles), positions3(numParticles);
+    for (int i = 0; i < (int) positions.size(); ++i) {
+        Vec3 p = positions[i];
+        Vec3 f = forces[i];
+        positions2[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
+        positions3[i] = Vec3(p[0]+f[0]*step, p[1]+f[1]*step, p[2]+f[2]*step);
+    }
+    context.setPositions(positions2);
+    State state2 = context.getState(State::Energy);
+    context.setPositions(positions3);
+    State state3 = context.getState(State::Energy);
+    ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state3.getPotentialEnergy())/stepSize, 1e-4);
+}
+
+void testNoCutoff() {
+    CustomManyParticleForce* force = new CustomManyParticleForce(3,
+        "C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
+        "theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
+        "r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
+    force->addGlobalParameter("C", 1.5);
+    vector<double> params;
+    force->addParticle(params);
+    force->addParticle(params);
+    force->addParticle(params);
+    force->addParticle(params);
+    vector<Vec3> positions;
+    positions.push_back(Vec3(0, 0, 0));
+    positions.push_back(Vec3(1, 0, 0));
+    positions.push_back(Vec3(0, 1.1, 0.3));
+    positions.push_back(Vec3(0.4, 0, -0.8));
+    int sets[4][3] = {{0,1,2}, {1,2,3}, {2,3,0}, {3,0,1}};
+    vector<const int*> expectedSets(&sets[0], &sets[4]);
+    validateAxilrodTeller(force, positions, expectedSets, 2.0);
+}
+
+void testCutoff() {
+    CustomManyParticleForce* force = new CustomManyParticleForce(3,
+        "C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
+        "theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
+        "r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
+    force->addGlobalParameter("C", 1.5);
+    force->setNonbondedMethod(CustomManyParticleForce::CutoffNonPeriodic);
+    force->setCutoffDistance(1.55);
+    vector<double> params;
+    force->addParticle(params);
+    force->addParticle(params);
+    force->addParticle(params);
+    force->addParticle(params);
+    force->addParticle(params);
+    vector<Vec3> positions;
+    positions.push_back(Vec3(0, 0, 0));
+    positions.push_back(Vec3(1, 0, 0));
+    positions.push_back(Vec3(0, 1.1, 0.3));
+    positions.push_back(Vec3(0.4, 0, -0.8));
+    positions.push_back(Vec3(0.2, 0.5, -0.1));
+    int sets[7][3] = {{0,1,2}, {0,1,3}, {0,1,4}, {0,2,4}, {0,3,4}, {1,2,4}, {1,3,4}};
+    vector<const int*> expectedSets(&sets[0], &sets[7]);
+    validateAxilrodTeller(force, positions, expectedSets, 2.0);
+}
+
+void testPeriodic() {
+    CustomManyParticleForce* force = new CustomManyParticleForce(3,
+        "C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
+        "theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
+        "r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
+    force->addGlobalParameter("C", 1.5);
+    force->setNonbondedMethod(CustomManyParticleForce::CutoffPeriodic);
+    force->setCutoffDistance(1.05);
+    vector<double> params;
+    force->addParticle(params);
+    force->addParticle(params);
+    force->addParticle(params);
+    force->addParticle(params);
+    force->addParticle(params);
+    vector<Vec3> positions;
+    positions.push_back(Vec3(0, 0, 0));
+    positions.push_back(Vec3(1, 0, 0));
+    positions.push_back(Vec3(0, 1.1, 0.3));
+    positions.push_back(Vec3(0.4, 0, -0.8));
+    positions.push_back(Vec3(0.2, 0.5, -0.1));
+    double boxSize = 2.1;
+    int sets[5][3] = {{0,1,3}, {0,1,4}, {0,2,4}, {0,3,4}, {1,3,4}};
+    vector<const int*> expectedSets(&sets[0], &sets[5]);
+    validateAxilrodTeller(force, positions, expectedSets, boxSize);
+}
+
+void testParameters() {
+    // Create a system.
+    
+    int numParticles = 5;
+    System system;
+    CustomManyParticleForce* force = new CustomManyParticleForce(3, "C*scale1*scale2*scale3*(distance(p1,p2)+distance(p2,p3)+distance(p1,p3))");
+    force->addGlobalParameter("C", 2.0);
+    force->addPerParticleParameter("scale");
+    vector<double> params(1);
+    vector<Vec3> positions;
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    for (int i = 0; i < numParticles; i++) {
+        params[0] = i+1;
+        force->addParticle(params);
+        positions.push_back(Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt)));
+        system.addParticle(1.0);
+    }
+    system.addForce(force);
+    VerletIntegrator integrator(0.001);
+    ReferencePlatform platform;
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    
+    // See if the energy is correct.
+
+    State state1 = context.getState(State::Energy);
+    double expectedEnergy = 0;
+    for (int i = 0; i < numParticles; i++)
+        for (int j = i+1; j < numParticles; j++)
+            for (int k = j+1; k < numParticles; k++) {
+                Vec3 d12 = positions[j]-positions[i];
+                Vec3 d13 = positions[k]-positions[i];
+                Vec3 d23 = positions[k]-positions[j];
+                double r12 = sqrt(d12.dot(d12));
+                double r13 = sqrt(d13.dot(d13));
+                double r23 = sqrt(d23.dot(d23));
+                expectedEnergy += 2.0*(i+1)*(j+1)*(k+1)*(r12+r13+r23);
+            }
+    ASSERT_EQUAL_TOL(expectedEnergy, state1.getPotentialEnergy(), 1e-5);
+}
+
+int main() {
+    try {
+        testNoCutoff();
+        testCutoff();
+        testPeriodic();
+        testParameters();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}