Merge pull request #1467 from peastman/amoebaperiodic

Periodic boundary conditions for AMOEBA bonded forces

plugins/amoeba/platforms/cuda/src/kernels/amoebaStretchBendForce.cu

 // compute the value of the bond angle
 
-real xab = pos1.x - pos2.x;
-real yab = pos1.y - pos2.y;
-real zab = pos1.z - pos2.z;
+real3 ab = make_real3(pos1.x-pos2.x, pos1.y-pos2.y, pos1.z-pos2.z);
+real3 cb = make_real3(pos3.x-pos2.x, pos3.y-pos2.y, pos3.z-pos2.z);
 
-real xcb = pos3.x - pos2.x;
-real ycb = pos3.y - pos2.y;
-real zcb = pos3.z - pos2.z;
+#if APPLY_PERIODIC
+APPLY_PERIODIC_TO_DELTA(ab)
+APPLY_PERIODIC_TO_DELTA(cb)
+#endif
 
-real rab = SQRT(xab*xab + yab*yab + zab*zab);
-real rcb = SQRT(xcb*xcb + ycb*ycb + zcb*zcb);
+real rab = SQRT(ab.x*ab.x + ab.y*ab.y + ab.z*ab.z);
+real rcb = SQRT(cb.x*cb.x + cb.y*cb.y + cb.z*cb.z);
 
-real xp = ycb*zab - zcb*yab;
-real yp = zcb*xab - xcb*zab;
-real zp = xcb*yab - ycb*xab;
+real xp = cb.y*ab.z - cb.z*ab.y;
+real yp = cb.z*ab.x - cb.x*ab.z;
+real zp = cb.x*ab.y - cb.y*ab.x;
 
 real rp = SQRT(xp*xp + yp*yp + zp*zp);
 
