Reference implementation of periodic boundary conditions for bonded forces

platforms/reference/src/SimTKReference/ReferenceCustomTorsionIxn.cpp

-/* Portions copyright (c) 2010-2013 Stanford University and Simbios.
+/* Portions copyright (c) 2010-2016 Stanford University and Simbios.
  * Contributors: Peter Eastman
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -39,7 +39,7 @@ using namespace OpenMM;
 
 ReferenceCustomTorsionIxn::ReferenceCustomTorsionIxn(const Lepton::CompiledExpression& energyExpression,
         const Lepton::CompiledExpression& forceExpression, const vector<string>& parameterNames, map<string, double> globalParameters) :
-        energyExpression(energyExpression), forceExpression(forceExpression) {
+        energyExpression(energyExpression), forceExpression(forceExpression), usePeriodic(false) {
 
     energyTheta = ReferenceForce::getVariablePointer(this->energyExpression, "theta");
     forceTheta = ReferenceForce::getVariablePointer(this->forceExpression, "theta");
@@ -61,13 +61,13 @@ ReferenceCustomTorsionIxn::ReferenceCustomTorsionIxn(const Lepton::CompiledExpre
    --------------------------------------------------------------------------------------- */
 
 ReferenceCustomTorsionIxn::~ReferenceCustomTorsionIxn() {
+}
 
-   // ---------------------------------------------------------------------------------------
-
-   // static const char* methodName = "\nReferenceCustomTorsionIxn::~ReferenceCustomTorsionIxn";
-
-   // ---------------------------------------------------------------------------------------
-
+void ReferenceCustomTorsionIxn::setPeriodic(OpenMM::RealVec* vectors) {
+    usePeriodic = true;
+    boxVectors[0] = vectors[0];
+    boxVectors[1] = vectors[1];
+    boxVectors[2] = vectors[2];
 }
 
 /**---------------------------------------------------------------------------------------
@@ -107,9 +107,16 @@ void ReferenceCustomTorsionIxn::calculateBondIxn(int* atomIndices,
    int atomBIndex = atomIndices[1];
    int atomCIndex = atomIndices[2];
    int atomDIndex = atomIndices[3];
-   ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], deltaR[0]);
-   ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomCIndex], deltaR[1]);
-   ReferenceForce::getDeltaR(atomCoordinates[atomDIndex], atomCoordinates[atomCIndex], deltaR[2]);
+   if (usePeriodic) {
+      ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], boxVectors, deltaR[0]);  
+      ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomBIndex], atomCoordinates[atomCIndex], boxVectors, deltaR[1]);  
+      ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomDIndex], atomCoordinates[atomCIndex], boxVectors, deltaR[2]);  
+   }
+   else {
+      ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], deltaR[0]);  
+      ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomCIndex], deltaR[1]);  
+      ReferenceForce::getDeltaR(atomCoordinates[atomDIndex], atomCoordinates[atomCIndex], deltaR[2]);  
+   }
 
    // Visual Studio complains if crossProduct declared as 'crossProduct[2][3]'
 

platforms/reference/src/SimTKReference/ReferenceHarmonicBondIxn.cpp

 
-/* Portions copyright (c) 2006 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2016 Stanford University and Simbios.
  * Contributors: Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -38,14 +38,7 @@ using namespace OpenMM;
 
    --------------------------------------------------------------------------------------- */
 
