A few optimizations. Also added a new test case.

platforms/opencl/src/OpenCLExpressionUtilities.cpp

@@ -177,8 +177,34 @@ void OpenCLExpressionUtilities::processExpression(stringstream& out, const Expre
             out << doubleToString(dynamic_cast<const Operation::MultiplyConstant*>(&node.getOperation())->getValue()) << "*" << getTempName(node.getChildren()[0], temps);
             break;
         case Operation::POWER_CONSTANT:
-            out << "pow(" << getTempName(node.getChildren()[0], temps) << ", " << doubleToString(dynamic_cast<const Operation::PowerConstant*>(&node.getOperation())->getValue()) << ")";
+        {
+            double exponent = dynamic_cast<const Operation::PowerConstant*>(&node.getOperation())->getValue();
+            if (exponent == 0.0)
+                out << "1.0f";
+            else if (exponent == (int) exponent) {
+                out << "0.0f;\n";
+                out << "{\n";
+                out << "float multiplier = " << (exponent < 0.0 ? "1.0f/" : "") << getTempName(node.getChildren()[0], temps) << ";\n";
+                int exp = (int) fabs(exponent);
+                bool hasAssigned = false;
+                while (exp != 0) {
+                    if (exp%2 == 1) {
+                        if (!hasAssigned)
+                            out << name << " = multiplier;\n";
+                        else
+                            out << name << " *= multiplier;\n";
+                        hasAssigned = true;
+                    }
+                    exp >>= 1;
+                    if (exp != 0)
+                        out << "multiplier *= multiplier;\n";
+                }
+                out << "}";
+            }
+            else
+                out << "pow(" << getTempName(node.getChildren()[0], temps) << ", " << doubleToString(exponent) << ")";
             break;
+        }
         default:
             throw OpenMMException("Internal error: Unknown operation in user-defined expression: "+node.getOperation().getName());
     }

platforms/opencl/src/OpenCLKernels.cpp

@@ -553,6 +553,8 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
     vector<mm_float2> sigmaEpsilonVector(numParticles);
     vector<vector<int> > exclusionList(numParticles);
     double sumSquaredCharges = 0.0;
