Merge pull request #1119 from peastman/centroid

Created CustomCentroidBondForce

platforms/reference/include/ReferenceCustomCentroidBondIxn.h

+
+/* Portions copyright (c) 2009-2015 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 __ReferenceCustomCentroidBondIxn_H__
+#define __ReferenceCustomCentroidBondIxn_H__
+
+#include "ReferenceBondIxn.h"
+#include "lepton/ExpressionProgram.h"
+#include "lepton/ParsedExpression.h"
+#include <map>
+#include <vector>
+
+namespace OpenMM {
+
+class ReferenceCustomCentroidBondIxn : public ReferenceBondIxn {
+
+   private:
+
+      class PositionTermInfo;
+      class DistanceTermInfo;
+      class AngleTermInfo;
+      class DihedralTermInfo;
+      std::vector<std::vector<int> > groupAtoms;
+      std::vector<std::vector<double> > normalizedWeights;
+      std::vector<std::vector<int> > bondGroups;
+      Lepton::ExpressionProgram energyExpression;
+      std::vector<std::string> bondParamNames;
+      std::vector<PositionTermInfo> positionTerms;
+      std::vector<DistanceTermInfo> distanceTerms;
+      std::vector<AngleTermInfo> angleTerms;
+      std::vector<DihedralTermInfo> dihedralTerms;
+
+
+      /**---------------------------------------------------------------------------------------
+
+         Calculate custom interaction for one bond
+
+         @param bond             the index of the bond
+         @param groupCenters     group center coordinates
+         @param variables        the values of variables that may appear in expressions
+         @param forces           force array (forces added)
+         @param totalEnergy      total energy
+
+         --------------------------------------------------------------------------------------- */
+
+      void calculateOneIxn(int bond, std::vector<OpenMM::RealVec>& groupCenters,
+                           std::map<std::string, double>& variables, std::vector<OpenMM::RealVec>& forces,
+                           RealOpenMM* totalEnergy) const;
+
+      void computeDelta(int group1, int group2, RealOpenMM* delta, std::vector<OpenMM::RealVec>& groupCenters) const;
+
+      static RealOpenMM computeAngle(RealOpenMM* vec1, RealOpenMM* vec2);
+
+
+   public:
+
+      /**---------------------------------------------------------------------------------------
+
+         Constructor
+
+         --------------------------------------------------------------------------------------- */
+
+       ReferenceCustomCentroidBondIxn(int numGroupsPerBond, const std::vector<std::vector<int> >& groupAtoms,
+                               const std::vector<std::vector<double> >& normalizedWeights, const std::vector<std::vector<int> >& bondGroups, const Lepton::ParsedExpression& energyExpression,
+                               const std::vector<std::string>& bondParameterNames, 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
+
+         --------------------------------------------------------------------------------------- */
+
+       ~ReferenceCustomCentroidBondIxn();
+
+      /**---------------------------------------------------------------------------------------
+
+         Get the list of groups in each bond.
+
+         --------------------------------------------------------------------------------------- */
+
+       const std::vector<std::vector<int> >& getBondGroups() const {
+           return bondGroups;
+       }
+
+      /**---------------------------------------------------------------------------------------
+
+         Calculate custom compound bond interaction
+
+         @param atomCoordinates    atom coordinates
+         @param bondParameters     bond parameters values       bondParameters[bondIndex][parameterIndex]
+         @param globalParameters   the values of global parameters
+         @param forces             force array (forces added)
+         @param totalEnergy        total energy
+
+         --------------------------------------------------------------------------------------- */
+
+      void calculatePairIxn(std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** bondParameters,
+                            const std::map<std::string, double>& globalParameters,
+                            std::vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy) const;
+
+// ---------------------------------------------------------------------------------------
+
+};
+
+class ReferenceCustomCentroidBondIxn::PositionTermInfo {
+public:
+    std::string name;
+    int group, component;
+    Lepton::ExpressionProgram forceExpression;
+    PositionTermInfo(const std::string& name, int group, int component, const Lepton::ExpressionProgram& forceExpression) :
+            name(name), group(group), component(component), forceExpression(forceExpression) {
+    }
+};
+
+class ReferenceCustomCentroidBondIxn::DistanceTermInfo {
+public:
+    std::string name;
+    int g1, g2;
+    Lepton::ExpressionProgram forceExpression;
+    mutable RealOpenMM delta[ReferenceForce::LastDeltaRIndex];
+    DistanceTermInfo(const std::string& name, const std::vector<int>& groups, const Lepton::ExpressionProgram& forceExpression) :
+            name(name), g1(groups[0]), g2(groups[1]), forceExpression(forceExpression) {
+    }
+};
+
+class ReferenceCustomCentroidBondIxn::AngleTermInfo {
+public:
+    std::string name;
+    int g1, g2, g3;
+    Lepton::ExpressionProgram forceExpression;
+    mutable RealOpenMM delta1[ReferenceForce::LastDeltaRIndex];
+    mutable RealOpenMM delta2[ReferenceForce::LastDeltaRIndex];
+    AngleTermInfo(const std::string& name, const std::vector<int>& groups, const Lepton::ExpressionProgram& forceExpression) :
+            name(name), g1(groups[0]), g2(groups[1]), g3(groups[2]), forceExpression(forceExpression) {
+    }
+};
+
+class ReferenceCustomCentroidBondIxn::DihedralTermInfo {
+public:
+    std::string name;
+    int g1, g2, g3, g4;
+    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>& groups, const Lepton::ExpressionProgram& forceExpression) :
+            name(name), g1(groups[0]), g2(groups[1]), g3(groups[2]), g4(groups[3]), forceExpression(forceExpression) {
+    }
+};
+
+} // namespace OpenMM
+
+#endif // __ReferenceCustomCentroidBondIxn_H__

