Created CustomBondForce, along with the reference implementation

olla/include/openmm/kernels.h

@@ -35,6 +35,7 @@
 #include "openmm/AndersenThermostat.h"
 #include "openmm/BrownianIntegrator.h"
 #include "openmm/CMMotionRemover.h"
+#include "openmm/CustomBondForce.h"
 #include "openmm/CustomNonbondedForce.h"
 #include "openmm/GBSAOBCForce.h"
 #include "openmm/GBVIForce.h"
@@ -200,6 +201,38 @@ public:
     virtual double executeEnergy(ContextImpl& context) = 0;
 };
 
+/**
+ * This kernel is invoked by CustomBondForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CalcCustomBondForceKernel : public KernelImpl {
+public:
+    static std::string Name() {
+        return "CalcCustomBondForce";
+    }
+    CalcCustomBondForceKernel(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 HarmonicBondForce this kernel will be used for
+     */
+    virtual void initialize(const System& system, const CustomBondForce& 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 HarmonicBondForce
+     */
+    virtual double executeEnergy(ContextImpl& context) = 0;
+};
+
 /**
  * This kernel is invoked by HarmonicAngleForce to calculate the forces acting on the system and the energy of the system.
  */

openmmapi/include/openmm/CustomBondForce.h

+#ifndef OPENMM_CUSTOMBONDEDFORCE_H_
+#define OPENMM_CUSTOMBONDEDFORCE_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-2009 Stanford University and the Authors.      *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "Force.h"
+#include "Vec3.h"
+#include <vector>
+#include "internal/windowsExport.h"
+
+namespace OpenMM {
+
+/**
+ * This class implements bonded interactions between pairs of particles.  Unlike HarmonicBondForce, the functional form
+ * of the interaction is completely customizable, and may involve arbitrary algebraic expressions.
+ * It may depend on the distance between particles, as well as on arbitrary global and
+ * per-bond parameters.
+ *
+ * To use this class, create a CustomBondForce object, passing an algebraic expression to the constructor
+ * that defines the interaction energy between each pair of bonded particles.  The expression may depend on r, the distance
+ * between the particles, as well as on any parameters you choose.  Then call addPerBondParameter() to define per-bond
+ * parameters, and addGlobalParameter() to define global parameters.  The values of per-bond parameters are specified as
+ * part of the system definition, while values of global parameters may be modified during a simulation by calling Context::setParameter().
+ * Finally, call addBond() once for each bond.  After a bond has been added, you can modify its parameters by calling setBondParameters().
+ *
+ * As an example, the following code creates a CustomBondForce that implements a harmonic potential:
+ *
+ * <tt>CustomBondForce* force = new CustomBondForce("0.5*k*(r-r0)^2");</tt>
+ *
+ * This force depends on two parameters: the spring constant k and equilibrium distance r0.  The following code defines these parameters:
+ *
+ * <tt><pre>
+ * force->addPerBondParameter("k");
+ * force->addPerBondParameter("r0");
+ * </pre></tt>
+ *
+ * Expressions may involve the operators + (add), - (subtract), * (multiply), / (divide), and ^ (power), and the following
+ * functions: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh.  All trigonometric functions
+ * are defined in radians, and log is the natural logarithm.
+ */
+
+class OPENMM_EXPORT CustomBondForce : public Force {
+public:
+    /**
+     * Create a CustomBondForce.
+     *
+     * @param energy    an algebraic expression giving the interaction energy between two bonded particles as a function
+     *                  of r, the distance between them
+     */
+    CustomBondForce(const std::string& energy);
+    /**
+     * Get the number of bonds for which force field parameters have been defined.
+     */
+    int getNumBonds() const {
+        return bonds.size();
+    }
+    /**
+     * Get the number of per-bond parameters that the interaction depends on.
+     */
+    int getNumPerBondParameters() const {
+        return parameters.size();
+    }
+    /**
+     * Get the number of global parameters that the interaction depends on.
+     */
+    int getNumGlobalParameters() const {
+        return globalParameters.size();
+    }
+    /**
+     * Get the algebraic expression that gives the interaction energy for each bond
+     */
+    const std::string& getEnergyFunction() const;
+    /**
+     * Set the algebraic expression that gives the interaction energy for each bond
+     */
+    void setEnergyFunction(const std::string& energy);
+    /**
+     * Add a new per-bond parmeter that the interaction may depend on.
+     *
+     * @param name             the name of the parameter
+     * @return the index of the parameter that was added
+     */
+    int addPerBondParameter(const std::string& name);
+    /**
+     * Get the name of a per-bond parameter.
+     *
+     * @param index     the index of the parameter for which to get the name
+     * @return the parameter name
+     */
+    const std::string& getPerBondParameterName(int index) const;
+    /**
+     * Set the name of a per-bond parameter.
+     *
+     * @param index          the index of the parameter for which to set the name
+     * @param name           the name of the parameter
+     */
+    void setPerBondParameterName(int index, const std::string& name);
+    /**
+     * Add a new global parmeter that the interaction may depend on.
+     *
+     * @param name             the name of the parameter
+     * @param defaultValue     the default value of the parameter
+     * @return the index of the parameter that was added
+     */
+    int addGlobalParameter(const std::string& name, double defaultValue);
+    /**
+     * Get the name of a global parameter.
+     *
+     * @param index     the index of the parameter for which to get the name
+     * @return the parameter name
+     */
+    const std::string& getGlobalParameterName(int index) const;
+    /**
+     * Set the name of a global parameter.
+     *
+     * @param index          the index of the parameter for which to set the name
+     * @param name           the name of the parameter
+     */
+    void setGlobalParameterName(int index, const std::string& name);
+    /**
+     * Get the default value of a global parameter.
+     *
+     * @param index     the index of the parameter for which to get the default value
+     * @return the parameter default value
+     */
+    double getGlobalParameterDefaultValue(int index) const;
+    /**
+     * Set the default value of a global parameter.
+     *
+     * @param index          the index of the parameter for which to set the default value
+     * @param name           the default value of the parameter
+     */
+    void setGlobalParameterDefaultValue(int index, double defaultValue);
+    /**
+     * Add a bond term to the force field.
+     *
+     * @param particle1     the index of the first particle connected by the bond
+     * @param particle2     the index of the second particle connected by the bond
+     * @param parameters    the list of parameters for the new bond
+     * @return the index of the bond that was added
+     */
+    int addBond(int particle1, int particle2, const std::vector<double>& parameters);
+    /**
+     * Get the force field parameters for a bond term.
+     *
+     * @param index         the index of the bond for which to get parameters
+     * @param particle1     the index of the first particle connected by the bond
+     * @param particle2     the index of the second particle connected by the bond
+     * @param parameters    the list of parameters for the bond
+     */
+    void getBondParameters(int index, int& particle1, int& particle2, std::vector<double>& parameters) const;
+    /**
+     * Set the force field parameters for a bond term.
+     *
+     * @param index         the index of the bond for which to set parameters
+     * @param particle1     the index of the first particle connected by the bond
+     * @param particle2     the index of the second particle connected by the bond
+     * @param parameters    the list of parameters for the bond
+     */
+    void setBondParameters(int index, int particle1, int particle2, const std::vector<double>& parameters);
+protected:
+    ForceImpl* createImpl();
+private:
+    class BondInfo;
+    class BondParameterInfo;
+    class GlobalParameterInfo;
+    std::string energyExpression;
+    std::vector<BondParameterInfo> parameters;
+    std::vector<GlobalParameterInfo> globalParameters;
+    std::vector<BondInfo> bonds;
+};
+
+class CustomBondForce::BondInfo {
+public:
+    int particle1, particle2;
+    std::vector<double> parameters;
+    BondInfo() : particle1(-1), particle2(-1) {
+    }
+    BondInfo(int particle1, int particle2, const std::vector<double>& parameters) :
+            particle1(particle1), particle2(particle2), parameters(parameters) {
+    }
+};
+
+class CustomBondForce::BondParameterInfo {
+public:
+    std::string name;
+    BondParameterInfo() {
+    }
+    BondParameterInfo(const std::string& name) : name(name) {
+    }
+};
+
+class CustomBondForce::GlobalParameterInfo {
+public:
+    std::string name;
+    double defaultValue;
+    GlobalParameterInfo() {
+    }
+    GlobalParameterInfo(const std::string& name, double defaultValue) : name(name), defaultValue(defaultValue) {
+    }
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_CUSTOMBONDEDFORCE_H_*/

openmmapi/include/openmm/internal/CustomBondForceImpl.h

+#ifndef OPENMM_CUSTOMBONDEDFORCEIMPL_H_
+#define OPENMM_CUSTOMBONDEDFORCEIMPL_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-2009 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/CustomBondForce.h"
+#include "openmm/Kernel.h"
+#include <utility>
+#include <map>
+#include <string>
+
+namespace OpenMM {
+
+/**
+ * This is the internal implementation of CustomBondForce.
+ */
+
+class CustomBondForceImpl : public ForceImpl {
+public:
+    CustomBondForceImpl(CustomBondForce& owner);
+    ~CustomBondForceImpl();
+    void initialize(ContextImpl& context);
+    CustomBondForce& 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:
+    CustomBondForce& owner;
+    Kernel kernel;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_CUSTOMBONDEDFORCEIMPL_H_*/

openmmapi/src/CustomBondForce.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-2009 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/CustomBondForce.h"
+#include "openmm/internal/CustomBondForceImpl.h"
+#include <cmath>
+#include <map>
+#include <sstream>
+#include <utility>
+
+using namespace OpenMM;
+using std::string;
+using std::stringstream;
+using std::vector;
+
+CustomBondForce::CustomBondForce(const string& energy) : energyExpression(energy) {
+}
+
+const string& CustomBondForce::getEnergyFunction() const {
+    return energyExpression;
+}
+
+void CustomBondForce::setEnergyFunction(const std::string& energy) {
+    energyExpression = energy;
+}
+
+int CustomBondForce::addPerBondParameter(const string& name) {
+    parameters.push_back(BondParameterInfo(name));
+    return parameters.size()-1;
+}
+
+const string& CustomBondForce::getPerBondParameterName(int index) const {
+    return parameters[index].name;
+}
+
+void CustomBondForce::setPerBondParameterName(int index, const string& name) {
+    parameters[index].name = name;
+}
+
+int CustomBondForce::addGlobalParameter(const string& name, double defaultValue) {
+    globalParameters.push_back(GlobalParameterInfo(name, defaultValue));
+    return globalParameters.size()-1;
+}
+
+const string& CustomBondForce::getGlobalParameterName(int index) const {
+    return globalParameters[index].name;
+}
+
+void CustomBondForce::setGlobalParameterName(int index, const string& name) {
+    globalParameters[index].name = name;
+}
+
+double CustomBondForce::getGlobalParameterDefaultValue(int index) const {
+    return globalParameters[index].defaultValue;
+}
+
+void CustomBondForce::setGlobalParameterDefaultValue(int index, double defaultValue) {
+    globalParameters[index].defaultValue = defaultValue;
+}
+
+int CustomBondForce::addBond(int particle1, int particle2, const vector<double>& parameters) {
+    bonds.push_back(BondInfo(particle1, particle2, parameters));
+    return bonds.size()-1;
+}
+
+void CustomBondForce::getBondParameters(int index, int& particle1, int& particle2, std::vector<double>& parameters) const {
+    particle1 = bonds[index].particle1;
+    particle2 = bonds[index].particle2;
+    parameters = bonds[index].parameters;
+}
+
+void CustomBondForce::setBondParameters(int index, int particle1, int particle2, const vector<double>& parameters) {
+    bonds[index].parameters = parameters;
+    bonds[index].particle1 = particle1;
+    bonds[index].particle2 = particle2;
+}
+
+ForceImpl* CustomBondForce::createImpl() {
+    return new CustomBondForceImpl(*this);
+}
+

openmmapi/src/CustomBondForceImpl.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-2009 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/CustomBondForceImpl.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;
+
+CustomBondForceImpl::CustomBondForceImpl(CustomBondForce& owner) : owner(owner) {
+}
+
+CustomBondForceImpl::~CustomBondForceImpl() {
+}
+
+void CustomBondForceImpl::initialize(ContextImpl& context) {
+    kernel = context.getPlatform().createKernel(CalcCustomBondForceKernel::Name(), context);
+
+    // Check for errors in the specification of bonds.
+
+    System& system = context.getSystem();
+    vector<double> parameters;
+    int numParameters = owner.getNumPerBondParameters();
+    for (int i = 0; i < owner.getNumBonds(); i++) {
+        int particle1, particle2;
+        owner.getBondParameters(i, particle1, particle2, parameters);
+        if (particle1 < 0 || particle1 >= system.getNumParticles()) {
+            stringstream msg;
+            msg << "CustomBondForce: Illegal particle index for a bond: ";
+            msg << particle1;
+            throw OpenMMException(msg.str());
+        }
+        if (particle2 < 0 || particle2 >= system.getNumParticles()) {
+            stringstream msg;
+            msg << "CustomBondForce: Illegal particle index for a bond: ";
+            msg << particle2;
+            throw OpenMMException(msg.str());
+        }
+        if (parameters.size() != numParameters) {
+            stringstream msg;
+            msg << "CustomBondForce: Wrong number of parameters for bond ";
+            msg << i;
+            throw OpenMMException(msg.str());
+        }
+    }
+    dynamic_cast<CalcCustomBondForceKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner);
+}
+
+void CustomBondForceImpl::calcForces(ContextImpl& context) {
+    dynamic_cast<CalcCustomBondForceKernel&>(kernel.getImpl()).executeForces(context);
+}
+
+double CustomBondForceImpl::calcEnergy(ContextImpl& context) {
+    return dynamic_cast<CalcCustomBondForceKernel&>(kernel.getImpl()).executeEnergy(context);
+}
+
+vector<string> CustomBondForceImpl::getKernelNames() {
+    vector<string> names;
+    names.push_back(CalcCustomBondForceKernel::Name());
+    return names;
+}
+
+map<string, double> CustomBondForceImpl::getDefaultParameters() {
+    map<string, double> parameters;
+    for (int i = 0; i < owner.getNumGlobalParameters(); i++)
+        parameters[owner.getGlobalParameterName(i)] = owner.getGlobalParameterDefaultValue(i);
+    return parameters;
+}
+

platforms/cuda/src/CudaKernels.cpp

@@ -682,6 +682,7 @@ void CudaIntegrateLangevinStepKernel::execute(ContextImpl& context, const Langev
     double stepSize = integrator.getStepSize();
     if (temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
         // Initialize the GPU parameters.
+        
         double tau = (friction == 0.0 ? 0.0 : 1.0/friction);
         gpuSetLangevinIntegrationParameters(gpu, (float) tau, (float) stepSize, (float) temperature, 0.0f);
         gpuSetConstants(gpu);

platforms/reference/src/ReferenceKernelFactory.cpp

@@ -48,6 +48,8 @@ KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Pla
         return new ReferenceCalcCustomNonbondedForceKernel(name, platform);
     if (name == CalcHarmonicBondForceKernel::Name())
         return new ReferenceCalcHarmonicBondForceKernel(name, platform);
+    if (name == CalcCustomBondForceKernel::Name())
+        return new ReferenceCalcCustomBondForceKernel(name, platform);
     if (name == CalcHarmonicAngleForceKernel::Name())
         return new ReferenceCalcHarmonicAngleForceKernel(name, platform);
     if (name == CalcHarmonicAngleForceKernel::Name())

platforms/reference/src/ReferenceKernels.cpp

@@ -37,6 +37,7 @@
 #include "SimTKReference/ReferenceBondForce.h"
 #include "SimTKReference/ReferenceBrownianDynamics.h"
 #include "SimTKReference/ReferenceCCMAAlgorithm.h"
+#include "SimTKReference/ReferenceCustomBondIxn.h"
 #include "SimTKReference/ReferenceCustomNonbondedIxn.h"
 #include "SimTKReference/ReferenceHarmonicBondIxn.h"
 #include "SimTKReference/ReferenceLJCoulomb14.h"
@@ -262,6 +263,69 @@ double ReferenceCalcHarmonicBondForceKernel::executeEnergy(ContextImpl& context)
     return energy;
 }
 
+ReferenceCalcCustomBondForceKernel::~ReferenceCalcCustomBondForceKernel() {
+    disposeIntArray(bondIndexArray, numBonds);
+    disposeRealArray(bondParamArray, numBonds);
+}
+
+void ReferenceCalcCustomBondForceKernel::initialize(const System& system, const CustomBondForce& force) {
+    numBonds = force.getNumBonds();
+    int numParameters = force.getNumPerBondParameters();
+
+    // Build the arrays.
+
+    bondIndexArray = allocateIntArray(numBonds, numParameters);
+    bondParamArray = allocateRealArray(numBonds, numParameters);
+    vector<double> params;
+    for (int i = 0; i < force.getNumBonds(); ++i) {
+        int particle1, particle2;
+        force.getBondParameters(i, particle1, particle2, params);
+        bondIndexArray[i][0] = particle1;
+        bondIndexArray[i][1] = particle2;
+        for (int j = 0; j < numParameters; j++)
+            bondParamArray[i][j] = (RealOpenMM) params[j];
+    }
+
+    // Parse the expression used to calculate the force.
+
+    Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction()).optimize();
+    energyExpression = expression.createProgram();
+    forceExpression = expression.differentiate("r").optimize().createProgram();
+    for (int i = 0; i < numParameters; i++)
+        parameterNames.push_back(force.getPerBondParameterName(i));
+    for (int i = 0; i < force.getNumGlobalParameters(); i++)
+        globalParameterNames.push_back(force.getGlobalParameterName(i));
+}
+
+void ReferenceCalcCustomBondForceKernel::executeForces(ContextImpl& context) {
+    RealOpenMM** posData = extractPositions(context);
+    RealOpenMM** forceData = extractForces(context);
+    map<string, double> globalParameters;
+    for (int i = 0; i < (int) globalParameterNames.size(); i++)
+        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
+    ReferenceBondForce refBondForce;
+    ReferenceCustomBondIxn harmonicBond(energyExpression, forceExpression, parameterNames, globalParameters);
+    refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, 0, 0, 0, harmonicBond);
+}
+
+double ReferenceCalcCustomBondForceKernel::executeEnergy(ContextImpl& context) {
+    RealOpenMM** posData = extractPositions(context);
+    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
+    RealOpenMM* energyArray = new RealOpenMM[numBonds];
+    RealOpenMM energy = 0;
+    map<string, double> globalParameters;
+    for (int i = 0; i < (int) globalParameterNames.size(); i++)
+        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
+    ReferenceBondForce refBondForce;
+    ReferenceCustomBondIxn harmonicBond(energyExpression, forceExpression, parameterNames, globalParameters);
+    for (int i = 0; i < numBonds; ++i)
+        energyArray[i] = 0;
+    refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, energyArray, 0, &energy, harmonicBond);
+    disposeRealArray(forceData, context.getSystem().getNumParticles());
+    delete[] energyArray;
+    return energy;
+}
+
 ReferenceCalcHarmonicAngleForceKernel::~ReferenceCalcHarmonicAngleForceKernel() {
     disposeIntArray(angleIndexArray, numAngles);
     disposeRealArray(angleParamArray, numAngles);

platforms/reference/src/ReferenceKernels.h

@@ -192,6 +192,42 @@ private:
     RealOpenMM **bondParamArray;
 };
 
+/**
+ * This kernel is invoked by CustomBondForce to calculate the forces acting on the system and the energy of the system.
+ */
+class ReferenceCalcCustomBondForceKernel : public CalcCustomBondForceKernel {
+public:
+    ReferenceCalcCustomBondForceKernel(std::string name, const Platform& platform) : CalcCustomBondForceKernel(name, platform) {
+    }
+    ~ReferenceCalcCustomBondForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CustomBondForce this kernel will be used for
+     */
+    void initialize(const System& system, const CustomBondForce& 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 CustomBondForce
+     */
+    double executeEnergy(ContextImpl& context);
+private:
+    int numBonds;
+    int **bondIndexArray;
+    RealOpenMM **bondParamArray;
+    Lepton::ExpressionProgram energyExpression, forceExpression;
+    std::vector<std::string> parameterNames, globalParameterNames;
+};
+
 /**
  * This kernel is invoked by HarmonicAngleForce to calculate the forces acting on the system and the energy of the system.
  */

platforms/reference/src/ReferencePlatform.cpp

@@ -42,6 +42,7 @@ ReferencePlatform::ReferencePlatform() {
     registerKernelFactory(CalcForcesAndEnergyKernel::Name(), factory);
     registerKernelFactory(UpdateStateDataKernel::Name(), factory);
     registerKernelFactory(CalcHarmonicBondForceKernel::Name(), factory);
+    registerKernelFactory(CalcCustomBondForceKernel::Name(), factory);
     registerKernelFactory(CalcHarmonicAngleForceKernel::Name(), factory);
     registerKernelFactory(CalcPeriodicTorsionForceKernel::Name(), factory);
     registerKernelFactory(CalcRBTorsionForceKernel::Name(), factory);

platforms/reference/src/SimTKReference/ReferenceCustomBondIxn.cpp

+
+/* Portions copyright (c) 2009 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 "../SimTKUtilities/SimTKOpenMMCommon.h"
+#include "../SimTKUtilities/SimTKOpenMMLog.h"
+#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
+#include "ReferenceCustomBondIxn.h"
+#include "ReferenceForce.h"
+
+using namespace std;
+
+/**---------------------------------------------------------------------------------------
+
+   ReferenceCustomBondIxn constructor
+
+   --------------------------------------------------------------------------------------- */
+
+ReferenceCustomBondIxn::ReferenceCustomBondIxn(const Lepton::ExpressionProgram& energyExpression,
+        const Lepton::ExpressionProgram& forceExpression, const vector<string>& parameterNames, map<string, double> globalParameters) :
+        energyExpression(energyExpression), forceExpression(forceExpression), paramNames(parameterNames), globalParameters(globalParameters) {
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nReferenceCustomBondIxn::ReferenceCustomBondIxn";
+
+   // ---------------------------------------------------------------------------------------
+
+}
+
+/**---------------------------------------------------------------------------------------
+
+   ReferenceCustomBondIxn destructor
+
+   --------------------------------------------------------------------------------------- */
+
+ReferenceCustomBondIxn::~ReferenceCustomBondIxn( ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nReferenceCustomBondIxn::~ReferenceCustomBondIxn";
+
+   // ---------------------------------------------------------------------------------------
+
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Calculate Custom Bond Ixn
+
+   @param atomIndices      atom indices of atom participating in bond
+   @param atomCoordinates  atom coordinates
+   @param parameters       parameters values
+   @param forces           force array (forces added to input values)
+   @param energiesByBond   energies by bond: energiesByBond[bondIndex]
+   @param energiesByAtom   energies by atom: energiesByAtom[atomIndex]
+
+   @return ReferenceForce::DefaultReturn;
+
+   --------------------------------------------------------------------------------------- */
+
+int ReferenceCustomBondIxn::calculateBondIxn( int* atomIndices,
+                                                RealOpenMM** atomCoordinates,
+                                                RealOpenMM* parameters,
+                                                RealOpenMM** forces,
+                                                RealOpenMM* energiesByBond,
+                                                RealOpenMM* energiesByAtom ) const {
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nReferenceCustomBondIxn::calculateBondIxn";
+
+   // ---------------------------------------------------------------------------------------
+
+   static const std::string methodName = "\nReferenceCustomBondIxn::calculateBondIxn";
+
+   static const int twoI               = 2;
+
+   static const RealOpenMM zero        = 0.0;
+   static const RealOpenMM two         = 2.0;
+   static const RealOpenMM half        = 0.5;
+
+   RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
+   map<string, double> variables = globalParameters;
+   for (int i = 0; i < (int) paramNames.size(); ++i)
+       variables[paramNames[i]] = parameters[i];
+
+   // ---------------------------------------------------------------------------------------
+
+   // get deltaR, R2, and R between 2 atoms
+
+   int atomAIndex = atomIndices[0];
+   int atomBIndex = atomIndices[1];
+   ReferenceForce::getDeltaR( atomCoordinates[atomAIndex], atomCoordinates[atomBIndex], deltaR );
+   variables["r"]             = deltaR[ReferenceForce::RIndex];
+   RealOpenMM dEdR            = forceExpression.evaluate(variables);
+   dEdR                       = deltaR[ReferenceForce::RIndex] > zero ? (dEdR/deltaR[ReferenceForce::RIndex]) : zero;
+
+   forces[atomAIndex][0]     += dEdR*deltaR[ReferenceForce::XIndex];
+   forces[atomAIndex][1]     += dEdR*deltaR[ReferenceForce::YIndex];
+   forces[atomAIndex][2]     += dEdR*deltaR[ReferenceForce::ZIndex];
+
+   forces[atomBIndex][0]     -= dEdR*deltaR[ReferenceForce::XIndex];
+   forces[atomBIndex][1]     -= dEdR*deltaR[ReferenceForce::YIndex];
+   forces[atomBIndex][2]     -= dEdR*deltaR[ReferenceForce::ZIndex];
+
+   RealOpenMM energy          = energyExpression.evaluate(variables);
+   updateEnergy( energy, energiesByBond, twoI, atomIndices, energiesByAtom );
+
+   return ReferenceForce::DefaultReturn;
+}

platforms/reference/src/SimTKReference/ReferenceCustomBondIxn.h

+
+/* Portions copyright (c) 2009 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 __ReferenceCustomBondIxn_H__
+#define __ReferenceCustomBondIxn_H__
+
+#include "ReferenceBondIxn.h"
+#include "lepton/ExpressionProgram.h"
+
+// ---------------------------------------------------------------------------------------
+
+class ReferenceCustomBondIxn : public ReferenceBondIxn {
+
+   private:
+      Lepton::ExpressionProgram energyExpression;
+      Lepton::ExpressionProgram forceExpression;
+      std::vector<std::string> paramNames;
+      std::map<std::string, double> globalParameters;
+
+   public:
+
+      /**---------------------------------------------------------------------------------------
+
+         Constructor
+
+         --------------------------------------------------------------------------------------- */
+
+       ReferenceCustomBondIxn(const Lepton::ExpressionProgram& energyExpression, const Lepton::ExpressionProgram& forceExpression,
+                              const std::vector<std::string>& parameterNames, std::map<std::string, double> globalParameters);
+
+      /**---------------------------------------------------------------------------------------
+
+         Destructor
+
+         --------------------------------------------------------------------------------------- */
+
+       ~ReferenceCustomBondIxn( );
+
+      /**---------------------------------------------------------------------------------------
+
+         Calculate Custom Bond Ixn
+
+         @param atomIndices      two bond indices
+         @param atomCoordinates  atom coordinates
+         @param parameters       parameter values
+         @param forces           force array (forces added)
+         @param energiesByBond   energies by bond: energiesByBond[bondIndex]
+         @param energiesByAtom   energies by atom: energiesByAtom[atomIndex]
+
+         --------------------------------------------------------------------------------------- */
+
+      int calculateBondIxn( int* atomIndices, RealOpenMM** atomCoordinates,
+                            RealOpenMM* parameters, RealOpenMM** forces,
+                            RealOpenMM* energiesByBond, RealOpenMM* energiesByAtom ) const;
+
+
+};
+
+// ---------------------------------------------------------------------------------------
+
+#endif // _ReferenceCustomBondIxn___

platforms/reference/tests/TestReferenceCustomBondForce.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-2009 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 all the different force terms in the reference implementation of HarmonicBondForce.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "ReferencePlatform.h"
+#include "openmm/CustomBondForce.h"
+#include "openmm/System.h"
+#include "openmm/VerletIntegrator.h"
+#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+const double TOL = 1e-5;
+
+void testBonds() {
+    ReferencePlatform platform;
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    CustomBondForce* forceField = new CustomBondForce("scale*k*(r-r0)^2");
+    forceField->addPerBondParameter("r0");
+    forceField->addPerBondParameter("k");
+    forceField->addGlobalParameter("scale", 0.5);
+    vector<double> parameters(2);
+    parameters[0] = 1.5;
+    parameters[1] = 0.8;
+    forceField->addBond(0, 1, parameters);
+    parameters[0] = 1.2;
+    parameters[1] = 0.7;
+    forceField->addBond(1, 2, parameters);
+    system.addForce(forceField);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(3);
+    positions[0] = Vec3(0, 2, 0);
+    positions[1] = Vec3(0, 0, 0);
+    positions[2] = Vec3(1, 0, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    ASSERT_EQUAL_VEC(Vec3(0, -0.8*0.5, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0.7*0.2, 0, 0), forces[2], TOL);
+    ASSERT_EQUAL_VEC(Vec3(-forces[0][0]-forces[2][0], -forces[0][1]-forces[2][1], -forces[0][2]-forces[2][2]), forces[1], TOL);
+    ASSERT_EQUAL_TOL(0.5*0.8*0.5*0.5 + 0.5*0.7*0.2*0.2, state.getPotentialEnergy(), TOL);
+}
+
+int main() {
+    try {
+        testBonds();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}
+