+    bool hasCoulomb = false;
+    bool hasLJ = false;
     for (int i = 0; i < numParticles; i++) {
         double charge, sigma, epsilon;
         force.getParticleParameters(i, charge, sigma, epsilon);
@@ -560,6 +562,10 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
         sigmaEpsilonVector[i] = (mm_float2) {(float) (0.5*sigma), (float) (2.0*sqrt(epsilon))};
         exclusionList[i].push_back(i);
         sumSquaredCharges += charge*charge;
+        if (charge != 0.0)
+            hasCoulomb = true;
+        if (epsilon != 0.0)
+            hasLJ = true;
     }
     for (int i = 0; i < (int) exclusions.size(); i++) {
         exclusionList[exclusions[i].first].push_back(exclusions[i].second);
@@ -572,6 +578,8 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
     Vec3 boxVectors[3];
     system.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
     map<string, string> defines;
+    defines["HAS_COULOMB"] = (hasCoulomb ? "1" : "0");
+    defines["HAS_LENNARD_JONES"] = (hasLJ ? "1" : "0");
     if (useCutoff) {
         // Compute the reaction field constants.
 

platforms/opencl/src/kernels/coulombLennardJones.cl

@@ -13,6 +13,7 @@ if (r2 < CUTOFF_SQUARED) {
             tempEnergy += -prefactor*(1.0f-erfcAlphaR);
         }
         else {
+#if HAS_LENNARD_JONES
             float sig = sigmaEpsilon1.x + sigmaEpsilon2.x;
             float sig2 = invR*sig;
             sig2 *= sig2;
@@ -20,6 +21,10 @@ if (r2 < CUTOFF_SQUARED) {
             float eps = sigmaEpsilon1.y*sigmaEpsilon2.y;
             tempForce = eps*(12.0f*sig6 - 6.0f)*sig6 + prefactor*(erfcAlphaR+alphaR*exp(-alphaR*alphaR)*TWO_OVER_SQRT_PI);
             tempEnergy += eps*(sig6 - 1.0f)*sig6 + prefactor*erfcAlphaR;
+#else
+            tempForce = prefactor*(erfcAlphaR+alphaR*exp(-alphaR*alphaR)*TWO_OVER_SQRT_PI);
+            tempEnergy += prefactor*erfcAlphaR;
+#endif
         }
         dEdR += tempForce*invR*invR;
     }
@@ -30,21 +35,26 @@ if (!isExcluded && r2 < CUTOFF_SQUARED) {
 #else
 if (!isExcluded) {
 #endif
+    float tempForce = 0.0f;
+  #if HAS_LENNARD_JONES
     float sig = sigmaEpsilon1.x + sigmaEpsilon2.x;
     float sig2 = invR*sig;
     sig2 *= sig2;
     float sig6 = sig2*sig2*sig2;
     float eps = sigmaEpsilon1.y*sigmaEpsilon2.y;
-    float tempForce = eps*(12.0f*sig6 - 6.0f)*sig6;
+    tempForce = eps*(12.0f*sig6 - 6.0f)*sig6;
     tempEnergy += eps*(sig6 - 1.0f)*sig6;
-#ifdef USE_CUTOFF
+  #endif
+#if HAS_COULOMB
+  #ifdef USE_CUTOFF
     const float prefactor = 138.935485f*posq1.w*posq2.w;
     tempForce += prefactor*(invR - 2.0f*REACTION_FIELD_K*r2);
     tempEnergy += prefactor*(invR + REACTION_FIELD_K*r2 - REACTION_FIELD_C);
-#else
+  #else
     const float prefactor = 138.935485f*posq1.w*posq2.w*invR;
     tempForce += prefactor;
     tempEnergy += prefactor;
+  #endif
 #endif
     dEdR += tempForce*invR*invR;
 }

platforms/opencl/tests/TestOpenCLCustomNonbondedForce.cpp

@@ -35,9 +35,11 @@
  */
 
 #include "../../../tests/AssertionUtilities.h"
+#include "../src/sfmt/SFMT.h"
 #include "openmm/Context.h"
 #include "OpenCLPlatform.h"
 #include "openmm/CustomNonbondedForce.h"
+#include "openmm/NonbondedForce.h"
 #include "openmm/System.h"
 #include "openmm/VerletIntegrator.h"
 #include <iostream>
@@ -234,6 +236,77 @@ void testTabulatedFunction(bool interpolating) {
     }
 }
 
+void testCoulombLennardJones() {
+    const int numMolecules = 300;
+    const int numParticles = numMolecules*2;
+    const double boxSize = 20.0;
+    OpenCLPlatform platform;
+
+    // Create two systems: one with a NonbondedForce, and one using a CustomNonbondedForce to implement the same interaction.
+
+    System standardSystem;
+    System customSystem;
+    for (int i = 0; i < numParticles; i++) {
+        standardSystem.addParticle(1.0);
+        customSystem.addParticle(1.0);
+    }
+    NonbondedForce* standardNonbonded = new NonbondedForce();
+    CustomNonbondedForce* customNonbonded = new CustomNonbondedForce("4*eps*((sigma/r)^12-(sigma/r)^6) + 138.935485*q/r");
+    customNonbonded->addParameter("q", "q1*q2");
+    customNonbonded->addParameter("sigma", "0.5*(sigma1+sigma2)");
+    customNonbonded->addParameter("eps", "sqrt(eps1*eps2)");
+    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) {
+            standardNonbonded->addParticle(1.0, 0.2, 0.1);
+            params[0] = 1.0;
+            params[1] = 0.2;
+            params[2] = 0.1;
+            customNonbonded->addParticle(params);
+            standardNonbonded->addParticle(1.0, 0.1, 0.1);
+            params[1] = 0.1;
+            customNonbonded->addParticle(params);
+        }
+        else {
+            standardNonbonded->addParticle(1.0, 0.2, 0.2);
+            params[0] = 1.0;
+            params[1] = 0.2;
+            params[2] = 0.2;
+            customNonbonded->addParticle(params);
+            standardNonbonded->addParticle(1.0, 0.1, 0.2);
+            params[1] = 0.1;
+            customNonbonded->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());
+        standardNonbonded->addException(2*i, 2*i+1, 0.0, 1.0, 0.0);
+        customNonbonded->addException(2*i, 2*i+1, vector<double>());
+    }
+    standardNonbonded->setNonbondedMethod(NonbondedForce::NoCutoff);
+    customNonbonded->setNonbondedMethod(CustomNonbondedForce::NoCutoff);
+    standardSystem.addForce(standardNonbonded);
+    customSystem.addForce(customNonbonded);
+    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);
+    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-4);
+    }
+}
+
 int main() {
     try {
         testSimpleExpression();
@@ -243,6 +316,7 @@ int main() {
         testPeriodic();
         testTabulatedFunction(true);
         testTabulatedFunction(false);
+        testCoulombLennardJones();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;

platforms/reference/tests/TestReferenceCustomNonbondedForce.cpp

@@ -7,7 +7,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Portions copyright (c) 2008-2009 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -35,9 +35,11 @@
  */
 
 #include "../../../tests/AssertionUtilities.h"
+#include "../src/sfmt/SFMT.h"
 #include "openmm/Context.h"
 #include "ReferencePlatform.h"
 #include "openmm/CustomNonbondedForce.h"
+#include "openmm/NonbondedForce.h"
 #include "openmm/System.h"
 #include "openmm/VerletIntegrator.h"
 #include <iostream>
@@ -234,6 +236,77 @@ void testTabulatedFunction(bool interpolating) {
     }
 }
 
+void testCoulombLennardJones() {
+    const int numMolecules = 300;
+    const int numParticles = numMolecules*2;
+    const double boxSize = 20.0;
+    ReferencePlatform platform;
+
+    // Create two systems: one with a NonbondedForce, and one using a CustomNonbondedForce to implement the same interaction.
+
+    System standardSystem;
+    System customSystem;
+    for (int i = 0; i < numParticles; i++) {
+        standardSystem.addParticle(1.0);
+        customSystem.addParticle(1.0);
+    }
+    NonbondedForce* standardNonbonded = new NonbondedForce();
+    CustomNonbondedForce* customNonbonded = new CustomNonbondedForce("4*eps*((sigma/r)^12-(sigma/r)^6) + 138.935485*q/r");
+    customNonbonded->addParameter("q", "q1*q2");
+    customNonbonded->addParameter("sigma", "0.5*(sigma1+sigma2)");
+    customNonbonded->addParameter("eps", "sqrt(eps1*eps2)");
+    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) {
+            standardNonbonded->addParticle(1.0, 0.2, 0.1);
+            params[0] = 1.0;
+            params[1] = 0.2;
+            params[2] = 0.1;
+            customNonbonded->addParticle(params);
+            standardNonbonded->addParticle(1.0, 0.1, 0.1);
+            params[1] = 0.1;
+            customNonbonded->addParticle(params);
+        }
+        else {
+            standardNonbonded->addParticle(1.0, 0.2, 0.2);
+            params[0] = 1.0;
+            params[1] = 0.2;
+            params[2] = 0.2;
+            customNonbonded->addParticle(params);
+            standardNonbonded->addParticle(1.0, 0.1, 0.2);
+            params[1] = 0.1;
+            customNonbonded->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());
+        standardNonbonded->addException(2*i, 2*i+1, 0.0, 1.0, 0.0);
+        customNonbonded->addException(2*i, 2*i+1, vector<double>());
+    }
+    standardNonbonded->setNonbondedMethod(NonbondedForce::NoCutoff);
+    customNonbonded->setNonbondedMethod(CustomNonbondedForce::NoCutoff);
+    standardSystem.addForce(standardNonbonded);
+    customSystem.addForce(customNonbonded);
+    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);
+    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-4);
+    }
+}
+
 int main() {
     try {
         testSimpleExpression();
@@ -243,6 +316,7 @@ int main() {
         testPeriodic();
         testTabulatedFunction(true);
         testTabulatedFunction(false);
+        testCoulombLennardJones();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;