Began CUDA implementation of periodicdistance() function for CustomExternalForce

platforms/cuda/include/CudaExpressionUtilities.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2009-2014 Stanford University and the Authors.      *
+ * Portions copyright (c) 2009-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -89,6 +89,10 @@ public:
      * @param function   the function for which to get a placeholder
      */
     Lepton::CustomFunction* getFunctionPlaceholder(const TabulatedFunction& function);
+    /**
+     * Get a Lepton::CustomFunction that can be used to represent the periodicdistance() function when parsing expressions.
+     */
+    Lepton::CustomFunction* getPeriodicDistancePlaceholder();
 private:
     class FunctionPlaceholder : public Lepton::CustomFunction {
         public:
@@ -114,13 +118,13 @@ private:
             const std::vector<const TabulatedFunction*>& functions, const std::vector<std::pair<std::string, std::string> >& functionNames,
             const std::string& prefix, const std::vector<std::vector<double> >& functionParams, const std::vector<Lepton::ParsedExpression>& allExpressions, const std::string& tempType);
     std::string getTempName(const Lepton::ExpressionTreeNode& node, const std::vector<std::pair<Lepton::ExpressionTreeNode, std::string> >& temps);
-    void findRelatedTabulatedFunctions(const Lepton::ExpressionTreeNode& node, const Lepton::ExpressionTreeNode& searchNode,
+    void findRelatedCustomFunctions(const Lepton::ExpressionTreeNode& node, const Lepton::ExpressionTreeNode& searchNode,
             std::vector<const Lepton::ExpressionTreeNode*>& nodes);
     void findRelatedPowers(const Lepton::ExpressionTreeNode& node, const Lepton::ExpressionTreeNode& searchNode,
             std::map<int, const Lepton::ExpressionTreeNode*>& powers);
     std::vector<std::vector<double> > computeFunctionParameters(const std::vector<const TabulatedFunction*>& functions);
     CudaContext& context;
-    FunctionPlaceholder fp1, fp2, fp3;
+    FunctionPlaceholder fp1, fp2, fp3, periodicDistance;
 };
 
 } // namespace OpenMM

platforms/cuda/src/CudaKernels.cpp

@@ -3652,7 +3652,9 @@ void CudaCalcCustomExternalForceKernel::initialize(const System& system, const C
         globalParamNames[i] = force.getGlobalParameterName(i);
         globalParamValues[i] = (float) force.getGlobalParameterDefaultValue(i);
     }
-    Lepton::ParsedExpression energyExpression = Lepton::Parser::parse(force.getEnergyFunction()).optimize();
+    map<string, Lepton::CustomFunction*> customFunctions;
+    customFunctions["periodicdistance"] = cu.getExpressionUtilities().getPeriodicDistancePlaceholder();
+    Lepton::ParsedExpression energyExpression = Lepton::Parser::parse(force.getEnergyFunction(), customFunctions).optimize();
     Lepton::ParsedExpression forceExpressionX = energyExpression.differentiate("x").optimize();
     Lepton::ParsedExpression forceExpressionY = energyExpression.differentiate("y").optimize();
     Lepton::ParsedExpression forceExpressionZ = energyExpression.differentiate("z").optimize();

platforms/cuda/tests/TestCudaCustomExternalForce.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -161,6 +161,47 @@ void testParallelComputation() {
         ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-5);
 }
 
+void testPeriodic() {
+    Vec3 vx(5, 0, 0);
+    Vec3 vy(0, 6, 0);
+    Vec3 vz(1, 2, 7);
+    double x0 = 51, y0 = -17, z0 = 11.2;
+    System system;
+    system.setDefaultPeriodicBoxVectors(vx, vy, vz);
+    system.addParticle(1.0);
+    CustomExternalForce* force = new CustomExternalForce("periodicdistance(x, y, z, x0, y0, z0)^2");
+    force->addPerParticleParameter("x0");
+    force->addPerParticleParameter("y0");
+    force->addPerParticleParameter("z0");
+    vector<double> params(3);
+    params[0] = x0;
+    params[1] = y0;
+    params[2] = z0;
+    force->addParticle(0, params);
+    system.addForce(force);
+    VerletIntegrator integrator(0.01);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(1);
+    positions[0] = Vec3(0, 2, 0);
+    context.setPositions(positions);
+    for (int i = 0; i < 100; i++) {
+        State state = context.getState(State::Positions | State::Forces | State::Energy);
+
+        // Apply periodic boundary conditions to the difference between the two positions.
+
+        Vec3 delta = Vec3(x0, y0, z0)-state.getPositions()[0];
+        delta -= vz*floor(delta[2]/vz[2]+0.5);
+        delta -= vy*floor(delta[1]/vy[1]+0.5);
+        delta -= vx*floor(delta[0]/vx[0]+0.5);
+
+        // Verify that the force and energy are correct.
+
+        ASSERT_EQUAL_VEC(delta*2, state.getForces()[0], 1e-6);
+        ASSERT_EQUAL_TOL(delta.dot(delta), state.getPotentialEnergy(), 1e-6);
+        integrator.step(1);
+    }
+}
+
 int main(int argc, char* argv[]) {
     try {
         if (argc > 1)
@@ -168,6 +209,7 @@ int main(int argc, char* argv[]) {
         testForce();
         testManyParameters();
         testParallelComputation();
+        testPeriodic();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;

platforms/reference/tests/TestReferenceCustomExternalForce.cpp

@@ -6,7 +6,7 @@
  * 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.      *
+ * Portions copyright (c) 2008-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *