Created CustomExternalForce, along with the reference implementation

olla/include/openmm/kernels.h

@@ -36,6 +36,7 @@
 #include "openmm/BrownianIntegrator.h"
 #include "openmm/CMMotionRemover.h"
 #include "openmm/CustomBondForce.h"
+#include "openmm/CustomExternalForce.h"
 #include "openmm/CustomNonbondedForce.h"
 #include "openmm/GBSAOBCForce.h"
 #include "openmm/GBVIForce.h"
@@ -215,7 +216,7 @@ public:
      * Initialize the kernel.
      *
      * @param system     the System this kernel will be applied to
-     * @param force      the HarmonicBondForce this kernel will be used for
+     * @param force      the CustomBondForce this kernel will be used for
      */
     virtual void initialize(const System& system, const CustomBondForce& force) = 0;
     /**
@@ -228,7 +229,7 @@ public:
      * 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
+     * @return the potential energy due to the CustomBondForce
      */
     virtual double executeEnergy(ContextImpl& context) = 0;
 };
@@ -470,6 +471,38 @@ public:
     virtual double executeEnergy(ContextImpl& context) = 0;
 };
 
+/**
+ * This kernel is invoked by CustomExternalForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CalcCustomExternalForceKernel : public KernelImpl {
+public:
+    static std::string Name() {
+        return "CalcCustomExternalForce";
+    }
+    CalcCustomExternalForceKernel(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 CustomExternalForce this kernel will be used for
+     */
+    virtual void initialize(const System& system, const CustomExternalForce& 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 CustomExternalForce
+     */
+    virtual double executeEnergy(ContextImpl& context) = 0;
+};
+
 /**
  * This kernel is invoked by VerletIntegrator to take one time step.
  */

openmmapi/include/openmm/CustomBondForce.h

