Merge pull request #1137 from peastman/periodic

Created periodicdistance() function for CustomExternalForce

openmmapi/include/openmm/CustomExternalForce.h

@@ -67,6 +67,14 @@ namespace OpenMM {
  * force->addPerParticleParameter("z0");
  * </pre></tt>
  *
+ * Special care is needed in systems that use periodic boundary conditions.  In that case, each particle really represents
+ * an infinite set of particles repeating through space.  The variables x, y, and z contain the coordinates of one of those
+ * periodic copies, but there is no guarantee about which.  It might even change from one time step to the next.  You can handle
+ * this situation by using the function periodicdistance(x1, y1, z1, x2, y2, z2), which returns the minimum distance between
+ * periodic copies of the points (x1, y1, z1) and (x2, y2, z2).  For example, the force given above would be rewritten as
+ *
+ * <tt>CustomExternalForce* force = new CustomExternalForce("k*periodicdistance(x, y, z, x0, y0, z0)^2");</tt>
+ *
  * Expressions may involve the operators + (add), - (subtract), * (multiply), / (divide), and ^ (power), and the following
  * functions: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, erf, erfc, min, max, abs, floor, ceil, step, delta, select.  All trigonometric functions
  * are defined in radians, and log is the natural logarithm.  step(x) = 0 if x is less than 0, 1 otherwise.  delta(x) = 1 if x is 0, 0 otherwise.

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/CudaBondedUtilities.cpp

@@ -99,7 +99,7 @@ void CudaBondedUtilities::initialize(const System& system) {
     s<<CudaKernelSources::vectorOps;
     for (int i = 0; i < (int) prefixCode.size(); i++)
         s<<prefixCode[i];
-    s<<"extern \"C\" __global__ void computeBondedForces(unsigned long long* __restrict__ forceBuffer, real* __restrict__ energyBuffer, const real4* __restrict__ posq, int groups";
+    s<<"extern \"C\" __global__ void computeBondedForces(unsigned long long* __restrict__ forceBuffer, real* __restrict__ energyBuffer, const real4* __restrict__ posq, int groups, real4 periodicBoxSize, real4 invPeriodicBoxSize, real4 periodicBoxVecX, real4 periodicBoxVecY, real4 periodicBoxVecZ";
     for (int force = 0; force < numForces; force++) {
         for (int i = 0; i < (int) atomIndices[force].size(); i++) {
             int indexWidth = atomIndices[force][i]->getElementSize()/4;
@@ -161,6 +161,11 @@ void CudaBondedUtilities::computeInteractions(int groups) {
         kernelArgs.push_back(&context.getEnergyBuffer().getDevicePointer());
         kernelArgs.push_back(&context.getPosq().getDevicePointer());
         kernelArgs.push_back(NULL);
+        kernelArgs.push_back(context.getPeriodicBoxSizePointer());
+        kernelArgs.push_back(context.getInvPeriodicBoxSizePointer());
+        kernelArgs.push_back(context.getPeriodicBoxVecXPointer());
+        kernelArgs.push_back(context.getPeriodicBoxVecYPointer());
+        kernelArgs.push_back(context.getPeriodicBoxVecZPointer());
         for (int i = 0; i < (int) atomIndices.size(); i++)
             for (int j = 0; j < (int) atomIndices[i].size(); j++)
                 kernelArgs.push_back(&atomIndices[i][j]->getDevicePointer());

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-5);
+        ASSERT_EQUAL_TOL(delta.dot(delta), state.getPotentialEnergy(), 1e-5);
+        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/opencl/include/OpenCLExpressionUtilities.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);
     OpenCLContext& context;
-    FunctionPlaceholder fp1, fp2, fp3;
+    FunctionPlaceholder fp1, fp2, fp3, periodicDistance;
 };
 
 } // namespace OpenMM

platforms/opencl/src/OpenCLBondedUtilities.cpp

@@ -181,7 +181,7 @@ void OpenCLBondedUtilities::initialize(const System& system) {
         for (int i = 0; i < (int) prefixCode.size(); i++)
             s<<prefixCode[i];
         string bufferType = (context.getSupports64BitGlobalAtomics() ? "long" : "real4");
-        s<<"__kernel void computeBondedForces(__global "<<bufferType<<"* restrict forceBuffers, __global real* restrict energyBuffer, __global const real4* restrict posq, int groups";
+        s<<"__kernel void computeBondedForces(__global "<<bufferType<<"* restrict forceBuffers, __global real* restrict energyBuffer, __global const real4* restrict posq, int groups, real4 periodicBoxSize, real4 invPeriodicBoxSize, real4 periodicBoxVecX, real4 periodicBoxVecY, real4 periodicBoxVecZ";
         for (int i = 0; i < setSize; i++) {
             int force = set[i];
             string indexType = "uint"+(indexWidth[force] == 1 ? "" : context.intToString(indexWidth[force]));
@@ -267,7 +267,7 @@ void OpenCLBondedUtilities::computeInteractions(int groups) {
                 kernel.setArg<cl::Buffer>(index++, context.getForceBuffers().getDeviceBuffer());
             kernel.setArg<cl::Buffer>(index++, context.getEnergyBuffer().getDeviceBuffer());
             kernel.setArg<cl::Buffer>(index++, context.getPosq().getDeviceBuffer());
-            index++;
+            index += 6;
             for (int j = 0; j < (int) forceSets[i].size(); j++) {
                 kernel.setArg<cl::Buffer>(index++, atomIndices[forceSets[i][j]]->getDeviceBuffer());
                 kernel.setArg<cl::Buffer>(index++, bufferIndices[forceSets[i][j]]->getDeviceBuffer());
@@ -277,7 +277,22 @@ void OpenCLBondedUtilities::computeInteractions(int groups) {
         }
     }
     for (int i = 0; i < (int) kernels.size(); i++) {
-        kernels[i].setArg<cl_int>(3, groups);
+        cl::Kernel& kernel = kernels[i];
+        kernel.setArg<cl_int>(3, groups);
+        if (context.getUseDoublePrecision()) {
+            kernel.setArg<mm_double4>(4, context.getPeriodicBoxSizeDouble());
+            kernel.setArg<mm_double4>(5, context.getInvPeriodicBoxSizeDouble());
+            kernel.setArg<mm_double4>(6, context.getPeriodicBoxVecXDouble());
+            kernel.setArg<mm_double4>(7, context.getPeriodicBoxVecYDouble());
+            kernel.setArg<mm_double4>(8, context.getPeriodicBoxVecZDouble());
+        }
+        else {
+            kernel.setArg<mm_float4>(4, context.getPeriodicBoxSize());
+            kernel.setArg<mm_float4>(5, context.getInvPeriodicBoxSize());
+            kernel.setArg<mm_float4>(6, context.getPeriodicBoxVecX());
+            kernel.setArg<mm_float4>(7, context.getPeriodicBoxVecY());
+            kernel.setArg<mm_float4>(8, context.getPeriodicBoxVecZ());
+        }
         context.executeKernel(kernels[i], maxBonds);
     }
 }

platforms/opencl/src/OpenCLKernels.cpp

@@ -3821,7 +3821,9 @@ void OpenCLCalcCustomExternalForceKernel::initialize(const System& system, const
         globalParamNames[i] = force.getGlobalParameterName(i);
         globalParamValues[i] = (cl_float) force.getGlobalParameterDefaultValue(i);
     }
-    Lepton::ParsedExpression energyExpression = Lepton::Parser::parse(force.getEnergyFunction()).optimize();
+    map<string, Lepton::CustomFunction*> customFunctions;
+    customFunctions["periodicdistance"] = cl.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/opencl/tests/TestOpenCLCustomExternalForce.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:                                                              *
  *                                                                            *
@@ -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-5);
+        ASSERT_EQUAL_TOL(delta.dot(delta), state.getPotentialEnergy(), 1e-5);
+        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/include/ReferenceKernels.h

@@ -37,6 +37,7 @@
 #include "SimTKOpenMMRealType.h"
 #include "ReferenceNeighborList.h"
 #include "lepton/CompiledExpression.h"
+#include "lepton/CustomFunction.h"
 #include "lepton/ExpressionProgram.h"
 
 namespace OpenMM {
@@ -795,11 +796,23 @@ public:
      */
     void copyParametersToContext(ContextImpl& context, const CustomExternalForce& force);
 private:
+    class PeriodicDistanceFunction;
     int numParticles;
     std::vector<int> particles;
     RealOpenMM **particleParamArray;
     Lepton::CompiledExpression energyExpression, forceExpressionX, forceExpressionY, forceExpressionZ;
     std::vector<std::string> parameterNames, globalParameterNames;
+    RealVec* boxVectors;
+};
+
+class ReferenceCalcCustomExternalForceKernel::PeriodicDistanceFunction : public Lepton::CustomFunction {
+public:
+    RealVec** boxVectorHandle;
+    PeriodicDistanceFunction(RealVec** boxVectorHandle);
+    int getNumArguments() const;
+    double evaluate(const double* arguments) const;
+    double evaluateDerivative(const double* arguments, const int* derivOrder) const;
+    Lepton::CustomFunction* clone() const;
 };
 
 /**

platforms/reference/src/ReferenceKernels.cpp

@@ -1402,6 +1402,46 @@ void ReferenceCalcCustomGBForceKernel::copyParametersToContext(ContextImpl& cont
     }
 }
 
+ReferenceCalcCustomExternalForceKernel::PeriodicDistanceFunction::PeriodicDistanceFunction(RealVec** boxVectorHandle) : boxVectorHandle(boxVectorHandle) {
+}
+
+int ReferenceCalcCustomExternalForceKernel::PeriodicDistanceFunction::getNumArguments() const {
+    return 6;
+}
+
+double ReferenceCalcCustomExternalForceKernel::PeriodicDistanceFunction::evaluate(const double* arguments) const {
+    RealVec* boxVectors = *boxVectorHandle;
+    RealVec delta = RealVec(arguments[0], arguments[1], arguments[2])-RealVec(arguments[3], arguments[4], arguments[5]);
+    delta -= boxVectors[2]*floor(delta[2]/boxVectors[2][2]+0.5);
+    delta -= boxVectors[1]*floor(delta[1]/boxVectors[1][1]+0.5);
+    delta -= boxVectors[0]*floor(delta[0]/boxVectors[0][0]+0.5);
+    return sqrt(delta.dot(delta));
+}
+
+double ReferenceCalcCustomExternalForceKernel::PeriodicDistanceFunction::evaluateDerivative(const double* arguments, const int* derivOrder) const {
+    int argIndex = -1;
+    for (int i = 0; i < 6; i++) {
+        if (derivOrder[i] > 0) {
+            if (derivOrder[i] > 1 || argIndex != -1)
+                throw OpenMMException("Unsupported derivative of periodicdistance"); // Should be impossible for this to happen.
+            argIndex = i;
+        }
+    }
+    RealVec* boxVectors = *boxVectorHandle;
+    RealVec delta = RealVec(arguments[0], arguments[1], arguments[2])-RealVec(arguments[3], arguments[4], arguments[5]);
+    delta -= boxVectors[2]*floor(delta[2]/boxVectors[2][2]+0.5);
+    delta -= boxVectors[1]*floor(delta[1]/boxVectors[1][1]+0.5);
+    delta -= boxVectors[0]*floor(delta[0]/boxVectors[0][0]+0.5);
+    double r = sqrt(delta.dot(delta));
+    if (argIndex < 3)
+        return delta[argIndex]/r;
+    return -delta[argIndex-3]/r;
+}
+
+Lepton::CustomFunction* ReferenceCalcCustomExternalForceKernel::PeriodicDistanceFunction::clone() const {
+    return new PeriodicDistanceFunction(boxVectorHandle);
+}
+
 ReferenceCalcCustomExternalForceKernel::~ReferenceCalcCustomExternalForceKernel() {
     disposeRealArray(particleParamArray, numParticles);
 }
@@ -1423,7 +1463,10 @@ void ReferenceCalcCustomExternalForceKernel::initialize(const System& system, co
 
     // Parse the expression used to calculate the force.
 
-    Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction()).optimize();
+    map<string, Lepton::CustomFunction*> functions;
+    PeriodicDistanceFunction periodicDistance(&boxVectors);
+    functions["periodicdistance"] = &periodicDistance;
+    Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction(), functions).optimize();
     energyExpression = expression.createCompiledExpression();
     forceExpressionX = expression.differentiate("x").createCompiledExpression();
     forceExpressionY = expression.differentiate("y").createCompiledExpression();
@@ -1437,6 +1480,7 @@ void ReferenceCalcCustomExternalForceKernel::initialize(const System& system, co
 double ReferenceCalcCustomExternalForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
     vector<RealVec>& posData = extractPositions(context);
     vector<RealVec>& forceData = extractForces(context);
+    boxVectors = extractBoxVectors(context);
     RealOpenMM energy = 0;
     map<string, double> globalParameters;
     for (int i = 0; i < (int) globalParameterNames.size(); i++)

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:                                                              *
  *                                                                            *
@@ -101,9 +101,51 @@ void testForce() {
     }
 }
 
+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() {
     try {
         testForce();
+        testPeriodic();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;