Implemented CustomHbondForce, including reference implementation.

olla/include/openmm/kernels.h

@@ -39,6 +39,7 @@
 #include "openmm/CustomBondForce.h"
 #include "openmm/CustomExternalForce.h"
 #include "openmm/CustomGBForce.h"
+#include "openmm/CustomHbondForce.h"
 #include "openmm/CustomNonbondedForce.h"
 #include "openmm/CustomTorsionForce.h"
 #include "openmm/GBSAOBCForce.h"
@@ -607,6 +608,43 @@ public:
     virtual double executeEnergy(ContextImpl& context) = 0;
 };
 
+/**
+ * This kernel is invoked by CustomHbondForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CalcCustomHbondForceKernel : public KernelImpl {
+public:
+    enum NonbondedMethod {
+        NoCutoff = 0,
+        CutoffNonPeriodic = 1,
+        CutoffPeriodic = 2
+    };
+    static std::string Name() {
+        return "CalcCustomHbondForce";
+    }
+    CalcCustomHbondForceKernel(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 CustomHbondForce this kernel will be used for
+     */
+    virtual void initialize(const System& system, const CustomHbondForce& force) = 0;
+    /**
+     * Execute the kernel to calculate the forces.
+     *
+     * @param context    the context in which to execute this kernel
+     */
+    virtual void executeForces(ContextImpl& context) = 0;
+    /**
+     * Execute the kernel to calculate the energy.
+     *
+     * @param context    the context in which to execute this kernel
+     * @return the potential energy due to the CustomHbondForce
+     */
+    virtual double executeEnergy(ContextImpl& context) = 0;
+};
+
 /**
  * This kernel is invoked by VerletIntegrator to take one time step.
  */

openmmapi/include/OpenMM.h

@@ -40,6 +40,7 @@
 #include "openmm/CustomTorsionForce.h"
 #include "openmm/CustomExternalForce.h"
 #include "openmm/CustomGBForce.h"
+#include "openmm/CustomHbondForce.h"
 #include "openmm/CustomNonbondedForce.h"
 #include "openmm/Force.h"
 #include "openmm/GBSAOBCForce.h"

openmmapi/include/openmm/CustomNonbondedForce.h