platforms/reference/include/ReferenceKernels.h

@@ -44,6 +44,7 @@ namespace OpenMM {
 class ReferenceObc;
 class ReferenceGBVI;
 class ReferenceAndersenThermostat;
+class ReferenceCustomCentroidBondIxn;
 class ReferenceCustomCompoundBondIxn;
 class ReferenceCustomHbondIxn;
 class ReferenceCustomManyParticleIxn;
@@ -833,6 +834,44 @@ private:
     std::vector<std::string> globalParameterNames;
 };
 
+/**
+ * This kernel is invoked by CustomCentroidBondForce to calculate the forces acting on the system.
+ */
+class ReferenceCalcCustomCentroidBondForceKernel : public CalcCustomCentroidBondForceKernel {
+public:
+    ReferenceCalcCustomCentroidBondForceKernel(std::string name, const Platform& platform) : CalcCustomCentroidBondForceKernel(name, platform), ixn(NULL) {
+    }
+    ~ReferenceCalcCustomCentroidBondForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CustomCentroidBondForce this kernel will be used for
+     */
+    void initialize(const System& system, const CustomCentroidBondForce& 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 CustomCentroidBondForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomCentroidBondForce& force);
+private:
+    int numBonds, numParticles;
+    RealOpenMM **bondParamArray;
+    ReferenceCustomCentroidBondIxn* ixn;
+    std::vector<std::string> globalParameterNames;
+};
+
 /**
  * This kernel is invoked by CustomCompoundBondForce to calculate the forces acting on the system.
  */
@@ -865,7 +904,7 @@ public:
      */
     void copyParametersToContext(ContextImpl& context, const CustomCompoundBondForce& force);
 private:
-    int numBonds, numParticles;
+    int numBonds;
     RealOpenMM **bondParamArray;
     ReferenceCustomCompoundBondIxn* ixn;
     std::vector<std::string> globalParameterNames;

platforms/reference/src/ReferenceKernelFactory.cpp