-ReferenceHarmonicBondIxn::ReferenceHarmonicBondIxn() {
-
-   // ---------------------------------------------------------------------------------------
-
-   // static const char* methodName = "\nReferenceHarmonicBondIxn::ReferenceHarmonicBondIxn";
-
-   // ---------------------------------------------------------------------------------------
-
+ReferenceHarmonicBondIxn::ReferenceHarmonicBondIxn() : usePeriodic(false) {
 }
 
 /**---------------------------------------------------------------------------------------
@@ -55,13 +48,13 @@ ReferenceHarmonicBondIxn::ReferenceHarmonicBondIxn() {
    --------------------------------------------------------------------------------------- */
 
 ReferenceHarmonicBondIxn::~ReferenceHarmonicBondIxn() {
+}
 
-   // ---------------------------------------------------------------------------------------
-
-   // static const char* methodName = "\nReferenceHarmonicBondIxn::~ReferenceHarmonicBondIxn";
-
-   // ---------------------------------------------------------------------------------------
-
+void ReferenceHarmonicBondIxn::setPeriodic(OpenMM::RealVec* vectors) {
+    usePeriodic = true;
+    boxVectors[0] = vectors[0];
+    boxVectors[1] = vectors[1];
+    boxVectors[2] = vectors[2];
 }
 
 /**---------------------------------------------------------------------------------------
@@ -99,7 +92,10 @@ void ReferenceHarmonicBondIxn::calculateBondIxn(int* atomIndices,
 
    int atomAIndex = atomIndices[0];
    int atomBIndex = atomIndices[1];
-   ReferenceForce::getDeltaR(atomCoordinates[atomAIndex], atomCoordinates[atomBIndex], deltaR);  
+   if (usePeriodic)
+       ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomAIndex], atomCoordinates[atomBIndex], boxVectors, deltaR);  
+   else
+       ReferenceForce::getDeltaR(atomCoordinates[atomAIndex], atomCoordinates[atomBIndex], deltaR);  
 
    // deltaIdeal = r - r_0
 

platforms/reference/src/SimTKReference/ReferenceProperDihedralBond.cpp

 
-/* Portions copyright (c) 2006 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2016 Stanford University and Simbios.
  * Contributors: Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -38,14 +38,7 @@ using namespace OpenMM;
 
    --------------------------------------------------------------------------------------- */
 
-ReferenceProperDihedralBond::ReferenceProperDihedralBond() {
-
-   // ---------------------------------------------------------------------------------------
-
-   // static const char* methodName = "\nReferenceProperDihedralBond::ReferenceProperDihedralBond";
-
-   // ---------------------------------------------------------------------------------------
-
+ReferenceProperDihedralBond::ReferenceProperDihedralBond() : usePeriodic(false) {
 }
 
 /**---------------------------------------------------------------------------------------
@@ -55,13 +48,13 @@ ReferenceProperDihedralBond::ReferenceProperDihedralBond() {
    --------------------------------------------------------------------------------------- */
 
 ReferenceProperDihedralBond::~ReferenceProperDihedralBond() {
+}
 
-   // ---------------------------------------------------------------------------------------
-
-   // static const char* methodName = "\nReferenceProperDihedralBond::~ReferenceProperDihedralBond";
-
-   // ---------------------------------------------------------------------------------------
-
+void ReferenceProperDihedralBond::setPeriodic(OpenMM::RealVec* vectors) {
+    usePeriodic = true;
+    boxVectors[0] = vectors[0];
+    boxVectors[1] = vectors[1];
+    boxVectors[2] = vectors[2];
 }
 
 /**---------------------------------------------------------------------------------------
@@ -110,9 +103,16 @@ void ReferenceProperDihedralBond::calculateBondIxn(int* atomIndices,
    int atomBIndex = atomIndices[1];
    int atomCIndex = atomIndices[2];
    int atomDIndex = atomIndices[3];
-   ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], deltaR[0]);  
-   ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomCIndex], deltaR[1]);  
-   ReferenceForce::getDeltaR(atomCoordinates[atomDIndex], atomCoordinates[atomCIndex], deltaR[2]);  
+   if (usePeriodic) {
+      ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], boxVectors, deltaR[0]);  
+      ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomBIndex], atomCoordinates[atomCIndex], boxVectors, deltaR[1]);  
+      ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomDIndex], atomCoordinates[atomCIndex], boxVectors, deltaR[2]);  
+   }
+   else {
+      ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], deltaR[0]);  
+      ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomCIndex], deltaR[1]);  
+      ReferenceForce::getDeltaR(atomCoordinates[atomDIndex], atomCoordinates[atomCIndex], deltaR[2]);  
+   }
 
    RealOpenMM dotDihedral;
    RealOpenMM signOfAngle;

platforms/reference/src/SimTKReference/ReferenceRbDihedralBond.cpp

 
-/* Portions copyright (c) 2006 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2016 Stanford University and Simbios.
  * Contributors: Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -38,14 +38,7 @@ using namespace OpenMM;
 
    --------------------------------------------------------------------------------------- */
 
-ReferenceRbDihedralBond::ReferenceRbDihedralBond() {
-
-   // ---------------------------------------------------------------------------------------
-
-   // static const char* methodName = "\nReferenceRbDihedralBond::ReferenceRbDihedralBond";
-
-   // ---------------------------------------------------------------------------------------
-
+ReferenceRbDihedralBond::ReferenceRbDihedralBond() : usePeriodic(false) {
 }
 
 /**---------------------------------------------------------------------------------------
@@ -55,13 +48,13 @@ ReferenceRbDihedralBond::ReferenceRbDihedralBond() {
    --------------------------------------------------------------------------------------- */
 
 ReferenceRbDihedralBond::~ReferenceRbDihedralBond() {
+}
 
-   // ---------------------------------------------------------------------------------------
-
-   // static const char* methodName = "\nReferenceRbDihedralBond::~ReferenceRbDihedralBond";
-
-   // ---------------------------------------------------------------------------------------
-
+void ReferenceRbDihedralBond::setPeriodic(OpenMM::RealVec* vectors) {
+    usePeriodic = true;
+    boxVectors[0] = vectors[0];
+    boxVectors[1] = vectors[1];
+    boxVectors[2] = vectors[2];
 }
 
 /**---------------------------------------------------------------------------------------
@@ -112,9 +105,16 @@ void ReferenceRbDihedralBond::calculateBondIxn(int* atomIndices,
    int atomBIndex = atomIndices[1];
    int atomCIndex = atomIndices[2];
    int atomDIndex = atomIndices[3];
-   ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], deltaR[0]);  
-   ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomCIndex], deltaR[1]);  
-   ReferenceForce::getDeltaR(atomCoordinates[atomDIndex], atomCoordinates[atomCIndex], deltaR[2]);  
+   if (usePeriodic) {
+      ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], boxVectors, deltaR[0]);  
+      ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomBIndex], atomCoordinates[atomCIndex], boxVectors, deltaR[1]);  
+      ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomDIndex], atomCoordinates[atomCIndex], boxVectors, deltaR[2]);  
+   }
+   else {
+      ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], deltaR[0]);  
+      ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomCIndex], deltaR[1]);  
+      ReferenceForce::getDeltaR(atomCoordinates[atomDIndex], atomCoordinates[atomCIndex], deltaR[2]);  
+   }
 
    RealOpenMM cosPhi;
    RealOpenMM signOfAngle;

tests/TestCMAPTorsionForce.h

@@ -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) 2010-2015 Stanford University and the Authors.      *
+ * Portions copyright (c) 2010-2016 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -158,6 +158,61 @@ void testChangingParameters() {
     ASSERT_EQUAL_TOL(2-20*sin(M_PI/4), energy, 1e-5);
 }
 
+void testPeriodic() {
+    const int mapSize = 36;
+
+    // Create two systems that use periodic boundary conditions: one with a pair
+    // of periodic torsions, and one with a CMAP torsion that approximates the same force.
+
+    System system1;
+    for (int i = 0; i < 5; i++)
+        system1.addParticle(1.0);
+    PeriodicTorsionForce* periodic = new PeriodicTorsionForce();
+    periodic->addTorsion(0, 1, 2, 3, 2, M_PI/4, 1.5);
+    periodic->addTorsion(1, 2, 3, 4, 3, M_PI/3, 2.0);
+    periodic->setUsesPeriodicBoundaryConditions(true);
+    system1.addForce(periodic);
+    System system2;
+    for (int i = 0; i < 5; i++)
+        system2.addParticle(1.0);
+    CMAPTorsionForce* cmap = new CMAPTorsionForce();
+    vector<double> mapEnergy(mapSize*mapSize);
+    for (int i = 0; i < mapSize; i++) {
+        double angle1 = i*2*M_PI/mapSize;
+        double energy1 = 1.5*(1+cos(2*angle1-M_PI/4));
+        for (int j = 0; j < mapSize; j++) {
+            double angle2 = j*2*M_PI/mapSize;
+            double energy2 = 2.0*(1+cos(3*angle2-M_PI/3));
+            mapEnergy[i+j*mapSize] = energy1+energy2;
+        }
+    }
+    cmap->addMap(mapSize, mapEnergy);
+    cmap->addTorsion(0, 0, 1, 2, 3, 1, 2, 3, 4);
+    cmap->setUsesPeriodicBoundaryConditions(true);
+    system2.addForce(cmap);
+
+    // Set the atoms in various positions, and verify that both systems give equal forces and energy.
+
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    vector<Vec3> positions(5);
+    VerletIntegrator integrator1(0.01);
+    VerletIntegrator integrator2(0.01);
+    Context c1(system1, integrator1, platform);
+    Context c2(system2, integrator2, platform);
+    for (int i = 0; i < 50; i++) {
+        for (int j = 0; j < (int) positions.size(); j++)
+            positions[j] = Vec3(5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt));
+        c1.setPositions(positions);
+        c2.setPositions(positions);
+        State s1 = c1.getState(State::Forces | State::Energy);
+        State s2 = c2.getState(State::Forces | State::Energy);
+        for (int i = 0; i < system1.getNumParticles(); i++)
+            ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], 0.05);
+        ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), 1e-3);
+    }
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
@@ -165,6 +220,7 @@ int main(int argc, char* argv[]) {
         initializeTests(argc, argv);
         testCMAPTorsions();
         testChangingParameters();
+        testPeriodic();
         runPlatformTests();
     }
     catch(const exception& e) {

tests/TestCustomAngleForce.h

@@ -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-2015 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -29,6 +29,9 @@
  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
  * -------------------------------------------------------------------------- */
 
+#ifdef WIN32
+  #define _USE_MATH_DEFINES // Needed to get M_PI
+#endif
 #include "openmm/internal/AssertionUtilities.h"
 #include "openmm/Context.h"
 #include "openmm/CustomAngleForce.h"
@@ -145,6 +148,39 @@ void testIllegalVariable() {
     ASSERT(threwException);
 }
 
+void testPeriodic() {
+    // Create a force that uses periodic boundary conditions.
+    
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 1.5, 0), Vec3(0, 0, 3));
+    VerletIntegrator integrator(0.01);
+    CustomAngleForce* angles = new CustomAngleForce("0.5*k*(theta-theta0)^2");
+    angles->addPerAngleParameter("theta0");
+    angles->addPerAngleParameter("k");
+    vector<double> parameters(2);
+    parameters[0] = M_PI/3;
+    parameters[1] = 1.1;
+    angles->addAngle(0, 1, 2, parameters);
+    angles->setUsesPeriodicBoundaryConditions(true);
+    system.addForce(angles);
+    angles->setUsesPeriodicBoundaryConditions(true);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(3);
+    positions[0] = Vec3(0, 1, 0);
+    positions[1] = Vec3(0, 0, 0);
+    positions[2] = Vec3(1, 0, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    double torque = 1.1*M_PI/6;
+    ASSERT_EQUAL_VEC(Vec3(2*torque, 0, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, -torque, 0), forces[2], TOL);
+    ASSERT_EQUAL_TOL(0.5*1.1*(M_PI/6)*(M_PI/6), state.getPotentialEnergy(), TOL);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
@@ -152,6 +188,7 @@ int main(int argc, char* argv[]) {
         initializeTests(argc, argv);
         testAngles();
         testIllegalVariable();
+        testPeriodic();
         runPlatformTests();
     }
     catch(const exception& e) {

tests/TestCustomBondForce.h

@@ -148,6 +148,36 @@ void testIllegalVariable() {
     ASSERT(threwException);
 }
 
+void testPeriodic() {
+    // Create a force that uses periodic boundary conditions.
+    
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
+    VerletIntegrator integrator(0.01);
+    CustomBondForce* forceField = new CustomBondForce("scale*k*(r-r0)^2");
+    forceField->addPerBondParameter("r0");
+    forceField->addPerBondParameter("k");
+    forceField->addGlobalParameter("scale", 0.5);
+    vector<double> parameters(2);
+    parameters[0] = 1.9;
+    parameters[1] = 0.8;
+    forceField->addBond(0, 1, parameters);
+    forceField->setUsesPeriodicBoundaryConditions(true);
+    system.addForce(forceField);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(0, 2, 0);
+    positions[1] = Vec3(0, 0, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    ASSERT_EQUAL_VEC(Vec3(0, -0.8*0.9, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, 0.8*0.9, 0), forces[1], TOL);
+    ASSERT_EQUAL_TOL(0.5*0.8*0.9*0.9, state.getPotentialEnergy(), TOL);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
@@ -156,6 +186,7 @@ int main(int argc, char* argv[]) {
         testBonds();
         testManyParameters();
         testIllegalVariable();
+        testPeriodic();
         runPlatformTests();
     }
     catch(const exception& e) {

tests/TestCustomCentroidBondForce.h

@@ -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) 2015 Stanford University and the Authors.           *
+ * Portions copyright (c) 2015-2016 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -277,6 +277,62 @@ void testIllegalVariable() {
     ASSERT(threwException);
 }
 
+void testPeriodic() {
+    // Create a force that uses periodic boundary conditions.
+    
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(2.0);
+    system.addParticle(3.0);
+    system.addParticle(4.0);
+    system.addParticle(5.0);
+    system.setDefaultPeriodicBoxVectors(Vec3(2, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
+    CustomCentroidBondForce* force = new CustomCentroidBondForce(2, "k*distance(g1,g2)^2");
+    force->addPerBondParameter("k");
+    vector<int> particles1;
+    particles1.push_back(0);
+    particles1.push_back(1);
+    vector<int> particles2;
+    particles2.push_back(2);
+    particles2.push_back(3);
+    particles2.push_back(4);
+    force->addGroup(particles1);
+    force->addGroup(particles2);
+    vector<int> groups;
+    groups.push_back(0);
+    groups.push_back(1);
+    vector<double> parameters;
+    parameters.push_back(1.0);
+    force->addBond(groups, parameters);
+    force->setUsesPeriodicBoundaryConditions(true);
+    system.addForce(force);
+
+    // The center of mass of group 0 is (1.5, 0, 0).
+
+    vector<Vec3> positions(5);
+    positions[0] = Vec3(2.5, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+
+    // The center of mass of group 1 is (-1, 0, 0).
+
+    positions[2] = Vec3(-6, 0, 0);
+    positions[3] = Vec3(-1, 0, 0);
+    positions[4] = Vec3(2, 0, 0);
+
+    // Check the forces and energy.
+
+    VerletIntegrator integrator(0.01);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    ASSERT_EQUAL_TOL(0.5*0.5, state.getPotentialEnergy(), TOL);
+    ASSERT_EQUAL_VEC(Vec3(-2*0.5*(1.0/3.0), 0, 0), state.getForces()[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(-2*0.5*(2.0/3.0), 0, 0), state.getForces()[1], TOL);
+    ASSERT_EQUAL_VEC(Vec3(2*0.5*(3.0/12.0), 0, 0), state.getForces()[2], TOL);
+    ASSERT_EQUAL_VEC(Vec3(2*0.5*(4.0/12.0), 0, 0), state.getForces()[3], TOL);
+    ASSERT_EQUAL_VEC(Vec3(2*0.5*(5.0/12.0), 0, 0), state.getForces()[4], TOL);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
@@ -286,6 +342,7 @@ int main(int argc, char* argv[]) {
         testComplexFunction();
         testCustomWeights();
         testIllegalVariable();
+        testPeriodic();
         runPlatformTests();
     }
     catch(const exception& e) {

tests/TestCustomCompoundBondForce.h

@@ -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) 2012-2015 Stanford University and the Authors.      *
+ * Portions copyright (c) 2012-2016 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -330,6 +330,82 @@ void testIllegalVariable() {
     ASSERT(threwException);
 }
 
+void testPeriodic() {
+    // Create a force that uses periodic boundary conditions.
+
+    System customSystem;
+    customSystem.addParticle(1.0);
+    customSystem.addParticle(1.0);
+    customSystem.addParticle(1.0);
+    customSystem.addParticle(1.0);
+    customSystem.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
+    CustomCompoundBondForce* custom = new CustomCompoundBondForce(4, "0.5*kb*((distance(p1,p2)-b0)^2+(distance(p2,p3)-b0)^2)+0.5*ka*(angle(p2,p3,p4)-a0)^2+kt*(1+cos(dihedral(p1,p2,p3,p4)-t0))");
+    custom->addPerBondParameter("kb");
+    custom->addPerBondParameter("ka");
+    custom->addPerBondParameter("kt");
+    custom->addPerBondParameter("b0");
+    custom->addPerBondParameter("a0");
+    custom->addPerBondParameter("t0");
+    vector<int> particles(4);
+    particles[0] = 0;
+    particles[1] = 1;
+    particles[2] = 3;
+    particles[3] = 2;
+    vector<double> parameters(6);
+    parameters[0] = 1.5;
+    parameters[1] = 0.8;
+    parameters[2] = 0.6;
+    parameters[3] = 1.1;
+    parameters[4] = 2.9;
+    parameters[5] = 1.3;
+    custom->addBond(particles, parameters);
+    custom->setUsesPeriodicBoundaryConditions(true);
+    customSystem.addForce(custom);
+
+    // Create an identical system using standard forces.
+
+    System standardSystem;
+    standardSystem.addParticle(1.0);
+    standardSystem.addParticle(1.0);
+    standardSystem.addParticle(1.0);
+    standardSystem.addParticle(1.0);
+    standardSystem.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
+    HarmonicBondForce* bonds = new HarmonicBondForce();
+    bonds->addBond(0, 1, 1.1, 1.5);
+    bonds->addBond(1, 3, 1.1, 1.5);
+    bonds->setUsesPeriodicBoundaryConditions(true);
+    standardSystem.addForce(bonds);
+    HarmonicAngleForce* angles = new HarmonicAngleForce();
+    angles->addAngle(1, 3, 2, 2.9, 0.8);
+    angles->setUsesPeriodicBoundaryConditions(true);
+    standardSystem.addForce(angles);
+    PeriodicTorsionForce* torsions = new PeriodicTorsionForce();
+    torsions->addTorsion(0, 1, 3, 2, 1, 1.3, 0.6);
+    torsions->setUsesPeriodicBoundaryConditions(true);
+    standardSystem.addForce(torsions);
+
+    // Set the atoms in various positions, and verify that both systems give identical forces and energy.
+
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    VerletIntegrator integrator1(0.01);
+    VerletIntegrator integrator2(0.01);
+    Context c1(customSystem, integrator1, platform);
+    Context c2(standardSystem, integrator2, platform);
+    vector<Vec3> positions(4);
+    for (int i = 0; i < 10; i++) {
+        for (int j = 0; j < (int) positions.size(); j++)
+            positions[j] = Vec3(5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt));
+        c1.setPositions(positions);
+        c2.setPositions(positions);
+        State s1 = c1.getState(State::Forces | State::Energy);
+        State s2 = c2.getState(State::Forces | State::Energy);
+        for (int i = 0; i < customSystem.getNumParticles(); i++)
+            ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], TOL);
+        ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), TOL);
+    }
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
@@ -341,6 +417,7 @@ int main(int argc, char* argv[]) {
         testContinuous3DFunction();
         testMultipleBonds();
         testIllegalVariable();
+        testPeriodic();
         runPlatformTests();
     }
     catch(const exception& e) {

tests/TestCustomTorsionForce.h

@@ -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-2015 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -186,6 +186,42 @@ void testIllegalVariable() {
     ASSERT(threwException);
 }
 
+void testPeriodic() {
+    // Create a force that uses periodic boundary conditions.
+    
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
+    VerletIntegrator integrator(0.01);
+    CustomTorsionForce* torsions = new CustomTorsionForce("k*(1+cos(n*theta-theta0))");
+    torsions->addPerTorsionParameter("theta0");
+    torsions->addPerTorsionParameter("n");
+    torsions->addGlobalParameter("k", 1.1);
+    vector<double> parameters(2);
+    parameters[0] = M_PI/3;
+    parameters[1] = 2;
+    torsions->addTorsion(0, 1, 2, 3, parameters);
+    torsions->setUsesPeriodicBoundaryConditions(true);
+    system.addForce(torsions);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(4);
+    positions[0] = Vec3(0, 1, 0);
+    positions[1] = Vec3(0, 0, 0);
+    positions[2] = Vec3(1, 0, 0);
+    positions[3] = Vec3(1, 0, 2);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    double torque = -2*1.1*std::sin(2*M_PI/3);
+    ASSERT_EQUAL_VEC(Vec3(0, 0, torque), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, -torque, 0), forces[3], TOL);
+    ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0], forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1], forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]), Vec3(0, 0, 0), TOL);
+    ASSERT_EQUAL_TOL(1.1*(1+std::cos(2*M_PI/3)), state.getPotentialEnergy(), TOL);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
@@ -194,6 +230,7 @@ int main(int argc, char* argv[]) {
         testTorsions();
         testRange();
         testIllegalVariable();
+        testPeriodic();
         runPlatformTests();
     }
     catch(const exception& e) {

tests/TestHarmonicAngleForce.h

@@ -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-2015 Stanford University and the Authors.s      *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.s      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -93,12 +93,40 @@ void testAngles() {
     }
 }
 
+void testPeriodic() {
+    // Create a force that uses periodic boundary conditions.
+    
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 1.5, 0), Vec3(0, 0, 3));
+    VerletIntegrator integrator(0.01);
+    HarmonicAngleForce* angles = new HarmonicAngleForce();
+    angles->addAngle(0, 1, 2, PI_M/3, 1.1);
+    system.addForce(angles);
+    angles->setUsesPeriodicBoundaryConditions(true);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(3);
+    positions[0] = Vec3(0, 1, 0);
+    positions[1] = Vec3(0, 0, 0);
+    positions[2] = Vec3(1, 0, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    double torque = 1.1*PI_M/6;
+    ASSERT_EQUAL_VEC(Vec3(2*torque, 0, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, -torque, 0), forces[2], TOL);
+    ASSERT_EQUAL_TOL(0.5*1.1*(PI_M/6)*(PI_M/6), state.getPotentialEnergy(), TOL);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
     try {
         initializeTests(argc, argv);
         testAngles();
+        testPeriodic();
         runPlatformTests();
     }
     catch(const exception& e) {

tests/TestHarmonicBondForce.h

@@ -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-2015 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -82,12 +82,37 @@ void testBonds() {
     }
 }
 
+void testPeriodic() {
+    // Create a force that uses periodic boundary conditions.
+    
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
+    VerletIntegrator integrator(0.01);
+    HarmonicBondForce* bonds = new HarmonicBondForce();
+    bonds->addBond(0, 1, 1.2, 0.8);
+    bonds->setUsesPeriodicBoundaryConditions(true);
+    system.addForce(bonds);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(0, 2, 0);
+    positions[1] = Vec3(0, 0, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    ASSERT_EQUAL_VEC(Vec3(0, -0.8*0.2, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, 0.8*0.2, 0), forces[1], TOL);
+    ASSERT_EQUAL_TOL(0.5*0.8*0.2*0.2, state.getPotentialEnergy(), TOL);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
     try {
         initializeTests(argc, argv);
         testBonds();
+        testPeriodic();
         runPlatformTests();
     }
     catch(const exception& e) {

tests/TestPeriodicTorsionForce.h

@@ -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-2015 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -88,12 +88,43 @@ void testPeriodicTorsions() {
     }
 }
 
+void testPeriodic() {
+    // Create a force that uses periodic boundary conditions.
+    
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
+    VerletIntegrator integrator(0.01);
+    PeriodicTorsionForce* torsions = new PeriodicTorsionForce();
+    torsions->addTorsion(0, 1, 2, 3, 2, PI_M/3, 1.1);
+    torsions->setUsesPeriodicBoundaryConditions(true);
+    system.addForce(torsions);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(4);
+    positions[0] = Vec3(0, 1, 0);
+    positions[1] = Vec3(0, 0, 0);
+    positions[2] = Vec3(1, 0, 0);
+    positions[3] = Vec3(1, 0, 2);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    double torque = -2*1.1*std::sin(2*PI_M/3);
+    ASSERT_EQUAL_VEC(Vec3(0, 0, torque), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, -torque, 0), forces[3], TOL);
+    ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0], forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1], forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]), Vec3(0, 0, 0), TOL);
+    ASSERT_EQUAL_TOL(1.1*(1+std::cos(2*PI_M/3)), state.getPotentialEnergy(), TOL);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
     try {
         initializeTests(argc, argv);
         testPeriodicTorsions();
+        testPeriodic();
         runPlatformTests();
     }
     catch(const exception& e) {

tests/TestRBTorsionForce.h

@@ -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-2015 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -107,12 +107,53 @@ void testRBTorsions() {
     }
 }
 