-real dotp = xab*xcb + yab*ycb + zab*zcb;
+real dotp = ab.x*cb.x + ab.y*cb.y + ab.z*cb.z;
 real cosine = rab*rcb > 0 ? (dotp / (rab*rcb)) : (real) 1;
 cosine = (cosine > 1 ? (real) 1 : cosine);
 cosine = (cosine < -1 ? -(real) 1 : cosine);
 real angle;
 if (cosine > 0.99f || cosine < -0.99f) {
     // Highly unlikely a stretch-bend angle will be near 0 or 180, but just in case...
-    real3 cross_prod = cross(make_real3(xab, yab, zab), make_real3(xcb, ycb, zcb));
+    real3 cross_prod = cross(make_real3(ab.x, ab.y, ab.z), make_real3(cb.x, cb.y, cb.z));
     angle = ASIN(SQRT(dot(cross_prod, cross_prod))/(rab*rcb))*RAD_TO_DEG;
     if (cosine < 0.0f)
         angle = 180-angle;
@@ -41,13 +41,13 @@ real dt = angle - RAD_TO_DEG*parameters.z;
 real terma = rab*rp != 0 ? (-RAD_TO_DEG/(rab*rab*rp)) : (real) 0;
 real termc = rcb*rp != 0 ? (RAD_TO_DEG/(rcb*rcb*rp)) : (real) 0;
 
-real ddtdxia = terma * (yab*zp-zab*yp);
-real ddtdyia = terma * (zab*xp-xab*zp);
-real ddtdzia = terma * (xab*yp-yab*xp);
+real ddtdxia = terma * (ab.y*zp-ab.z*yp);
+real ddtdyia = terma * (ab.z*xp-ab.x*zp);
+real ddtdzia = terma * (ab.x*yp-ab.y*xp);
 
-real ddtdxic = termc * (ycb*zp-zcb*yp);
-real ddtdyic = termc * (zcb*xp-xcb*zp);
-real ddtdzic = termc * (xcb*yp-ycb*xp);
+real ddtdxic = termc * (cb.y*zp-cb.z*yp);
+real ddtdyic = termc * (cb.z*xp-cb.x*zp);
+real ddtdzic = termc * (cb.x*yp-cb.y*xp);
 
 // find chain rule terms for the bond length deviations
 
@@ -61,13 +61,13 @@ real frc2 = ((rp != 0) ? force_constants.y : (real) 0);
 
 real drkk = dr1*frc1 + dr2*frc2;
 
-real ddrdxia = terma * xab;
-real ddrdyia = terma * yab;
-real ddrdzia = terma * zab;
+real ddrdxia = terma * ab.x;
+real ddrdyia = terma * ab.y;
+real ddrdzia = terma * ab.z;
 
-real ddrdxic = termc * xcb;
-real ddrdyic = termc * ycb;
-real ddrdzic = termc * zcb;
+real ddrdxic = termc * cb.x;
+real ddrdyic = termc * cb.y;
+real ddrdzic = termc * cb.z;
 
 // get the energy and master chain rule terms for derivatives
 

plugins/amoeba/platforms/cuda/src/kernels/amoebaTorsionTorsionForce.cu

 int2 torsionParams = TORSION_PARAMS[index];
 
-real xba = pos2.x - pos1.x;
-real yba = pos2.y - pos1.y;
-real zba = pos2.z - pos1.z;
-
-real xcb = pos3.x - pos2.x;
-real ycb = pos3.y - pos2.y;
-real zcb = pos3.z - pos2.z;
-
-real xdc = pos4.x - pos3.x;
-real ydc = pos4.y - pos3.y;
-real zdc = pos4.z - pos3.z;
-
-real xed = pos5.x - pos4.x;
-real yed = pos5.y - pos4.y;
-real zed = pos5.z - pos4.z;
-
-real xt = yba*zcb - ycb*zba;
-real yt = zba*xcb - zcb*xba;
-real zt = xba*ycb - xcb*yba;
-
-real xu = ycb*zdc - ydc*zcb;
-real yu = zcb*xdc - zdc*xcb;
-real zu = xcb*ydc - xdc*ycb;
+real3 ba = make_real3(pos2.x-pos1.x, pos2.y-pos1.y, pos2.z-pos1.z);
+real3 cb = make_real3(pos3.x-pos2.x, pos3.y-pos2.y, pos3.z-pos2.z);
+real3 dc = make_real3(pos4.x-pos3.x, pos4.y-pos3.y, pos4.z-pos3.z);
+real3 ed = make_real3(pos5.x-pos4.x, pos5.y-pos4.y, pos5.z-pos4.z);
+
+#if APPLY_PERIODIC
+APPLY_PERIODIC_TO_DELTA(ba)
+APPLY_PERIODIC_TO_DELTA(cb)
+APPLY_PERIODIC_TO_DELTA(dc)
+APPLY_PERIODIC_TO_DELTA(ed)
+#endif
+
+real xt = ba.y*cb.z - cb.y*ba.z;
+real yt = ba.z*cb.x - cb.z*ba.x;
+real zt = ba.x*cb.y - cb.x*ba.y;
+
+real xu = cb.y*dc.z - dc.y*cb.z;
+real yu = cb.z*dc.x - dc.z*cb.x;
+real zu = cb.x*dc.y - dc.x*cb.y;
 
 real rt2 = xt*xt + yt*yt + zt*zt;
 real ru2 = xu*xu + yu*yu + zu*zu;
 
 real rtru = SQRT(rt2 * ru2);
 
-real xv = ydc*zed - yed*zdc;
-real yv = zdc*xed - zed*xdc;
-real zv = xdc*yed - xed*ydc;
+real xv = dc.y*ed.z - ed.y*dc.z;
+real yv = dc.z*ed.x - ed.z*dc.x;
+real zv = dc.x*ed.y - ed.x*dc.y;
 
 real rv2 = xv*xv + yv*yv + zv*zv;
 real rurv = SQRT(ru2 * rv2);
 
-real rcb = SQRT(xcb*xcb + ycb*ycb + zcb*zcb);
+real rcb = SQRT(cb.x*cb.x + cb.y*cb.y + cb.z*cb.z);
 real cosine1 = (rtru != 0 ? (xt*xu+yt*yu+zt*zu)/rtru : (real) 0);
 cosine1 = (cosine1 > 1 ? (real) 1 : cosine1);
 cosine1 = (cosine1 < -1 ? (real) -1 : cosine1);
@@ -51,11 +47,11 @@ if (cosine1 > 0.99f || cosine1 < -0.99f) {
 }
 else
    angle1 = RAD_TO_DEG*ACOS(cosine1);
-real sign = xba*xu + yba*yu + zba*zu;
+real sign = ba.x*xu + ba.y*yu + ba.z*zu;
 angle1 = (sign < 0 ? -angle1 : angle1);
 real value1 = angle1;
 
-real rdc = SQRT(xdc*xdc + ydc*ydc + zdc*zdc);
+real rdc = SQRT(dc.x*dc.x + dc.y*dc.y + dc.z*dc.z);
 real cosine2 = (xu*xv + yu*yv + zu*zv) / rurv;
 cosine2 = (cosine2 > 1 ? (real) 1 : cosine2);
 cosine2 = (cosine2 < -1 ? (real) -1 : cosine2);
@@ -70,7 +66,7 @@ if (cosine2 > 0.99f || cosine2 < -0.99f) {
 }
 else
    angle2 = RAD_TO_DEG*ACOS(cosine2);
-sign = xcb*xv + ycb*yv + zcb*zv;
+sign = cb.x*xv + cb.y*yv + cb.z*zv;
 angle2 = (sign < 0 ? -angle2 : angle2);
 real value2 = angle2;
 
@@ -83,24 +79,20 @@ real value2 = angle2;
 int chiralAtomIndex = (torsionParams.x > -1 ? torsionParams.x : atom5);
 real4 pos6 = posq[chiralAtomIndex];
 
-real xac = pos6.x - pos3.x;
-real yac = pos6.y - pos3.y;
-real zac = pos6.z - pos3.z;
-
-real xbc = pos2.x - pos3.x;
-real ybc = pos2.y - pos3.y;
-real zbc = pos2.z - pos3.z;
+real3 ac = make_real3(pos6.x-pos3.x, pos6.y-pos3.y, pos6.z-pos3.z);
+real3 bc = make_real3(pos2.x-pos3.x, pos2.y-pos3.y, pos2.z-pos3.z);
+real3 dc1 = make_real3(pos4.x-pos3.x, pos4.y-pos3.y, pos4.z-pos3.z);
 
-// xdc, ydc, zdc appear above
+#if APPLY_PERIODIC
+APPLY_PERIODIC_TO_DELTA(ac)
+APPLY_PERIODIC_TO_DELTA(bc)
+APPLY_PERIODIC_TO_DELTA(dc1)
+#endif
 
-real xdc1 = pos4.x - pos3.x;
-real ydc1 = pos4.y - pos3.y;
-real zdc1 = pos4.z - pos3.z;
-
-real c1 = ybc*zdc1 - zbc*ydc1;
-real c2 = ydc1*zac - zdc1*yac;
-real c3 = yac*zbc - zac*ybc;
-real vol = xac*c1 + xbc*c2 + xdc1*c3;
+real c1 = bc.y*dc1.z - bc.z*dc1.y;
+real c2 = dc1.y*ac.z - dc1.z*ac.y;
+real c3 = ac.y*bc.z - ac.z*bc.y;
+real vol = ac.x*c1 + bc.x*c2 + dc1.x*c3;
 sign = (vol > 0 ? (real) 1 : (real) -1);
 sign = (torsionParams.x < 0 ? (real) 1 : sign);
 value1 *= sign;
@@ -170,66 +162,68 @@ dedang2 *= sign * RAD_TO_DEG;
 
 // chain rule terms for first angle derivative components
 
-real xca = pos3.x - pos1.x;
-real yca = pos3.y - pos1.y;
-real zca = pos3.z - pos1.z;
+real3 ca = make_real3(pos3.x-pos1.x, pos3.y-pos1.y, pos3.z-pos1.z);
+real3 db = make_real3(pos4.x-pos2.x, pos4.y-pos2.y, pos4.z-pos2.z);
 
-real xdb = pos4.x - pos2.x;
-real ydb = pos4.y - pos2.y;
-real zdb = pos4.z - pos2.z;
+#if APPLY_PERIODIC
+APPLY_PERIODIC_TO_DELTA(ca)
+APPLY_PERIODIC_TO_DELTA(db)
+#endif
 
-real dedxt = dedang1 * (yt*zcb - ycb*zt) / (rt2*rcb);
-real dedyt = dedang1 * (zt*xcb - zcb*xt) / (rt2*rcb);
-real dedzt = dedang1 * (xt*ycb - xcb*yt) / (rt2*rcb);
-real dedxu = -dedang1 * (yu*zcb - ycb*zu) / (ru2*rcb);
-real dedyu = -dedang1 * (zu*xcb - zcb*xu) / (ru2*rcb);
-real dedzu = -dedang1 * (xu*ycb - xcb*yu) / (ru2*rcb);
+real dedxt = dedang1 * (yt*cb.z - cb.y*zt) / (rt2*rcb);
+real dedyt = dedang1 * (zt*cb.x - cb.z*xt) / (rt2*rcb);
+real dedzt = dedang1 * (xt*cb.y - cb.x*yt) / (rt2*rcb);
+real dedxu = -dedang1 * (yu*cb.z - cb.y*zu) / (ru2*rcb);
+real dedyu = -dedang1 * (zu*cb.x - cb.z*xu) / (ru2*rcb);
+real dedzu = -dedang1 * (xu*cb.y - cb.x*yu) / (ru2*rcb);
 
 // compute first derivative components for first angle
 
-real dedxia = zcb*dedyt - ycb*dedzt;
-real dedyia = xcb*dedzt - zcb*dedxt;
-real dedzia = ycb*dedxt - xcb*dedyt;
+real dedxia = cb.z*dedyt - cb.y*dedzt;
+real dedyia = cb.x*dedzt - cb.z*dedxt;
+real dedzia = cb.y*dedxt - cb.x*dedyt;
 
-real dedxib = yca*dedzt - zca*dedyt + zdc*dedyu - ydc*dedzu;
-real dedyib = zca*dedxt - xca*dedzt + xdc*dedzu - zdc*dedxu;
-real dedzib = xca*dedyt - yca*dedxt + ydc*dedxu - xdc*dedyu;
+real dedxib = ca.y*dedzt - ca.z*dedyt + dc.z*dedyu - dc.y*dedzu;
+real dedyib = ca.z*dedxt - ca.x*dedzt + dc.x*dedzu - dc.z*dedxu;
+real dedzib = ca.x*dedyt - ca.y*dedxt + dc.y*dedxu - dc.x*dedyu;
 
-real dedxic = zba*dedyt - yba*dedzt + ydb*dedzu - zdb*dedyu;
-real dedyic = xba*dedzt - zba*dedxt + zdb*dedxu - xdb*dedzu;
-real dedzic = yba*dedxt - xba*dedyt + xdb*dedyu - ydb*dedxu;
+real dedxic = ba.z*dedyt - ba.y*dedzt + db.y*dedzu - db.z*dedyu;
+real dedyic = ba.x*dedzt - ba.z*dedxt + db.z*dedxu - db.x*dedzu;
+real dedzic = ba.y*dedxt - ba.x*dedyt + db.x*dedyu - db.y*dedxu;
 
-real dedxid = zcb*dedyu - ycb*dedzu;
-real dedyid = xcb*dedzu - zcb*dedxu;
-real dedzid = ycb*dedxu - xcb*dedyu;
+real dedxid = cb.z*dedyu - cb.y*dedzu;
+real dedyid = cb.x*dedzu - cb.z*dedxu;
+real dedzid = cb.y*dedxu - cb.x*dedyu;
 
 // chain rule terms for second angle derivative components
 
-real xec = pos5.x - pos3.x;
-real yec = pos5.y - pos3.y;
-real zec = pos5.z - pos3.z;
+real3 ec = make_real3(pos5.x-pos3.x, pos5.y-pos3.y, pos5.z-pos3.z);
+
+#if APPLY_PERIODIC
+APPLY_PERIODIC_TO_DELTA(ec)
+#endif
 
-real dedxu2 = dedang2 * (yu*zdc - ydc*zu) / (ru2*rdc);
-real dedyu2 = dedang2 * (zu*xdc - zdc*xu) / (ru2*rdc);
-real dedzu2 = dedang2 * (xu*ydc - xdc*yu) / (ru2*rdc);
-real dedxv2 = -dedang2 * (yv*zdc - ydc*zv) / (rv2*rdc);
-real dedyv2 = -dedang2 * (zv*xdc - zdc*xv) / (rv2*rdc);
-real dedzv2 = -dedang2 * (xv*ydc - xdc*yv) / (rv2*rdc);
+real dedxu2 = dedang2 * (yu*dc.z - dc.y*zu) / (ru2*rdc);
+real dedyu2 = dedang2 * (zu*dc.x - dc.z*xu) / (ru2*rdc);
+real dedzu2 = dedang2 * (xu*dc.y - dc.x*yu) / (ru2*rdc);
+real dedxv2 = -dedang2 * (yv*dc.z - dc.y*zv) / (rv2*rdc);
+real dedyv2 = -dedang2 * (zv*dc.x - dc.z*xv) / (rv2*rdc);
+real dedzv2 = -dedang2 * (xv*dc.y - dc.x*yv) / (rv2*rdc);
 
 // compute first derivative components for second angle
 
-real dedxib2 = zdc*dedyu2 - ydc*dedzu2;
-real dedyib2 = xdc*dedzu2 - zdc*dedxu2;
-real dedzib2 = ydc*dedxu2 - xdc*dedyu2;
-real dedxic2 = ydb*dedzu2 - zdb*dedyu2 + zed*dedyv2 - yed*dedzv2;
-real dedyic2 = zdb*dedxu2 - xdb*dedzu2 + xed*dedzv2 - zed*dedxv2;
-real dedzic2 = xdb*dedyu2 - ydb*dedxu2 + yed*dedxv2 - xed*dedyv2;
-real dedxid2 = zcb*dedyu2 - ycb*dedzu2 + yec*dedzv2 - zec*dedyv2;
-real dedyid2 = xcb*dedzu2 - zcb*dedxu2 + zec*dedxv2 - xec*dedzv2;
-real dedzid2 = ycb*dedxu2 - xcb*dedyu2 + xec*dedyv2 - yec*dedxv2;
-real dedxie2 = zdc*dedyv2 - ydc*dedzv2;
-real dedyie2 = xdc*dedzv2 - zdc*dedxv2;
-real dedzie2 = ydc*dedxv2 - xdc*dedyv2;
+real dedxib2 = dc.z*dedyu2 - dc.y*dedzu2;
+real dedyib2 = dc.x*dedzu2 - dc.z*dedxu2;
+real dedzib2 = dc.y*dedxu2 - dc.x*dedyu2;
+real dedxic2 = db.y*dedzu2 - db.z*dedyu2 + ed.z*dedyv2 - ed.y*dedzv2;
+real dedyic2 = db.z*dedxu2 - db.x*dedzu2 + ed.x*dedzv2 - ed.z*dedxv2;
+real dedzic2 = db.x*dedyu2 - db.y*dedxu2 + ed.y*dedxv2 - ed.x*dedyv2;
+real dedxid2 = cb.z*dedyu2 - cb.y*dedzu2 + ec.y*dedzv2 - ec.z*dedyv2;
+real dedyid2 = cb.x*dedzu2 - cb.z*dedxu2 + ec.z*dedxv2 - ec.x*dedzv2;
+real dedzid2 = cb.y*dedxu2 - cb.x*dedyu2 + ec.x*dedyv2 - ec.y*dedxv2;
+real dedxie2 = dc.z*dedyv2 - dc.y*dedzv2;
+real dedyie2 = dc.x*dedzv2 - dc.z*dedxv2;
+real dedzie2 = dc.y*dedxv2 - dc.x*dedyv2;
 
 real3 force1 = make_real3(-dedxia, -dedyia, -dedzia);
 real3 force2 = make_real3(-dedxib-dedxib2, -dedyib-dedyib2, -dedzib-dedzib2);

plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaAngleForce.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 Stanford University and the Authors.           *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
  * Authors: Mark Friedrichs                                                   *
  * Contributors:                                                              *
  *                                                                            *
@@ -35,6 +35,7 @@
 
 #include "openmm/internal/AssertionUtilities.h"
 #include "openmm/Context.h"
+#include "openmm/CustomAngleForce.h"
 #include "OpenMMAmoeba.h"
 #include "openmm/System.h"
 #include "openmm/LangevinIntegrator.h"
@@ -273,6 +274,56 @@ void testOneAngle() {
     compareWithExpectedForceAndEnergy(context, *amoebaAngleForce, TOL, "testOneAngle");
 }
 
+void testPeriodic() {
+    // Create a force that uses periodic boundary conditions, then compare to an identical custom force.
+    
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
+    int numParticles = 3;
+    for (int ii = 0; ii < numParticles; ii++)
+        system.addParticle(1.0);
+    LangevinIntegrator integrator(0.0, 0.1, 0.01);
+    AmoebaAngleForce* amoebaAngleForce = new AmoebaAngleForce();
+    double angle      = 100.0;
+    double quadraticK = 1.0;
+    double cubicK     = 1.0e-01;
+    double quarticK   = 1.0e-02;
+    double penticK    = 1.0e-03;
+    double sexticK    = 1.0e-04;
+    amoebaAngleForce->addAngle(0, 1, 2, angle, quadraticK);
+    amoebaAngleForce->setAmoebaGlobalAngleCubic(cubicK);
+    amoebaAngleForce->setAmoebaGlobalAngleQuartic(quarticK);
+    amoebaAngleForce->setAmoebaGlobalAnglePentic(penticK);
+    amoebaAngleForce->setAmoebaGlobalAngleSextic(sexticK);
+    amoebaAngleForce->setUsesPeriodicBoundaryConditions(true);
+    system.addForce(amoebaAngleForce);
+    CustomAngleForce* customForce = new CustomAngleForce("k2*delta^2 + k3*delta^3 + k4*delta^4 + k5*delta^5 + k6*delta^6; delta=theta-theta0");
+    customForce->addGlobalParameter("theta0", angle*M_PI/180);
+    customForce->addGlobalParameter("k2", quadraticK*pow(180/M_PI, 2.0));
+    customForce->addGlobalParameter("k3", cubicK*pow(180/M_PI, 3.0));
+    customForce->addGlobalParameter("k4", quarticK*pow(180/M_PI, 4.0));
+    customForce->addGlobalParameter("k5", penticK*pow(180/M_PI, 5.0));
+    customForce->addGlobalParameter("k6", sexticK*pow(180/M_PI, 6.0));
+    customForce->addAngle(0, 1, 2);
+    customForce->setUsesPeriodicBoundaryConditions(true);
+    customForce->setForceGroup(1);
+    system.addForce(customForce);
+    Context context(system, integrator, Platform::getPlatformByName("CUDA"));
+
+    std::vector<Vec3> positions(numParticles);
+
+    positions[0] = Vec3(0, 1, 0);
+    positions[1] = Vec3(0, 0, 0);
+    positions[2] = Vec3(0, 0, 2);
+
+    context.setPositions(positions);
+    State s1 = context.getState(State::Forces | State::Energy, true, 1);
+    State s2 = context.getState(State::Forces | State::Energy, true, 2);
+    ASSERT_EQUAL_TOL(s2.getPotentialEnergy(), s1.getPotentialEnergy(), 1e-5);
+    for (int i = 0; i < numParticles; i++)
+        ASSERT_EQUAL_VEC(s2.getForces()[i], s1.getForces()[i], 1e-5);
+}
+
 int main(int argc, char* argv[]) {
     try {
         std::cout << "TestCudaAmoebaAngleForce running test..." << std::endl;
@@ -280,6 +331,7 @@ int main(int argc, char* argv[]) {
         if (argc > 1)
             Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1]));
         testOneAngle();
+        testPeriodic();
 
     } catch(const std::exception& e) {
         std::cout << "exception: " << e.what() << std::endl;

plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaBondForce.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 Stanford University and the Authors.           *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.           *
  * Authors: Mark Friedrichs                                                   *
  * Contributors:                                                              *
  *                                                                            *
@@ -36,6 +36,7 @@
 #include "openmm/internal/AssertionUtilities.h"
 #include "CudaPlatform.h"
 #include "openmm/Context.h"
+#include "openmm/CustomBondForce.h"
 #include "OpenMMAmoeba.h"
 #include "openmm/System.h"
 #include "openmm/LangevinIntegrator.h"
@@ -199,6 +200,49 @@ void testTwoBond() {
     compareWithExpectedForceAndEnergy(context, *amoebaBondForce, TOL, "testTwoBond");
 }
 
+void testPeriodic() {
+    // Create a force that uses periodic boundary conditions, then compare to an identical custom force.
+    
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
+    int numParticles = 2;
+    for (int ii = 0; ii < numParticles; ii++)
+        system.addParticle(1.0);
+    LangevinIntegrator integrator(0.0, 0.1, 0.01);
+    AmoebaBondForce* amoebaBondForce = new AmoebaBondForce();
+    double bondLength = 1.5;
+    double quadraticK = 1.0;
+    double cubicK     = 2.0;
+    double quarticK   = 3.0;
+    amoebaBondForce->setAmoebaGlobalBondCubic(cubicK);
+    amoebaBondForce->setAmoebaGlobalBondQuartic(quarticK);
+    amoebaBondForce->addBond(0, 1, bondLength, quadraticK);
+    amoebaBondForce->setUsesPeriodicBoundaryConditions(true);
+    system.addForce(amoebaBondForce);
+    CustomBondForce* customForce = new CustomBondForce("k2*delta^2 + k3*delta^3 + k4*delta^4; delta=r-r0");
+    customForce->addGlobalParameter("r0", bondLength);
+    customForce->addGlobalParameter("k2", quadraticK);
+    customForce->addGlobalParameter("k3", cubicK);
+    customForce->addGlobalParameter("k4", quarticK);
+    customForce->addBond(0, 1);
+    customForce->setUsesPeriodicBoundaryConditions(true);
+    customForce->setForceGroup(1);
+    system.addForce(customForce);
+    Context context(system, integrator, Platform::getPlatformByName("CUDA"));
+
+    std::vector<Vec3> positions(numParticles);
+
+    positions[0] = Vec3(0, 2, 0);
+    positions[1] = Vec3(0, 0, 0);
+
+    context.setPositions(positions);
+    State s1 = context.getState(State::Forces | State::Energy, true, 1);
+    State s2 = context.getState(State::Forces | State::Energy, true, 2);
+    ASSERT_EQUAL_TOL(s2.getPotentialEnergy(), s1.getPotentialEnergy(), 1e-5);
+    for (int i = 0; i < numParticles; i++)
+        ASSERT_EQUAL_VEC(s2.getForces()[i], s1.getForces()[i], 1e-5);
+}
+
 int main(int argc, char* argv[]) {
     try {
         std::cout << "TestCudaAmoebaBondForce running test..." << std::endl;
@@ -206,6 +250,7 @@ int main(int argc, char* argv[]) {
         if (argc > 1)
             Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1]));
         testTwoBond();
+        testPeriodic();
     } catch(const std::exception& e) {
         std::cout << "exception: " << e.what() << std::endl;
         std::cout << "FAIL - ERROR.  Test failed." << std::endl;

plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaInPlaneAngleForce.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 Stanford University and the Authors.           *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
  * Authors: Mark Friedrichs                                                   *
  * Contributors:                                                              *
  *                                                                            *
@@ -346,6 +346,50 @@ void testOneAngle() {
     compareWithExpectedForceAndEnergy(context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle");
 }
 
+void testPeriodic() {
+    // Create a force that uses periodic boundary conditions.
+
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
+    int numberOfParticles = 4;
+    for (int ii = 0; ii < numberOfParticles; ii++)
+        system.addParticle(1.0);
+    LangevinIntegrator integrator(0.0, 0.1, 0.01);
+    AmoebaInPlaneAngleForce* amoebaInPlaneAngleForce = new AmoebaInPlaneAngleForce();
+    double angle      = 65.0;
+    double quadraticK = 1.0;
+    double cubicK     = 0.0e-01;
+    double quarticK   = 0.0e-02;
+    double penticK    = 0.0e-03;
+    double sexticK    = 0.0e-04;
+    amoebaInPlaneAngleForce->addAngle(0, 1, 2, 3, angle, quadraticK);
+    amoebaInPlaneAngleForce->setAmoebaGlobalInPlaneAngleCubic(cubicK);
+    amoebaInPlaneAngleForce->setAmoebaGlobalInPlaneAngleQuartic(quarticK);
+    amoebaInPlaneAngleForce->setAmoebaGlobalInPlaneAnglePentic(penticK);
+    amoebaInPlaneAngleForce->setAmoebaGlobalInPlaneAngleSextic(sexticK);
+    amoebaInPlaneAngleForce->setUsesPeriodicBoundaryConditions(true);
+    system.addForce(amoebaInPlaneAngleForce);
+    Context context(system, integrator, Platform::getPlatformByName("CUDA"));
+
+    std::vector<Vec3> positions(numberOfParticles);
+    positions[0] = Vec3(0, 1, 0);
+    positions[1] = Vec3(0, 0, 0);
+    positions[2] = Vec3(0, 0, 1);
+    positions[3] = Vec3(1, 1, 1);
+
+    context.setPositions(positions);
+    State s1 = context.getState(State::Forces | State::Energy);
+    
+    // Move one atom to a position that should give identical results.
+
+    positions[2] = Vec3(0, 0, -2);
+    context.setPositions(positions);
+    State s2 = context.getState(State::Forces | State::Energy);
+    ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), 1e-5);
+    for (int i = 0; i < numberOfParticles; i++)
+        ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], 1e-5);
+}
+
 int main(int argc, char* argv[]) {
     try {
         std::cout << "TestCudaAmoebaInPlaneAngleForce running test..." << std::endl;
@@ -353,6 +397,7 @@ int main(int argc, char* argv[]) {
         if (argc > 1)
             Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1]));
         testOneAngle();