@@ -78,6 +78,8 @@ KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Pla
         return new ReferenceCalcCustomExternalForceKernel(name, platform);
     if (name == CalcCustomHbondForceKernel::Name())
         return new ReferenceCalcCustomHbondForceKernel(name, platform);
+    if (name == CalcCustomCentroidBondForceKernel::Name())
+        return new ReferenceCalcCustomCentroidBondForceKernel(name, platform);
     if (name == CalcCustomCompoundBondForceKernel::Name())
         return new ReferenceCalcCustomCompoundBondForceKernel(name, platform);
     if (name == CalcCustomManyParticleForceKernel::Name())

platforms/reference/src/ReferenceKernels.cpp

@@ -41,6 +41,7 @@
 #include "ReferenceConstraints.h"
 #include "ReferenceCustomAngleIxn.h"
 #include "ReferenceCustomBondIxn.h"
+#include "ReferenceCustomCentroidBondIxn.h"
 #include "ReferenceCustomCompoundBondIxn.h"
 #include "ReferenceCustomDynamics.h"
 #include "ReferenceCustomExternalIxn.h"
@@ -66,6 +67,7 @@
 #include "openmm/System.h"
 #include "openmm/internal/AndersenThermostatImpl.h"
 #include "openmm/internal/ContextImpl.h"
+#include "openmm/internal/CustomCentroidBondForceImpl.h"
 #include "openmm/internal/CustomCompoundBondForceImpl.h"
 #include "openmm/internal/CustomHbondForceImpl.h"
 #include "openmm/internal/CustomNonbondedForceImpl.h"
@@ -1582,6 +1584,90 @@ void ReferenceCalcCustomHbondForceKernel::copyParametersToContext(ContextImpl& c
     }
 }
 
