Preliminary implementation of CustomGBForce

olla/include/openmm/kernels.h

@@ -37,6 +37,7 @@
 #include "openmm/CMMotionRemover.h"
 #include "openmm/CustomBondForce.h"
 #include "openmm/CustomExternalForce.h"
+#include "openmm/CustomGBForce.h"
 #include "openmm/CustomNonbondedForce.h"
 #include "openmm/GBSAOBCForce.h"
 #include "openmm/GBVIForce.h"
@@ -471,6 +472,43 @@ public:
     virtual double executeEnergy(ContextImpl& context) = 0;
 };
 
+/**
+ * This kernel is invoked by CustomGBForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CalcCustomGBForceKernel : public KernelImpl {
+public:
+    enum NonbondedMethod {
+        NoCutoff = 0,
+        CutoffNonPeriodic = 1,
+        CutoffPeriodic = 2
+    };
+    static std::string Name() {
+        return "CalcCustomGBForce";
+    }
+    CalcCustomGBForceKernel(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 CustomGBForce this kernel will be used for
+     */
+    virtual void initialize(const System& system, const CustomGBForce& 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 CustomGBForce
+     */
+    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.
  */

openmmapi/include/OpenMM.h

@@ -37,6 +37,7 @@
 #include "openmm/CMMotionRemover.h"
 #include "openmm/CustomBondForce.h"
 #include "openmm/CustomExternalForce.h"
+#include "openmm/CustomGBForce.h"
 #include "openmm/CustomNonbondedForce.h"
 #include "openmm/Force.h"
 #include "openmm/GBSAOBCForce.h"

openmmapi/include/openmm/CustomGBForce.h

+#ifndef OPENMM_CUSTOMGBFORCE_H_
+#define OPENMM_CUSTOMGBFORCE_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 <map>
+#include <set>
+#include <utility>
+#include <vector>
+#include "internal/windowsExport.h"
+
+namespace OpenMM {
+
+/**
+ */
+
+class OPENMM_EXPORT CustomGBForce : public Force {
+public:
+    /**
+     * This is an enumeration of the different methods that may be used for handling long range nonbonded forces.
+     */
+    enum NonbondedMethod {
+        /**
+         * No cutoff is applied to nonbonded interactions.  The full set of N^2 interactions is computed exactly.
+         * This necessarily means that periodic boundary conditions cannot be used.  This is the default.
+         */
+        NoCutoff = 0,
+        /**
+         * Interactions beyond the cutoff distance are ignored.
+         */
+        CutoffNonPeriodic = 1,
+        /**
+         * Periodic boundary conditions are used, so that each particle interacts only with the nearest periodic copy of
+         * each other particle.  Interactions beyond the cutoff distance are ignored.
+         */
+        CutoffPeriodic = 2,
+    };
+    enum ComputationType {
+        SingleParticle = 0,
+        ParticlePair = 1,
+        ParticlePairNoExclusions = 2
+    };
+    /**
+     * Create a CustomGBForce.
+     */
+    CustomGBForce();
+    /**
+     * Get the number of particles for which force field parameters have been defined.
+     */
+    int getNumParticles() const {
+        return particles.size();
+    }
+    /**
+     * Get the number of particle pairs whose interactions should be excluded.
+     */
+    int getNumExclusions() const {
+        return exclusions.size();
+    }
+    /**
+     * Get the number of per-particle parameters that the interaction depends on.
+     */
+    int getNumPerParticleParameters() const {
+        return parameters.size();
+    }
+    /**
+     * Get the number of global parameters that the interaction depends on.
+     */
+    int getNumGlobalParameters() const {
+        return globalParameters.size();
+    }
+    /**
+     * Get the number of tabulated functions that have been defined.
+     */
+    int getNumFunctions() const {
+        return functions.size();
+    }
+    int getNumComputedValues() const {
+        return computedValues.size();
+    }
+    int getNumEnergyTerms() const {
+        return energyTerms.size();
+    }
+    /**
+     * Get the method used for handling long range nonbonded interactions.
+     */
+    NonbondedMethod getNonbondedMethod() const;
+    /**
+     * Set the method used for handling long range nonbonded interactions.
+     */
+    void setNonbondedMethod(NonbondedMethod method);
+    /**
+     * Get the cutoff distance (in nm) being used for nonbonded interactions.  If the NonbondedMethod in use
+     * is NoCutoff, this value will have no effect.
+     */
+    double getCutoffDistance() const;
+    /**
+     * Set the cutoff distance (in nm) being used for nonbonded interactions.  If the NonbondedMethod in use
+     * is NoCutoff, this value will have no effect.
+     */
+    void setCutoffDistance(double distance);
+    /**
+     * Add a new per-particle parameter that the interaction may depend on.
+     *
+     * @param name     the name of the parameter
+     * @return the index of the parameter that was added
+     */
+    int addPerParticleParameter(const std::string& name);
+    /**
+     * Get the name of a per-particle parameter.
+     *
+     * @param index     the index of the parameter for which to get the name
+     * @return the parameter name
+     */
+    const std::string& getPerParticleParameterName(int index) const;
+    /**
+     * Set the name of a per-particle parameter.
+     *
+     * @param index          the index of the parameter for which to set the name
+     * @param name           the name of the parameter
+     */
+    void setPerParticleParameterName(int index, const std::string& name);
+    /**
+     * Add a new global parameter that the interaction may depend on.
+     *
+     * @param name             the name of the parameter
+     * @param defaultValue     the default value of the parameter
+     * @return the index of the parameter that was added
+     */
+    int addGlobalParameter(const std::string& name, double defaultValue);
+    /**
+     * Get the name of a global parameter.
+     *
+     * @param index     the index of the parameter for which to get the name
+     * @return the parameter name
+     */
+    const std::string& getGlobalParameterName(int index) const;
+    /**
+     * Set the name of a global parameter.
+     *
+     * @param index          the index of the parameter for which to set the name
+     * @param name           the name of the parameter
+     */
+    void setGlobalParameterName(int index, const std::string& name);
+    /**
+     * Get the default value of a global parameter.
+     *
+     * @param index     the index of the parameter for which to get the default value
+     * @return the parameter default value
+     */
+    double getGlobalParameterDefaultValue(int index) const;
+    /**
+     * Set the default value of a global parameter.
+     *
+     * @param index          the index of the parameter for which to set the default value
+     * @param name           the default value of the parameter
+     */
+    void setGlobalParameterDefaultValue(int index, double defaultValue);
+    /**
+     * Add the nonbonded force parameters for a particle.  This should be called once for each particle
+     * in the System.  When it is called for the i'th time, it specifies the parameters for the i'th particle.
+     *
+     * @param parameters    the list of parameters for the new particle
+     * @return the index of the particle that was added
+     */
+    int addParticle(const std::vector<double>& parameters);
+    /**
+     * Get the nonbonded force parameters for a particle.
+     *
+     * @param index       the index of the particle for which to get parameters
+     * @param parameters  the list of parameters for the specified particle
+     */
+    void getParticleParameters(int index, std::vector<double>& parameters) const;
+    /**
+     * Set the nonbonded force parameters for a particle.
+     *
+     * @param index       the index of the particle for which to set parameters
+     * @param parameters  the list of parameters for the specified particle
+     */
+    void setParticleParameters(int index, const std::vector<double>& parameters);
+    int addComputedValue(const std::string& name, const std::string& expression, ComputationType type);
+    void getComputedValueParameters(int index, std::string& name, std::string& expression, ComputationType& type) const;
+    void setComputedValueParameters(int index, const std::string& name, const std::string& expression, ComputationType type);
+    int addEnergyTerm(const std::string& expression, ComputationType type);
+    void getEnergyTermParameters(int index, std::string& expression, ComputationType& type) const;
+    void setEnergyTermParameters(int index, const std::string& expression, ComputationType type);
+    /**
+     * Add a particle pair to the list of interactions that should be excluded.
+     *
+     * @param particle1  the index of the first particle in the pair
+     * @param particle2  the index of the second particle in the pair
+     * @return the index of the exclusion that was added
+     */
+    int addExclusion(int particle1, int particle2);
+    /**
+     * Get the particles in a pair whose interaction should be excluded.
+     *
+     * @param index      the index of the exclusion for which to get particle indices
+     * @param particle1  the index of the first particle in the pair
+     * @param particle2  the index of the second particle in the pair
+     */
+    void getExclusionParticles(int index, int& particle1, int& particle2) const;
+    /**
+     * Set the particles in a pair whose interaction should be excluded.
+     *
+     * @param index      the index of the exclusion for which to set particle indices
+     * @param particle1  the index of the first particle in the pair
+     * @param particle2  the index of the second particle in the pair
+     */
+    void setExclusionParticles(int index, int particle1, int particle2);
+    /**
+     * Add a tabulated function that may appear in the energy expression.
+     *
+     * @param name           the name of the function as it appears in expressions
+     * @param values         the tabulated values of the function f(x) at uniformly spaced values of x between min and max.
+     *                       The function is assumed to be zero for x &lt; min or x &gt; max.
+     * @param min            the value of the independent variable corresponding to the first element of values
+     * @param max            the value of the independent variable corresponding to the last element of values
+     * @param interpolating  if true, an interpolating (Catmull-Rom) spline will be used to represent the function.
+     *                       If false, an approximating spline (B-spline) will be used.
+     * @return the index of the function that was added
+     */
+    int addFunction(const std::string& name, const std::vector<double>& values, double min, double max, bool interpolating);
+    /**
+     * Get the parameters for a tabulated function that may appear in the energy expression.
+     *
+     * @param index          the index of the function for which to get parameters
+     * @param name           the name of the function as it appears in expressions
+     * @param values         the tabulated values of the function f(x) at uniformly spaced values of x between min and max.
+     *                       The function is assumed to be zero for x &lt; min or x &gt; max.
+     * @param min            the value of the independent variable corresponding to the first element of values
+     * @param max            the value of the independent variable corresponding to the last element of values
+     * @param interpolating  if true, an interpolating (Catmull-Rom) spline will be used to represent the function.
+     *                       If false, an approximating spline (B-spline) will be used.
+     */
+    void getFunctionParameters(int index, std::string& name, std::vector<double>& values, double& min, double& max, bool& interpolating) const;
+    /**
+     * Set the parameters for a tabulated function that may appear in algebraic expressions.
+     *
+     * @param index          the index of the function for which to set parameters
+     * @param name           the name of the function as it appears in expressions
+     * @param values         the tabulated values of the function f(x) at uniformly spaced values of x between min and max.
+     *                       The function is assumed to be zero for x &lt; min or x &gt; max.
+     * @param min            the value of the independent variable corresponding to the first element of values
+     * @param max            the value of the independent variable corresponding to the last element of values
+     * @param interpolating  if true, an interpolating (Catmull-Rom) spline will be used to represent the function.
+     *                       If false, an approximating spline (B-spline) will be used.
+     */
+    void setFunctionParameters(int index, const std::string& name, const std::vector<double>& values, double min, double max, bool interpolating);
+protected:
+    ForceImpl* createImpl();
+private:
+    class ParticleInfo;
+    class PerParticleParameterInfo;
+    class GlobalParameterInfo;
+    class ExclusionInfo;
+    class FunctionInfo;
+    class ComputationInfo;
+    NonbondedMethod nonbondedMethod;
+    double cutoffDistance;
+    std::string energyExpression;
+    std::vector<PerParticleParameterInfo> parameters;
+    std::vector<GlobalParameterInfo> globalParameters;
+    std::vector<ParticleInfo> particles;
+    std::vector<ExclusionInfo> exclusions;
+    std::vector<FunctionInfo> functions;
+    std::vector<ComputationInfo> computedValues;
+    std::vector<ComputationInfo> energyTerms;
+};
+
+class CustomGBForce::ParticleInfo {
+public:
+    std::vector<double> parameters;
+    ParticleInfo() {
+    }
+    ParticleInfo(const std::vector<double>& parameters) : parameters(parameters) {
+    }
+};
+
+class CustomGBForce::PerParticleParameterInfo {
+public:
+    std::string name;
+    PerParticleParameterInfo() {
+    }
+    PerParticleParameterInfo(const std::string& name) : name(name) {
+    }
+};
+
+class CustomGBForce::GlobalParameterInfo {
+public:
+    std::string name;
+    double defaultValue;
+    GlobalParameterInfo() {
+    }
+    GlobalParameterInfo(const std::string& name, double defaultValue) : name(name), defaultValue(defaultValue) {
+    }
+};
+
+class CustomGBForce::ExclusionInfo {
+public:
+    int particle1, particle2;
+    ExclusionInfo() {
+        particle1 = particle2 = -1;
+    }
+    ExclusionInfo(int particle1, int particle2) :
+        particle1(particle1), particle2(particle2) {
+    }
+};
+
+class CustomGBForce::FunctionInfo {
+public:
+    std::string name;
+    std::vector<double> values;
+    double min, max;
+    bool interpolating;
+    FunctionInfo() {
+    }
+    FunctionInfo(const std::string& name, const std::vector<double>& values, double min, double max, bool interpolating) :
+        name(name), values(values), min(min), max(max), interpolating(interpolating) {
+    }
+};
+
+class CustomGBForce::ComputationInfo {
+public:
+    std::string name;
+    std::string expression;
+    CustomGBForce::ComputationType type;
+    ComputationInfo() {
+    }
+    ComputationInfo(const std::string& name, const std::string& expression, CustomGBForce::ComputationType type) :
+        name(name), expression(expression), type(type) {
+    }
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_CUSTOMGBFORCE_H_*/

openmmapi/include/openmm/CustomNonbondedForce.h

@@ -326,7 +326,6 @@ private:
     std::vector<ParticleInfo> particles;
     std::vector<ExclusionInfo> exclusions;
     std::vector<FunctionInfo> functions;
-    std::map<std::pair<int, int>, int> exclusionMap;
 };
 
 class CustomNonbondedForce::ParticleInfo {

openmmapi/include/openmm/internal/CustomGBForceImpl.h

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

openmmapi/src/CustomGBForce.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/CustomGBForce.h"
+#include "openmm/internal/CustomGBForceImpl.h"
+#include <cmath>
+#include <map>
+#include <sstream>
+#include <utility>
+
+using namespace OpenMM;
+using std::map;
+using std::pair;
+using std::set;
+using std::string;
+using std::stringstream;
+using std::vector;
+
+CustomGBForce::CustomGBForce() : nonbondedMethod(NoCutoff), cutoffDistance(1.0) {
+}
+
+CustomGBForce::NonbondedMethod CustomGBForce::getNonbondedMethod() const {
+    return nonbondedMethod;
+}
+
+void CustomGBForce::setNonbondedMethod(NonbondedMethod method) {
+    nonbondedMethod = method;
+}
+
+double CustomGBForce::getCutoffDistance() const {
+    return cutoffDistance;
+}
+
+void CustomGBForce::setCutoffDistance(double distance) {
+    cutoffDistance = distance;
+}
+
+int CustomGBForce::addPerParticleParameter(const string& name) {
+    parameters.push_back(PerParticleParameterInfo(name));
+    return parameters.size()-1;
+}
+
+const string& CustomGBForce::getPerParticleParameterName(int index) const {
+    return parameters[index].name;
+}
+
+void CustomGBForce::setPerParticleParameterName(int index, const string& name) {
+    parameters[index].name = name;
+}
+
+int CustomGBForce::addGlobalParameter(const string& name, double defaultValue) {
+    globalParameters.push_back(GlobalParameterInfo(name, defaultValue));
+    return globalParameters.size()-1;
+}
+
+const string& CustomGBForce::getGlobalParameterName(int index) const {
+    return globalParameters[index].name;
+}
+
+void CustomGBForce::setGlobalParameterName(int index, const string& name) {
+    globalParameters[index].name = name;
+}
+
+double CustomGBForce::getGlobalParameterDefaultValue(int index) const {
+    return globalParameters[index].defaultValue;
+}
+
+void CustomGBForce::setGlobalParameterDefaultValue(int index, double defaultValue) {
+    globalParameters[index].defaultValue = defaultValue;
+}
+
+int CustomGBForce::addParticle(const vector<double>& parameters) {
+    particles.push_back(ParticleInfo(parameters));
+    return particles.size()-1;
+}
+
+void CustomGBForce::getParticleParameters(int index, std::vector<double>& parameters) const {
+    parameters = particles[index].parameters;
+}
+
+void CustomGBForce::setParticleParameters(int index, const vector<double>& parameters) {
+    particles[index].parameters = parameters;
+}
+
+int CustomGBForce::addComputedValue(const std::string& name, const std::string& expression, ComputationType type) {
+    computedValues.push_back(CustomGBForce::ComputationInfo(name, expression, type));
+    return computedValues.size()-1;
+}
+
+void CustomGBForce::getComputedValueParameters(int index, std::string& name, std::string& expression, ComputationType& type) const {
+    name = computedValues[index].name;
+    expression = computedValues[index].expression;
+    type = computedValues[index].type;
+}
+
+void CustomGBForce::setComputedValueParameters(int index, const std::string& name, const std::string& expression, ComputationType type) {
+    computedValues[index].name = name;
+    computedValues[index].expression = expression;
+    computedValues[index].type = type;
+}
+
+int CustomGBForce::addEnergyTerm(const std::string& expression, ComputationType type) {
+    energyTerms.push_back(CustomGBForce::ComputationInfo("", expression, type));
+    return energyTerms.size()-1;
+}
+
+void CustomGBForce::getEnergyTermParameters(int index, std::string& expression, ComputationType& type) const {
+    expression = energyTerms[index].expression;
+    type = energyTerms[index].type;
+}
+
+void CustomGBForce::setEnergyTermParameters(int index, const std::string& expression, ComputationType type) {
+    energyTerms[index].expression = expression;
+    energyTerms[index].type = type;
+}
+
+int CustomGBForce::addExclusion(int particle1, int particle2) {
+    exclusions.push_back(ExclusionInfo(particle1, particle2));
+    return exclusions.size()-1;
+}
+void CustomGBForce::getExclusionParticles(int index, int& particle1, int& particle2) const {
+    particle1 = exclusions[index].particle1;
+    particle2 = exclusions[index].particle2;
+}
+
+void CustomGBForce::setExclusionParticles(int index, int particle1, int particle2) {
+    exclusions[index].particle1 = particle1;
+    exclusions[index].particle2 = particle2;
+}
+
+int CustomGBForce::addFunction(const std::string& name, const std::vector<double>& values, double min, double max, bool interpolating) {
+    if (max <= min)
+        throw OpenMMException("CustomGBForce: max <= min for a tabulated function.");
+    if (values.size() < 2)
+        throw OpenMMException("CustomGBForce: a tabulated function must have at least two points");
+    functions.push_back(FunctionInfo(name, values, min, max, interpolating));
+    return functions.size()-1;
+}
+
+void CustomGBForce::getFunctionParameters(int index, std::string& name, std::vector<double>& values, double& min, double& max, bool& interpolating) const {
+    name = functions[index].name;
+    values = functions[index].values;
+    min = functions[index].min;
+    max = functions[index].max;
+    interpolating = functions[index].interpolating;
+}
+
+void CustomGBForce::setFunctionParameters(int index, const std::string& name, const std::vector<double>& values, double min, double max, bool interpolating) {
+    if (max <= min)
+        throw OpenMMException("CustomGBForce: max <= min for a tabulated function.");
+    if (values.size() < 2)
+        throw OpenMMException("CustomGBForce: a tabulated function must have at least two points");
+    functions[index].name = name;
+    functions[index].values = values;
+    functions[index].min = min;
+    functions[index].max = max;
+    functions[index].interpolating = interpolating;
+}
+
+ForceImpl* CustomGBForce::createImpl() {
+    return new CustomGBForceImpl(*this);
+}
+

openmmapi/src/CustomGBForceImpl.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/CustomGBForceImpl.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;
+
+CustomGBForceImpl::CustomGBForceImpl(CustomGBForce& owner) : owner(owner) {
+}
+
+CustomGBForceImpl::~CustomGBForceImpl() {
+}
+
+void CustomGBForceImpl::initialize(ContextImpl& context) {
+    kernel = context.getPlatform().createKernel(CalcCustomGBForceKernel::Name(), context);
+
+    // Check for errors in the specification of parameters and exclusions.
+
+    System& system = context.getSystem();
+    if (owner.getNumParticles() != system.getNumParticles())
+        throw OpenMMException("CustomGBForce must have exactly as many particles as the System it belongs to.");
+    vector<set<int> > exclusions(owner.getNumParticles());
+    vector<double> parameters;
+    int numParameters = owner.getNumPerParticleParameters();
+    for (int i = 0; i < owner.getNumParticles(); i++) {
+        owner.getParticleParameters(i, parameters);
+        if (parameters.size() != numParameters) {
+            stringstream msg;
+            msg << "CustomGBForce: Wrong number of parameters for particle ";
+            msg << i;
+            throw OpenMMException(msg.str());
+        }
+    }
+    for (int i = 0; i < owner.getNumExclusions(); i++) {
+        int particle1, particle2;
+        owner.getExclusionParticles(i, particle1, particle2);
+        if (particle1 < 0 || particle1 >= owner.getNumParticles()) {
+            stringstream msg;
+            msg << "CustomGBForce: Illegal particle index for an exclusion: ";
+            msg << particle1;
+            throw OpenMMException(msg.str());
+        }
+        if (particle2 < 0 || particle2 >= owner.getNumParticles()) {
+            stringstream msg;
+            msg << "CustomGBForce: Illegal particle index for an exclusion: ";
+            msg << particle2;
+            throw OpenMMException(msg.str());
+        }
+        if (exclusions[particle1].count(particle2) > 0 || exclusions[particle2].count(particle1) > 0) {
+            stringstream msg;
+            msg << "CustomGBForce: Multiple exclusions are specified for particles ";
+            msg << particle1;
+            msg << " and ";
+            msg << particle2;
+            throw OpenMMException(msg.str());
+        }
+        exclusions[particle1].insert(particle2);
+        exclusions[particle2].insert(particle1);
+    }
+    dynamic_cast<CalcCustomGBForceKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner);
+}
+
+void CustomGBForceImpl::calcForces(ContextImpl& context) {
+    dynamic_cast<CalcCustomGBForceKernel&>(kernel.getImpl()).executeForces(context);
+}
+
+double CustomGBForceImpl::calcEnergy(ContextImpl& context) {
+    return dynamic_cast<CalcCustomGBForceKernel&>(kernel.getImpl()).executeEnergy(context);
+}
+
+vector<string> CustomGBForceImpl::getKernelNames() {
+    vector<string> names;
+    names.push_back(CalcCustomGBForceKernel::Name());
+    return names;
+}
+
+map<string, double> CustomGBForceImpl::getDefaultParameters() {
+    map<string, double> parameters;
+    for (int i = 0; i < owner.getNumGlobalParameters(); i++)
+        parameters[owner.getGlobalParameterName(i)] = owner.getGlobalParameterDefaultValue(i);
+    return parameters;
+}
+

platforms/opencl/src/OpenCLKernels.cpp

@@ -1109,7 +1109,7 @@ void OpenCLCalcGBSAOBCForceKernel::initialize(const System& system, const GBSAOB
     }
     posq.upload();
     params->upload(paramsVector);
-    prefactor = 2.0*-166.02691*0.4184*((1.0/force.getSoluteDielectric())-(1.0/force.getSolventDielectric()));
+    prefactor = -138.935485*((1.0/force.getSoluteDielectric())-(1.0/force.getSolventDielectric()));
     bool useCutoff = (force.getNonbondedMethod() != GBSAOBCForce::NoCutoff);
     bool usePeriodic = (force.getNonbondedMethod() != GBSAOBCForce::NoCutoff && force.getNonbondedMethod() != GBSAOBCForce::CutoffNonPeriodic);
     string source = cl.loadSourceFromFile("gbsaObc2.cl");

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 == CalcCustomGBForceKernel::Name())
+        return new ReferenceCalcCustomGBForceKernel(name, platform);
     if (name == CalcCustomExternalForceKernel::Name())
         return new ReferenceCalcCustomExternalForceKernel(name, platform);
     if (name == IntegrateVerletStepKernel::Name())

platforms/reference/src/ReferenceKernels.cpp

@@ -39,6 +39,7 @@
 #include "SimTKReference/ReferenceCCMAAlgorithm.h"
 #include "SimTKReference/ReferenceCustomBondIxn.h"
 #include "SimTKReference/ReferenceCustomExternalIxn.h"
+#include "SimTKReference/ReferenceCustomGBIxn.h"
 #include "SimTKReference/ReferenceCustomNonbondedIxn.h"
 #include "SimTKReference/ReferenceHarmonicBondIxn.h"
 #include "SimTKReference/ReferenceLJCoulomb14.h"
@@ -891,6 +892,151 @@ double ReferenceCalcGBVIForceKernel::executeEnergy(ContextImpl& context) {
     return static_cast<double>(energy);
 }
 
+ReferenceCalcCustomGBForceKernel::~ReferenceCalcCustomGBForceKernel() {
+    disposeRealArray(particleParamArray, numParticles);
+    if (neighborList != NULL)
+        delete neighborList;
+}
+
+void ReferenceCalcCustomGBForceKernel::initialize(const System& system, const CustomGBForce& force) {
+
+    // Record the exclusions.
+
+    numParticles = force.getNumParticles();
+    exclusions.resize(numParticles);
+    for (int i = 0; i < force.getNumExclusions(); i++) {
+        int particle1, particle2;
+        force.getExclusionParticles(i, particle1, particle2);
+        exclusions[particle1].insert(particle2);
+        exclusions[particle2].insert(particle1);
+    }
+
+    // Build the arrays.
+
+    int numPerParticleParameters = force.getNumPerParticleParameters();
+    particleParamArray = allocateRealArray(numParticles, numPerParticleParameters);
+    for (int i = 0; i < numParticles; ++i) {
+        vector<double> parameters;
+        force.getParticleParameters(i, parameters);
+        for (int j = 0; j < numPerParticleParameters; j++)
+            particleParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
+    }
+    for (int i = 0; i < numPerParticleParameters; i++)
+        particleParameterNames.push_back(force.getPerParticleParameterName(i));
+    for (int i = 0; i < force.getNumGlobalParameters(); i++)
+        globalParameterNames.push_back(force.getGlobalParameterName(i));
+    nonbondedMethod = CalcCustomGBForceKernel::NonbondedMethod(force.getNonbondedMethod());
+    nonbondedCutoff = (RealOpenMM) force.getCutoffDistance();
+    Vec3 boxVectors[3];
+    system.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
+    periodicBoxSize[0] = (RealOpenMM) boxVectors[0][0];
+    periodicBoxSize[1] = (RealOpenMM) boxVectors[1][1];
+    periodicBoxSize[2] = (RealOpenMM) boxVectors[2][2];
+    if (nonbondedMethod == NoCutoff)
+        neighborList = NULL;
+    else
+        neighborList = new NeighborList();
+
+    // Create custom functions for the tabulated functions.
+
+    map<string, Lepton::CustomFunction*> functions;
+//    for (int i = 0; i < force.getNumFunctions(); i++) {
+//        string name;
+//        vector<double> values;
+//        double min, max;
+//        bool interpolating;
+//        force.getFunctionParameters(i, name, values, min, max, interpolating);
+//        functions[name] = new TabulatedFunction(min, max, values, interpolating);
+//    }
+
+    // Parse the expressions for computed values.
+
+    valueDerivExpressions.resize(force.getNumComputedValues());
+    for (int i = 0; i < force.getNumComputedValues(); i++) {
+        string name, expression;
+        CustomGBForce::ComputationType type;
+        force.getComputedValueParameters(i, name, expression, type);
+        Lepton::ParsedExpression ex = Lepton::Parser::parse(expression, functions).optimize();
+        valueExpressions.push_back(ex.createProgram());
+        valueTypes.push_back(type);
+        valueNames.push_back(name);
+        if (type != CustomGBForce::SingleParticle)
+            valueDerivExpressions[i].push_back(ex.differentiate("r").optimize().createProgram());
+        for (int j = 0; j < i; j++) {
+            if (type == CustomGBForce::SingleParticle)
+                valueDerivExpressions[i].push_back(ex.differentiate(valueNames[j]).optimize().createProgram());
+            else
+                valueDerivExpressions[i].push_back(ex.differentiate(valueNames[j]+"1").optimize().createProgram());
+        }
+    }
+
+    // Parse the various expressions used to calculate the force.
+
+    energyDerivExpressions.resize(force.getNumEnergyTerms());
+    for (int i = 0; i < force.getNumEnergyTerms(); i++) {
+        string expression;
+        CustomGBForce::ComputationType type;
+        force.getEnergyTermParameters(i, expression, type);
+        Lepton::ParsedExpression ex = Lepton::Parser::parse(expression, functions).optimize();
+        energyExpressions.push_back(ex.createProgram());
+        energyTypes.push_back(type);
+        if (type != CustomGBForce::SingleParticle)
+            energyDerivExpressions[i].push_back(ex.differentiate("r").optimize().createProgram());
+        for (int j = 0; j < force.getNumComputedValues(); j++) {
+            if (type == CustomGBForce::SingleParticle)
+                energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]).optimize().createProgram());
+            else {
+                energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]+"1").optimize().createProgram());
+                energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]+"2").optimize().createProgram());
+            }
+        }
+    }
+
+    // Delete the custom functions.
+
+    for (map<string, Lepton::CustomFunction*>::iterator iter = functions.begin(); iter != functions.end(); iter++)
+        delete iter->second;
+}
+
+void ReferenceCalcCustomGBForceKernel::executeForces(ContextImpl& context) {
+    RealOpenMM** posData = extractPositions(context);
+    RealOpenMM** forceData = extractForces(context);
+    ReferenceCustomGBIxn ixn(valueExpressions, valueDerivExpressions, valueNames, valueTypes, energyExpressions,
+        energyDerivExpressions, energyTypes, particleParameterNames);
+    bool periodic = (nonbondedMethod == CutoffPeriodic);
+    if (nonbondedMethod != NoCutoff) {
+        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, periodic ? periodicBoxSize : NULL, nonbondedCutoff, 0.0);
+        ixn.setUseCutoff(nonbondedCutoff, *neighborList);
+    }
+    if (periodic)
+        ixn.setPeriodic(periodicBoxSize);
+    map<string, double> globalParameters;
+    for (int i = 0; i < (int) globalParameterNames.size(); i++)
+        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
+    ixn.calculateIxn(numParticles, posData, particleParamArray, exclusions, globalParameters, forceData, 0, 0);
+}
+
+double ReferenceCalcCustomGBForceKernel::executeEnergy(ContextImpl& context) {
+    RealOpenMM** posData = extractPositions(context);
+    RealOpenMM** forceData = allocateRealArray(numParticles, 3);
+    RealOpenMM energy = 0;
+    ReferenceCustomGBIxn ixn(valueExpressions, valueDerivExpressions, valueNames, valueTypes, energyExpressions,
+        energyDerivExpressions, energyTypes, particleParameterNames);
+    bool periodic = (nonbondedMethod == CutoffPeriodic);
+    if (nonbondedMethod != NoCutoff) {
+        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, periodic ? periodicBoxSize : NULL, nonbondedCutoff, 0.0);
+        ixn.setUseCutoff(nonbondedCutoff, *neighborList);
+    }
+    if (periodic)
+        ixn.setPeriodic(periodicBoxSize);
+    map<string, double> globalParameters;
+    for (int i = 0; i < (int) globalParameterNames.size(); i++)
+        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
+    ixn.calculateIxn(numParticles, posData, particleParamArray, exclusions, globalParameters, forceData, 0, &energy);
+    disposeRealArray(forceData, numParticles);
+    return energy;
+}
+
 ReferenceCalcCustomExternalForceKernel::~ReferenceCalcCustomExternalForceKernel() {
     disposeRealArray(particleParamArray, numParticles);
 }