+void testPeriodic() {
+    // Create a force that uses periodic boundary conditions.
+    
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
+    VerletIntegrator integrator(0.01);
+    RBTorsionForce* torsions = new RBTorsionForce();
+    torsions->addTorsion(0, 1, 2, 3, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6);
+    torsions->setUsesPeriodicBoundaryConditions(true);
+    system.addForce(torsions);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(4);
+    positions[0] = Vec3(0, 1, 0);
+    positions[1] = Vec3(0, 0, 0);
+    positions[2] = Vec3(1, 0, 0);
+    positions[3] = Vec3(0, 0, 2);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    double psi = 0.5*PI_M;
+    double torque = 0.0;
+    for (int i = 1; i < 6; ++i) {
+        double c = 0.1*(i+1);
+        torque += -c*i*std::pow(std::cos(psi), i-1)*std::sin(psi);
+    }
+    ASSERT_EQUAL_VEC(Vec3(0, 0, torque), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, -torque, 0), forces[3], TOL);
+    ASSERT_EQUAL_VEC(Vec3(forces[0][0]+forces[1][0]+forces[2][0]+forces[3][0], forces[0][1]+forces[1][1]+forces[2][1]+forces[3][1], forces[0][2]+forces[1][2]+forces[2][2]+forces[3][2]), Vec3(0, 0, 0), TOL);
+    double energy = 0.0;
+    for (int i = 0; i < 6; ++i) {
+        double c = 0.1*(i+1);
+        energy += c*std::pow(std::cos(psi), i);
+    }
+    ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
     try {
         initializeTests(argc, argv);
         testRBTorsions();
+        testPeriodic();
         runPlatformTests();
     }
     catch(const exception& e) {