@@ -66,8 +66,7 @@ namespace OpenMM {
  *
  * <tt>CustomNonbondedForce* force = new CustomNonbondedForce("4*epsilon*((sigma/r)^12-(sigma/r)^6); sigma=0.5*(sigma1*sigma2); epsilon=sqrt(epsilon1*epsilon2)");</tt>
  *
- * This force depends on two parameters: sigma and epsilon.  The following code defines these parameters, and
- * specifies combining rules for them which correspond to the standard Lorentz-Bertelot combining rules:
+ * This force depends on two parameters: sigma and epsilon.  The following code defines these as per-particle parameters:
  *
  * <tt><pre>
  * force->addPerParticleParameter("sigma");

openmmapi/include/openmm/internal/CustomHbondForceImpl.h

+#ifndef OPENMM_CUSTOMHBONDFORCEIMPL_H_
+#define OPENMM_CUSTOMHBONDFORCEIMPL_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-2010 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/CustomHbondForce.h"
+#include "openmm/Kernel.h"
+#include <utility>
+#include <map>
+#include <string>
+
+namespace OpenMM {
+
+/**
+ * This is the internal implementation of CustomNonbondedForce.
+ */
+
+class CustomHbondForceImpl : public ForceImpl {
+public:
+    CustomHbondForceImpl(CustomHbondForce& owner);
+    ~CustomHbondForceImpl();
+    void initialize(ContextImpl& context);
+    CustomHbondForce& getOwner() {
+        return owner;
+    }
+    void updateContextState(ContextImpl& context) {
+        // This force field doesn't update the state directly.
+    }
+    void calcForces(ContextImpl& context);
+    double calcEnergy(ContextImpl& context);
+    std::map<std::string, double> getDefaultParameters();
+    std::vector<std::string> getKernelNames();
+private:
+    CustomHbondForce& owner;
+    Kernel kernel;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_CUSTOMHBONDFORCEIMPL_H_*/

openmmapi/src/CustomHbondForce.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-2010 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/CustomHbondForce.h"
+#include "openmm/internal/CustomHbondForceImpl.h"
+#include <cmath>
+#include <map>
+#include <sstream>
+#include <utility>
+
+using namespace OpenMM;
+using std::map;
+using std::pair;
+using std::set;
+using std::string;
+using std::stringstream;
+using std::vector;
+
+CustomHbondForce::CustomHbondForce(const string& energy) : energyExpression(energy), nonbondedMethod(NoCutoff), cutoffDistance(1.0) {
+}
+
+const string& CustomHbondForce::getEnergyFunction() const {
+    return energyExpression;
+}
+
+void CustomHbondForce::setEnergyFunction(const std::string& energy) {
+    energyExpression = energy;
+}
+
+CustomHbondForce::NonbondedMethod CustomHbondForce::getNonbondedMethod() const {
+    return nonbondedMethod;
+}
+
+void CustomHbondForce::setNonbondedMethod(NonbondedMethod method) {
+    nonbondedMethod = method;
+}
+
+double CustomHbondForce::getCutoffDistance() const {
+    return cutoffDistance;
+}
+
+void CustomHbondForce::setCutoffDistance(double distance) {
+    cutoffDistance = distance;
+}
+
+int CustomHbondForce::addPerDonorParameter(const string& name) {
+    donorParameters.push_back(PerPairParameterInfo(name));
+    return donorParameters.size()-1;
+}
+
+const string& CustomHbondForce::getPerDonorParameterName(int index) const {
+    return donorParameters[index].name;
+}
+
+void CustomHbondForce::setPerDonorParameterName(int index, const string& name) {
+    donorParameters[index].name = name;
+}
+
+int CustomHbondForce::addPerAcceptorParameter(const string& name) {
+    acceptorParameters.push_back(PerPairParameterInfo(name));
+    return acceptorParameters.size()-1;
+}
+
+const string& CustomHbondForce::getPerAcceptorParameterName(int index) const {
+    return acceptorParameters[index].name;
+}
+
+void CustomHbondForce::setPerAcceptorParameterName(int index, const string& name) {
+    acceptorParameters[index].name = name;
+}
+
+int CustomHbondForce::addGlobalParameter(const string& name, double defaultValue) {
+    globalParameters.push_back(GlobalParameterInfo(name, defaultValue));
+    return globalParameters.size()-1;
+}
+
+const string& CustomHbondForce::getGlobalParameterName(int index) const {
+    return globalParameters[index].name;
+}
+
+void CustomHbondForce::setGlobalParameterName(int index, const string& name) {
+    globalParameters[index].name = name;
+}
+
+double CustomHbondForce::getGlobalParameterDefaultValue(int index) const {
+    return globalParameters[index].defaultValue;
+}
+
+void CustomHbondForce::setGlobalParameterDefaultValue(int index, double defaultValue) {
+    globalParameters[index].defaultValue = defaultValue;
+}
+
+int CustomHbondForce::addDonor(int primaryParticle, int secondaryParticle, const vector<double>& parameters) {
+    donors.push_back(PairInfo(primaryParticle, secondaryParticle, parameters));
+    return donors.size()-1;
+}
+
+void CustomHbondForce::getDonorParameters(int index, int& primaryParticle, int& secondaryParticle, std::vector<double>& parameters) const {
+    primaryParticle = donors[index].primaryParticle;
+    secondaryParticle = donors[index].secondaryParticle;
+    parameters = donors[index].parameters;
+}
+
+void CustomHbondForce::setDonorParameters(int index, int primaryParticle, int secondaryParticle, const vector<double>& parameters) {
+    donors[index].primaryParticle = primaryParticle;
+    donors[index].secondaryParticle = secondaryParticle;
+    donors[index].parameters = parameters;
+}
+
+int CustomHbondForce::addAcceptor(int primaryParticle, int secondaryParticle, const vector<double>& parameters) {
+    acceptors.push_back(PairInfo(primaryParticle, secondaryParticle, parameters));
+    return acceptors.size()-1;
+}
+
+void CustomHbondForce::getAcceptorParameters(int index, int& primaryParticle, int& secondaryParticle, std::vector<double>& parameters) const {
+    primaryParticle = acceptors[index].primaryParticle;
+    secondaryParticle = acceptors[index].secondaryParticle;
+    parameters = acceptors[index].parameters;
+}
+
+void CustomHbondForce::setAcceptorParameters(int index, int primaryParticle, int secondaryParticle, const vector<double>& parameters) {
+    acceptors[index].primaryParticle = primaryParticle;
+    acceptors[index].secondaryParticle = secondaryParticle;
+    acceptors[index].parameters = parameters;
+}
+
+int CustomHbondForce::addExclusion(int donor, int acceptor) {
+    exclusions.push_back(ExclusionInfo(donor, acceptor));
+    return exclusions.size()-1;
+}
+void CustomHbondForce::getExclusionParticles(int index, int& donor, int& acceptor) const {
+    donor = exclusions[index].donor;
+    acceptor = exclusions[index].acceptor;
+}
+
+void CustomHbondForce::setExclusionParticles(int index, int donor, int acceptor) {
+    exclusions[index].donor = donor;
+    exclusions[index].acceptor = acceptor;
+}
+
+int CustomHbondForce::addFunction(const std::string& name, const std::vector<double>& values, double min, double max, bool interpolating) {
+    if (max <= min)
+        throw OpenMMException("CustomHbondForce: max <= min for a tabulated function.");
+    if (values.size() < 2)
+        throw OpenMMException("CustomHbondForce: a tabulated function must have at least two points");
+    functions.push_back(FunctionInfo(name, values, min, max, interpolating));
+    return functions.size()-1;
+}
+
+void CustomHbondForce::getFunctionParameters(int index, std::string& name, std::vector<double>& values, double& min, double& max, bool& interpolating) const {
+    name = functions[index].name;
+    values = functions[index].values;
+    min = functions[index].min;
+    max = functions[index].max;
+    interpolating = functions[index].interpolating;
+}
+
+void CustomHbondForce::setFunctionParameters(int index, const std::string& name, const std::vector<double>& values, double min, double max, bool interpolating) {
+    if (max <= min)
+        throw OpenMMException("CustomHbondForce: max <= min for a tabulated function.");
+    if (values.size() < 2)
+        throw OpenMMException("CustomHbondForce: a tabulated function must have at least two points");
+    functions[index].name = name;
+    functions[index].values = values;
+    functions[index].min = min;
+    functions[index].max = max;
+    functions[index].interpolating = interpolating;
+}
+
+ForceImpl* CustomHbondForce::createImpl() {
+    return new CustomHbondForceImpl(*this);
+}
+

openmmapi/src/CustomHbondForceImpl.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-2010 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/CustomHbondForceImpl.h"
+#include "openmm/kernels.h"
+#include <sstream>
+
+using namespace OpenMM;
+using std::map;
+using std::pair;
+using std::vector;
+using std::set;
+using std::string;
+using std::stringstream;
+
+CustomHbondForceImpl::CustomHbondForceImpl(CustomHbondForce& owner) : owner(owner) {
+}
+
+CustomHbondForceImpl::~CustomHbondForceImpl() {
+}
+
+void CustomHbondForceImpl::initialize(ContextImpl& context) {
+    kernel = context.getPlatform().createKernel(CalcCustomHbondForceKernel::Name(), context);
+
+    // Check for errors in the specification of parameters and exclusions.
+
+    System& system = context.getSystem();
+    vector<set<int> > exclusions(owner.getNumDonors());
+    vector<double> parameters;
+    int numDonorParameters = owner.getNumPerDonorParameters();
+    for (int i = 0; i < owner.getNumDonors(); i++) {
+        int primary, secondary;
+        owner.getDonorParameters(i, primary, secondary, parameters);
+        if (primary < 0 || primary >= system.getNumParticles()) {
+            stringstream msg;
+            msg << "CustomHbondForce: Illegal particle index for a donor: ";
+            msg << primary;
+            throw OpenMMException(msg.str());
+        }
+        if (secondary < 0 || secondary >= system.getNumParticles()) {
+            stringstream msg;
+            msg << "CustomHbondForce: Illegal particle index for a donor: ";
+            msg << secondary;
+            throw OpenMMException(msg.str());
+        }
+        if (parameters.size() != numDonorParameters) {
+            stringstream msg;
+            msg << "CustomHbondForce: Wrong number of parameters for donor ";
+            msg << i;
+            throw OpenMMException(msg.str());
+        }
+    }
+    int numAcceptorParameters = owner.getNumPerAcceptorParameters();
+    for (int i = 0; i < owner.getNumAcceptors(); i++) {
+        int primary, secondary;
+        owner.getAcceptorParameters(i, primary, secondary, parameters);
+        if (primary < 0 || primary >= system.getNumParticles()) {
+            stringstream msg;
+            msg << "CustomHbondForce: Illegal particle index for an acceptor: ";
+            msg << primary;
+            throw OpenMMException(msg.str());
+        }
+        if (secondary < 0 || secondary >= system.getNumParticles()) {
+            stringstream msg;
+            msg << "CustomHbondForce: Illegal particle index for an acceptor: ";
+            msg << secondary;
+            throw OpenMMException(msg.str());
+        }
+        if (parameters.size() != numAcceptorParameters) {
+            stringstream msg;
+            msg << "CustomHbondForce: Wrong number of parameters for acceptor ";
+            msg << i;
+            throw OpenMMException(msg.str());
+        }
+    }
+    for (int i = 0; i < owner.getNumExclusions(); i++) {
+        int donor, acceptor;
+        owner.getExclusionParticles(i, donor, acceptor);
+        if (donor < 0 || donor >= owner.getNumDonors()) {
+            stringstream msg;
+            msg << "CustomHbondForce: Illegal donor index for an exclusion: ";
+            msg << donor;
+            throw OpenMMException(msg.str());
+        }
+        if (acceptor < 0 || acceptor >= owner.getNumAcceptors()) {
+            stringstream msg;
+            msg << "CustomHbondForce: Illegal acceptor index for an exclusion: ";
+            msg << acceptor;
+            throw OpenMMException(msg.str());
+        }
+        if (exclusions[donor].count(acceptor) > 0) {
+            stringstream msg;
+            msg << "CustomHbondForce: Multiple exclusions are specified for donor ";
+            msg << donor;
+            msg << " and acceptor ";
+            msg << acceptor;
+            throw OpenMMException(msg.str());
+        }
+        exclusions[donor].insert(acceptor);
+    }
+    if (owner.getNonbondedMethod() == CustomHbondForce::CutoffPeriodic) {
+        Vec3 boxVectors[3];
+        system.getPeriodicBoxVectors(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("CustomHbondForce: The cutoff distance cannot be greater than half the periodic box size.");
+    }
+    dynamic_cast<CalcCustomHbondForceKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner);
+}
+
+void CustomHbondForceImpl::calcForces(ContextImpl& context) {
+    dynamic_cast<CalcCustomHbondForceKernel&>(kernel.getImpl()).executeForces(context);
+}
+
+double CustomHbondForceImpl::calcEnergy(ContextImpl& context) {
+    return dynamic_cast<CalcCustomHbondForceKernel&>(kernel.getImpl()).executeEnergy(context);
+}
+
+vector<string> CustomHbondForceImpl::getKernelNames() {
+    vector<string> names;
+    names.push_back(CalcCustomHbondForceKernel::Name());
+    return names;
+}
+
+map<string, double> CustomHbondForceImpl::getDefaultParameters() {
+    map<string, double> parameters;
+    for (int i = 0; i < owner.getNumGlobalParameters(); i++)
+        parameters[owner.getGlobalParameterName(i)] = owner.getGlobalParameterDefaultValue(i);
+    return parameters;
+}

platforms/reference/src/ReferenceKernelFactory.cpp

@@ -70,6 +70,8 @@ KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Pla
         return new ReferenceCalcCustomGBForceKernel(name, platform);
     if (name == CalcCustomExternalForceKernel::Name())
         return new ReferenceCalcCustomExternalForceKernel(name, platform);
+    if (name == CalcCustomHbondForceKernel::Name())
+        return new ReferenceCalcCustomHbondForceKernel(name, platform);
     if (name == IntegrateVerletStepKernel::Name())
         return new ReferenceIntegrateVerletStepKernel(name, platform, data);
     if (name == IntegrateLangevinStepKernel::Name())

platforms/reference/src/ReferenceKernels.cpp

@@ -41,6 +41,7 @@
 #include "SimTKReference/ReferenceCustomBondIxn.h"
 #include "SimTKReference/ReferenceCustomExternalIxn.h"
 #include "SimTKReference/ReferenceCustomGBIxn.h"
+#include "SimTKReference/ReferenceCustomHbondIxn.h"
 #include "SimTKReference/ReferenceCustomNonbondedIxn.h"
 #include "SimTKReference/ReferenceCustomTorsionIxn.h"
 #include "SimTKReference/ReferenceHarmonicBondIxn.h"
@@ -1233,6 +1234,130 @@ double ReferenceCalcCustomExternalForceKernel::executeEnergy(ContextImpl& contex
     return energy;
 }
 
+ReferenceCalcCustomHbondForceKernel::~ReferenceCalcCustomHbondForceKernel() {
+    disposeRealArray(donorParamArray, numDonors);
+    disposeRealArray(acceptorParamArray, numAcceptors);
+    disposeIntArray(exclusionArray, numDonors);
+}
+
+void ReferenceCalcCustomHbondForceKernel::initialize(const System& system, const CustomHbondForce& force) {
+
+    // Record the exclusions.
+
+    numDonors = force.getNumDonors();
+    numAcceptors = force.getNumAcceptors();
+    numParticles = system.getNumParticles();
+    exclusions.resize(numDonors);
+    for (int i = 0; i < force.getNumExclusions(); i++) {
+        int donor, acceptor;
+        force.getExclusionParticles(i, donor, acceptor);
+        exclusions[donor].insert(acceptor);
+    }
+
+    // Build the arrays.
+
+    donorParticles.resize(numDonors);
+    int numDonorParameters = force.getNumPerDonorParameters();
+    donorParamArray = allocateRealArray(numDonors, numDonorParameters);
+    for (int i = 0; i < numDonors; ++i) {
+        vector<double> parameters;
+        force.getDonorParameters(i, donorParticles[i].first, donorParticles[i].second, parameters);
+        for (int j = 0; j < numDonorParameters; j++)
+            donorParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
+    }
+    acceptorParticles.resize(numAcceptors);
+    int numAcceptorParameters = force.getNumPerAcceptorParameters();
+    acceptorParamArray = allocateRealArray(numAcceptors, numAcceptorParameters);
+    for (int i = 0; i < numAcceptors; ++i) {
+        vector<double> parameters;
+        force.getAcceptorParameters(i, acceptorParticles[i].first, acceptorParticles[i].second, parameters);
+        for (int j = 0; j < numAcceptorParameters; j++)
+            acceptorParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
+    }
+    exclusionArray = new int*[numDonors];
+    for (int i = 0; i < numDonors; ++i) {
+        exclusionArray[i] = new int[exclusions[i].size()+1];
+        exclusionArray[i][0] = exclusions[i].size();
+        int index = 0;
+        for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter)
+            exclusionArray[i][++index] = *iter;
+    }
+    nonbondedMethod = CalcCustomHbondForceKernel::NonbondedMethod(force.getNonbondedMethod());
+    nonbondedCutoff = (RealOpenMM) force.getCutoffDistance();
+    Vec3 boxVectors[3];
+    system.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
+    periodicBoxSize[0] = (RealOpenMM) boxVectors[0][0];
+    periodicBoxSize[1] = (RealOpenMM) boxVectors[1][1];
+    periodicBoxSize[2] = (RealOpenMM) boxVectors[2][2];
+
+    // Create custom functions for the tabulated functions.
+
+    map<string, Lepton::CustomFunction*> functions;
+    for (int i = 0; i < force.getNumFunctions(); i++) {
+        string name;
+        vector<double> values;
+        double min, max;
+        bool interpolating;
+        force.getFunctionParameters(i, name, values, min, max, interpolating);
+        functions[name] = new ReferenceTabulatedFunction(min, max, values, interpolating);
+    }
+
+    // Parse the various expressions used to calculate the force.
+
+    Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction(), functions).optimize();
+    energyExpression = expression.createProgram();
+    rForceExpression = expression.differentiate("r").optimize().createProgram();
+    thetaForceExpression = expression.differentiate("theta").optimize().createProgram();
+    psiForceExpression = expression.differentiate("psi").optimize().createProgram();
+    chiForceExpression = expression.differentiate("chi").optimize().createProgram();
+    for (int i = 0; i < numDonorParameters; i++)
+        donorParameterNames.push_back(force.getPerDonorParameterName(i));
+    for (int i = 0; i < numAcceptorParameters; i++)
+        acceptorParameterNames.push_back(force.getPerAcceptorParameterName(i));
+    for (int i = 0; i < force.getNumGlobalParameters(); i++)
+        globalParameterNames.push_back(force.getGlobalParameterName(i));
+
+    // Delete the custom functions.
+
+    for (map<string, Lepton::CustomFunction*>::iterator iter = functions.begin(); iter != functions.end(); iter++)
+        delete iter->second;
+}
+
+void ReferenceCalcCustomHbondForceKernel::executeForces(ContextImpl& context) {
+    RealOpenMM** posData = extractPositions(context);
+    RealOpenMM** forceData = extractForces(context);
+    ReferenceCustomHbondIxn ixn(donorParticles, acceptorParticles, energyExpression, rForceExpression, thetaForceExpression, psiForceExpression,
+            chiForceExpression, donorParameterNames, acceptorParameterNames);
+    bool periodic = (nonbondedMethod == CutoffPeriodic);
+    if (nonbondedMethod != NoCutoff)
+        ixn.setUseCutoff(nonbondedCutoff);
+    if (periodic)
+        ixn.setPeriodic(periodicBoxSize);
+    map<string, double> globalParameters;
+    for (int i = 0; i < (int) globalParameterNames.size(); i++)
+        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
+    ixn.calculatePairIxn(posData, donorParamArray, acceptorParamArray, exclusionArray, globalParameters, forceData, 0);
+}
+
+double ReferenceCalcCustomHbondForceKernel::executeEnergy(ContextImpl& context) {
+    RealOpenMM** posData = extractPositions(context);
+    RealOpenMM** forceData = allocateRealArray(numParticles, 3);
+    RealOpenMM energy = 0;
+    ReferenceCustomHbondIxn ixn(donorParticles, acceptorParticles, energyExpression, rForceExpression, thetaForceExpression, psiForceExpression,
+            chiForceExpression, donorParameterNames, acceptorParameterNames);
+    bool periodic = (nonbondedMethod == CutoffPeriodic);
+    if (nonbondedMethod != NoCutoff)
+        ixn.setUseCutoff(nonbondedCutoff);
+    if (periodic)
+        ixn.setPeriodic(periodicBoxSize);
+    map<string, double> globalParameters;
+    for (int i = 0; i < (int) globalParameterNames.size(); i++)
+        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
+    ixn.calculatePairIxn(posData, donorParamArray, acceptorParamArray, exclusionArray, globalParameters, forceData, &energy);
+    disposeRealArray(forceData, numParticles);
+    return energy;
+}
+
 ReferenceIntegrateVerletStepKernel::~ReferenceIntegrateVerletStepKernel() {
     if (dynamics)
         delete dynamics;

platforms/reference/src/ReferenceKernels.h

@@ -628,6 +628,46 @@ private:
     std::vector<std::string> parameterNames, globalParameterNames;
 };
 
+/**
+ * This kernel is invoked by CustomHbondForce to calculate the forces acting on the system.
+ */
+class ReferenceCalcCustomHbondForceKernel : public CalcCustomHbondForceKernel {
+public:
+    ReferenceCalcCustomHbondForceKernel(std::string name, const Platform& platform) : CalcCustomHbondForceKernel(name, platform) {
+    }
+    ~ReferenceCalcCustomHbondForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CustomHbondForce this kernel will be used for
+     */
+    void initialize(const System& system, const CustomHbondForce& force);
+    /**
+     * Execute the kernel to calculate the forces.
+     *
+     * @param context    the context in which to execute this kernel
+     */
+    void executeForces(ContextImpl& context);
+    /**
+     * Execute the kernel to calculate the energy.
+     *
+     * @param context    the context in which to execute this kernel
+     * @return the potential energy due to the CustomHbondForce
+     */
+    double executeEnergy(ContextImpl& context);
+private:
+    int numDonors, numAcceptors, numParticles;
+    int **exclusionArray;
+    RealOpenMM **donorParamArray, **acceptorParamArray;
+    RealOpenMM nonbondedCutoff, periodicBoxSize[3];
+    std::vector<std::set<int> > exclusions;
+    std::vector<std::pair<int, int> > donorParticles, acceptorParticles;
+    Lepton::ExpressionProgram energyExpression, rForceExpression, thetaForceExpression, psiForceExpression, chiForceExpression;
+    std::vector<std::string> donorParameterNames, acceptorParameterNames, globalParameterNames;
+    NonbondedMethod nonbondedMethod;
+};
+
 /**
  * This kernel is invoked by VerletIntegrator to take one time step.
  */

platforms/reference/src/ReferencePlatform.cpp

@@ -55,6 +55,7 @@ ReferencePlatform::ReferencePlatform() {
     registerKernelFactory(CalcGBVIForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomGBForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomExternalForceKernel::Name(), factory);
+    registerKernelFactory(CalcCustomHbondForceKernel::Name(), factory);
     registerKernelFactory(IntegrateVerletStepKernel::Name(), factory);
     registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);
     registerKernelFactory(IntegrateBrownianStepKernel::Name(), factory);

platforms/reference/src/SimTKReference/ReferenceCustomHbondIxn.cpp

+
+/* Portions copyright (c) 2009-2010 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 "../SimTKUtilities/SimTKOpenMMCommon.h"
+#include "../SimTKUtilities/SimTKOpenMMLog.h"
+#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
+#include "ReferenceForce.h"
+#include "ReferenceCustomHbondIxn.h"
+
+using std::map;
+using std::pair;
+using std::string;
+using std::stringstream;
+using std::vector;
+
+/**---------------------------------------------------------------------------------------
+
+   ReferenceCustomHbondIxn constructor
+
+   --------------------------------------------------------------------------------------- */
+
+ReferenceCustomHbondIxn::ReferenceCustomHbondIxn(const vector<pair<int, int> >& donorAtoms, const vector<pair<int, int> >& acceptorAtoms,
+            const Lepton::ExpressionProgram& energyExpression,
+            const Lepton::ExpressionProgram& rForceExpression, const Lepton::ExpressionProgram& thetaForceExpression,
+            const Lepton::ExpressionProgram& psiForceExpression, const Lepton::ExpressionProgram& chiForceExpression,
+            const vector<string>& donorParameterNames, const vector<string>& acceptorParameterNames) :
+            cutoff(false), periodic(false), donorAtoms(donorAtoms), acceptorAtoms(acceptorAtoms), energyExpression(energyExpression),
+            rForceExpression(rForceExpression), thetaForceExpression(thetaForceExpression), psiForceExpression(psiForceExpression),
+            chiForceExpression(chiForceExpression), donorParamNames(donorParameterNames), acceptorParamNames(acceptorParameterNames) {
+}
+
+/**---------------------------------------------------------------------------------------
+
+   ReferenceCustomHbondIxn destructor
+
+   --------------------------------------------------------------------------------------- */
+
+ReferenceCustomHbondIxn::~ReferenceCustomHbondIxn( ){
+}
+
+  /**---------------------------------------------------------------------------------------
+
+     Set the force to use a cutoff.
+
+     @param distance            the cutoff distance
+
+     --------------------------------------------------------------------------------------- */
+
+void ReferenceCustomHbondIxn::setUseCutoff(RealOpenMM distance) {
+    cutoff = true;
+    cutoffDistance = distance;
+}
+
+  /**---------------------------------------------------------------------------------------
+
+     Set the force to use periodic boundary conditions.  This requires that a cutoff has
+     also 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 ReferenceCustomHbondIxn::setPeriodic(RealOpenMM* boxSize) {
+    assert(cutoff);
+    assert(boxSize[0] >= 2.0*cutoffDistance);
+    assert(boxSize[1] >= 2.0*cutoffDistance);
+    assert(boxSize[2] >= 2.0*cutoffDistance);
+    periodic = true;
+    periodicBoxSize[0] = boxSize[0];
+    periodicBoxSize[1] = boxSize[1];
+    periodicBoxSize[2] = boxSize[2];
+  }
+
+
+/**---------------------------------------------------------------------------------------
+
+   Calculate custom hbond interaction
+
+   @param atomCoordinates    atom coordinates
+   @param donorParameters    donor parameters values       donorParameters[donorIndex][parameterIndex]
+   @param acceptorParameters acceptor parameters values    acceptorParameters[acceptorIndex][parameterIndex]
+   @param exclusions         exclusion indices             exclusions[donorIndex][acceptorToExcludeIndex]
+                             exclusions[donorIndex][0] = number of exclusions
+                             exclusions[donorIndex][no.-1] = indices of acceptors to excluded from
+                             interacting w/ donor donorIndex
+   @param globalParameters   the values of global parameters
+   @param forces             force array (forces added)
+   @param totalEnergy        total energy
+
+   --------------------------------------------------------------------------------------- */
+
+void ReferenceCustomHbondIxn::calculatePairIxn(RealOpenMM** atomCoordinates, RealOpenMM** donorParameters, RealOpenMM** acceptorParameters,
+                                             int** exclusions, const map<string, double>& globalParameters, RealOpenMM** forces,
+                                             RealOpenMM* totalEnergy) const {
+
+   map<string, double> variables = globalParameters;
+
+   // allocate and initialize exclusion array
+
+   int numDonors = donorAtoms.size();
+   int numAcceptors = acceptorAtoms.size();
+   int* exclusionIndices = new int[numAcceptors];
+   for( int ii = 0; ii < numAcceptors; ii++ ){
+      exclusionIndices[ii] = -1;
+   }
+
+   for( int donor = 0; donor < numDonors; donor++ ){
+
+      // set exclusions
+
+      for (int j = 1; j <= exclusions[donor][0]; j++)
+         exclusionIndices[exclusions[donor][j]] = donor;
+
+      // Initialize per-donor parameters.
+
+      for (int j = 0; j < (int) donorParamNames.size(); j++)
+          variables[donorParamNames[j]] = donorParameters[donor][j];
+
+      // loop over atom pairs
+
+      for( int acceptor = 0; acceptor < numAcceptors; acceptor++ ){
+
+         if( exclusionIndices[acceptor] != donor ){
+             for (int j = 0; j < (int) acceptorParamNames.size(); j++)
+                 variables[acceptorParamNames[j]] = acceptorParameters[acceptor][j];
+             calculateOneIxn(donor, acceptor, atomCoordinates, variables, forces, totalEnergy);
+         }
+      }
+   }
+
+   delete[] exclusionIndices;
+}
+
+  /**---------------------------------------------------------------------------------------
+
+     Calculate custom interaction between a donor and an acceptor
+
+     @param donor            the index of the donor
+     @param acceptor         the index of the acceptor
+     @param atomCoordinates  atom coordinates
+     @param variables        the values of variables that may appear in expressions
+     @param forces           force array (forces added)
+     @param energyByAtom     atom energy
+     @param totalEnergy      total energy
+
+     --------------------------------------------------------------------------------------- */
+
+void ReferenceCustomHbondIxn::calculateOneIxn(int donor, int acceptor, RealOpenMM** atomCoordinates,
+                        map<string, double>& variables, RealOpenMM** forces, RealOpenMM* totalEnergy) const {
+
+    // ---------------------------------------------------------------------------------------
+
+    static const std::string methodName = "\nReferenceCustomHbondIxn::calculateOneIxn";
+
+    // ---------------------------------------------------------------------------------------
+
+    // Visual Studio complains if crossProduct declared as 'crossProduct[2][3]'
+
+    RealOpenMM crossProductMemory[6];
+    RealOpenMM* crossProduct[2];
+    crossProduct[0] = crossProductMemory;
+    crossProduct[1] = crossProductMemory + 3;
+
+    // Compute the distance between the primary donor and acceptor atoms, and compare to the cutoff.
+
+    int d1 = donorAtoms[donor].first;
+    int a1 = acceptorAtoms[acceptor].first;
+    RealOpenMM deltaD1A1[ReferenceForce::LastDeltaRIndex];
+    if (periodic)
+        ReferenceForce::getDeltaRPeriodic(atomCoordinates[a1], atomCoordinates[d1], periodicBoxSize, deltaD1A1);
+    else
+        ReferenceForce::getDeltaR(atomCoordinates[a1], atomCoordinates[d1], deltaD1A1);
+    if (cutoff && deltaD1A1[ReferenceForce::RIndex] >= cutoffDistance)
+        return;
+
+    // Compute all of the variables the energy can depend on.
+    
+    int d2 = donorAtoms[donor].second;
+    int a2 = acceptorAtoms[acceptor].second;
+    RealOpenMM deltaD1D2[ReferenceForce::LastDeltaRIndex];
+    RealOpenMM deltaA2A1[ReferenceForce::LastDeltaRIndex];
+    if (periodic) {
+        ReferenceForce::getDeltaRPeriodic(atomCoordinates[d2], atomCoordinates[d1], periodicBoxSize, deltaD1D2);
+        ReferenceForce::getDeltaRPeriodic(atomCoordinates[a1], atomCoordinates[a2], periodicBoxSize, deltaA2A1);
+    }
+    else {
+        ReferenceForce::getDeltaR(atomCoordinates[d2], atomCoordinates[d1], deltaD1D2);
+        ReferenceForce::getDeltaR(atomCoordinates[a1], atomCoordinates[a2], deltaA2A1);
+    }
+    variables["r"] = deltaD1A1[ReferenceForce::RIndex];
+    variables["theta"] = computeAngle(deltaD1A1, deltaD1D2, 1);
+    variables["psi"] = computeAngle(deltaD1A1, deltaA2A1, 1);
+    RealOpenMM dotDihedral, signOfDihedral;
+    variables["chi"] =  getDihedralAngleBetweenThreeVectors(deltaA2A1, deltaD1A1, deltaD1D2,
+                                  crossProduct, &dotDihedral, deltaA2A1, &signOfDihedral, 1);
+
+    // Apply forces based on r.
+
+    RealOpenMM dEdR = (RealOpenMM) (rForceExpression.evaluate(variables)/(deltaD1A1[ReferenceForce::RIndex]));
+    if (dEdR != 0) {
+        for (int i = 0; i < 3; i++) {
+           RealOpenMM force  = -dEdR*deltaD1A1[i];
+           forces[d1][i] += force;
+           forces[a1][i] -= force;
+        }
+    }
+
+    // Apply forces based on theta.
+
+    RealOpenMM dEdTheta = (RealOpenMM) thetaForceExpression.evaluate(variables);
+    if (dEdTheta != 0) {
+        RealOpenMM thetaCross[ReferenceForce::LastDeltaRIndex];
+        SimTKOpenMMUtilities::crossProductVector3(deltaD1D2, deltaD1A1, thetaCross);
+        RealOpenMM lengthThetaCross = SQRT(DOT3(thetaCross, thetaCross));
+        if (lengthThetaCross < 1.0e-06)
+            lengthThetaCross = (RealOpenMM) 1.0e-06;
+        RealOpenMM termA = dEdTheta/(deltaD1D2[ReferenceForce::R2Index]*lengthThetaCross);
+        RealOpenMM termC = -dEdTheta/(deltaD1A1[ReferenceForce::R2Index]*lengthThetaCross);
+        RealOpenMM deltaCrossP[3][3];
+        SimTKOpenMMUtilities::crossProductVector3(deltaD1D2, thetaCross, deltaCrossP[0]);
+        SimTKOpenMMUtilities::crossProductVector3(deltaD1A1, 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[d2][i] += deltaCrossP[0][i];
+            forces[d1][i] += deltaCrossP[1][i];
+            forces[a1][i] += deltaCrossP[2][i];
+        }
+    }
+
+    // Apply forces based on psi.
+
+    RealOpenMM dEdPsi = (RealOpenMM) psiForceExpression.evaluate(variables);
+    if (dEdPsi != 0) {
+        RealOpenMM psiCross[ReferenceForce::LastDeltaRIndex];
+        SimTKOpenMMUtilities::crossProductVector3(deltaA2A1, deltaD1A1, psiCross);
+        RealOpenMM lengthPsiCross = SQRT(DOT3(psiCross, psiCross));
+        if (lengthPsiCross < 1.0e-06)
+            lengthPsiCross = (RealOpenMM) 1.0e-06;
+        RealOpenMM termA =  dEdPsi/(deltaD1A1[ReferenceForce::R2Index]*lengthPsiCross);
+        RealOpenMM termC = -dEdPsi/(deltaA2A1[ReferenceForce::R2Index]*lengthPsiCross);
+        RealOpenMM deltaCrossP[3][3];
+        SimTKOpenMMUtilities::crossProductVector3(deltaD1A1, psiCross, deltaCrossP[0]);
+        SimTKOpenMMUtilities::crossProductVector3(deltaA2A1, psiCross, 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[d1][i] += deltaCrossP[0][i];
+            forces[a1][i] += deltaCrossP[1][i];
+            forces[a2][i] += deltaCrossP[2][i];
+        }
+    }
+
+    // Apply forces based on chi.
+
+    RealOpenMM dEdChi = (RealOpenMM) chiForceExpression.evaluate(variables);
+    if (dEdChi != 0) {
+        RealOpenMM internalF[4][3];
+        RealOpenMM forceFactors[4];
+        RealOpenMM normCross1 = DOT3(crossProduct[0], crossProduct[0]);
+        RealOpenMM normBC = deltaD1A1[ReferenceForce::RIndex];
+        forceFactors[0] = (-dEdChi*normBC)/normCross1;
+        RealOpenMM normCross2 = DOT3(crossProduct[1], crossProduct[1]);
+                   forceFactors[3] = (dEdChi*normBC)/normCross2;
+                   forceFactors[1] = DOT3(deltaA2A1, deltaD1A1);
+                   forceFactors[1] /= deltaD1A1[ReferenceForce::R2Index];
+                   forceFactors[2] = DOT3(deltaD1D2, deltaD1A1);
+                   forceFactors[2] /= deltaD1A1[ReferenceForce::R2Index];
+        for (int i = 0; i < 3; i++) {
+            internalF[0][i] = forceFactors[0]*crossProduct[0][i];
+            internalF[3][i] = forceFactors[3]*crossProduct[1][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[a2][i] += internalF[0][i];
+            forces[a1][i] -= internalF[1][i];
+            forces[d1][i] -= internalF[2][i];
+            forces[d2][i] += internalF[3][i];
+        }
+    }
+
+    // Add the energy
+
+    if (totalEnergy)
+        *totalEnergy += (RealOpenMM) energyExpression.evaluate(variables);
+}
+
+RealOpenMM ReferenceCustomHbondIxn::computeAngle(RealOpenMM* vec1, RealOpenMM* vec2, RealOpenMM sign) {
+    RealOpenMM dot = sign*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/src/SimTKReference/ReferenceCustomHbondIxn.h

+
+/* Portions copyright (c) 2009-2010 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 __ReferenceCustomHbondIxn_H__
+#define __ReferenceCustomHbondIxn_H__
+
+#include "ReferenceBondIxn.h"
+#include "lepton/ExpressionProgram.h"
+#include <map>
+#include <vector>
+
+// ---------------------------------------------------------------------------------------
+
+class ReferenceCustomHbondIxn : public ReferenceBondIxn {
+
+   private:
+
+      bool cutoff;
+      bool periodic;
+      RealOpenMM periodicBoxSize[3];
+      RealOpenMM cutoffDistance;
+      Lepton::ExpressionProgram energyExpression;
+      Lepton::ExpressionProgram rForceExpression;
+      Lepton::ExpressionProgram thetaForceExpression;
+      Lepton::ExpressionProgram psiForceExpression;
+      Lepton::ExpressionProgram chiForceExpression;
+      std::vector<std::string> donorParamNames, acceptorParamNames;
+      std::vector<std::pair<int, int> > donorAtoms, acceptorAtoms;
+
+      /**---------------------------------------------------------------------------------------
+
+         Calculate custom interaction between a donor and an acceptor
+
+         @param donor            the index of the donor
+         @param acceptor         the index of the acceptor
+         @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(int donor, int acceptor, RealOpenMM** atomCoordinates,
+                           std::map<std::string, double>& variables, RealOpenMM** forces,
+                           RealOpenMM* totalEnergy) const;
+
+      static RealOpenMM computeAngle(RealOpenMM* vec1, RealOpenMM* vec2, RealOpenMM sign);
+
+
+   public:
+
+      /**---------------------------------------------------------------------------------------
+
+         Constructor
+
+         --------------------------------------------------------------------------------------- */
+
+       ReferenceCustomHbondIxn(const std::vector<std::pair<int, int> >& donorAtoms, const std::vector<std::pair<int, int> >& acceptorAtoms,
+                               const Lepton::ExpressionProgram& energyExpression, const Lepton::ExpressionProgram& rForceExpression,
+                               const Lepton::ExpressionProgram& thetaForceExpression, const Lepton::ExpressionProgram& psiForceExpression,
+                               const Lepton::ExpressionProgram& chiForceExpression, const std::vector<std::string>& donorParameterNames,
+                               const std::vector<std::string>& acceptorParameterNames);
+
+      /**---------------------------------------------------------------------------------------
+
+         Destructor
+
+         --------------------------------------------------------------------------------------- */
+
+       ~ReferenceCustomHbondIxn();
+
+      /**---------------------------------------------------------------------------------------
+
+         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(RealOpenMM* boxSize);
+
+      /**---------------------------------------------------------------------------------------
+
+         Calculate custom hbond interaction
+
+         @param atomCoordinates    atom coordinates
+         @param donorParameters    donor parameters values       donorParameters[donorIndex][parameterIndex]
+         @param acceptorParameters acceptor parameters values    acceptorParameters[acceptorIndex][parameterIndex]
+         @param exclusions         exclusion indices             exclusions[donorIndex][acceptorToExcludeIndex]
+                                   exclusions[donorIndex][0] = number of exclusions
+                                   exclusions[donorIndex][no.-1] = indices of acceptors to excluded from
+                                   interacting w/ donor donorIndex
+         @param globalParameters   the values of global parameters
+         @param forces             force array (forces added)
+         @param totalEnergy        total energy
+
+         --------------------------------------------------------------------------------------- */
+
+      void calculatePairIxn(RealOpenMM** atomCoordinates, RealOpenMM** donorParameters, RealOpenMM** acceptorParameters,
+                            int** exclusions, const std::map<std::string, double>& globalParameters,
+                            RealOpenMM** forces, RealOpenMM* totalEnergy) const;
+
+// ---------------------------------------------------------------------------------------
+
+};
+
+#endif // __ReferenceCustomHbondIxn_H__

platforms/reference/tests/TestReferenceCustomHbondForce.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-2010 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 CustomTorsionForce.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "ReferencePlatform.h"
+#include "openmm/CustomHbondForce.h"
+#include "openmm/HarmonicAngleForce.h"
+#include "openmm/HarmonicBondForce.h"
+#include "openmm/PeriodicTorsionForce.h"
+#include "openmm/System.h"
+#include "openmm/VerletIntegrator.h"
+#include "../src/sfmt/SFMT.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+const double TOL = 1e-5;
+
+void testHbond() {
+    ReferencePlatform platform;
+
+    // Create a system using a CustomHbondForce.
+
+    System customSystem;
+    customSystem.addParticle(1.0);
+    customSystem.addParticle(1.0);
+    customSystem.addParticle(1.0);
+    customSystem.addParticle(1.0);
+    CustomHbondForce* custom = new CustomHbondForce("0.5*kr*(r-r0)^2 + 0.5*ktheta*(theta-theta0)^2 + 0.5*kpsi*(psi-psi0)^2 + kchi*(1+cos(n*chi-chi0))");
+    custom->addPerDonorParameter("r0");
+    custom->addPerDonorParameter("theta0");
+    custom->addPerDonorParameter("psi0");
+    custom->addPerAcceptorParameter("chi0");
+    custom->addPerAcceptorParameter("n");
+    custom->addGlobalParameter("kr", 0.4);
+    custom->addGlobalParameter("ktheta", 0.5);
+    custom->addGlobalParameter("kpsi", 0.6);
+    custom->addGlobalParameter("kchi", 0.7);
+    vector<double> parameters(3);
+    parameters[0] = 1.5;
+    parameters[1] = 1.7;
+    parameters[2] = 1.9;
+    custom->addDonor(1, 0, parameters);
+    parameters.resize(2);
+    parameters[0] = 2.1;
+    parameters[1] = 2;
+    custom->addAcceptor(2, 3, parameters);
+    custom->setCutoffDistance(10.0);
+    customSystem.addForce(custom);
+
+    // Create an identical system using HarmonicBondForce, HarmonicAngleForce, and PeriodicTorsionForce.
+
+    System standardSystem;
+    standardSystem.addParticle(1.0);
+    standardSystem.addParticle(1.0);
+    standardSystem.addParticle(1.0);
+    standardSystem.addParticle(1.0);
+    HarmonicBondForce* bond = new HarmonicBondForce();
+    bond->addBond(1, 2, 1.5, 0.4);
+    standardSystem.addForce(bond);
+    HarmonicAngleForce* angle = new HarmonicAngleForce();
+    angle->addAngle(0, 1, 2, 1.7, 0.5);
+    angle->addAngle(1, 2, 3, 1.9, 0.6);
+    standardSystem.addForce(angle);
+    PeriodicTorsionForce* torsion = new PeriodicTorsionForce();
+    torsion->addTorsion(0, 1, 2, 3, 2, 2.1, 0.7);;
+    standardSystem.addForce(torsion);
+
+    // Set the atoms in various positions, and verify that both systems give identical forces and energy.
+
+    init_gen_rand(0);
+    vector<Vec3> positions(4);
+    VerletIntegrator integrator1(0.01);
+    VerletIntegrator integrator2(0.01);
+    for (int i = 0; i < 10; i++) {
+        Context c1(customSystem, integrator1, platform);
+        Context c2(standardSystem, integrator2, platform);
+        for (int j = 0; j < (int) positions.size(); j++)
+            positions[j] = Vec3(2.0*genrand_real2(), 2.0*genrand_real2(), 2.0*genrand_real2());
+        c1.setPositions(positions);
+        c2.setPositions(positions);
+        State s1 = c1.getState(State::Forces | State::Energy);
+        State s2 = c2.getState(State::Forces | State::Energy);
+        for (int i = 0; i < customSystem.getNumParticles(); i++)
+            ASSERT_EQUAL_VEC(s2.getForces()[i], s1.getForces()[i], TOL);
+        ASSERT_EQUAL_TOL(s2.getPotentialEnergy(), s1.getPotentialEnergy(), TOL);
+    }
+}
+
+void testExclusions() {
+    ReferencePlatform platform;
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    CustomHbondForce* custom = new CustomHbondForce("(r-1)^2");
+    custom->addDonor(0, 1, vector<double>());
+    custom->addDonor(1, 0, vector<double>());
+    custom->addAcceptor(2, 0, vector<double>());
+    custom->addExclusion(1, 0);
+    system.addForce(custom);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(3);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(0, 2, 0);
+    positions[2] = Vec3(2, 0, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    ASSERT_EQUAL_VEC(Vec3(2, 0, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[1], TOL);
+    ASSERT_EQUAL_VEC(Vec3(-2, 0, 0), forces[2], TOL);
+    ASSERT_EQUAL_TOL(1, state.getPotentialEnergy(), TOL);
+}
+
+void testCutoff() {
+    ReferencePlatform platform;
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    CustomHbondForce* custom = new CustomHbondForce("(r-1)^2");
+    custom->addDonor(0, 1, vector<double>());
+    custom->addDonor(1, 0, vector<double>());
+    custom->addAcceptor(2, 0, vector<double>());
+    custom->setNonbondedMethod(CustomHbondForce::CutoffNonPeriodic);
+    custom->setCutoffDistance(2.5);
+    system.addForce(custom);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(3);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(0, 3, 0);
+    positions[2] = Vec3(2, 0, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    ASSERT_EQUAL_VEC(Vec3(2, 0, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[1], TOL);
+    ASSERT_EQUAL_VEC(Vec3(-2, 0, 0), forces[2], TOL);
+    ASSERT_EQUAL_TOL(1, state.getPotentialEnergy(), TOL);
+}
+
+int main() {
+    try {
+        testHbond();
+        testExclusions();
+        testCutoff();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}
+