Finished reference implementation of CustomGBForce

openmmapi/include/openmm/CustomGBForce.h

@@ -52,12 +52,12 @@ namespace OpenMM {
  * The computation consists of calculating some number of per-particle <i>computed values</i>, followed by one or more
  * <i>energy terms</i>.  A computed value is a scalar value that is computed for each particle in the system.  It may
  * depend on an arbitrary set of global and per-particle parameters, and well as on other computed values that have
- * been calculating before it.  Once all computed values have been calculated, the energy terms and their derivatives
+ * been calculated before it.  Once all computed values have been calculated, the energy terms and their derivatives
  * are evaluated to determine the system energy and particle forces.  The energy terms may depend on global parameters,
  * per-particle parameters, and per-particle computed values.
  *
  * When specifying a computed value or energy term, you provide an algebraic expression to evaluate and a <i>computation type</i>
- * describing how is the expression is to be evaluated.  There are two main types of computations:
+ * describing how the expression is to be evaluated.  There are two main types of computations:
  *
  * <ul>
  * <li><b>Single Particle</b>: The expression is evaluated once for each particle in the System.  In the case of a computed
@@ -68,14 +68,14 @@ namespace OpenMM {
  * value, the value for a particular particle is calculated by pairing it with every other particle in the system, evaluating
  * the expression for each pair, and summing them.  For an energy term, each particle pair makes an independent contribution to
  * the System energy.  (Note that energy terms are assumed to be symmetric with respect to the two interacting particles, and
- * therefore are evaluated only once per pair.  In contrast, computed values need not be symmetric and therefore are calculated
+ * therefore are evaluated only once per pair.  In contrast, expressions for computed values need not be symmetric and therefore are calculated
  * twice for each pair: once when calculating the value for the first particle, and again when calculating the value for the
  * second particle.)</li>
  * </ul>
  *
  * Be aware that, although this class is extremely general in the computations it can define, particular Platforms may only support
  * more restricted types of computations.  In particular, all currently existing Platforms require that the first computed value
- * <i>must</i> be a particle pair computation, and all computed values that follow it <i>must</i> be single particle computations.
+ * <i>must</i> be a particle pair computation, and all computed values after the first <i>must</i> be single particle computations.
  * This is sufficient for most Generalized Born models, but might not permit some other types of calculations to be implemented.
  *
  * This is a complicated class to use, and an example may help to clarify it.  The following code implements the OBC variant
@@ -107,7 +107,7 @@ namespace OpenMM {
  * (the dielectric constants for the solute and solvent).  It then defines a computed value "I" of type ParticlePair.  The
  * expression for evaluating it is a complicated function of the distance between each pair of particles (r), their atomic
  * radii (radius1 and radius2), and their scale factors (scale1 and scale2).  Very roughly speaking, it is a measure of the
- * overlap between each particle and other nearby particles.
+ * distance between each particle and other nearby particles.
  *
  * Next a computation is defined for the Born Radius (B).  It is computed independently for each particle, and is a function of
  * that particle's atomic radius and the intermediate value I defined above.
@@ -128,12 +128,13 @@ namespace OpenMM {
  *
  * 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, step.  All trigonometric functions
- * are defined in radians, and log is the natural logarithm.  step(x) = 0 if x is less than 0, 1 otherwise.  The names of per-particle parameters
+ * are defined in radians, and log is the natural logarithm.  step(x) = 0 if x is less than 0, 1 otherwise.  In expressions for
+ * particle pair calculations, the names of per-particle parameters and computed values
  * have the suffix "1" or "2" appended to them to indicate the values for the two interacting particles.  As seen in the above example,
- * the expression may also involve intermediate quantities that are defined following the main expression, using ";" as a separator.
+ * an expression may also involve intermediate quantities that are defined following the main expression, using ";" as a separator.
  *
  * In addition, you can call addFunction() to define a new function based on tabulated values.  You specify a vector of
- * values, and an interpolating or approximating spline is created from them.  That function can then appear in the expression.
+ * values, and an interpolating or approximating spline is created from them.  That function can then appear in expressions.
  */
 
 class OPENMM_EXPORT CustomGBForce : public Force {

platforms/reference/src/ReferenceKernels.cpp

@@ -607,9 +607,9 @@ double ReferenceCalcNonbondedForceKernel::executeEnergy(ContextImpl& context) {
     return energy;
 }
 
-class ReferenceCalcCustomNonbondedForceKernel::TabulatedFunction : public Lepton::CustomFunction {
+class ReferenceTabulatedFunction : public Lepton::CustomFunction {
 public:
-    TabulatedFunction(double min, double max, const vector<double>& values, bool interpolating) :
+    ReferenceTabulatedFunction(double min, double max, const vector<double>& values, bool interpolating) :
             min(min), max(max), values(values), interpolating(interpolating) {
     }
     int getNumArguments() const {
@@ -659,7 +659,7 @@ public:
         return scale*(coeff[1]+x*(2.0*coeff[2]+x*3.0*coeff[3])); // We assume a first derivative, because that's the only order ever used by CustomNonbondedForce.
     }
     CustomFunction* clone() const {
-        return new TabulatedFunction(min, max, values, interpolating);
+        return new ReferenceTabulatedFunction(min, max, values, interpolating);
     }
     double min, max;
     vector<double> values;
@@ -725,7 +725,7 @@ void ReferenceCalcCustomNonbondedForceKernel::initialize(const System& system, c
         double min, max;
         bool interpolating;
         force.getFunctionParameters(i, name, values, min, max, interpolating);
-        functions[name] = new TabulatedFunction(min, max, values, interpolating);
+        functions[name] = new ReferenceTabulatedFunction(min, max, values, interpolating);
     }
 
     // Parse the various expressions used to calculate the force.
@@ -953,14 +953,14 @@ void ReferenceCalcCustomGBForceKernel::initialize(const System& system, const Cu
     // 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);
-//    }
+    for (int i = 0; i < force.getNumFunctions(); i++) {
+        string name;
+        vector<double> values;
+        double min, max;
+        bool interpolating;
+        force.getFunctionParameters(i, name, values, min, max, interpolating);
+        functions[name] = new ReferenceTabulatedFunction(min, max, values, interpolating);
+    }
 
     // Parse the expressions for computed values.
 
@@ -978,7 +978,7 @@ void ReferenceCalcCustomGBForceKernel::initialize(const System& system, const Cu
             valueDerivExpressions.push_back(ex.differentiate("r").optimize().createProgram());
     }
 
-    // Parse the various expressions used to calculate the force.
+    // Parse the expressions for energy terms.
 
     energyDerivExpressions.resize(force.getNumEnergyTerms());
     for (int i = 0; i < force.getNumEnergyTerms(); i++) {

platforms/reference/src/ReferenceKernels.h

@@ -407,7 +407,6 @@ private:
     std::vector<std::string> parameterNames, globalParameterNames;
     NonbondedMethod nonbondedMethod;
     NeighborList* neighborList;
-    class TabulatedFunction;
 };
 
 /**

platforms/reference/tests/TestReferenceCustomGBForce.cpp

@@ -139,11 +139,46 @@ void testOBC(GBSAOBCForce::NonbondedMethod obcMethod, CustomGBForce::NonbondedMe
     }
 }
 
+void testTabulatedFunction(bool interpolating) {
+    ReferencePlatform platform;
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    CustomGBForce* force = new CustomGBForce();
+    force->addComputedValue("a", "0", CustomGBForce::ParticlePair);
+    force->addEnergyTerm("fn(r)+1", CustomGBForce::ParticlePair);
+    force->addParticle(vector<double>());
+    force->addParticle(vector<double>());
+    vector<double> table;
+    for (int i = 0; i < 21; i++)
+        table.push_back(std::sin(0.25*i));
+    force->addFunction("fn", table, 1.0, 6.0, interpolating);
+    system.addForce(force);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(0, 0, 0);
+    for (int i = 1; i < 30; i++) {
+        double x = (7.0/30.0)*i;
+        positions[1] = Vec3(x, 0, 0);
+        context.setPositions(positions);
+        State state = context.getState(State::Forces | State::Energy);
+        const vector<Vec3>& forces = state.getForces();
+        double force = (x < 1.0 || x > 6.0 ? 0.0 : -std::cos(x-1.0));
+        double energy = (x < 1.0 || x > 6.0 ? 0.0 : std::sin(x-1.0))+1.0;
+        ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], 0.1);
+        ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], 0.1);
+        ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), 0.02);
+    }
+}
+
 int main() {
     try {
         testOBC(GBSAOBCForce::NoCutoff, CustomGBForce::NoCutoff);
         testOBC(GBSAOBCForce::CutoffNonPeriodic, CustomGBForce::CutoffNonPeriodic);
         testOBC(GBSAOBCForce::CutoffPeriodic, CustomGBForce::CutoffPeriodic);
+        testTabulatedFunction(true);
+        testTabulatedFunction(false);
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;