@@ -104,7 +104,7 @@ public:
      */
     void setEnergyFunction(const std::string& energy);
     /**
-     * Add a new per-bond parmeter that the interaction may depend on.
+     * Add a new per-bond parameter that the interaction may depend on.
      *
      * @param name             the name of the parameter
      * @return the index of the parameter that was added
@@ -125,7 +125,7 @@ public:
      */
     void setPerBondParameterName(int index, const std::string& name);
     /**
-     * Add a new global parmeter that the interaction may depend on.
+     * Add a new global parameter that the interaction may depend on.
      *
      * @param name             the name of the parameter
      * @param defaultValue     the default value of the parameter

openmmapi/include/openmm/CustomExternalForce.h

+#ifndef OPENMM_CUSTOMEXTERNALFORCE_H_
+#define OPENMM_CUSTOMEXTERNALFORCE_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 an "external" force on particles.  The force may be applied to any subset of the particles
+ * in the System.  The force on each particle is specified by an arbitrary algebraic expression, which may depend
+ * on the current position of the particle as well as on arbitrary global and per-particle parameters.
+ *
+ * To use this class, create a CustomExternalForce object, passing an algebraic expression to the constructor
+ * that defines the potential energy of each affected particle.  The expression may depend on the particle's x, y, and
+ * z coordinates, as well as on any parameters you choose.  Then call addPerParticleParameter() to define per-particle
+ * parameters, and addGlobalParameter() to define global parameters.  The values of per-particle parameters are specified as
+ * part of the system definition, while values of global parameters may be modified during a simulation by calling Context::setParameter().
+ * Finally, call addParticle() once for each particle that should be affected by the force.  After a particle has been added,
+ * you can modify its parameters by calling setParticleParameters().
+ *
+ * As an example, the following code creates a CustomExternalForce that attracts each particle to a target position (x0, y0, z0)
+ * via a harmonic potential:
+ *
+ * <tt>CustomExternalForce* force = new CustomExternalForce("k*((x-x0)^2+(y-y0)^2+(z-z0)^2");</tt>
+ *
+ * This force depends on four parameters: the spring constant k and equilibrium coordinates x0, y0, and z0.  The following code defines these parameters:
+ *
+ * <tt><pre>
+ * force->addGlobalParameter("k");
+ * force->addPerParticleParameter("x0");
+ * force->addPerParticleParameter("y0");
+ * force->addPerParticleParameter("z0");
+ * </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 CustomExternalForce : public Force {
+public:
+    /**
+     * Create a CustomExternalForce.
+     *
+     * @param energy    an algebraic expression giving the potential energy of each particle as a function
+     *                  of its x, y, and z coordinates
+     */
+    CustomExternalForce(const std::string& energy);
+    /**
+     * Get the number of particles for which force field parameters have been defined.
+     */
+    int getNumParticles() const {
+        return particles.size();
+    }
+    /**
+     * Get the number of per-particle parameters that the force depends on
+     */
+    int getNumPerParticleParameters() const {
+        return parameters.size();
+    }
+    /**
+     * Get the number of global parameters that the force depends on.
+     */
+    int getNumGlobalParameters() const {
+        return globalParameters.size();
+    }
+    /**
+     * Get the algebraic expression that gives the potential energy of each particle
+     */
+    const std::string& getEnergyFunction() const;
+    /**
+     * Set the algebraic expression that gives the potential energy of each particle
+     */
+    void setEnergyFunction(const std::string& energy);
+    /**
+     * Add a new per-particle parameter that the force may depend on.
+     *
+     * @param name             the name of the parameter
+     * @return the index of the parameter that was added
+     */
+    int addPerParticleParameter(const std::string& name);
+    /**
+     * Get the name of a per-particle parameter.
+     *
+     * @param index     the index of the parameter for which to get the name
+     * @return the parameter name
+     */
+    const std::string& getPerParticleParameterName(int index) const;
+    /**
+     * Set the name of a per-particle parameter.
+     *
+     * @param index          the index of the parameter for which to set the name
+     * @param name           the name of the parameter
+     */
+    void setPerParticleParameterName(int index, const std::string& name);
+    /**
+     * Add a new global parameter that the force 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 particle term to the force field.
+     *
+     * @param particle     the index of the particle this term is applied to
+     * @param parameters   the list of parameters for the new force term
+     * @return the index of the particle term that was added
+     */
+    int addParticle(int particle, const std::vector<double>& parameters);
+    /**
+     * Get the force field parameters for a force field term.
+     *
+     * @param index         the index of the particle term for which to get parameters
+     * @param particle      the index of the particle this term is applied to
+     * @param parameters    the list of parameters for the force field term
+     */
+    void getParticleParameters(int index, int& particle, std::vector<double>& parameters) const;
+    /**
+     * Set the force field parameters for a force field term.
+     *
+     * @param index         the index of the particle term for which to set parameters
+     * @param particle      the index of the particle this term is applied to
+     * @param parameters    the list of parameters for the force field term
+     */
+    void setParticleParameters(int index, int particle, const std::vector<double>& parameters);
+protected:
+    ForceImpl* createImpl();
+private:
+    class ParticleInfo;
+    class ParticleParameterInfo;
+    class GlobalParameterInfo;
+    std::string energyExpression;
+    std::vector<ParticleParameterInfo> parameters;
+    std::vector<GlobalParameterInfo> globalParameters;
+    std::vector<ParticleInfo> particles;
+};
+
+class CustomExternalForce::ParticleInfo {
+public:
+    int particle;
+    std::vector<double> parameters;
+    ParticleInfo() : particle(-1) {
+    }
+    ParticleInfo(int particle, const std::vector<double>& parameters) : particle(particle), parameters(parameters) {
+    }
+};
+
+class CustomExternalForce::ParticleParameterInfo {
+public:
+    std::string name;
+    ParticleParameterInfo() {
+    }
+    ParticleParameterInfo(const std::string& name) : name(name) {
+    }
+};
+
+class CustomExternalForce::GlobalParameterInfo {
+public:
+    std::string name;
+    double defaultValue;
+    GlobalParameterInfo() {
+    }
+    GlobalParameterInfo(const std::string& name, double defaultValue) : name(name), defaultValue(defaultValue) {
+    }
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_CUSTOMEXTERNALFORCE_H_*/

openmmapi/include/openmm/CustomNonbondedForce.h

@@ -171,7 +171,7 @@ public:
      */
     void setCutoffDistance(double distance);
     /**
-     * Add a new per-particle parmeter that the interaction may depend on.
+     * Add a new per-particle parameter that the interaction may depend on.
      *
      * @param name             the name of the parameter
      * @param combiningRule    an algebraic expression giving the combining rule for this parameter
@@ -207,7 +207,7 @@ public:
      */
     void setParameterCombiningRule(int index, const std::string& combiningRule);
     /**
-     * Add a new global parmeter that the interaction may depend on.
+     * Add a new global parameter that the interaction may depend on.
      *
      * @param name             the name of the parameter
      * @param defaultValue     the default value of the parameter

openmmapi/include/openmm/internal/CustomExternalForceImpl.h

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

openmmapi/src/CustomExternalForce.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/CustomExternalForce.h"
+#include "openmm/internal/CustomExternalForceImpl.h"
+#include <cmath>
+#include <map>
+#include <sstream>
+#include <utility>
+
+using namespace OpenMM;
+using std::string;
+using std::stringstream;
+using std::vector;
+
+CustomExternalForce::CustomExternalForce(const string& energy) : energyExpression(energy) {
+}
+
+const string& CustomExternalForce::getEnergyFunction() const {
+    return energyExpression;
+}
+
+void CustomExternalForce::setEnergyFunction(const std::string& energy) {
+    energyExpression = energy;
+}
+
+int CustomExternalForce::addPerParticleParameter(const string& name) {
+    parameters.push_back(ParticleParameterInfo(name));
+    return parameters.size()-1;
+}
+
+const string& CustomExternalForce::getPerParticleParameterName(int index) const {
+    return parameters[index].name;
+}
+
+void CustomExternalForce::setPerParticleParameterName(int index, const string& name) {
+    parameters[index].name = name;
+}
+
+int CustomExternalForce::addGlobalParameter(const string& name, double defaultValue) {
+    globalParameters.push_back(GlobalParameterInfo(name, defaultValue));
+    return globalParameters.size()-1;
+}
+
+const string& CustomExternalForce::getGlobalParameterName(int index) const {
+    return globalParameters[index].name;
+}
+
+void CustomExternalForce::setGlobalParameterName(int index, const string& name) {
+    globalParameters[index].name = name;
+}
+
+double CustomExternalForce::getGlobalParameterDefaultValue(int index) const {
+    return globalParameters[index].defaultValue;
+}
+
+void CustomExternalForce::setGlobalParameterDefaultValue(int index, double defaultValue) {
+    globalParameters[index].defaultValue = defaultValue;
+}
+
+int CustomExternalForce::addParticle(int particle, const vector<double>& parameters) {
+    particles.push_back(ParticleInfo(particle, parameters));
+    return particles.size()-1;
+}
+
+void CustomExternalForce::getParticleParameters(int index, int& particle, std::vector<double>& parameters) const {
+    particle = particles[index].particle;
+    parameters = particles[index].parameters;
+}
+
+void CustomExternalForce::setParticleParameters(int index, int particle, const vector<double>& parameters) {
+    particles[index].parameters = parameters;
+    particles[index].particle = particle;
+}
+
+ForceImpl* CustomExternalForce::createImpl() {
+    return new CustomExternalForceImpl(*this);
+}

openmmapi/src/CustomExternalForceImpl.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/CustomExternalForceImpl.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;
+
+CustomExternalForceImpl::CustomExternalForceImpl(CustomExternalForce& owner) : owner(owner) {
+}
+
+CustomExternalForceImpl::~CustomExternalForceImpl() {
+}
+
+void CustomExternalForceImpl::initialize(ContextImpl& context) {
+    kernel = context.getPlatform().createKernel(CalcCustomExternalForceKernel::Name(), context);
+
+    // Check for errors in the specification of bonds.
+
+    System& system = context.getSystem();
+    vector<double> parameters;
+    int numParameters = owner.getNumPerParticleParameters();
+    for (int i = 0; i < owner.getNumParticles(); i++) {
+        int particle;
+        owner.getParticleParameters(i, particle, parameters);
+        if (particle < 0 || particle >= system.getNumParticles()) {
+            stringstream msg;
+            msg << "CustomExternalForce: Illegal particle index: ";
+            msg << particle;
+            throw OpenMMException(msg.str());
+        }
+        if (parameters.size() != numParameters) {
+            stringstream msg;
+            msg << "CustomExternalForce: Wrong number of parameters for particle ";
+            msg << i;
+            throw OpenMMException(msg.str());
+        }
+    }
+    dynamic_cast<CalcCustomExternalForceKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner);
+}
+
+void CustomExternalForceImpl::calcForces(ContextImpl& context) {
+    dynamic_cast<CalcCustomExternalForceKernel&>(kernel.getImpl()).executeForces(context);
+}
+
+double CustomExternalForceImpl::calcEnergy(ContextImpl& context) {
+    return dynamic_cast<CalcCustomExternalForceKernel&>(kernel.getImpl()).executeEnergy(context);
+}
+
+vector<string> CustomExternalForceImpl::getKernelNames() {
+    vector<string> names;
+    names.push_back(CalcCustomExternalForceKernel::Name());
+    return names;
+}
+
+map<string, double> CustomExternalForceImpl::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

@@ -62,6 +62,8 @@ KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Pla
         return new ReferenceCalcGBSAOBCForceKernel(name, platform);
     if (name == CalcGBVIForceKernel::Name())
         return new ReferenceCalcGBVIForceKernel(name, platform);
+    if (name == CalcCustomExternalForceKernel::Name())
+        return new ReferenceCalcCustomExternalForceKernel(name, platform);
     if (name == IntegrateVerletStepKernel::Name())
         return new ReferenceIntegrateVerletStepKernel(name, platform, data);
     if (name == IntegrateLangevinStepKernel::Name())

platforms/reference/src/ReferenceKernels.cpp

@@ -38,6 +38,7 @@
 #include "SimTKReference/ReferenceBrownianDynamics.h"
 #include "SimTKReference/ReferenceCCMAAlgorithm.h"
 #include "SimTKReference/ReferenceCustomBondIxn.h"
+#include "SimTKReference/ReferenceCustomExternalIxn.h"
 #include "SimTKReference/ReferenceCustomNonbondedIxn.h"
 #include "SimTKReference/ReferenceHarmonicBondIxn.h"
 #include "SimTKReference/ReferenceLJCoulomb14.h"
@@ -912,6 +913,63 @@ double ReferenceCalcGBVIForceKernel::executeEnergy(ContextImpl& context) {
     return static_cast<double>(energy);
 }
 
+ReferenceCalcCustomExternalForceKernel::~ReferenceCalcCustomExternalForceKernel() {
+    disposeRealArray(particleParamArray, numParticles);
+}
+
+void ReferenceCalcCustomExternalForceKernel::initialize(const System& system, const CustomExternalForce& force) {
+    numParticles = force.getNumParticles();
+    int numParameters = force.getNumPerParticleParameters();
+
+    // Build the arrays.
+
+    particles.resize(numParticles);
+    particleParamArray = allocateRealArray(numParticles, numParameters);
+    vector<double> params;
+    for (int i = 0; i < numParticles; ++i) {
+        force.getParticleParameters(i, particles[i], params);
+        for (int j = 0; j < numParameters; j++)
+            particleParamArray[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();
+    forceExpressionX = expression.differentiate("x").optimize().createProgram();
+    forceExpressionY = expression.differentiate("y").optimize().createProgram();
+    forceExpressionZ = expression.differentiate("z").optimize().createProgram();
+    for (int i = 0; i < numParameters; i++)
+        parameterNames.push_back(force.getPerParticleParameterName(i));
+    for (int i = 0; i < force.getNumGlobalParameters(); i++)
+        globalParameterNames.push_back(force.getGlobalParameterName(i));
+}
+
+void ReferenceCalcCustomExternalForceKernel::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]);
+    ReferenceCustomExternalIxn force(energyExpression, forceExpressionX, forceExpressionY, forceExpressionZ, parameterNames, globalParameters);
+    for (int i = 0; i < numParticles; ++i)
+        force.calculateForce(particles[i], posData, particleParamArray[i], forceData, 0);
+}
+
+double ReferenceCalcCustomExternalForceKernel::executeEnergy(ContextImpl& context) {
+    RealOpenMM** posData = extractPositions(context);
+    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
+    RealOpenMM energy = 0;
+    map<string, double> globalParameters;
+    for (int i = 0; i < (int) globalParameterNames.size(); i++)
+        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
+    ReferenceCustomExternalIxn force(energyExpression, forceExpressionX, forceExpressionY, forceExpressionZ, parameterNames, globalParameters);
+    for (int i = 0; i < numParticles; ++i)
+        force.calculateForce(particles[i], posData, particleParamArray[i], forceData, &energy);
+    disposeRealArray(forceData, context.getSystem().getNumParticles());
+    return energy;
+}
+
 ReferenceIntegrateVerletStepKernel::~ReferenceIntegrateVerletStepKernel() {
     if (dynamics)
         delete dynamics;

platforms/reference/src/ReferenceKernels.h

@@ -478,6 +478,42 @@ private:
     std::vector<RealOpenMM> charges;
 };
 
+/**
+ * This kernel is invoked by CustomExternalForce to calculate the forces acting on the system and the energy of the system.
+ */
+class ReferenceCalcCustomExternalForceKernel : public CalcCustomExternalForceKernel {
+public:
+    ReferenceCalcCustomExternalForceKernel(std::string name, const Platform& platform) : CalcCustomExternalForceKernel(name, platform) {
+    }
+    ~ReferenceCalcCustomExternalForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CustomExternalForce this kernel will be used for
+     */
+    void initialize(const System& system, const CustomExternalForce& 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 CustomExternalForce
+     */
+    double executeEnergy(ContextImpl& context);
+private:
+    int numParticles;
+    std::vector<int> particles;
+    RealOpenMM **particleParamArray;
+    Lepton::ExpressionProgram energyExpression, forceExpressionX, forceExpressionY, forceExpressionZ;
+    std::vector<std::string> parameterNames, globalParameterNames;
+};
+
 /**
  * This kernel is invoked by VerletIntegrator to take one time step.
  */

platforms/reference/src/ReferencePlatform.cpp

@@ -50,6 +50,7 @@ ReferencePlatform::ReferencePlatform() {
     registerKernelFactory(CalcCustomNonbondedForceKernel::Name(), factory);
     registerKernelFactory(CalcGBSAOBCForceKernel::Name(), factory);
     registerKernelFactory(CalcGBVIForceKernel::Name(), factory);
+    registerKernelFactory(CalcCustomExternalForceKernel::Name(), factory);
     registerKernelFactory(IntegrateVerletStepKernel::Name(), factory);
     registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);
     registerKernelFactory(IntegrateBrownianStepKernel::Name(), factory);

platforms/reference/src/SimTKReference/ReferenceCustomExternalIxn.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 "ReferenceCustomExternalIxn.h"
+#include "ReferenceForce.h"
+
+using namespace std;
+
+/**---------------------------------------------------------------------------------------
+
+   ReferenceCustomExternalIxn constructor
+
+   --------------------------------------------------------------------------------------- */
+
+ReferenceCustomExternalIxn::ReferenceCustomExternalIxn(const Lepton::ExpressionProgram& energyExpression,
+        const Lepton::ExpressionProgram& forceExpressionX, const Lepton::ExpressionProgram& forceExpressionY,
+        const Lepton::ExpressionProgram& forceExpressionZ, const vector<string>& parameterNames, map<string, double> globalParameters) :
+        energyExpression(energyExpression), forceExpressionX(forceExpressionX), forceExpressionY(forceExpressionY),
+        forceExpressionZ(forceExpressionZ), paramNames(parameterNames), globalParameters(globalParameters) {
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nReferenceCustomExternalIxn::ReferenceCustomExternalIxn";
+
+   // ---------------------------------------------------------------------------------------
+
+}
+
+/**---------------------------------------------------------------------------------------
+
+   ReferenceCustomExternalIxn destructor
+
+   --------------------------------------------------------------------------------------- */
+
+ReferenceCustomExternalIxn::~ReferenceCustomExternalIxn( ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nReferenceCustomExternalIxn::~ReferenceCustomExternalIxn";
+
+   // ---------------------------------------------------------------------------------------
+
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Calculate Custom External Ixn
+
+   @param atomIndex        the index of the atom to apply the force to
+   @param atomCoordinates  atom coordinates
+   @param parameters       parameters values
+   @param forces           force array (forces added to input values)
+   @param energy           energy is added to this
+
+   @return ReferenceForce::DefaultReturn;
+
+   --------------------------------------------------------------------------------------- */
+
+int ReferenceCustomExternalIxn::calculateForce( int atomIndex,
+                                                RealOpenMM** atomCoordinates,
+                                                RealOpenMM* parameters,
+                                                RealOpenMM** forces,
+                                                RealOpenMM* energy ) const {
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nReferenceCustomExternalIxn::calculateBondIxn";
+
+   // ---------------------------------------------------------------------------------------
+
+   static const std::string methodName = "\nReferenceCustomExternalIxn::calculateBondIxn";
+
+   map<string, double> variables = globalParameters;
+   for (int i = 0; i < (int) paramNames.size(); ++i)
+       variables[paramNames[i]] = parameters[i];
+   variables["x"] = atomCoordinates[atomIndex][0];
+   variables["y"] = atomCoordinates[atomIndex][1];
+   variables["z"] = atomCoordinates[atomIndex][2];
+
+   // ---------------------------------------------------------------------------------------
+
+   forces[atomIndex][0] += forceExpressionX.evaluate(variables);
+   forces[atomIndex][1] += forceExpressionY.evaluate(variables);
+   forces[atomIndex][2] += forceExpressionZ.evaluate(variables);
+   if (energy != NULL)
+       *energy += energyExpression.evaluate(variables);
+
+   return ReferenceForce::DefaultReturn;
+}

platforms/reference/src/SimTKReference/ReferenceCustomExternalIxn.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 __ReferenceCustomExternalIxn_H__
+#define __ReferenceCustomExternalIxn_H__
+
+#include "ReferenceCustomExternalIxn.h"
+#include "lepton/ExpressionProgram.h"
+
+// ---------------------------------------------------------------------------------------
+
+class ReferenceCustomExternalIxn {
+
+   private:
+      Lepton::ExpressionProgram energyExpression;
+      Lepton::ExpressionProgram forceExpressionX;
+      Lepton::ExpressionProgram forceExpressionY;
+      Lepton::ExpressionProgram forceExpressionZ;
+      std::vector<std::string> paramNames;
+      std::map<std::string, double> globalParameters;
+
+   public:
+
+      /**---------------------------------------------------------------------------------------
+
+         Constructor
+
+         --------------------------------------------------------------------------------------- */
+
+       ReferenceCustomExternalIxn(const Lepton::ExpressionProgram& energyExpression, const Lepton::ExpressionProgram& forceExpressionX,
+                              const Lepton::ExpressionProgram& forceExpressionY, const Lepton::ExpressionProgram& forceExpressionZ,
+                              const std::vector<std::string>& parameterNames, std::map<std::string, double> globalParameters);
+
+      /**---------------------------------------------------------------------------------------
+
+         Destructor
+
+         --------------------------------------------------------------------------------------- */
+
+       ~ReferenceCustomExternalIxn( );
+
+      /**---------------------------------------------------------------------------------------
+
+         Calculate Custom External Force
+
+         @param atomIndex        the index of the atom to apply the force to
+         @param atomCoordinates  atom coordinates
+         @param parameters       parameter values
+         @param forces           force array (forces added)
+         @param energy           energy is added to this
+
+         --------------------------------------------------------------------------------------- */
+
+      int calculateForce( int atomIndex, RealOpenMM** atomCoordinates,
+                            RealOpenMM* parameters, RealOpenMM** forces, RealOpenMM* energy ) const;
+
+
+};
+
+// ---------------------------------------------------------------------------------------
+
+#endif // _ReferenceCustomBondIxn___

platforms/reference/tests/TestReferenceCustomExternalForce.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 the reference implementation of CustomExternalForce.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "ReferencePlatform.h"
+#include "openmm/CustomExternalForce.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 testForce() {
+    ReferencePlatform platform;
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    CustomExternalForce* forceField = new CustomExternalForce("scale*(x+yscale*(y-y0)^2)");
+    forceField->addPerParticleParameter("y0");
+    forceField->addPerParticleParameter("yscale");
+    forceField->addGlobalParameter("scale", 0.5);
+    vector<double> parameters(2);
+    parameters[0] = 0.5;
+    parameters[1] = 2.0;
+    forceField->addParticle(0, parameters);
+    parameters[0] = 1.5;
+    parameters[1] = 3.0;
+    forceField->addParticle(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, 1);
+    positions[2] = Vec3(1, 0, 1);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    ASSERT_EQUAL_VEC(Vec3(0.5, 0.5*2.0*2.0*1.5, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[1], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0.5, -0.5*3.0*2.0*1.5, 0), forces[2], TOL);
+    ASSERT_EQUAL_TOL(0.5*(1.0 + 2.0*1.5*1.5 + 3.0*1.5*1.5), state.getPotentialEnergy(), TOL);
+}
+
+int main() {
+    try {
+        testForce();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}
+
+