+ReferenceCalcCustomCentroidBondForceKernel::~ReferenceCalcCustomCentroidBondForceKernel() {
+    disposeRealArray(bondParamArray, numBonds);
+    if (ixn != NULL)
+        delete ixn;
+}
+
+void ReferenceCalcCustomCentroidBondForceKernel::initialize(const System& system, const CustomCentroidBondForce& force) {
+
+    // Build the arrays.
+
+    int numGroups = force.getNumGroups();
+    vector<vector<int> > groupAtoms(numGroups);
+    vector<double> ignored;
+    for (int i = 0; i < numGroups; i++)
+        force.getGroupParameters(i, groupAtoms[i], ignored);
+    vector<vector<double> > normalizedWeights;
+    CustomCentroidBondForceImpl::computeNormalizedWeights(force, system, normalizedWeights);
+    numBonds = force.getNumBonds();
+    vector<vector<int> > bondGroups(numBonds);
+    int numBondParameters = force.getNumPerBondParameters();
+    bondParamArray = allocateRealArray(numBonds, numBondParameters);
+    for (int i = 0; i < numBonds; ++i) {
+        vector<double> parameters;
+        force.getBondParameters(i, bondGroups[i], parameters);
+        for (int j = 0; j < numBondParameters; j++)
+            bondParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
+    }
+
+    // Create custom functions for the tabulated functions.
+
+    map<string, Lepton::CustomFunction*> functions;
+    for (int i = 0; i < force.getNumFunctions(); 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 = CustomCentroidBondForceImpl::prepareExpression(force, functions, distances, angles, dihedrals);
+    vector<string> bondParameterNames;
+    for (int i = 0; i < numBondParameters; i++)
+        bondParameterNames.push_back(force.getPerBondParameterName(i));
+    for (int i = 0; i < force.getNumGlobalParameters(); i++)
+        globalParameterNames.push_back(force.getGlobalParameterName(i));
+    ixn = new ReferenceCustomCentroidBondIxn(force.getNumGroupsPerBond(), groupAtoms, normalizedWeights, bondGroups, energyExpression, bondParameterNames, distances, angles, dihedrals);
+
+    // Delete the custom functions.
+
+    for (map<string, Lepton::CustomFunction*>::iterator iter = functions.begin(); iter != functions.end(); iter++)
+        delete iter->second;
+}
+
+double ReferenceCalcCustomCentroidBondForceKernel::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]);
+    ixn->calculatePairIxn(posData, bondParamArray, globalParameters, forceData, includeEnergy ? &energy : NULL);
+    return energy;
+}
+
+void ReferenceCalcCustomCentroidBondForceKernel::copyParametersToContext(ContextImpl& context, const CustomCentroidBondForce& force) {
+    if (numBonds != force.getNumBonds())
+        throw OpenMMException("updateParametersInContext: The number of bonds has changed");
+
+    // Record the values.
+
+    int numParameters = force.getNumPerBondParameters();
+    const vector<vector<int> >& bondGroups = ixn->getBondGroups();
+    vector<int> groups;
+    vector<double> params;
+    for (int i = 0; i < numBonds; ++i) {
+        force.getBondParameters(i, groups, params);
+        for (int j = 0; j < groups.size(); j++)
+            if (groups[j] != bondGroups[i][j])
+                throw OpenMMException("updateParametersInContext: The set of groups in a bond has changed");
+        for (int j = 0; j < numParameters; j++)
+            bondParamArray[i][j] = (RealOpenMM) params[j];
+    }
+}
+
 ReferenceCalcCustomCompoundBondForceKernel::~ReferenceCalcCustomCompoundBondForceKernel() {
     disposeRealArray(bondParamArray, numBonds);
     if (ixn != NULL)
@@ -1593,7 +1679,6 @@ void ReferenceCalcCustomCompoundBondForceKernel::initialize(const System& system
     // Build the arrays.
 
     numBonds = force.getNumBonds();
-    numParticles = system.getNumParticles();
     vector<vector<int> > bondParticles(numBonds);
     int numBondParameters = force.getNumPerBondParameters();
     bondParamArray = allocateRealArray(numBonds, numBondParameters);
@@ -1652,7 +1737,7 @@ void ReferenceCalcCustomCompoundBondForceKernel::copyParametersToContext(Context
     vector<double> params;
     for (int i = 0; i < numBonds; ++i) {
         force.getBondParameters(i, particles, params);
-        for (int j = 0; j < numParticles; j++)
+        for (int j = 0; j < particles.size(); j++)
             if (particles[j] != bondAtoms[i][j])
                 throw OpenMMException("updateParametersInContext: The set of particles in a bond has changed");
         for (int j = 0; j < numParameters; j++)

platforms/reference/src/ReferencePlatform.cpp

@@ -62,6 +62,7 @@ ReferencePlatform::ReferencePlatform() {
     registerKernelFactory(CalcCustomGBForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomExternalForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomHbondForceKernel::Name(), factory);
+    registerKernelFactory(CalcCustomCentroidBondForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomCompoundBondForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomManyParticleForceKernel::Name(), factory);
     registerKernelFactory(IntegrateVerletStepKernel::Name(), factory);

platforms/reference/src/SimTKReference/ReferenceCustomCentroidBondIxn.cpp

+
+/* Portions copyright (c) 2009-2015 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 "SimTKOpenMMUtilities.h"
+#include "ReferenceForce.h"
+#include "ReferenceCustomCentroidBondIxn.h"
+
+using std::map;
+using std::pair;
+using std::string;
+using std::stringstream;
+using std::vector;
+using namespace OpenMM;
+
+ReferenceCustomCentroidBondIxn::ReferenceCustomCentroidBondIxn(int numGroupsPerBond, const vector<vector<int> >& groupAtoms,
+            const vector<vector<double> >& normalizedWeights, const vector<vector<int> >& bondGroups,
+            const Lepton::ParsedExpression& energyExpression, const vector<string>& bondParameterNames,
+            const map<string, vector<int> >& distances, const map<string, vector<int> >& angles, const map<string, vector<int> >& dihedrals) :
+            groupAtoms(groupAtoms), normalizedWeights(normalizedWeights), bondGroups(bondGroups), energyExpression(energyExpression.createProgram()), bondParamNames(bondParameterNames) {
+    for (int i = 0; i < numGroupsPerBond; i++) {
+        stringstream xname, yname, zname;
+        xname << 'x' << (i+1);
+        yname << 'y' << (i+1);
+        zname << 'z' << (i+1);
+        positionTerms.push_back(ReferenceCustomCentroidBondIxn::PositionTermInfo(xname.str(), i, 0, energyExpression.differentiate(xname.str()).optimize().createProgram()));
+        positionTerms.push_back(ReferenceCustomCentroidBondIxn::PositionTermInfo(yname.str(), i, 1, energyExpression.differentiate(yname.str()).optimize().createProgram()));
+        positionTerms.push_back(ReferenceCustomCentroidBondIxn::PositionTermInfo(zname.str(), i, 2, energyExpression.differentiate(zname.str()).optimize().createProgram()));
+    }
+    for (map<string, vector<int> >::const_iterator iter = distances.begin(); iter != distances.end(); ++iter)
+        distanceTerms.push_back(ReferenceCustomCentroidBondIxn::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(ReferenceCustomCentroidBondIxn::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(ReferenceCustomCentroidBondIxn::DihedralTermInfo(iter->first, iter->second, energyExpression.differentiate(iter->first).optimize().createProgram()));
+}
+
+ReferenceCustomCentroidBondIxn::~ReferenceCustomCentroidBondIxn() {
+}
+
+void ReferenceCustomCentroidBondIxn::calculatePairIxn(vector<RealVec>& atomCoordinates, RealOpenMM** bondParameters,
+                                             const map<string, double>& globalParameters, vector<RealVec>& forces,
+                                             RealOpenMM* totalEnergy) const {
+
+    // First compute the center of each group.
+
+    int numGroups = groupAtoms.size();
+    vector<RealVec> groupCenters(numGroups);
+    for (int group = 0; group < numGroups; group++) {
+        for (int i = 0; i < groupAtoms[group].size(); i++)
+            groupCenters[group] += atomCoordinates[groupAtoms[group][i]]*normalizedWeights[group][i];
+    }
+
+    // Compute the forces on groups.
+
+    map<string, double> variables = globalParameters;
+    vector<RealVec> groupForces(numGroups);
+    int numBonds = bondGroups.size();
+    for (int bond = 0; bond < numBonds; bond++) {
+        for (int j = 0; j < (int) bondParamNames.size(); j++)
+            variables[bondParamNames[j]] = bondParameters[bond][j];
+        calculateOneIxn(bond, groupCenters, variables, groupForces, totalEnergy);
+    }
+
+    // Apply the forces to the individual atoms.
+
+    for (int group = 0; group < numGroups; group++) {
+        for (int i = 0; i < groupAtoms[group].size(); i++)
+            forces[groupAtoms[group][i]] += groupForces[group]*normalizedWeights[group][i];
+    }
+}
+
+void ReferenceCustomCentroidBondIxn::calculateOneIxn(int bond, vector<RealVec>& groupCenters,
+                        map<string, double>& variables, vector<RealVec>& forces, RealOpenMM* totalEnergy) const {
+
+    // ---------------------------------------------------------------------------------------
+
+    static const std::string methodName = "\nReferenceCustomCentroidBondIxn::calculateOneIxn";
+
+    // ---------------------------------------------------------------------------------------
+
+    // Compute all of the variables the energy can depend on.
+
+    const vector<int>& groups = bondGroups[bond];
+    for (int i = 0; i < (int) positionTerms.size(); i++) {
+        const PositionTermInfo& term = positionTerms[i];
+        variables[term.name] = groupCenters[groups[term.group]][term.component];
+    }
+    for (int i = 0; i < (int) distanceTerms.size(); i++) {
+        const DistanceTermInfo& term = distanceTerms[i];
+        computeDelta(groups[term.g1], groups[term.g2], term.delta, groupCenters);
+        variables[term.name] = term.delta[ReferenceForce::RIndex];
+    }
+    for (int i = 0; i < (int) angleTerms.size(); i++) {
+        const AngleTermInfo& term = angleTerms[i];
+        computeDelta(groups[term.g1], groups[term.g2], term.delta1, groupCenters);
+        computeDelta(groups[term.g3], groups[term.g2], term.delta2, groupCenters);
+        variables[term.name] = computeAngle(term.delta1, term.delta2);
+    }
+    for (int i = 0; i < (int) dihedralTerms.size(); i++) {
+        const DihedralTermInfo& term = dihedralTerms[i];
+        computeDelta(groups[term.g2], groups[term.g1], term.delta1, groupCenters);
+        computeDelta(groups[term.g2], groups[term.g3], term.delta2, groupCenters);
+        computeDelta(groups[term.g4], groups[term.g3], term.delta3, groupCenters);
+        RealOpenMM dotDihedral, signOfDihedral;
+        RealOpenMM* crossProduct[] = {term.cross1, term.cross2};
+        variables[term.name] = 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) positionTerms.size(); i++) {
+        const PositionTermInfo& term = positionTerms[i];
+        forces[groups[term.group]][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[groups[term.g1]][i] -= force;
+           forces[groups[term.g2]][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[groups[term.g1]][i] += deltaCrossP[0][i];
+            forces[groups[term.g2]][i] += deltaCrossP[1][i];
+            forces[groups[term.g3]][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[groups[term.g1]][i] += internalF[0][i];
+            forces[groups[term.g2]][i] -= internalF[1][i];
+            forces[groups[term.g3]][i] -= internalF[2][i];
+            forces[groups[term.g4]][i] += internalF[3][i];
+        }
+    }
+
+    // Add the energy
+
+    if (totalEnergy)
+        *totalEnergy += (RealOpenMM) energyExpression.evaluate(variables);
+}
+
+void ReferenceCustomCentroidBondIxn::computeDelta(int group1, int group2, RealOpenMM* delta, vector<RealVec>& groupCenters) const {
+    ReferenceForce::getDeltaR(groupCenters[group1], groupCenters[group2], delta);
+}
+
+RealOpenMM ReferenceCustomCentroidBondIxn::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/TestReferenceCustomCentroidBondForce.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) 2015 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 CustomCompoundBondForce.
+ */
+
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "ReferencePlatform.h"
+#include "openmm/CustomCentroidBondForce.h"
+#include "openmm/CustomCompoundBondForce.h"
+#include "openmm/System.h"
+#include "openmm/TabulatedFunction.h"
+#include "openmm/VerletIntegrator.h"
+#include "sfmt/SFMT.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+ReferencePlatform platform;
+
+const double TOL = 1e-5;
+
+void testHarmonicBond() {
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(2.0);
+    system.addParticle(3.0);
+    system.addParticle(4.0);
+    system.addParticle(5.0);
+    CustomCentroidBondForce* force = new CustomCentroidBondForce(2, "k*distance(g1,g2)^2");
+    force->addPerBondParameter("k");
+    vector<int> particles1;
+    particles1.push_back(0);
+    particles1.push_back(1);
+    vector<int> particles2;
+    particles2.push_back(2);
+    particles2.push_back(3);
+    particles2.push_back(4);
+    force->addGroup(particles1);
+    force->addGroup(particles2);
+    vector<int> groups;
+    groups.push_back(0);
+    groups.push_back(1);
+    vector<double> parameters;
+    parameters.push_back(1.0);
+    force->addBond(groups, parameters);
+    system.addForce(force);
+    ASSERT(!system.usesPeriodicBoundaryConditions());
+
+    // The center of mass of group 0 is (1.5, 0, 0).
+
+    vector<Vec3> positions(5);
+    positions[0] = Vec3(2.5, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+
+    // The center of mass of group 1 is (-1, 0, 0).
+
+    positions[2] = Vec3(-6, 0, 0);
+    positions[3] = Vec3(-1, 0, 0);
+    positions[4] = Vec3(2, 0, 0);
+
+    // Check the forces and energy.
+
+    VerletIntegrator integrator(0.01);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    ASSERT_EQUAL_TOL(2.5*2.5, state.getPotentialEnergy(), TOL);
+    ASSERT_EQUAL_VEC(Vec3(-2*2.5*(1.0/3.0), 0, 0), state.getForces()[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(-2*2.5*(2.0/3.0), 0, 0), state.getForces()[1], TOL);
+    ASSERT_EQUAL_VEC(Vec3(2*2.5*(3.0/12.0), 0, 0), state.getForces()[2], TOL);
+    ASSERT_EQUAL_VEC(Vec3(2*2.5*(4.0/12.0), 0, 0), state.getForces()[3], TOL);
+    ASSERT_EQUAL_VEC(Vec3(2*2.5*(5.0/12.0), 0, 0), state.getForces()[4], TOL);
+
+    // Update the per-bond parameter and see if the results change.
+
+    parameters[0] = 2.0;
+    force->setBondParameters(0, groups, parameters);
+    force->updateParametersInContext(context);
+    state = context.getState(State::Forces | State::Energy);
+    ASSERT_EQUAL_TOL(2*2.5*2.5, state.getPotentialEnergy(), TOL);
+    ASSERT_EQUAL_VEC(Vec3(-4*2.5*(1.0/3.0), 0, 0), state.getForces()[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(-4*2.5*(2.0/3.0), 0, 0), state.getForces()[1], TOL);
+    ASSERT_EQUAL_VEC(Vec3(4*2.5*(3.0/12.0), 0, 0), state.getForces()[2], TOL);
+    ASSERT_EQUAL_VEC(Vec3(4*2.5*(4.0/12.0), 0, 0), state.getForces()[3], TOL);
+    ASSERT_EQUAL_VEC(Vec3(4*2.5*(5.0/12.0), 0, 0), state.getForces()[4], TOL);
+    
+    // All the particles should be treated as a single molecule.
+    
+    vector<std::vector<int> > molecules = context.getMolecules();
+    ASSERT_EQUAL(1, molecules.size());
+    ASSERT_EQUAL(5, molecules[0].size());
+}
+
+void testComplexFunction() {
+    int numParticles = 5;
+    System system;
+    for (int i = 0; i < numParticles; i++)
+        system.addParticle(2.0);
+    vector<double> table(20);
+    for (int i = 0; i < 20; i++)
+        table[i] = sin(0.11*i);
+
+    // When every group contains only one particle, a CustomCentroidBondForce is identical to a
+    // CustomCompoundBondForce.  Use that to test a complicated energy function with lots of terms.
+
+    CustomCompoundBondForce* compound = new CustomCompoundBondForce(4, "x1+y2+z4+fn(distance(p1,p2))*angle(p3,p2,p4)+scale*dihedral(p2,p1,p4,p3)");
+    CustomCentroidBondForce* centroid = new CustomCentroidBondForce(4, "x1+y2+z4+fn(distance(g1,g2))*angle(g3,g2,g4)+scale*dihedral(g2,g1,g4,g3)");
+    compound->addGlobalParameter("scale", 0.5);
+    centroid->addGlobalParameter("scale", 0.5);
+    compound->addTabulatedFunction("fn", new Continuous1DFunction(table, -1, 10));
+    centroid->addTabulatedFunction("fn", new Continuous1DFunction(table, -1, 10));
+
+    // Add two bonds to the CustomCompoundBondForce.
+
+    vector<int> particles(4);
+    vector<double> parameters;
+    particles[0] = 0;
+    particles[1] = 1;
+    particles[2] = 2;
+    particles[3] = 3;
+    compound->addBond(particles, parameters);
+    particles[0] = 2;
+    particles[1] = 4;
+    particles[2] = 3;
+    particles[3] = 1;
+    compound->addBond(particles, parameters);
+
+    // Add identical bonds to the CustomCentroidBondForce.  As a stronger test, make sure that
+    // group number is different from particle number.
+
+    vector<int> groupMembers(1);
+    groupMembers[0] = 3;
+    centroid->addGroup(groupMembers);
+    groupMembers[0] = 0;
+    centroid->addGroup(groupMembers);
+    groupMembers[0] = 1;
+    centroid->addGroup(groupMembers);
+    groupMembers[0] = 2;
+    centroid->addGroup(groupMembers);
+    groupMembers[0] = 4;
+    centroid->addGroup(groupMembers);
+    vector<int> groups(4);
+    groups[0] = 1;
+    groups[1] = 2;
+    groups[2] = 3;
+    groups[3] = 0;
+    centroid->addBond(groups, parameters);
+    groups[0] = 3;
+    groups[1] = 4;
+    groups[2] = 0;
+    groups[3] = 2;
+    centroid->addBond(groups, parameters);
+
+    // Add both forces as different force groups, and create a context.
+
+    centroid->setForceGroup(1);
+    system.addForce(compound);
+    system.addForce(centroid);
+    VerletIntegrator integrator(0.01);
+    Context context(system, integrator, platform);
+
+    // Evaluate the force and energy for various positions and see if they match.
+
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    vector<Vec3> positions(numParticles);
+    for (int i = 0; i < 10; i++) {
+        for (int j = 0; j < numParticles; j++)
+            positions[j] = Vec3(5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt));
+        context.setPositions(positions);
+        State state1 = context.getState(State::Forces | State::Energy, false, 1<<0);
+        State state2 = context.getState(State::Forces | State::Energy, false, 1<<1);
+        ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), TOL);
+        for (int i = 0; i < numParticles; i++)
+            ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], TOL);
+    }
+}
+
+void testCustomWeights() {
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(2.0);
+    system.addParticle(3.0);
+    system.addParticle(4.0);
+    CustomCentroidBondForce* force = new CustomCentroidBondForce(2, "distance(g1,g2)^2");
+    vector<int> particles(2);
+    vector<double> weights(2);
+    particles[0] = 0;
+    particles[1] = 1;
+    weights[0] = 0.5;
+    weights[1] = 1.5;
+    force->addGroup(particles, weights);
+    particles[0] = 2;
+    particles[1] = 3;
+    weights[0] = 2.0;
+    weights[1] = 1.0;
+    force->addGroup(particles, weights);
+    vector<int> groups;
+    groups.push_back(0);
+    groups.push_back(1);
+    vector<double> parameters;
+    force->addBond(groups, parameters);
+    system.addForce(force);
+
+    // The center of mass of group 0 is (0, 1, 0).
+
+    vector<Vec3> positions(4);
+    positions[0] = Vec3(0, 4, 0);
+    positions[1] = Vec3(0, 0, 0);
+
+    // The center of mass of group 1 is (0, 10, 0).
+
+    positions[2] = Vec3(0, 9, 0);
+    positions[3] = Vec3(0, 12, 0);
+
+    // Check the forces and energy.
+
+    VerletIntegrator integrator(0.01);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    ASSERT_EQUAL_TOL(9*9, state.getPotentialEnergy(), TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, 2*9*(0.5/2.0), 0), state.getForces()[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, 2*9*(1.5/2.0), 0), state.getForces()[1], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, -2*9*(2.0/3.0), 0), state.getForces()[2], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, -2*9*(1.0/3.0), 0), state.getForces()[3], TOL);
+}
+
+int main() {
+    try {
+        testHarmonicBond();
+        testComplexFunction();
+        testCustomWeights();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}