+        testPeriodic();
     } catch(const std::exception& e) {
         std::cout << "exception: " << e.what() << std::endl;
         std::cout << "FAIL - ERROR.  Test failed." << std::endl;

plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaOutOfPlaneBendForce.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 Stanford University and the Authors.           *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
  * Authors: Mark Friedrichs                                                   *
  * Contributors:                                                              *
  *                                                                            *
@@ -448,6 +448,43 @@ void testOneOutOfPlaneBend2(int setId) {
     }
 }
 
+void testPeriodic() {
+    // Create a force that uses periodic boundary conditions.
+
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
+    int numberOfParticles = 4;
+    for (int ii = 0; ii < numberOfParticles; ii++)
+        system.addParticle(1.0);
+    LangevinIntegrator integrator(0.0, 0.1, 0.01);
+    AmoebaOutOfPlaneBendForce* amoebaOutOfPlaneBendForce = new AmoebaOutOfPlaneBendForce();
+    amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendCubic( -0.1400000E-01);
+    amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendQuartic(0.5600000E-04);
+    amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendPentic(-0.7000000E-06);
+    amoebaOutOfPlaneBendForce->setAmoebaGlobalOutOfPlaneBendSextic( 0.2200000E-07);
+    double kOutOfPlaneBend = 0.328682196E-01;
+    amoebaOutOfPlaneBendForce->addOutOfPlaneBend(0, 1, 2, 3, kOutOfPlaneBend);
+    amoebaOutOfPlaneBendForce->setUsesPeriodicBoundaryConditions(true);
+    system.addForce(amoebaOutOfPlaneBendForce);
+    Context context(system, integrator, Platform::getPlatformByName("CUDA"));
+    std::vector<Vec3> positions(numberOfParticles);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    positions[2] = Vec3(0, 1, 0);
+    positions[3] = Vec3(0, 0, 1);
+    context.setPositions(positions);
+    State s1 = context.getState(State::Forces | State::Energy);
+    
+    // Move one atom to a position that should give identical results.
+
+    positions[3] = Vec3(0, 0, -2);
+    context.setPositions(positions);
+    State s2 = context.getState(State::Forces | State::Energy);
+    ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), 1e-5);
+    for (int i = 0; i < numberOfParticles; i++)
+        ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], 1e-5);
+}
+
 int main(int argc, char* argv[]) {
     try {
         std::cout << "TestCudaAmoebaOutOfPlaneBendForce running test..." << std::endl;
@@ -456,6 +493,7 @@ int main(int argc, char* argv[]) {
             Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1]));
 
         testOneOutOfPlaneBend();