platforms/reference/src/ReferenceKernels.h

@@ -477,6 +477,50 @@ private:
     std::vector<RealOpenMM> charges;
 };
 
+/**
+ * This kernel is invoked by CustomGBForce to calculate the forces acting on the system.
+ */
+class ReferenceCalcCustomGBForceKernel : public CalcCustomGBForceKernel {
+public:
+    ReferenceCalcCustomGBForceKernel(std::string name, const Platform& platform) : CalcCustomGBForceKernel(name, platform) {
+    }
+    ~ReferenceCalcCustomGBForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CustomGBForce this kernel will be used for
+     */
+    void initialize(const System& system, const CustomGBForce& 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 CustomGBForce
+     */
+    double executeEnergy(ContextImpl& context);
+private:
+    int numParticles;
+    RealOpenMM **particleParamArray;
+    RealOpenMM nonbondedCutoff, periodicBoxSize[3];
+    std::vector<std::set<int> > exclusions;
+    std::vector<std::string> particleParameterNames, globalParameterNames, valueNames;
+    std::vector<Lepton::ExpressionProgram> valueExpressions;
+    std::vector<std::vector<Lepton::ExpressionProgram> > valueDerivExpressions;
+    std::vector<OpenMM::CustomGBForce::ComputationType> valueTypes;
+    std::vector<Lepton::ExpressionProgram> energyExpressions;
+    std::vector<std::vector<Lepton::ExpressionProgram> > energyDerivExpressions;
+    std::vector<OpenMM::CustomGBForce::ComputationType> energyTypes;
+    NonbondedMethod nonbondedMethod;
+    NeighborList* neighborList;
+};
+
 /**
  * This kernel is invoked by CustomExternalForce to calculate the forces acting on the system and the energy of the system.
  */

platforms/reference/src/ReferencePlatform.cpp

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

platforms/reference/src/SimTKReference/ReferenceCustomGBIxn.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 "ReferenceForce.h"
+#include "ReferenceCustomGBIxn.h"
+
+using std::map;
+using std::set;
+using std::string;
+using std::stringstream;
+using std::vector;
+
+/**---------------------------------------------------------------------------------------
+
+   ReferenceCustomGBIxn constructor
+
+   --------------------------------------------------------------------------------------- */
+
+ReferenceCustomGBIxn::ReferenceCustomGBIxn(const vector<Lepton::ExpressionProgram>& valueExpressions,
+                     const vector<vector<Lepton::ExpressionProgram> >& valueDerivExpressions,
+                     const vector<string>& valueNames,
+                     const vector<OpenMM::CustomGBForce::ComputationType>& valueTypes,
+                     const vector<Lepton::ExpressionProgram>& energyExpressions,
+                     const vector<vector<Lepton::ExpressionProgram> > energyDerivExpressions,
+                     const vector<OpenMM::CustomGBForce::ComputationType>& energyTypes,
+                     const vector<string>& parameterNames) :
+            cutoff(false), periodic(false), valueExpressions(valueExpressions), valueDerivExpressions(valueDerivExpressions), valueNames(valueNames),
+            valueTypes(valueTypes), energyExpressions(energyExpressions), energyDerivExpressions(energyDerivExpressions),
+            energyTypes(energyTypes), paramNames(parameterNames) {
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nReferenceCustomGBIxn::ReferenceCustomGBIxn";
+
+   // ---------------------------------------------------------------------------------------
+
+    for (int i = 0; i < (int) paramNames.size(); i++) {
+        for (int j = 1; j < 3; j++) {
+            stringstream name;
+            name << paramNames[i] << j;
+            particleParamNames.push_back(name.str());
+        }
+    }
+    for (int i = 0; i < (int) valueNames.size(); i++) {
+        for (int j = 1; j < 3; j++) {
+            stringstream name;
+            name << valueNames[i] << j;
+            particleValueNames.push_back(name.str());
+        }
+    }
+}
+
+/**---------------------------------------------------------------------------------------
+
+   ReferenceCustomGBIxn destructor
+
+   --------------------------------------------------------------------------------------- */
+
+ReferenceCustomGBIxn::~ReferenceCustomGBIxn( ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nReferenceCustomGBIxn::~ReferenceCustomGBIxn";
+
+   // ---------------------------------------------------------------------------------------
+
+}
+
+  /**---------------------------------------------------------------------------------------
+
+     Set the force to use a cutoff.
+
+     @param distance            the cutoff distance
+     @param neighbors           the neighbor list to use
+
+     @return ReferenceForce::DefaultReturn
+
+     --------------------------------------------------------------------------------------- */
+
+  int ReferenceCustomGBIxn::setUseCutoff( RealOpenMM distance, const OpenMM::NeighborList& neighbors ) {
+
+    cutoff = true;
+    cutoffDistance = distance;
+    neighborList = &neighbors;
+
+    return ReferenceForce::DefaultReturn;
+  }
+
+  /**---------------------------------------------------------------------------------------
+
+     Set the force to use periodic boundary conditions.  This requires that a cutoff has
+     also been set, and the smallest side of the periodic box is at least twice the cutoff
+     distance.
+
+     @param boxSize             the X, Y, and Z widths of the periodic box
+
+     @return ReferenceForce::DefaultReturn
+
+     --------------------------------------------------------------------------------------- */
+
+  int ReferenceCustomGBIxn::setPeriodic( RealOpenMM* boxSize ) {
+
+    assert(cutoff);
+    assert(boxSize[0] >= 2.0*cutoffDistance);
+    assert(boxSize[1] >= 2.0*cutoffDistance);
+    assert(boxSize[2] >= 2.0*cutoffDistance);
+    periodic = true;
+    periodicBoxSize[0] = boxSize[0];
+    periodicBoxSize[1] = boxSize[1];
+    periodicBoxSize[2] = boxSize[2];
+    return ReferenceForce::DefaultReturn;
+
+  }
+
+
+/**---------------------------------------------------------------------------------------
+
+   Calculate the custom pair ixn
+
+   @param numberOfAtoms    number of atoms
+   @param atomCoordinates  atom coordinates
+   @param atomParameters   atom parameters                             atomParameters[atomIndex][paramterIndex]
+   @param exclusions       atom exclusion indices                      exclusions[atomIndex][atomToExcludeIndex]
+                           exclusions[atomIndex][0] = number of exclusions
+                           exclusions[atomIndex][1-no.] = atom indices of atoms to excluded from
+                           interacting w/ atom atomIndex
+   @param fixedParameters  non atom parameters (not currently used)
+   @param globalParameters the values of global parameters
+   @param forces           force array (forces added)
+   @param energyByAtom     atom energy
+   @param totalEnergy      total energy
+
+   @return ReferenceForce::DefaultReturn
+
+   --------------------------------------------------------------------------------------- */
+
+int ReferenceCustomGBIxn::calculateIxn(int numberOfAtoms, RealOpenMM** atomCoordinates, RealOpenMM** atomParameters,
+                                           const vector<set<int> >& exclusions, map<string, double>& globalParameters, RealOpenMM** forces,
+                                           RealOpenMM* energyByAtom, RealOpenMM* totalEnergy) const {
+    int numValues = valueTypes.size();
+    vector<vector<ComputedValue> > values(numValues);
+    for (int valueIndex = 0; valueIndex < numValues; valueIndex++) {
+        if (valueTypes[valueIndex] == OpenMM::CustomGBForce::SingleParticle)
+            calculateSingleParticleValue(valueIndex, numberOfAtoms, values, globalParameters, atomParameters);
+        else if (valueTypes[valueIndex] == OpenMM::CustomGBForce::ParticlePair)
+            calculateParticlePairValue(valueIndex, numberOfAtoms, atomCoordinates, atomParameters, values, globalParameters, exclusions, true);
+        else
+            calculateParticlePairValue(valueIndex, numberOfAtoms, atomCoordinates, atomParameters, values, globalParameters, exclusions, false);
+    }
+    for (int termIndex = 0; termIndex < (int) energyExpressions.size(); termIndex++) {
+        if (energyTypes[termIndex] == OpenMM::CustomGBForce::SingleParticle)
+            calculateSingleParticleEnergyTerm(termIndex, numberOfAtoms, values, globalParameters, atomParameters, forces, totalEnergy);
+        else if (energyTypes[termIndex] == OpenMM::CustomGBForce::ParticlePair)
+            calculateParticlePairEnergyTerm(termIndex, numberOfAtoms, atomCoordinates, atomParameters, values, globalParameters, exclusions, true, forces, totalEnergy);
+        else
+            calculateParticlePairEnergyTerm(termIndex, numberOfAtoms, atomCoordinates, atomParameters, values, globalParameters, exclusions, false, forces, totalEnergy);
+    }
+   return ReferenceForce::DefaultReturn;
+}
+
+void ReferenceCustomGBIxn::calculateSingleParticleValue(int index, int numAtoms, vector<vector<ComputedValue> >& values,
+        const map<string, double>& globalParameters, RealOpenMM** atomParameters) const {
+    values[index].resize(numAtoms);
+    map<string, double> variables = globalParameters;
+    for (int i = 0; i < numAtoms; i++) {
+        for (int j = 0; j < (int) paramNames.size(); j++)
+            variables[paramNames[j]] = atomParameters[i][j];
+        for (int j = 0; j < index; j++)
+            variables[valueNames[j]] = values[j][i].value;
+        values[index][i].value = (RealOpenMM) valueExpressions[index].evaluate(variables);
+        for (int j = 0; j < index; j++) {
+            RealOpenMM scale = (RealOpenMM) valueDerivExpressions[index][j].evaluate(variables);
+            values[index][i].gradient[0] += scale*values[j][i].gradient[0];
+            values[index][i].gradient[1] += scale*values[j][i].gradient[1];
+            values[index][i].gradient[2] += scale*values[j][i].gradient[2];
+        }
+    }
+}
+
+void ReferenceCustomGBIxn::calculateParticlePairValue(int index, int numAtoms, RealOpenMM** atomCoordinates, RealOpenMM** atomParameters,
+        vector<vector<ComputedValue> >& values, const map<string, double>& globalParameters, const vector<set<int> >& exclusions, bool useExclusions) const {
+    values[index].resize(numAtoms);
+    for (int i = 0; i < numAtoms; i++)
+        values[index][i].value = (RealOpenMM) 0.0;
+    if (cutoff) {
+        for (int i = 0; i < (int) neighborList->size(); i++) {
+            OpenMM::AtomPair pair = (*neighborList)[i];
+            if (useExclusions && exclusions[pair.first].find(pair.second) != exclusions[pair.first].end())
+                continue;
+            calculateOnePairValue(index, pair.first, pair.second, atomCoordinates, atomParameters, globalParameters, values);
+            calculateOnePairValue(index, pair.second, pair.first, atomCoordinates, atomParameters, globalParameters, values);
+        }
+    }
+    else {
+        for (int i = 0; i < numAtoms; i++){
+            for (int j = i+1; j < numAtoms; j++ ){
+                if (useExclusions && exclusions[i].find(j) != exclusions[i].end())
+                    continue;
+                calculateOnePairValue(index, i, j, atomCoordinates, atomParameters, globalParameters, values);
+                calculateOnePairValue(index, j, i, atomCoordinates, atomParameters, globalParameters, values);
+           }
+        }
+    }
+}
+
+void ReferenceCustomGBIxn::calculateOnePairValue(int index, int atom1, int atom2, RealOpenMM** atomCoordinates, RealOpenMM** atomParameters,
+        const map<string, double>& globalParameters, vector<vector<ComputedValue> >& values) const {
+    RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
+    if (periodic)
+        ReferenceForce::getDeltaRPeriodic(atomCoordinates[atom2], atomCoordinates[atom1], periodicBoxSize, deltaR);
+    else
+        ReferenceForce::getDeltaR(atomCoordinates[atom2], atomCoordinates[atom1], deltaR);
+    RealOpenMM r = deltaR[ReferenceForce::RIndex];
+    if (cutoff && r >= cutoffDistance)
+        return;
+    map<string, double> variables = globalParameters;
+    for (int i = 0; i < paramNames.size(); i++) {
+        variables[particleParamNames[i*2]] = atomParameters[atom1][i];
+        variables[particleParamNames[i*2+1]] = atomParameters[atom2][i];
+    }
+    variables["r"] = r;
+    for (int i = 0; i < index; i++) {
+        variables[particleValueNames[i*2]] = values[i][atom1].value;
+        variables[particleValueNames[i*2+1]] = values[i][atom2].value;
+    }
+    values[index][atom1].value += (RealOpenMM) valueExpressions[index].evaluate(variables);
+    RealOpenMM dVdR = (RealOpenMM) valueDerivExpressions[index][0].evaluate(variables);
+    RealOpenMM rinv = 1/r;
+    for (int i = 0; i < 3; i++)
+       values[index][atom1].gradient[i] += dVdR*deltaR[i]*rinv;
+    for (int i = 0; i < index; i++) {
+        RealOpenMM scale = (RealOpenMM) valueDerivExpressions[index][i+1].evaluate(variables);
+        values[index][atom1].gradient[0] += scale*values[i][atom1].gradient[0];
+        values[index][atom1].gradient[1] += scale*values[i][atom1].gradient[1];
+        values[index][atom1].gradient[2] += scale*values[i][atom1].gradient[2];
+    }
+}
+
+void ReferenceCustomGBIxn::calculateSingleParticleEnergyTerm(int index, int numAtoms, const vector<vector<ComputedValue> >& values,
+        const map<string, double>& globalParameters, RealOpenMM** atomParameters, RealOpenMM** forces, RealOpenMM* totalEnergy) const {
+    map<string, double> variables = globalParameters;
+    for (int i = 0; i < numAtoms; i++) {
+        for (int j = 0; j < (int) paramNames.size(); j++)
+            variables[paramNames[j]] = atomParameters[i][j];
+        for (int j = 0; j < (int) valueNames.size(); j++)
+            variables[valueNames[j]] = values[j][i].value;
+        if (totalEnergy != NULL)
+            *totalEnergy += (RealOpenMM) energyExpressions[index].evaluate(variables);
+        for (int j = 0; j < (int) valueNames.size(); j++) {
+            RealOpenMM scale = (RealOpenMM) energyDerivExpressions[index][j].evaluate(variables);
+            forces[i][0] -= scale*values[j][i].gradient[0];
+            forces[i][1] -= scale*values[j][i].gradient[1];
+            forces[i][2] -= scale*values[j][i].gradient[2];
+        }
+    }
+}
+
+void ReferenceCustomGBIxn::calculateParticlePairEnergyTerm(int index, int numAtoms, RealOpenMM** atomCoordinates, RealOpenMM** atomParameters,
+        const vector<vector<ComputedValue> >& values, const map<string, double>& globalParameters, const vector<set<int> >& exclusions, bool useExclusions,
+        RealOpenMM** forces, RealOpenMM* totalEnergy) const {
+    if (cutoff) {
+        for (int i = 0; i < (int) neighborList->size(); i++) {
+            OpenMM::AtomPair pair = (*neighborList)[i];
+            if (useExclusions && exclusions[pair.first].find(pair.second) != exclusions[pair.first].end())
+                continue;
+            calculateOnePairEnergyTerm(index, pair.first, pair.second, atomCoordinates, atomParameters, globalParameters, values, forces, totalEnergy);
+        }
+    }
+    else {
+        for (int i = 0; i < numAtoms; i++){
+            for (int j = i+1; j < numAtoms; j++ ){
+                if (useExclusions && exclusions[i].find(j) != exclusions[i].end())
+                    continue;
+                calculateOnePairEnergyTerm(index, i, j, atomCoordinates, atomParameters, globalParameters, values, forces, totalEnergy);
+           }
+        }
+    }
+}
+
+void ReferenceCustomGBIxn::calculateOnePairEnergyTerm(int index, int atom1, int atom2, RealOpenMM** atomCoordinates, RealOpenMM** atomParameters,
+        const map<string, double>& globalParameters, const vector<vector<ComputedValue> >& values, RealOpenMM** forces, RealOpenMM* totalEnergy) const {
+    RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
+    if (periodic)
+        ReferenceForce::getDeltaRPeriodic(atomCoordinates[atom2], atomCoordinates[atom1], periodicBoxSize, deltaR);
+    else
+        ReferenceForce::getDeltaR(atomCoordinates[atom2], atomCoordinates[atom1], deltaR);
+    RealOpenMM r = deltaR[ReferenceForce::RIndex];
+    if (cutoff && r >= cutoffDistance)
+        return;
+    map<string, double> variables = globalParameters;
+    for (int i = 0; i < paramNames.size(); i++) {
+        variables[particleParamNames[i*2]] = atomParameters[atom1][i];
+        variables[particleParamNames[i*2+1]] = atomParameters[atom2][i];
+    }
+    variables["r"] = r;
+    for (int i = 0; i < (int) valueNames.size(); i++) {
+        variables[particleValueNames[i*2]] = values[i][atom1].value;
+        variables[particleValueNames[i*2+1]] = values[i][atom2].value;
+    }
+    if (totalEnergy != NULL)
+        *totalEnergy += (RealOpenMM) energyExpressions[index].evaluate(variables);
+    RealOpenMM dEdR = (RealOpenMM) energyDerivExpressions[index][0].evaluate(variables);
+    dEdR *= 1/r;
+    for (int i = 0; i < 3; i++) {
+       forces[atom1][i] -= dEdR*deltaR[i];
+       forces[atom2][i] += dEdR*deltaR[i];
+    }
+    for (int i = 0; i < (int) valueNames.size(); i++) {
+        RealOpenMM scale1 = (RealOpenMM) energyDerivExpressions[index][2*i+1].evaluate(variables);
+        RealOpenMM scale2 = (RealOpenMM) energyDerivExpressions[index][2*i+2].evaluate(variables);
+        forces[atom1][0] -= scale1*values[i][atom1].gradient[0];
+        forces[atom1][1] -= scale1*values[i][atom1].gradient[1];
+        forces[atom1][2] -= scale1*values[i][atom1].gradient[2];
+        forces[atom2][0] += scale2*values[i][atom2].gradient[0];
+        forces[atom2][1] += scale2*values[i][atom2].gradient[1];
+        forces[atom2][2] += scale2*values[i][atom2].gradient[2];
+    }
+}

platforms/reference/src/SimTKReference/ReferenceCustomGBIxn.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 __ReferenceCustomGBIxn_H__
+#define __ReferenceCustomGBIxn_H__
+
+#include "ReferenceNeighborList.h"
+#include "lepton/ExpressionProgram.h"
+#include "openmm/CustomGBForce.h"
+#include <map>
+#include <set>
+#include <vector>
+
+// ---------------------------------------------------------------------------------------
+
+class ReferenceCustomGBIxn {
+
+   private:
+
+      bool cutoff;
+      bool periodic;
+      const OpenMM::NeighborList* neighborList;
+      RealOpenMM periodicBoxSize[3];
+      RealOpenMM cutoffDistance;
+      std::vector<Lepton::ExpressionProgram> valueExpressions;
+      std::vector<std::vector<Lepton::ExpressionProgram> > valueDerivExpressions;
+      std::vector<std::string> valueNames;
+      std::vector<OpenMM::CustomGBForce::ComputationType> valueTypes;
+      std::vector<Lepton::ExpressionProgram> energyExpressions;
+      std::vector<std::vector<Lepton::ExpressionProgram> > energyDerivExpressions;
+      std::vector<std::string> paramNames;
+      std::vector<OpenMM::CustomGBForce::ComputationType> energyTypes;
+      std::vector<std::string> particleParamNames;
+      std::vector<std::string> particleValueNames;
+      struct ComputedValue {
+          RealOpenMM value;
+          RealOpenMM gradient[3];
+          ComputedValue() {
+              value = (RealOpenMM) 0;
+              gradient[0] = (RealOpenMM) 0;
+              gradient[1] = (RealOpenMM) 0;
+              gradient[2] = (RealOpenMM) 0;
+          }
+      };
+
+      /**---------------------------------------------------------------------------------------
+
+         Calculate custom pair ixn between two atoms
+
+         @param atom1            the index of the first atom
+         @param atom2            the index of the second atom
+         @param atomCoordinates  atom coordinates
+         @param atomParameters   atomParameters[atomIndex][parameterIndex]
+         @param forces           force array (forces added)
+         @param energyByAtom     atom energy
+         @param totalEnergy      total energy
+
+         --------------------------------------------------------------------------------------- */
+
+      void calculateSingleParticleValue(int index, int numAtoms, std::vector<std::vector<ComputedValue> >& values,
+                                        const std::map<std::string, double>& globalParameters, RealOpenMM** atomParameters) const;
+
+      void calculateParticlePairValue(int index, int numAtoms, RealOpenMM** atomCoordinates, RealOpenMM** atomParameters,
+                                      std::vector<std::vector<ComputedValue> >& values,
+                                      const std::map<std::string, double>& globalParameters,
+                                      const std::vector<std::set<int> >& exclusions, bool useExclusions) const;
+
+      void calculateOnePairValue(int index, int atom1, int atom2, RealOpenMM** atomCoordinates, RealOpenMM** atomParameters,
+                                 const std::map<std::string, double>& globalParameters,
+                                 std::vector<std::vector<ComputedValue> >& values) const;
+
+      void calculateSingleParticleEnergyTerm(int index, int numAtoms, const std::vector<std::vector<ComputedValue> >& values,
+                                        const std::map<std::string, double>& globalParameters, RealOpenMM** atomParameters,
+                                        RealOpenMM** forces, RealOpenMM* totalEnergy) const;
+
+      void calculateParticlePairEnergyTerm(int index, int numAtoms, RealOpenMM** atomCoordinates, RealOpenMM** atomParameters,
+                                      const std::vector<std::vector<ComputedValue> >& values,
+                                      const std::map<std::string, double>& globalParameters,
+                                      const std::vector<std::set<int> >& exclusions, bool useExclusions,
+                                      RealOpenMM** forces, RealOpenMM* totalEnergy) const;
+
+      void calculateOnePairEnergyTerm(int index, int atom1, int atom2, RealOpenMM** atomCoordinates, RealOpenMM** atomParameters,
+                                 const std::map<std::string, double>& globalParameters,
+                                 const std::vector<std::vector<ComputedValue> >& values,
+                                 RealOpenMM** forces, RealOpenMM* totalEnergy) const;
+
+   public:
+
+      /**---------------------------------------------------------------------------------------
+
+         Constructor
+
+         --------------------------------------------------------------------------------------- */
+
+       ReferenceCustomGBIxn(const std::vector<Lepton::ExpressionProgram>& valueExpressions,
+                            const std::vector<std::vector<Lepton::ExpressionProgram> >& valueDerivExpressions,
+                            const std::vector<std::string>& valueNames,
+                            const std::vector<OpenMM::CustomGBForce::ComputationType>& valueTypes,
+                            const std::vector<Lepton::ExpressionProgram>& energyExpressions,
+                            const std::vector<std::vector<Lepton::ExpressionProgram> > energyDerivExpressions,
+                            const std::vector<OpenMM::CustomGBForce::ComputationType>& energyTypes,
+                            const std::vector<std::string>& parameterNames);
+
+      /**---------------------------------------------------------------------------------------
+
+         Destructor
+
+         --------------------------------------------------------------------------------------- */
+
+       ~ReferenceCustomGBIxn( );
+
+      /**---------------------------------------------------------------------------------------
+
+         Set the force to use a cutoff.
+
+         @param distance            the cutoff distance
+         @param neighbors           the neighbor list to use
+
+         @return ReferenceForce::DefaultReturn
+
+         --------------------------------------------------------------------------------------- */
+
+      int setUseCutoff( RealOpenMM distance, const OpenMM::NeighborList& neighbors );
+
+      /**---------------------------------------------------------------------------------------
+
+         Set the force to use periodic boundary conditions.  This requires that a cutoff has
+         already been set, and the smallest side of the periodic box is at least twice the cutoff
+         distance.
+
+         @param boxSize             the X, Y, and Z widths of the periodic box
+
+         @return ReferenceForce::DefaultReturn
+
+         --------------------------------------------------------------------------------------- */
+
+      int setPeriodic( RealOpenMM* boxSize );
+
+      /**---------------------------------------------------------------------------------------
+
+         Calculate custom pair ixn
+
+         @param numberOfAtoms    number of atoms
+         @param atomCoordinates  atom coordinates
+         @param atomParameters   atom parameters (charges, c6, c12, ...)     atomParameters[atomIndex][paramterIndex]
+         @param exclusions       atom exclusion indices                      exclusions[atomIndex][atomToExcludeIndex]
+                                 exclusions[atomIndex][0] = number of exclusions
+                                 exclusions[atomIndex][1-no.] = atom indices of atoms to excluded from
+                                 interacting w/ atom atomIndex
+         @param fixedParameters  non atom parameters (not currently used)
+         @param globalParameters the values of global parameters
+         @param forces           force array (forces added)
+         @param energyByAtom     atom energy
+         @param totalEnergy      total energy
+
+         @return ReferenceForce::DefaultReturn
+
+         --------------------------------------------------------------------------------------- */
+
+      int calculateIxn(int numberOfAtoms, RealOpenMM** atomCoordinates, RealOpenMM** atomParameters, const std::vector<std::set<int> >& exclusions,
+                       std::map<std::string, double>& globalParameters, RealOpenMM** forces,
+                       RealOpenMM* energyByAtom, RealOpenMM* totalEnergy) const;
+
+// ---------------------------------------------------------------------------------------
+
+};
+
+#endif // __ReferenceCustomGBIxn_H__

platforms/reference/src/SimTKReference/ReferenceCustomNonbondedIxn.cpp

@@ -147,7 +147,7 @@ ReferenceCustomNonbondedIxn::~ReferenceCustomNonbondedIxn( ){
 
 int ReferenceCustomNonbondedIxn::calculatePairIxn( int numberOfAtoms, RealOpenMM** atomCoordinates,
                                              RealOpenMM** atomParameters, int** exclusions,
-                                             RealOpenMM* fixedParameters, map<string, double> globalParameters, RealOpenMM** forces,
+                                             RealOpenMM* fixedParameters, const map<string, double>& globalParameters, RealOpenMM** forces,
                                              RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const {
 
    map<string, double> variables = globalParameters;

platforms/reference/src/SimTKReference/ReferenceCustomNonbondedIxn.h

@@ -135,7 +135,7 @@ class ReferenceCustomNonbondedIxn {
 
       int calculatePairIxn( int numberOfAtoms, RealOpenMM** atomCoordinates,
                             RealOpenMM** atomParameters, int** exclusions,
-                            RealOpenMM* fixedParameters, std::map<std::string, double> globalParameters,
+                            RealOpenMM* fixedParameters, const std::map<std::string, double>& globalParameters,
                             RealOpenMM** forces, RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const;
 
 // ---------------------------------------------------------------------------------------

platforms/reference/src/gbsa/CpuObc.cpp

@@ -624,7 +624,7 @@ int CpuObc::computeBornEnergyForces( RealOpenMM* bornRadii, RealOpenMM** atomCoo
       RealOpenMM sum3            = bornSum*sum2;
       RealOpenMM tanhSum         = TANH( alphaObc*bornSum - betaObc*sum2 + gammaObc*sum3 );
       
-      bornRadiiTemp[atomI]       = one/( one/offsetRadiusI - tanhSum/radiusI ); 
+      bornRadiiTemp[atomI]       = one/( one/offsetRadiusI - tanhSum/radiusI );
  
       obcChainTemp[atomI]        = offsetRadiusI*( alphaObc - two*betaObc*bornSum + three*gammaObc*sum2 );
       obcChainTemp[atomI]        = (one - tanhSum*tanhSum)*obcChainTemp[atomI]/radiusI;

platforms/reference/tests/TestReferenceCustomGBForce.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 CustomGBForce.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "../src/sfmt/SFMT.h"
+#include "openmm/Context.h"
+#include "ReferencePlatform.h"
+#include "openmm/CustomGBForce.h"
+#include "openmm/GBSAOBCForce.h"
+#include "openmm/System.h"
+#include "openmm/VerletIntegrator.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+const double TOL = 1e-5;
+
+void testOBC() {
+    const int numMolecules = 100;
+    const int numParticles = numMolecules*2;
+    const double boxSize = 10.0;
+    ReferencePlatform platform;
+
+    // Create two systems: one with a GBSAOBCForce, and one using a CustomGBForce to implement the same interaction.
+
+    System standardSystem;
+    System customSystem;
+    for (int i = 0; i < numParticles; i++) {
+        standardSystem.addParticle(1.0);
+        customSystem.addParticle(1.0);
+    }
+    GBSAOBCForce* obc = new GBSAOBCForce();
+    CustomGBForce* custom = new CustomGBForce();
+    custom->addPerParticleParameter("q");
+    custom->addPerParticleParameter("radius");
+    custom->addPerParticleParameter("scale");
+    custom->addGlobalParameter("solventDielectric", obc->getSolventDielectric());
+    custom->addGlobalParameter("soluteDielectric", obc->getSoluteDielectric());
+    custom->addComputedValue("I", "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r*(1/U^2-1/L^2))+0.5*log(L/U)/r+0.25*(sr2*sr2/r)*(1/L^2-1/U^2)+C);"
+                                  "U=r+sr2;"
+                                  "C=2*(1/or1-1/L)*step(sr2-r-or1);"
+                                  "L=step(or1-D)*or1+step(D-or1)*D;"
+                                  "D=step(r-sr2)*(r-sr2)+step(sr2-r)*(sr2-r);"
+                                  "sr1 = scale1*or1; sr2 = scale2*or2;"
+                                  "or1 = radius1-0.009; or2 = radius2-0.009", CustomGBForce::ParticlePairNoExclusions);
+    custom->addComputedValue("B", "1/(1/or-tanh(1*psi-0.8*psi^2+4.85*psi^3)/radius);"
+                                  "psi=I*or; or=radius-0.009", CustomGBForce::SingleParticle);
+    custom->addEnergyTerm("28.3919551*(radius+0.14)^2*(radius/B)^6-0.5*138.935485*(1/soluteDielectric-1/solventDielectric)*q^2/B", CustomGBForce::SingleParticle);
+    custom->addEnergyTerm("-138.935485*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;"
+                          "f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce::ParticlePairNoExclusions);
+    vector<Vec3> positions(numParticles);
+    vector<Vec3> velocities(numParticles);
+    init_gen_rand(0);
+    vector<double> params(3);
+    for (int i = 0; i < numMolecules; i++) {
+        if (i < numMolecules/2) {
+            obc->addParticle(1.0, 0.2, 0.5);
+            params[0] = 1.0;
+            params[1] = 0.2;
+            params[2] = 0.5;
+            custom->addParticle(params);
+            obc->addParticle(-1.0, 0.1, 0.5);
+            params[0] = -1.0;
+            params[1] = 0.1;
+            custom->addParticle(params);
+        }
+        else {
+            obc->addParticle(1.0, 0.2, 0.8);
+            params[0] = 1.0;
+            params[1] = 0.2;
+            params[2] = 0.8;
+            custom->addParticle(params);
+            obc->addParticle(-1.0, 0.1, 0.8);
+            params[0] = -1.0;
+            params[1] = 0.1;
+            custom->addParticle(params);
+        }
+        positions[2*i] = Vec3(boxSize*genrand_real2(), boxSize*genrand_real2(), boxSize*genrand_real2());
+        positions[2*i+1] = Vec3(positions[2*i][0]+1.0, positions[2*i][1], positions[2*i][2]);
+        velocities[2*i] = Vec3(genrand_real2(), genrand_real2(), genrand_real2());
+        velocities[2*i+1] = Vec3(genrand_real2(), genrand_real2(), genrand_real2());
+    }
+    obc->setNonbondedMethod(GBSAOBCForce::NoCutoff);
+    custom->setNonbondedMethod(CustomGBForce::NoCutoff);
+    standardSystem.addForce(obc);
+    customSystem.addForce(custom);
+    VerletIntegrator integrator1(0.01);
+    VerletIntegrator integrator2(0.01);
+    Context context1(standardSystem, integrator1, platform);
+    Context context2(customSystem, integrator2, platform);
+    context1.setPositions(positions);
+    context2.setPositions(positions);
+    context1.setVelocities(velocities);
+    context2.setVelocities(velocities);
+    State state1 = context1.getState(State::Forces | State::Energy);
+    State state2 = context2.getState(State::Forces | State::Energy);
+    for (int i = 0; i < numParticles; i++)
+        std::cout << state1.getForces()[i]<< state2.getForces()[i]<< std::endl;
+    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-4);
+    for (int i = 0; i < numParticles; i++) {
+        ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-3);
+    }
+}
+
+int main() {
+    try {
+        testOBC();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}
+