+        testPeriodic();
 
     } catch(const std::exception& e) {
         std::cout << "exception: " << e.what() << std::endl;

plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaPiTorsionForce.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 Stanford University and the Authors.           *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
  * Authors: Mark Friedrichs                                                   *
  * Contributors:                                                              *
  *                                                                            *
@@ -288,6 +288,41 @@ void testOnePiTorsion() {
     compareWithExpectedForceAndEnergy(context, *amoebaPiTorsionForce, TOL, "testOnePiTorsion");
 }
 
+void testPeriodic() {
+    // Create a force that uses periodic boundary conditions.
+
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
+    int numberOfParticles = 6;
+    for (int ii = 0; ii < numberOfParticles; ii++)
+        system.addParticle(1.0);
+    LangevinIntegrator integrator(0.0, 0.1, 0.01);
+    AmoebaPiTorsionForce* amoebaPiTorsionForce = new AmoebaPiTorsionForce();
+    double kTorsion = 6.85;
+    amoebaPiTorsionForce->addPiTorsion(0, 1, 2, 3, 4, 5, kTorsion);
+    amoebaPiTorsionForce->setUsesPeriodicBoundaryConditions(true);
+    system.addForce(amoebaPiTorsionForce);
+    Context context(system, integrator, Platform::getPlatformByName("CUDA"));
+    std::vector<Vec3> positions(numberOfParticles);
+    positions[0] = Vec3(0, 1, 0);
+    positions[1] = Vec3(0, 0, 0);
+    positions[2] = Vec3(0, 0, 0.5);
+    positions[3] = Vec3(0.4, 0.4, 0.4);
+    positions[4] = Vec3(1, 0, 1);
+    positions[5] = Vec3(1, 1, 0);
+    context.setPositions(positions);
+    State s1 = context.getState(State::Forces | State::Energy);
+    
+    // Move one atom to a position that should give identical results.
+
+    positions[0] = Vec3(0, -2, 0);
+    context.setPositions(positions);
+    State s2 = context.getState(State::Forces | State::Energy);
+    ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), 1e-5);
+    for (int i = 0; i < numberOfParticles; i++)
+        ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], 1e-5);
+}
+
 int main(int argc, char* argv[]) {
     try {
         std::cout << "TestCudaAmoebaPiTorsionForce running test..." << std::endl;
@@ -295,6 +330,7 @@ int main(int argc, char* argv[]) {
         if (argc > 1)
             Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1]));
         testOnePiTorsion();
+        testPeriodic();
     } catch(const std::exception& e) {
         std::cout << "exception: " << e.what() << std::endl;
         std::cout << "FAIL - ERROR.  Test failed." << std::endl;

plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaStretchBendForce.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 Stanford University and the Authors.           *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
  * Authors: Mark Friedrichs                                                   *
  * Contributors:                                                              *
  *                                                                            *
@@ -266,6 +266,41 @@ void testOneStretchBend() {
     compareWithExpectedForceAndEnergy(context, *amoebaStretchBendForce, TOL, "testOneStretchBend");
 }
 
+void testPeriodic() {
+    // Create a force that uses periodic boundary conditions.
+
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
+    int numberOfParticles = 3;
+    for (int ii = 0; ii < numberOfParticles; ii++)
+        system.addParticle(1.0);
+    LangevinIntegrator integrator(0.0, 0.1, 0.01);
+    AmoebaStretchBendForce* amoebaStretchBendForce = new AmoebaStretchBendForce();
+    double abLength         = 0.144800000E+01;
+    double cbLength         = 0.101500000E+01;
+    double angleStretchBend = 0.108500000E+03*DegreesToRadians;
+    double kStretchBend     = 1.0;
+    amoebaStretchBendForce->addStretchBend(0, 1, 2, abLength, cbLength, angleStretchBend, kStretchBend, kStretchBend);
+    amoebaStretchBendForce->setUsesPeriodicBoundaryConditions(true);
+    system.addForce(amoebaStretchBendForce);
+    Context context(system, integrator, Platform::getPlatformByName("CUDA"));
+    std::vector<Vec3> positions(numberOfParticles);
+    positions[0] = Vec3(0, 1, 0);
+    positions[1] = Vec3(0, 0, 0);
+    positions[2] = Vec3(0, 0, 1);
+    context.setPositions(positions);
+    State s1 = context.getState(State::Forces | State::Energy);
+    
+    // Move one atom to a position that should give identical results.
+
+    positions[2] = Vec3(0, 0, -2);
+    context.setPositions(positions);
+    State s2 = context.getState(State::Forces | State::Energy);
+    ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), 1e-5);
+    for (int i = 0; i < numberOfParticles; i++)
+        ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], 1e-5);
+}
+
 int main(int argc, char* argv[]) {
     try {
         std::cout << "TestCudaAmoebaStretchBendForce running test..." << std::endl;
@@ -273,6 +308,7 @@ int main(int argc, char* argv[]) {
         if (argc > 1)
             Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1]));
         testOneStretchBend();
+        testPeriodic();
     } catch(const std::exception& e) {
         std::cout << "exception: " << e.what() << std::endl;
         std::cout << "FAIL - ERROR.  Test failed." << std::endl;

plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaTorsionTorsionForce.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 Stanford University and the Authors.           *
+ * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
  * Authors: Mark Friedrichs                                                   *
  * Contributors:                                                              *
  *                                                                            *
@@ -2667,6 +2667,44 @@ void testTorsionTorsion(int systemId) {
     ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), tolerance);
 }
 
+void testPeriodic() {
+    // Create a force that uses periodic boundary conditions.
+
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
+    int numberOfParticles = 6;
+    for (int ii = 0; ii < numberOfParticles; ii++)
+        system.addParticle(1.0);
+    LangevinIntegrator integrator(0.0, 0.1, 0.01);
+    AmoebaTorsionTorsionForce* amoebaTorsionTorsionForce = new AmoebaTorsionTorsionForce();
+    int chiralCheckAtomIndex;
+    int gridIndex;
+    chiralCheckAtomIndex  = 5;
+    gridIndex             = 2;
+    amoebaTorsionTorsionForce->addTorsionTorsion(0, 1, 2, 3, 4, chiralCheckAtomIndex, 0);
+    amoebaTorsionTorsionForce->setTorsionTorsionGrid(0, getTorsionGrid(gridIndex));
+    amoebaTorsionTorsionForce->setUsesPeriodicBoundaryConditions(true);
+    system.addForce(amoebaTorsionTorsionForce);
+    Context context(system, integrator, Platform::getPlatformByName("CUDA"));
+    std::vector<Vec3> positions(numberOfParticles);
+    positions[0] = Vec3(0, 1, 0);
+    positions[1] = Vec3(0, 0, 0);
+    positions[2] = Vec3(0, 0, 0.5);
+    positions[3] = Vec3(0.4, 0.4, 0.4);
+    positions[4] = Vec3(1, 0, 1);
+    positions[5] = Vec3(1, 1, 0);
+    context.setPositions(positions);
+    State s1 = context.getState(State::Forces | State::Energy);
+    
+    // Move one atom to a position that should give identical results.
+
+    positions[0] = Vec3(0, -2, 0);
+    context.setPositions(positions);
+    State s2 = context.getState(State::Forces | State::Energy);
+    ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), 1e-5);
+    for (int i = 0; i < numberOfParticles; i++)
+        ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], 1e-5);
+}
 
 int main(int argc, char* argv[]) {
     try {
@@ -2675,6 +2713,7 @@ int main(int argc, char* argv[]) {
         if (argc > 1)
             Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", std::string(argv[1]));
         testTorsionTorsion(1);
+        testPeriodic();
     } catch(const std::exception& e) {
         std::cout << "exception: " << e.what() << std::endl;
         std::cout << "FAIL - ERROR.  Test failed." << std::endl;

plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.cpp

@@ -104,12 +104,15 @@ void ReferenceCalcAmoebaBondForceKernel::initialize(const System& system, const
     } 
     globalBondCubic   = static_cast<RealOpenMM>(force.getAmoebaGlobalBondCubic());
     globalBondQuartic = static_cast<RealOpenMM>(force.getAmoebaGlobalBondQuartic());
+    usePeriodic = force.usesPeriodicBoundaryConditions();
 }
 
 double ReferenceCalcAmoebaBondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
     vector<RealVec>& posData   = extractPositions(context);
     vector<RealVec>& forceData = extractForces(context);
     AmoebaReferenceBondForce amoebaReferenceBondForce;
+    if (usePeriodic)
+        amoebaReferenceBondForce.setPeriodic(extractBoxVectors(context));
     RealOpenMM energy      = amoebaReferenceBondForce.calculateForceAndEnergy(numBonds, posData, particle1, particle2, length, kQuadratic,
                                                                                        globalBondCubic, globalBondQuartic,
                                                                                        forceData);
@@ -160,12 +163,15 @@ void ReferenceCalcAmoebaAngleForceKernel::initialize(const System& system, const
     globalAngleQuartic  = static_cast<RealOpenMM>(force.getAmoebaGlobalAngleQuartic());
     globalAnglePentic   = static_cast<RealOpenMM>(force.getAmoebaGlobalAnglePentic());
     globalAngleSextic   = static_cast<RealOpenMM>(force.getAmoebaGlobalAngleSextic());
+    usePeriodic = force.usesPeriodicBoundaryConditions();
 }
 
 double ReferenceCalcAmoebaAngleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
     vector<RealVec>& posData   = extractPositions(context);
     vector<RealVec>& forceData = extractForces(context);
     AmoebaReferenceAngleForce amoebaReferenceAngleForce;
+    if (usePeriodic)
+        amoebaReferenceAngleForce.setPeriodic(extractBoxVectors(context));
     RealOpenMM energy      = amoebaReferenceAngleForce.calculateForceAndEnergy(numAngles, 
                                        posData, particle1, particle2, particle3, angle, kQuadratic, globalAngleCubic, globalAngleQuartic, globalAnglePentic, globalAngleSextic, forceData);
     return static_cast<double>(energy);
@@ -213,6 +219,7 @@ void ReferenceCalcAmoebaInPlaneAngleForceKernel::initialize(const System& system
     globalInPlaneAngleQuartic  = static_cast<RealOpenMM>(force.getAmoebaGlobalInPlaneAngleQuartic());
     globalInPlaneAnglePentic   = static_cast<RealOpenMM>(force.getAmoebaGlobalInPlaneAnglePentic());
     globalInPlaneAngleSextic   = static_cast<RealOpenMM>(force.getAmoebaGlobalInPlaneAngleSextic());
+    usePeriodic = force.usesPeriodicBoundaryConditions();
 }
 
 double ReferenceCalcAmoebaInPlaneAngleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
@@ -220,6 +227,8 @@ double ReferenceCalcAmoebaInPlaneAngleForceKernel::execute(ContextImpl& context,
     vector<RealVec>& posData   = extractPositions(context);
     vector<RealVec>& forceData = extractForces(context);
     AmoebaReferenceInPlaneAngleForce amoebaReferenceInPlaneAngleForce;
+    if (usePeriodic)
+        amoebaReferenceInPlaneAngleForce.setPeriodic(extractBoxVectors(context));
     RealOpenMM energy      = amoebaReferenceInPlaneAngleForce.calculateForceAndEnergy(numAngles, posData, particle1, particle2, particle3, particle4, 
                                                                                                angle, kQuadratic, globalInPlaneAngleCubic, globalInPlaneAngleQuartic,
                                                                                                globalInPlaneAnglePentic, globalInPlaneAngleSextic, forceData);
@@ -266,12 +275,15 @@ void ReferenceCalcAmoebaPiTorsionForceKernel::initialize(const System& system, c
         particle6.push_back(particle6Index); 
         kTorsion.push_back(static_cast<RealOpenMM>(kTorsionParameter));
     }
+    usePeriodic = force.usesPeriodicBoundaryConditions();
 }
 
 double ReferenceCalcAmoebaPiTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
     vector<RealVec>& posData   = extractPositions(context);
     vector<RealVec>& forceData = extractForces(context);
     AmoebaReferencePiTorsionForce amoebaReferencePiTorsionForce;
+    if (usePeriodic)
+        amoebaReferencePiTorsionForce.setPeriodic(extractBoxVectors(context));
     RealOpenMM energy      = amoebaReferencePiTorsionForce.calculateForceAndEnergy(numPiTorsions, posData, particle1, particle2,
                                                                                     particle3, particle4, particle5, particle6,
                                                                                     kTorsion, forceData);
@@ -318,12 +330,15 @@ void ReferenceCalcAmoebaStretchBendForceKernel::initialize(const System& system,
         k1Parameters.push_back(static_cast<RealOpenMM>(k1));
         k2Parameters.push_back(static_cast<RealOpenMM>(k2));
     }
+    usePeriodic = force.usesPeriodicBoundaryConditions();
 }
 
 double ReferenceCalcAmoebaStretchBendForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
     vector<RealVec>& posData   = extractPositions(context);
     vector<RealVec>& forceData = extractForces(context);
     AmoebaReferenceStretchBendForce amoebaReferenceStretchBendForce;
+    if (usePeriodic)
+        amoebaReferenceStretchBendForce.setPeriodic(extractBoxVectors(context));
     RealOpenMM energy      = amoebaReferenceStretchBendForce.calculateForceAndEnergy(numStretchBends, posData, particle1, particle2, particle3,
                                                                                       lengthABParameters, lengthCBParameters, angleParameters, k1Parameters,
                                                                                       k2Parameters, forceData);
@@ -376,13 +391,15 @@ void ReferenceCalcAmoebaOutOfPlaneBendForceKernel::initialize(const System& syst
     globalOutOfPlaneBendAngleQuartic    = static_cast<RealOpenMM>(force.getAmoebaGlobalOutOfPlaneBendQuartic());
     globalOutOfPlaneBendAnglePentic     = static_cast<RealOpenMM>(force.getAmoebaGlobalOutOfPlaneBendPentic());
     globalOutOfPlaneBendAngleSextic     = static_cast<RealOpenMM>(force.getAmoebaGlobalOutOfPlaneBendSextic());
-
+    usePeriodic = force.usesPeriodicBoundaryConditions();
 }
 
 double ReferenceCalcAmoebaOutOfPlaneBendForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
     vector<RealVec>& posData   = extractPositions(context);
     vector<RealVec>& forceData = extractForces(context);
     AmoebaReferenceOutOfPlaneBendForce amoebaReferenceOutOfPlaneBendForce;
+    if (usePeriodic)
+        amoebaReferenceOutOfPlaneBendForce.setPeriodic(extractBoxVectors(context));
     RealOpenMM energy      = amoebaReferenceOutOfPlaneBendForce.calculateForceAndEnergy(numOutOfPlaneBends, posData,
                                                                                          particle1, particle2, particle3, particle4,
                                                                                          kParameters, 
@@ -434,6 +451,7 @@ void ReferenceCalcAmoebaTorsionTorsionForceKernel::initialize(const System& syst
         chiralCheckAtom.push_back(chiralCheckAtomIndex); 
         gridIndices.push_back(gridIndex); 
     }
+    usePeriodic = force.usesPeriodicBoundaryConditions();
 
     // torsion-torsion grids
 
@@ -478,6 +496,8 @@ double ReferenceCalcAmoebaTorsionTorsionForceKernel::execute(ContextImpl& contex
     vector<RealVec>& posData   = extractPositions(context);
     vector<RealVec>& forceData = extractForces(context);
     AmoebaReferenceTorsionTorsionForce amoebaReferenceTorsionTorsionForce;
+    if (usePeriodic)
+        amoebaReferenceTorsionTorsionForce.setPeriodic(extractBoxVectors(context));
     RealOpenMM energy      = amoebaReferenceTorsionTorsionForce.calculateForceAndEnergy(numTorsionTorsions, posData,
                                                                                          particle1, particle2, particle3, particle4, particle5,
                                                                                          chiralCheckAtom, gridIndices, torsionTorsionGrids, forceData);

plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.h

@@ -77,6 +77,7 @@ private:
     RealOpenMM globalBondCubic;
     RealOpenMM globalBondQuartic;
     const System& system;
+    bool usePeriodic;
 };
 
 /**
@@ -121,6 +122,7 @@ private:
     RealOpenMM globalAnglePentic;
     RealOpenMM globalAngleSextic;
     const System& system;
+    bool usePeriodic;
 };
 
 /**
@@ -166,6 +168,7 @@ private:
     RealOpenMM globalInPlaneAnglePentic;
     RealOpenMM globalInPlaneAngleSextic;
     const System& system;
+    bool usePeriodic;
 };
 
 /**
@@ -208,6 +211,7 @@ private:
     std::vector<int>   particle6;
     std::vector<RealOpenMM> kTorsion;
     const System& system;
+    bool usePeriodic;
 };
 
 /**
@@ -251,6 +255,7 @@ private:
     std::vector<RealOpenMM> k1Parameters;
     std::vector<RealOpenMM> k2Parameters;
     const System& system;
+    bool usePeriodic;
 };
 
 /**
@@ -295,6 +300,7 @@ private:
     RealOpenMM globalOutOfPlaneBendAnglePentic;
     RealOpenMM globalOutOfPlaneBendAngleSextic;
     const System& system;
+    bool usePeriodic;
 };
 
 /**
@@ -334,6 +340,7 @@ private:
     std::vector< std::vector< std::vector< std::vector<RealOpenMM> > > > torsionTorsionGrids;
 
     const System& system;
+    bool usePeriodic;
 };
 
 /**

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceAngleForce.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
@@ -28,6 +28,13 @@
 using std::vector;
 using namespace OpenMM;
 
+void AmoebaReferenceAngleForce::setPeriodic(OpenMM::RealVec* vectors) {
+    usePeriodic = true;
+    boxVectors[0] = vectors[0];
+    boxVectors[1] = vectors[1];
+    boxVectors[2] = vectors[2];
+}
+
 /**---------------------------------------------------------------------------------------
 
    Get dEdT and energy prefactor given cosine of angle :: the calculation for different
@@ -132,11 +139,17 @@ RealOpenMM AmoebaReferenceAngleForce::calculateAngleIxn(const RealVec& positionA
    // ---------------------------------------------------------------------------------------
 
     std::vector<RealOpenMM> deltaR[2];
-    AmoebaReferenceForce::loadDeltaR(positionAtomA, positionAtomB, deltaR[0]);
+    if (usePeriodic)
+        AmoebaReferenceForce::loadDeltaRPeriodic(positionAtomA, positionAtomB, deltaR[0], boxVectors);
+    else
+        AmoebaReferenceForce::loadDeltaR(positionAtomA, positionAtomB, deltaR[0]);
     RealOpenMM rAB2      = AmoebaReferenceForce::getNormSquared3(deltaR[0]);
     RealOpenMM rAB       = SQRT(rAB2);
  
-    AmoebaReferenceForce::loadDeltaR(positionAtomC, positionAtomB, deltaR[1]);
+    if (usePeriodic)
+        AmoebaReferenceForce::loadDeltaRPeriodic(positionAtomC, positionAtomB, deltaR[1], boxVectors);
+    else
+        AmoebaReferenceForce::loadDeltaR(positionAtomC, positionAtomB, deltaR[1]);
     RealOpenMM rCB2      = AmoebaReferenceForce::getNormSquared3(deltaR[1]);
     RealOpenMM rCB       = SQRT(rCB2);
  

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceAngleForce.h

 
-/* 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
@@ -40,7 +40,7 @@ public:
        
        --------------------------------------------------------------------------------------- */
  
-    AmoebaReferenceAngleForce() {};
+    AmoebaReferenceAngleForce() : usePeriodic(false) {};
  
     /**---------------------------------------------------------------------------------------
        
@@ -50,6 +50,16 @@ public:
  
     ~AmoebaReferenceAngleForce() {};
  
+    /**---------------------------------------------------------------------------------------
+
+       Set the force to use periodic boundary conditions.
+      
+       @param vectors    the vectors defining the periodic box
+      
+       --------------------------------------------------------------------------------------- */
+      
+    void setPeriodic(OpenMM::RealVec* vectors);
+
      /**---------------------------------------------------------------------------------------
      
         Calculate Amoeba angle ixns (force and energy)
@@ -86,6 +96,9 @@ public:
 
 private:
 
+    bool usePeriodic;
+    RealVec boxVectors[3];
+
     /**---------------------------------------------------------------------------------------
     
        Get dEdT and energy prefactor given cosine of angle :: the calculation for different

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceBondForce.cpp

@@ -28,6 +28,13 @@
 using std::vector;
 using namespace OpenMM;
 
+void AmoebaReferenceBondForce::setPeriodic(OpenMM::RealVec* vectors) {
+    usePeriodic = true;
+    boxVectors[0] = vectors[0];
+    boxVectors[1] = vectors[1];
+    boxVectors[2] = vectors[2];
+}
+
 /**---------------------------------------------------------------------------------------
 
    Calculate Amoeba bond ixn (force and energy)
@@ -63,7 +70,10 @@ RealOpenMM AmoebaReferenceBondForce::calculateBondIxn(const RealVec& positionAto
    // get deltaR, R2, and R between 2 atoms
 
    std::vector<RealOpenMM> deltaR;
-   AmoebaReferenceForce::loadDeltaR(positionAtomA, positionAtomB, deltaR);
+   if (usePeriodic)
+       AmoebaReferenceForce::loadDeltaRPeriodic(positionAtomA, positionAtomB, deltaR, boxVectors);
+   else
+       AmoebaReferenceForce::loadDeltaR(positionAtomA, positionAtomB, deltaR);
    RealOpenMM r               = AmoebaReferenceForce::getNorm3(deltaR);
 
    // deltaIdeal = r - r_0

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceBondForce.h

 
-/* 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
@@ -40,7 +40,7 @@ public:
        
         --------------------------------------------------------------------------------------- */
  
-    AmoebaReferenceBondForce() {};
+    AmoebaReferenceBondForce() : usePeriodic(false) {};
  
     /**---------------------------------------------------------------------------------------
        
@@ -49,7 +49,16 @@ public:
         --------------------------------------------------------------------------------------- */
  
     ~AmoebaReferenceBondForce() {};
- 
+
+    /**---------------------------------------------------------------------------------------
+
+       Set the force to use periodic boundary conditions.
+      
+       @param vectors    the vectors defining the periodic box
+      
+       --------------------------------------------------------------------------------------- */
+      
+    void setPeriodic(OpenMM::RealVec* vectors);
  
      /**---------------------------------------------------------------------------------------
      
@@ -79,6 +88,9 @@ public:
 
 private:
 
+    bool usePeriodic;
+    RealVec boxVectors[3];
+
      /**---------------------------------------------------------------------------------------
      
         Calculate Amoeba bond ixns (force and energy)

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceForce.cpp

@@ -39,18 +39,33 @@ using namespace OpenMM;
 
 void AmoebaReferenceForce::loadDeltaR(const RealVec& xVector, const RealVec& yVector,
                                       std::vector<RealOpenMM>& deltaR) {
+    deltaR.resize(0);
+    deltaR.push_back(yVector[0] - xVector[0]);
+    deltaR.push_back(yVector[1] - xVector[1]);
+    deltaR.push_back(yVector[2] - xVector[2]);
+}
 
+/**---------------------------------------------------------------------------------------
 
-   // ---------------------------------------------------------------------------------------
+   Load delta of two vectors, applying periodic boundary conditions
 
-   //static const std::string methodName = "AmoebaReferenceForce::loadDeltaR";
+   @param xVector      first vector
+   @param yVector      second vector
+   @param deltaR       output vector: y - x
+   @param boxVectors   periodic box vectors
 
-   // ---------------------------------------------------------------------------------------
+   --------------------------------------------------------------------------------------- */
 
+void AmoebaReferenceForce::loadDeltaRPeriodic(const RealVec& xVector, const RealVec& yVector,
+                                      std::vector<RealOpenMM>& deltaR, const RealVec* boxVectors) {
+    RealVec diff = yVector-xVector;
+    diff -= boxVectors[2]*floor(diff[2]/boxVectors[2][2]+0.5);
+    diff -= boxVectors[1]*floor(diff[1]/boxVectors[1][1]+0.5);
+    diff -= boxVectors[0]*floor(diff[0]/boxVectors[0][0]+0.5);
     deltaR.resize(0);
-    deltaR.push_back(yVector[0] - xVector[0]);
-    deltaR.push_back(yVector[1] - xVector[1]);
-    deltaR.push_back(yVector[2] - xVector[2]);
+    deltaR.push_back(diff[0]);
+    deltaR.push_back(diff[1]);
+    deltaR.push_back(diff[2]);
 }
 
 /**---------------------------------------------------------------------------------------
@@ -64,13 +79,6 @@ void AmoebaReferenceForce::loadDeltaR(const RealVec& xVector, const RealVec& yVe
    --------------------------------------------------------------------------------------- */
 
 RealOpenMM AmoebaReferenceForce::getNormSquared3(const std::vector<RealOpenMM>& inputVector) {
-
-   // ---------------------------------------------------------------------------------------
-
-   //static const std::string methodName = "AmoebaReferenceForce::getNorm3";
-
-   // ---------------------------------------------------------------------------------------
-
    // get 3 norm
 
    return (inputVector[0]*inputVector[0] + inputVector[1]*inputVector[1] + inputVector[2]*inputVector[2]);
@@ -87,13 +95,6 @@ RealOpenMM AmoebaReferenceForce::getNormSquared3(const std::vector<RealOpenMM>&
    --------------------------------------------------------------------------------------- */
 
 RealOpenMM AmoebaReferenceForce::getNormSquared3(const RealOpenMM* inputVector) {
-
-   // ---------------------------------------------------------------------------------------
-
-   //static const std::string methodName = "AmoebaReferenceForce::getNorm3";
-
-   // ---------------------------------------------------------------------------------------
-
    // get 3 norm
 
    return (inputVector[0]*inputVector[0] + inputVector[1]*inputVector[1] + inputVector[2]*inputVector[2]);
@@ -110,39 +111,18 @@ RealOpenMM AmoebaReferenceForce::getNormSquared3(const RealOpenMM* inputVector)
    --------------------------------------------------------------------------------------- */
 
 RealOpenMM AmoebaReferenceForce::getNorm3(const std::vector<RealOpenMM>& inputVector) {
-
-   // ---------------------------------------------------------------------------------------
-
-   //static const std::string methodName = "AmoebaReferenceForce::getNorm3";
-
-   // ---------------------------------------------------------------------------------------
-
    // get 3 norm
 
    return SQRT(inputVector[0]*inputVector[0] + inputVector[1]*inputVector[1] + inputVector[2]*inputVector[2]);
 }
 
 RealOpenMM AmoebaReferenceForce::getNorm3(const RealOpenMM* inputVector) {
-
-   // ---------------------------------------------------------------------------------------
-
-   //static const std::string methodName = "AmoebaReferenceForce::getNorm3";
-
-   // ---------------------------------------------------------------------------------------
-
    // get 3 norm
 
    return SQRT(inputVector[0]*inputVector[0] + inputVector[1]*inputVector[1] + inputVector[2]*inputVector[2]);
 }
 
 RealOpenMM AmoebaReferenceForce::normalizeVector3(RealOpenMM* inputVector) {
-
-   // ---------------------------------------------------------------------------------------
-
-   //static const std::string methodName = "AmoebaReferenceForce::normalizeVector3";
-
-   // ---------------------------------------------------------------------------------------
-
     RealOpenMM norm   = SQRT(inputVector[0]*inputVector[0] + inputVector[1]*inputVector[1] + inputVector[2]*inputVector[2]);
     if (norm > 0.0) {
         RealOpenMM normI  = 1.0/norm;
@@ -166,13 +146,6 @@ RealOpenMM AmoebaReferenceForce::normalizeVector3(RealOpenMM* inputVector) {
    --------------------------------------------------------------------------------------- */
 
 RealOpenMM AmoebaReferenceForce::getDotProduct3(const std::vector<RealOpenMM>& xVector, const std::vector<RealOpenMM>& yVector) {
-
-   // ---------------------------------------------------------------------------------------
-
-   //static const std::string methodName = "AmoebaReferenceForce::getDotProduct3";
-
-   // ---------------------------------------------------------------------------------------
-
    // get dot product
 
    return xVector[0]*yVector[0] + xVector[1]*yVector[1] + xVector[2]*yVector[2];
@@ -190,26 +163,12 @@ RealOpenMM AmoebaReferenceForce::getDotProduct3(const std::vector<RealOpenMM>& x
    --------------------------------------------------------------------------------------- */
 
 RealOpenMM AmoebaReferenceForce::getDotProduct3(const RealOpenMM* xVector, const RealOpenMM* yVector) {
-
-   // ---------------------------------------------------------------------------------------
-
-   //static const std::string methodName = "AmoebaReferenceForce::getDotProduct3";
-
-   // ---------------------------------------------------------------------------------------
-
    // get dot product
 
    return xVector[0]*yVector[0] + xVector[1]*yVector[1] + xVector[2]*yVector[2];
 }
 
 RealOpenMM AmoebaReferenceForce::getDotProduct3(const RealOpenMM* xVector, const OpenMM::Vec3& yVector) {
-
-   // ---------------------------------------------------------------------------------------
-
-   //static const std::string methodName = "AmoebaReferenceForce::getDotProduct3";
-
-   // ---------------------------------------------------------------------------------------
-
    // get dot product
 
    return xVector[0]*yVector[0] + xVector[1]*yVector[1] + xVector[2]*yVector[2];
@@ -228,13 +187,6 @@ RealOpenMM AmoebaReferenceForce::getDotProduct3(const RealOpenMM* xVector, const
    --------------------------------------------------------------------------------------- */
 
 RealOpenMM AmoebaReferenceForce::getDotProduct3(unsigned int vectorOffset, const std::vector<RealOpenMM>& xVector, const RealOpenMM* yVector) {
-
-   // ---------------------------------------------------------------------------------------
-
-   //static const std::string methodName = "AmoebaReferenceForce::getDotProduct3";
-
-   // ---------------------------------------------------------------------------------------
-
    // get dot product
 
    return xVector[vectorOffset+0]*yVector[0] + xVector[vectorOffset+1]*yVector[1] + xVector[vectorOffset+2]*yVector[2];
@@ -253,18 +205,9 @@ RealOpenMM AmoebaReferenceForce::getDotProduct3(unsigned int vectorOffset, const
 void AmoebaReferenceForce::getCrossProduct(const std::vector<RealOpenMM>& xVector,
                                            const std::vector<RealOpenMM>& yVector,
                                            std::vector<RealOpenMM>& zVector) {
-
-   // ---------------------------------------------------------------------------------------
-
-   //static const std::string methodName = "AmoebaReferenceForce::getCrossProduct";
-
-   // ---------------------------------------------------------------------------------------
-
    zVector[0]  = xVector[1]*yVector[2] - xVector[2]*yVector[1];
    zVector[1]  = xVector[2]*yVector[0] - xVector[0]*yVector[2];
    zVector[2]  = xVector[0]*yVector[1] - xVector[1]*yVector[0];
-
-   return;
 }
 
 /**---------------------------------------------------------------------------------------
@@ -280,17 +223,8 @@ void AmoebaReferenceForce::getCrossProduct(const std::vector<RealOpenMM>& xVecto
 void AmoebaReferenceForce::getCrossProduct(const RealOpenMM* xVector,
                                            const RealOpenMM* yVector,
                                            RealOpenMM* zVector) {
-
-   // ---------------------------------------------------------------------------------------
-
-   //static const std::string methodName = "AmoebaReferenceForce::getCrossProduct";
-
-   // ---------------------------------------------------------------------------------------
-
    zVector[0]  = xVector[1]*yVector[2] - xVector[2]*yVector[1];
    zVector[1]  = xVector[2]*yVector[0] - xVector[0]*yVector[2];
    zVector[2]  = xVector[0]*yVector[1] - xVector[1]*yVector[0];
-
-   return;
 }
 

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceForce.h

@@ -64,6 +64,18 @@ public:
     
     static void loadDeltaR(const OpenMM::RealVec& xVector, const OpenMM::RealVec& yVector,
                            std::vector<RealOpenMM>& deltaR);
+    /**---------------------------------------------------------------------------------------
+
+       Load delta of two vectors, applying periodic boundary conditions
+
+       @param xVector      first vector
+       @param yVector      second vector
+       @param deltaR       output vector: y - x
+       @param boxVectors   periodic box vectors
+
+       --------------------------------------------------------------------------------------- */
+
+    static void loadDeltaRPeriodic(const RealVec& xVector, const RealVec& yVector, std::vector<RealOpenMM>& deltaR, const RealVec* boxVectors);
     
     /**---------------------------------------------------------------------------------------
     

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceInPlaneAngleForce.cpp

@@ -24,10 +24,18 @@
 
 #include "AmoebaReferenceForce.h"
 #include "AmoebaReferenceInPlaneAngleForce.h"
+#include "ReferenceForce.h"
 
 using std::vector;
 using namespace OpenMM;
 
+void AmoebaReferenceInPlaneAngleForce::setPeriodic(OpenMM::RealVec* vectors) {
+    usePeriodic = true;
+    boxVectors[0] = vectors[0];
+    boxVectors[1] = vectors[1];
+    boxVectors[2] = vectors[2];
+}
+
 /**---------------------------------------------------------------------------------------
 
    Get dEdT and energy prefactor given cosine of angle :: the calculation for different
@@ -120,17 +128,6 @@ RealOpenMM AmoebaReferenceInPlaneAngleForce::calculateAngleIxn(const RealVec& po
                                                                         RealOpenMM angleCubic,                 RealOpenMM angleQuartic,
                                                                         RealOpenMM anglePentic,                RealOpenMM angleSextic,
                                                                         RealVec* forces) const {
-
-   // ---------------------------------------------------------------------------------------
-
-    //static const std::string methodName = "AmoebaReferenceInPlaneAngleForce::calculateAngleIxn";
- 
-    static const RealOpenMM zero          = 0.0;
-    static const RealOpenMM one           = 1.0;
-    static const RealOpenMM two           = 2.0;
-
-   // ---------------------------------------------------------------------------------------
-
     // T   = AD x CD
     // P   = B + T*delta
     // AP  = A - P
@@ -139,40 +136,48 @@ RealOpenMM AmoebaReferenceInPlaneAngleForce::calculateAngleIxn(const RealVec& po
 
     enum { AD, BD, CD, T, AP, P, CP, M, APxM, CPxM, ADxBD, BDxCD, TxCD, ADxT, dBxAD, CDxdB, LastDeltaAtomIndex };
 
-    std::vector<RealOpenMM> deltaR[LastDeltaAtomIndex];
-    for (int ii = 0; ii < LastDeltaAtomIndex; ii++) {
-        deltaR[ii].resize(3);
+    RealVec deltaR[LastDeltaAtomIndex];
+    if (usePeriodic) {
+        deltaR[AD] = ReferenceForce::getDeltaRPeriodic(positionAtomD, positionAtomA, boxVectors);
+        deltaR[BD] = ReferenceForce::getDeltaRPeriodic(positionAtomD, positionAtomB, boxVectors);
+        deltaR[CD] = ReferenceForce::getDeltaRPeriodic(positionAtomD, positionAtomC, boxVectors);
+    }
+    else {
+        deltaR[AD] = ReferenceForce::getDeltaR(positionAtomD, positionAtomA);
+        deltaR[BD] = ReferenceForce::getDeltaR(positionAtomD, positionAtomB);
+        deltaR[CD] = ReferenceForce::getDeltaR(positionAtomD, positionAtomC);
     }
-    AmoebaReferenceForce::loadDeltaR(positionAtomD, positionAtomA, deltaR[AD]);
-    AmoebaReferenceForce::loadDeltaR(positionAtomD, positionAtomB, deltaR[BD]);
-    AmoebaReferenceForce::loadDeltaR(positionAtomD, positionAtomC, deltaR[CD]);
  
-    AmoebaReferenceForce::getCrossProduct(deltaR[AD], deltaR[CD], deltaR[T]);
+    deltaR[T] = deltaR[AD].cross(deltaR[CD]);
 
-    RealOpenMM rT2     = AmoebaReferenceForce::getNormSquared3(deltaR[T]);
-    RealOpenMM delta   = AmoebaReferenceForce::getDotProduct3(deltaR[T], deltaR[BD])/rT2;
-         delta        *= -one;
- 
-    for (int ii = 0; ii < 3; ii++) {
-       deltaR[P][ii]  = positionAtomB[ii] + deltaR[T][ii]*delta;
-       deltaR[AP][ii] = positionAtomA[ii] - deltaR[P][ii];
-       deltaR[CP][ii] = positionAtomC[ii] - deltaR[P][ii];
+    RealOpenMM rT2     = deltaR[T].dot(deltaR[T]);
+    RealOpenMM delta   = deltaR[T].dot(deltaR[BD])/rT2;
+         delta        *= -1;
+
+    deltaR[P]  = positionAtomB + deltaR[T]*delta;
+    if (usePeriodic) {
+        deltaR[AP] = ReferenceForce::getDeltaRPeriodic(deltaR[P], positionAtomA, boxVectors);
+        deltaR[CP] = ReferenceForce::getDeltaRPeriodic(deltaR[P], positionAtomC, boxVectors);
+    }
+    else {
+        deltaR[AP] = ReferenceForce::getDeltaR(deltaR[P], positionAtomA);
+        deltaR[CP] = ReferenceForce::getDeltaR(deltaR[P], positionAtomC);
     }
  
-    RealOpenMM rAp2 = AmoebaReferenceForce::getNormSquared3(deltaR[AP]);
-    RealOpenMM rCp2 = AmoebaReferenceForce::getNormSquared3(deltaR[CP]);
-    if (rAp2 <= zero && rCp2 <= zero) {
-       return zero;
+    RealOpenMM rAp2 = deltaR[AP].dot(deltaR[AP]);
+    RealOpenMM rCp2 = deltaR[CP].dot(deltaR[CP]);
+    if (rAp2 <= 0 && rCp2 <= 0) {
+       return 0;
     }
  
-    AmoebaReferenceForce::getCrossProduct(deltaR[CP], deltaR[AP], deltaR[M]);
+    deltaR[M] = deltaR[CP].cross(deltaR[AP]);
  
-    RealOpenMM rm = AmoebaReferenceForce::getNorm3(deltaR[M]);
+    RealOpenMM rm = SQRT(deltaR[M].dot(deltaR[M]));
     if (rm < 1.0e-06) {
        rm = 1.0e-06;
     }
  
-    RealOpenMM dot     = AmoebaReferenceForce::getDotProduct3(deltaR[AP], deltaR[CP]);
+    RealOpenMM dot     = deltaR[AP].dot(deltaR[CP]);
     RealOpenMM cosine  = dot/SQRT(rAp2*rCp2);
  
     RealOpenMM dEdR;
@@ -182,64 +187,47 @@ RealOpenMM AmoebaReferenceInPlaneAngleForce::calculateAngleIxn(const RealVec& po
     RealOpenMM termA   = -dEdR/(rAp2*rm);
     RealOpenMM termC   =  dEdR/(rCp2*rm);
  
-    AmoebaReferenceForce::getCrossProduct(deltaR[AP], deltaR[M], deltaR[APxM]);
-    AmoebaReferenceForce::getCrossProduct(deltaR[CP], deltaR[M], deltaR[CPxM]);
+    deltaR[APxM] = deltaR[AP].cross(deltaR[M]);
+    deltaR[CPxM] = deltaR[CP].cross(deltaR[M]);
  
     // forces will be gathered here
  
     enum { dA, dB, dC, dD, LastDIndex };
-    std::vector<RealOpenMM> forceTerm[LastDIndex];
-    for (int ii = 0; ii < LastDIndex; ii++) {
-        forceTerm[ii].resize(3);
-    }
+    RealVec forceTerm[LastDIndex];
  
-    for (int ii = 0; ii < 3; ii++) {
-       forceTerm[dA][ii] = deltaR[APxM][ii]*termA;
-       forceTerm[dC][ii] = deltaR[CPxM][ii]*termC;
-       forceTerm[dB][ii] = -one*(forceTerm[dA][ii] + forceTerm[dC][ii]);
-    }
- 
-    RealOpenMM pTrT2  = AmoebaReferenceForce::getDotProduct3(forceTerm[dB], deltaR[T]);
+    forceTerm[dA] = deltaR[APxM]*termA;
+    forceTerm[dC] = deltaR[CPxM]*termC;
+    forceTerm[dB] = -(forceTerm[dA] + forceTerm[dC]);
+
+    RealOpenMM pTrT2  = forceTerm[dB].dot(deltaR[T]);
                pTrT2 /= rT2;
  
-    AmoebaReferenceForce::getCrossProduct(deltaR[CD], forceTerm[dB], deltaR[CDxdB]);
-    AmoebaReferenceForce::getCrossProduct(forceTerm[dB], deltaR[AD], deltaR[dBxAD]);
+    deltaR[CDxdB] = deltaR[CD].cross(forceTerm[dB]);
+    deltaR[dBxAD] = forceTerm[dB].cross(deltaR[AD]);
  
     if (FABS(pTrT2) > 1.0e-08) {
  
-       RealOpenMM delta2 = delta*two;
+        RealOpenMM delta2 = delta*2;
  
-       AmoebaReferenceForce::getCrossProduct(deltaR[BD], deltaR[CD], deltaR[BDxCD]);
-       AmoebaReferenceForce::getCrossProduct(deltaR[T],  deltaR[CD], deltaR[TxCD]);
-       AmoebaReferenceForce::getCrossProduct(deltaR[AD], deltaR[BD], deltaR[ADxBD]);
-       AmoebaReferenceForce::getCrossProduct(deltaR[AD], deltaR[T],  deltaR[ADxT]);
-       for (int ii = 0; ii < 3; ii++) {
-    
-          RealOpenMM term     = deltaR[BDxCD][ii] + delta2*deltaR[TxCD][ii];
-          forceTerm[dA][ii]  += delta*deltaR[CDxdB][ii] + term*pTrT2;
-    
-               term           = deltaR[ADxBD][ii] + delta2*deltaR[ADxT][ii];
-          forceTerm[dC][ii]  += delta*deltaR[dBxAD][ii] + term*pTrT2;
-    
-          forceTerm[dD][ii]  = -(forceTerm[dA][ii] + forceTerm[dB][ii] + forceTerm[dC][ii]);
-       }
+        deltaR[BDxCD] = forceTerm[dB].cross(deltaR[CD]);
+        deltaR[TxCD] = forceTerm[T].cross(deltaR[CD]);
+        deltaR[ADxBD] = forceTerm[AD].cross(deltaR[BD]);
+        deltaR[ADxT] = forceTerm[AD].cross(deltaR[T]);
+        RealVec term = deltaR[BDxCD] + deltaR[TxCD]*delta2;
+        forceTerm[dA] += deltaR[CDxdB]*delta + term*pTrT2;
+        term = deltaR[ADxBD] + deltaR[ADxT]*delta2;
+        forceTerm[dC] += deltaR[dBxAD]*delta + term*pTrT2;
+        forceTerm[dD] = -(forceTerm[dA] + forceTerm[dB] + forceTerm[dC]);
     } else {
-       for (int ii = 0; ii < 3; ii++) {
-    
-          forceTerm[dA][ii] += delta*deltaR[CDxdB][ii];
-          forceTerm[dC][ii] += delta*deltaR[dBxAD][ii];
- 
-          forceTerm[dD][ii]  = -(forceTerm[dA][ii] + forceTerm[dB][ii] + forceTerm[dC][ii]);
-       }
+        forceTerm[dA] += deltaR[CDxdB]*delta;
+        forceTerm[dC] += deltaR[dBxAD]*delta;
+        forceTerm[dD]  = -(forceTerm[dA] + forceTerm[dB] + forceTerm[dC]);
     }
  
     // accumulate forces
  
-    for (int jj = 0; jj < 4; jj++) {
-        forces[jj][0] = forceTerm[jj][0];
-        forces[jj][1] = forceTerm[jj][1];
-        forces[jj][2] = forceTerm[jj][2];
-    }
+    for (int jj = 0; jj < 4; jj++)
+        forces[jj] = forceTerm[jj];
  
     return energy;
 

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceInPlaneAngleForce.h

 
-/* 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
@@ -40,7 +40,7 @@ public:
        
        --------------------------------------------------------------------------------------- */
  
-    AmoebaReferenceInPlaneAngleForce() {};
+    AmoebaReferenceInPlaneAngleForce() : usePeriodic(false) {};
  
     /**---------------------------------------------------------------------------------------
        
@@ -50,6 +50,16 @@ public:
  
     ~AmoebaReferenceInPlaneAngleForce() {};
  
+    /**---------------------------------------------------------------------------------------
+
+       Set the force to use periodic boundary conditions.
+      
+       @param vectors    the vectors defining the periodic box
+      
+       --------------------------------------------------------------------------------------- */
+      
+    void setPeriodic(OpenMM::RealVec* vectors);
+
      /**---------------------------------------------------------------------------------------
      
         Calculate Amoeba in-plane angle ixns (force and energy)
@@ -88,6 +98,9 @@ public:
 
 private:
 
+    bool usePeriodic;
+    RealVec boxVectors[3];
+
     /**---------------------------------------------------------------------------------------
     
        Get dEdT and energy prefactor given cosine of angle :: the calculation for different

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceOutOfPlaneBendForce.cpp

@@ -28,6 +28,13 @@
 using std::vector;
 using namespace OpenMM;
 
+void AmoebaReferenceOutOfPlaneBendForce::setPeriodic(OpenMM::RealVec* vectors) {
+    usePeriodic = true;
+    boxVectors[0] = vectors[0];
+    boxVectors[1] = vectors[1];
+    boxVectors[2] = vectors[2];
+}
+
 /**---------------------------------------------------------------------------------------
 
    Calculate Amoeba Out-Of-Plane-Bend  ixn (force and energy)
@@ -78,12 +85,21 @@ RealOpenMM AmoebaReferenceOutOfPlaneBendForce::calculateOutOfPlaneBendIxn(const
     std::vector<RealOpenMM> deltaR[LastDeltaIndex];
     for (int ii = 0; ii < LastDeltaIndex; ii++) {
         deltaR[ii].resize(3);
-    }   
-    AmoebaReferenceForce::loadDeltaR(positionAtomB, positionAtomA, deltaR[AB]);
-    AmoebaReferenceForce::loadDeltaR(positionAtomB, positionAtomC, deltaR[CB]);
-    AmoebaReferenceForce::loadDeltaR(positionAtomB, positionAtomD, deltaR[DB]);
-    AmoebaReferenceForce::loadDeltaR(positionAtomD, positionAtomA, deltaR[AD]);
-    AmoebaReferenceForce::loadDeltaR(positionAtomD, positionAtomC, deltaR[CD]);
+    }
+    if (usePeriodic) {
+        AmoebaReferenceForce::loadDeltaRPeriodic(positionAtomB, positionAtomA, deltaR[AB], boxVectors);
+        AmoebaReferenceForce::loadDeltaRPeriodic(positionAtomB, positionAtomC, deltaR[CB], boxVectors);
+        AmoebaReferenceForce::loadDeltaRPeriodic(positionAtomB, positionAtomD, deltaR[DB], boxVectors);
+        AmoebaReferenceForce::loadDeltaRPeriodic(positionAtomD, positionAtomA, deltaR[AD], boxVectors);
+        AmoebaReferenceForce::loadDeltaRPeriodic(positionAtomD, positionAtomC, deltaR[CD], boxVectors);
+    }
+    else {
+        AmoebaReferenceForce::loadDeltaR(positionAtomB, positionAtomA, deltaR[AB]);
+        AmoebaReferenceForce::loadDeltaR(positionAtomB, positionAtomC, deltaR[CB]);
+        AmoebaReferenceForce::loadDeltaR(positionAtomB, positionAtomD, deltaR[DB]);
+        AmoebaReferenceForce::loadDeltaR(positionAtomD, positionAtomA, deltaR[AD]);
+        AmoebaReferenceForce::loadDeltaR(positionAtomD, positionAtomC, deltaR[CD]);
+    }
 
     RealOpenMM rDB2  = AmoebaReferenceForce::getNormSquared3(deltaR[DB]);
     RealOpenMM rAD2  = AmoebaReferenceForce::getNormSquared3(deltaR[AD]);