Continuing to implement triclinic boxes in CPU platform

platforms/cpu/include/CpuCustomManyParticleForce.h

@@ -51,9 +51,11 @@ private:
     class ComputeForceTask;
     class ThreadData;
     int numParticles, numParticlesPerSet, numPerParticleParameters, numTypes;
-    bool useCutoff, usePeriodic, centralParticleMode;
+    bool useCutoff, usePeriodic, triclinic, centralParticleMode;
     RealOpenMM cutoffDistance;
-    RealOpenMM periodicBoxSize[3];
+    float recipBoxSize[3];
+    RealVec periodicBoxVectors[3];
+    AlignedArray<fvec4> periodicBoxVec4;
     CpuNeighborList* neighborList;
     ThreadPool& threads;
     std::vector<std::set<int> > exclusions;
@@ -138,9 +140,9 @@ public:
      * already been set, and the smallest side of the periodic box is at least twice the cutoff
      * distance.
      * 
-     * @param boxSize    the X, Y, and Z widths of the periodic box
+     * @param periodicBoxVectors    the vectors defining the periodic box
      */
-    void setPeriodic(OpenMM::RealVec& boxSize);
+    void setPeriodic(RealVec* periodicBoxVectors);
 
     /**
      * Calculate the interaction.

platforms/cpu/include/CpuCustomNonbondedForce.h

@@ -95,11 +95,11 @@ class CpuCustomNonbondedForce {
          already been set, and the smallest side of the periodic box is at least twice the cutoff
          distance.
 
-         @param boxSize             the X, Y, and Z widths of the periodic box
+         @param periodicBoxVectors    the vectors defining the periodic box
 
          --------------------------------------------------------------------------------------- */
 
-      void setPeriodic(OpenMM::RealVec& boxSize);
+      void setPeriodic(RealVec* periodicBoxVectors);
 
       /**---------------------------------------------------------------------------------------
 
@@ -127,8 +127,11 @@ private:
     bool cutoff;
     bool useSwitch;
     bool periodic;
+    bool triclinic;
     const CpuNeighborList* neighborList;
-    RealOpenMM periodicBoxSize[3];
+    float recipBoxSize[3];
+    RealVec periodicBoxVectors[3];
+    AlignedArray<fvec4> periodicBoxVec4;
     RealOpenMM cutoffDistance, switchingDistance;
     ThreadPool& threads;
     const std::vector<std::set<int> > exclusions;

platforms/cpu/include/CpuNeighborList.h

@@ -49,8 +49,6 @@ public:
     CpuNeighborList(int blockSize);
     void computeNeighborList(int numAtoms, const AlignedArray<float>& atomLocations, const std::vector<std::set<int> >& exclusions,
             const RealVec* periodicBoxVectors, bool usePeriodic, float maxDistance, ThreadPool& threads);
-    void computeNeighborList(int numAtoms, const AlignedArray<float>& atomLocations, const std::vector<std::set<int> >& exclusions,
-            const float* periodicBoxSize, bool usePeriodic, float maxDistance, ThreadPool& threads);
     int getNumBlocks() const;
     const std::vector<int>& getSortedAtoms() const;
     const std::vector<int>& getBlockNeighbors(int blockIndex) const;

platforms/cpu/src/CpuCustomManyParticleForce.cpp

@@ -122,8 +122,7 @@ void CpuCustomManyParticleForce::calculateIxn(AlignedArray<float>& posq, RealOpe
         particleNeighbors.resize(numParticles);
         for (int i = 0; i < numParticles; i++)
             particleNeighbors[i].clear();
-        float boxSizeFloat[] = {(float) periodicBoxSize[0], (float) periodicBoxSize[1], (float) periodicBoxSize[2]};
-        neighborList->computeNeighborList(numParticles, posq, exclusions, boxSizeFloat, usePeriodic, cutoffDistance, threads);
+        neighborList->computeNeighborList(numParticles, posq, exclusions, periodicBoxVectors, usePeriodic, cutoffDistance, threads);
         for (int blockIndex = 0; blockIndex < neighborList->getNumBlocks(); blockIndex++) {
             const vector<int>& neighbors = neighborList->getBlockNeighbors(blockIndex);
             const vector<char>& exclusions = neighborList->getBlockExclusions(blockIndex);
@@ -159,8 +158,8 @@ void CpuCustomManyParticleForce::calculateIxn(AlignedArray<float>& posq, RealOpe
 
 void CpuCustomManyParticleForce::threadComputeForce(ThreadPool& threads, int threadIndex) {
     vector<int> particleIndices(numParticlesPerSet);
-    fvec4 boxSize(periodicBoxSize[0], periodicBoxSize[1], periodicBoxSize[2], 0);
-    fvec4 invBoxSize((1/periodicBoxSize[0]), (1/periodicBoxSize[1]), (1/periodicBoxSize[2]), 0);
+    fvec4 boxSize(periodicBoxVectors[0][0], periodicBoxVectors[1][1], periodicBoxVectors[2][2], 0);
+    fvec4 invBoxSize(recipBoxSize[0], recipBoxSize[1], recipBoxSize[2], 0);
     float* forces = &(*threadForce)[threadIndex][0];
     ThreadData& data = *threadData[threadIndex];
     data.energy = 0;
@@ -201,15 +200,25 @@ void CpuCustomManyParticleForce::setUseCutoff(RealOpenMM distance) {
         neighborList = new CpuNeighborList(4);
 }
 
-void CpuCustomManyParticleForce::setPeriodic(RealVec& boxSize) {
-    assert(useCutoff);
-    assert(boxSize[0] >= 2.0*cutoffDistance);
-    assert(boxSize[1] >= 2.0*cutoffDistance);
-    assert(boxSize[2] >= 2.0*cutoffDistance);
+void CpuCustomManyParticleForce::setPeriodic(RealVec* periodicBoxVectors) {
+    assert(cutoff);
+    assert(periodicBoxVectors[0][0] >= 2.0*cutoffDistance);
+    assert(periodicBoxVectors[1][1] >= 2.0*cutoffDistance);
+    assert(periodicBoxVectors[2][2] >= 2.0*cutoffDistance);
     usePeriodic = true;
-    periodicBoxSize[0] = boxSize[0];
-    periodicBoxSize[1] = boxSize[1];
-    periodicBoxSize[2] = boxSize[2];
+    this->periodicBoxVectors[0] = periodicBoxVectors[0];
+    this->periodicBoxVectors[1] = periodicBoxVectors[1];
+    this->periodicBoxVectors[2] = periodicBoxVectors[2];
+    recipBoxSize[0] = (float) (1.0/periodicBoxVectors[0][0]);
+    recipBoxSize[1] = (float) (1.0/periodicBoxVectors[1][1]);
+    recipBoxSize[2] = (float) (1.0/periodicBoxVectors[2][2]);
+    periodicBoxVec4.resize(3);
+    periodicBoxVec4[0] = fvec4(periodicBoxVectors[0][0], periodicBoxVectors[0][1], periodicBoxVectors[0][2], 0);
+    periodicBoxVec4[1] = fvec4(periodicBoxVectors[1][0], periodicBoxVectors[1][1], periodicBoxVectors[1][2], 0);
+    periodicBoxVec4[2] = fvec4(periodicBoxVectors[2][0], periodicBoxVectors[2][1], periodicBoxVectors[2][2], 0);
+    triclinic = (periodicBoxVectors[0][1] != 0.0 || periodicBoxVectors[0][2] != 0.0 ||
+                 periodicBoxVectors[1][0] != 0.0 || periodicBoxVectors[1][2] != 0.0 ||
+                 periodicBoxVectors[2][0] != 0.0 || periodicBoxVectors[2][1] != 0.0);
 }
 
 void CpuCustomManyParticleForce::loopOverInteractions(vector<int>& availableParticles, vector<int>& particleSet, int loopIndex, int startIndex,
@@ -394,9 +403,16 @@ void CpuCustomManyParticleForce::calculateOneIxn(vector<int>& particleSet, RealO
 void CpuCustomManyParticleForce::computeDelta(const fvec4& posI, const fvec4& posJ, fvec4& deltaR, float& r2, const fvec4& boxSize, const fvec4& invBoxSize) const {
     deltaR = posJ-posI;
     if (usePeriodic) {
+        if (triclinic) {
+            deltaR -= periodicBoxVec4[2]*floorf(deltaR[2]*recipBoxSize[2]+0.5f);
+            deltaR -= periodicBoxVec4[1]*floorf(deltaR[1]*recipBoxSize[1]+0.5f);
+            deltaR -= periodicBoxVec4[0]*floorf(deltaR[0]*recipBoxSize[0]+0.5f);
+        }
+        else {
             fvec4 base = round(deltaR*invBoxSize)*boxSize;
             deltaR = deltaR-base;
         }
+    }
     r2 = dot3(deltaR, deltaR);
 }
 

platforms/cpu/src/CpuCustomNonbondedForce.cpp

@@ -98,15 +98,25 @@ void CpuCustomNonbondedForce::setUseSwitchingFunction(RealOpenMM distance) {
     switchingDistance = distance;
 }
 
-void CpuCustomNonbondedForce::setPeriodic(RealVec& boxSize) {
+void CpuCustomNonbondedForce::setPeriodic(RealVec* periodicBoxVectors) {
     assert(cutoff);
-    assert(boxSize[0] >= 2.0*cutoffDistance);
-    assert(boxSize[1] >= 2.0*cutoffDistance);
-    assert(boxSize[2] >= 2.0*cutoffDistance);
+    assert(periodicBoxVectors[0][0] >= 2.0*cutoffDistance);
+    assert(periodicBoxVectors[1][1] >= 2.0*cutoffDistance);
+    assert(periodicBoxVectors[2][2] >= 2.0*cutoffDistance);
     periodic = true;
-    periodicBoxSize[0] = boxSize[0];
-    periodicBoxSize[1] = boxSize[1];
-    periodicBoxSize[2] = boxSize[2];
+    this->periodicBoxVectors[0] = periodicBoxVectors[0];
+    this->periodicBoxVectors[1] = periodicBoxVectors[1];
+    this->periodicBoxVectors[2] = periodicBoxVectors[2];
+    recipBoxSize[0] = (float) (1.0/periodicBoxVectors[0][0]);
+    recipBoxSize[1] = (float) (1.0/periodicBoxVectors[1][1]);
+    recipBoxSize[2] = (float) (1.0/periodicBoxVectors[2][2]);
+    periodicBoxVec4.resize(3);
+    periodicBoxVec4[0] = fvec4(periodicBoxVectors[0][0], periodicBoxVectors[0][1], periodicBoxVectors[0][2], 0);
+    periodicBoxVec4[1] = fvec4(periodicBoxVectors[1][0], periodicBoxVectors[1][1], periodicBoxVectors[1][2], 0);
+    periodicBoxVec4[2] = fvec4(periodicBoxVectors[2][0], periodicBoxVectors[2][1], periodicBoxVectors[2][2], 0);
+    triclinic = (periodicBoxVectors[0][1] != 0.0 || periodicBoxVectors[0][2] != 0.0 ||
+                 periodicBoxVectors[1][0] != 0.0 || periodicBoxVectors[1][2] != 0.0 ||
+                 periodicBoxVectors[2][0] != 0.0 || periodicBoxVectors[2][1] != 0.0);
 }
 
 
@@ -155,8 +165,8 @@ void CpuCustomNonbondedForce::threadComputeForce(ThreadPool& threads, int thread
         ReferenceForce::setVariable(ReferenceForce::getVariablePointer(data.energyExpression, iter->first), iter->second);
         ReferenceForce::setVariable(ReferenceForce::getVariablePointer(data.forceExpression, iter->first), iter->second);
     }
-    fvec4 boxSize(periodicBoxSize[0], periodicBoxSize[1], periodicBoxSize[2], 0);
-    fvec4 invBoxSize((1/periodicBoxSize[0]), (1/periodicBoxSize[1]), (1/periodicBoxSize[2]), 0);
+    fvec4 boxSize(periodicBoxVectors[0][0], periodicBoxVectors[1][1], periodicBoxVectors[2][2], 0);
+    fvec4 invBoxSize(recipBoxSize[0], recipBoxSize[1], recipBoxSize[2], 0);
     if (groupInteractions.size() > 0) {
         // The user has specified interaction groups, so compute only the requested interactions.
         
@@ -266,8 +276,15 @@ void CpuCustomNonbondedForce::calculateOneIxn(int ii, int jj, ThreadData& data,
 void CpuCustomNonbondedForce::getDeltaR(const fvec4& posI, const fvec4& posJ, fvec4& deltaR, float& r2, const fvec4& boxSize, const fvec4& invBoxSize) const {
     deltaR = posJ-posI;
     if (periodic) {
+        if (triclinic) {
+            deltaR -= periodicBoxVec4[2]*floorf(deltaR[2]*recipBoxSize[2]+0.5f);
+            deltaR -= periodicBoxVec4[1]*floorf(deltaR[1]*recipBoxSize[1]+0.5f);
+            deltaR -= periodicBoxVec4[0]*floorf(deltaR[0]*recipBoxSize[0]+0.5f);
+        }
+        else {
             fvec4 base = round(deltaR*invBoxSize)*boxSize;
             deltaR = deltaR-base;
         }
+    }
     r2 = dot3(deltaR, deltaR);
 }

platforms/cpu/src/CpuKernels.cpp

@@ -512,8 +512,7 @@ double CpuCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeFo
     AlignedArray<float>& posq = data.posq;
     vector<RealVec>& posData = extractPositions(context);
     vector<RealVec>& forceData = extractForces(context);
-    RealVec boxSize = extractBoxSize(context);
-    float floatBoxSize[3] = {(float) boxSize[0], (float) boxSize[1], (float) boxSize[2]};
+    RealVec* boxVectors = extractBoxVectors(context);
     double energy = (includeReciprocal ? ewaldSelfEnergy : 0.0);
     bool ewald  = (nonbondedMethod == Ewald);
     bool pme  = (nonbondedMethod == PME);
@@ -559,7 +558,7 @@ double CpuCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeFo
                 }
         }
         if (needRecompute) {
-            neighborList->computeNeighborList(numParticles, posq, exclusions, floatBoxSize, data.isPeriodic, nonbondedCutoff+padding, data.threads);
+            neighborList->computeNeighborList(numParticles, posq, exclusions, boxVectors, data.isPeriodic, nonbondedCutoff+padding, data.threads);
             lastPositions = posData;
         }
         nonbonded->setUseCutoff(nonbondedCutoff, *neighborList, rfDielectric);
@@ -583,7 +582,7 @@ double CpuCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeFo
     if (includeReciprocal) {
         if (useOptimizedPme) {
             PmeIO io(&posq[0], &data.threadForce[0][0], numParticles);
-            Vec3 periodicBoxSize(boxSize[0], boxSize[1], boxSize[2]);
+            Vec3 periodicBoxSize(boxVectors[0][0], boxVectors[1][1], boxVectors[2][2]);
             optimizedPme.getAs<CalcPmeReciprocalForceKernel>().beginComputation(io, periodicBoxSize, includeEnergy);
             nonbondedEnergy += optimizedPme.getAs<CalcPmeReciprocalForceKernel>().finishComputation(io);
         }
@@ -596,7 +595,7 @@ double CpuCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeFo
         ReferenceLJCoulomb14 nonbonded14;
         refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, includeEnergy ? &energy : NULL, nonbonded14);
         if (data.isPeriodic)
-            energy += dispersionCoefficient/(boxSize[0]*boxSize[1]*boxSize[2]);
+            energy += dispersionCoefficient/(boxVectors[0][0]*boxVectors[1][1]*boxVectors[2][2]);
     }
     return energy;
 }
@@ -750,19 +749,18 @@ void CpuCalcCustomNonbondedForceKernel::initialize(const System& system, const C
 double CpuCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
     vector<RealVec>& posData = extractPositions(context);
     vector<RealVec>& forceData = extractForces(context);
-    RealVec& box = extractBoxSize(context);
-    float floatBoxSize[3] = {(float) box[0], (float) box[1], (float) box[2]};
+    RealVec* boxVectors = extractBoxVectors(context);
     double energy = 0;
     bool periodic = (nonbondedMethod == CutoffPeriodic);
     if (nonbondedMethod != NoCutoff) {
-        neighborList->computeNeighborList(numParticles, data.posq, exclusions, floatBoxSize, data.isPeriodic, nonbondedCutoff, data.threads);
+        neighborList->computeNeighborList(numParticles, data.posq, exclusions, boxVectors, data.isPeriodic, nonbondedCutoff, data.threads);
         nonbonded->setUseCutoff(nonbondedCutoff, *neighborList);
     }
     if (periodic) {
         double minAllowedSize = 2*nonbondedCutoff;
-        if (box[0] < minAllowedSize || box[1] < minAllowedSize || box[2] < minAllowedSize)
+        if (boxVectors[0][0] < minAllowedSize || boxVectors[1][1] < minAllowedSize || boxVectors[2][2] < minAllowedSize)
             throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
-        nonbonded->setPeriodic(box);
+        nonbonded->setPeriodic(boxVectors);
     }
     bool globalParamsChanged = false;
     for (int i = 0; i < (int) globalParameterNames.size(); i++) {
@@ -781,7 +779,7 @@ double CpuCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bool inc
         longRangeCoefficient = CustomNonbondedForceImpl::calcLongRangeCorrection(*forceCopy, context.getOwner());
         hasInitializedLongRangeCorrection = true;
     }
-    energy += longRangeCoefficient/(box[0]*box[1]*box[2]);
+    energy += longRangeCoefficient/(boxVectors[0][0]*boxVectors[1][1]*boxVectors[2][2]);
     return energy;
 }
 
@@ -989,13 +987,12 @@ void CpuCalcCustomGBForceKernel::initialize(const System& system, const CustomGB
 double CpuCalcCustomGBForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
     vector<RealVec>& forceData = extractForces(context);
     RealOpenMM energy = 0;
-    RealVec& box = extractBoxSize(context);
-    float floatBoxSize[3] = {(float) box[0], (float) box[1], (float) box[2]};
+    RealVec* boxVectors = extractBoxVectors(context);
     if (data.isPeriodic)
         ixn->setPeriodic(extractBoxSize(context));
     if (nonbondedMethod != NoCutoff) {
         vector<set<int> > noExclusions(numParticles);
-        neighborList->computeNeighborList(numParticles, data.posq, exclusions, floatBoxSize, data.isPeriodic, nonbondedCutoff, data.threads);
+        neighborList->computeNeighborList(numParticles, data.posq, exclusions, boxVectors, data.isPeriodic, nonbondedCutoff, data.threads);
         ixn->setUseCutoff(nonbondedCutoff, *neighborList);
     }
     map<string, double> globalParameters;
@@ -1052,6 +1049,7 @@ void CpuCalcCustomManyParticleForceKernel::initialize(const System& system, cons
     ixn = new CpuCustomManyParticleForce(force, data.threads);
     nonbondedMethod = CalcCustomManyParticleForceKernel::NonbondedMethod(force.getNonbondedMethod());
     cutoffDistance = force.getCutoffDistance();
+    data.isPeriodic = (nonbondedMethod == CutoffPeriodic);
 }
 
 double CpuCalcCustomManyParticleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
@@ -1059,11 +1057,11 @@ double CpuCalcCustomManyParticleForceKernel::execute(ContextImpl& context, bool
     for (int i = 0; i < (int) globalParameterNames.size(); i++)
         globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
     if (nonbondedMethod == CutoffPeriodic) {
-        RealVec& box = extractBoxSize(context);
+        RealVec* boxVectors = extractBoxVectors(context);
         double minAllowedSize = 2*cutoffDistance;
-        if (box[0] < minAllowedSize || box[1] < minAllowedSize || box[2] < minAllowedSize)
+        if (boxVectors[0][0] < minAllowedSize || boxVectors[1][1] < minAllowedSize || boxVectors[2][2] < minAllowedSize)
             throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
-        ixn->setPeriodic(box);
+        ixn->setPeriodic(boxVectors);
     }
     double energy = 0;
     ixn->calculateIxn(data.posq, particleParamArray, globalParameters, data.threadForce, includeForces, includeEnergy, energy);

platforms/cpu/src/CpuNeighborList.cpp

@@ -373,15 +373,6 @@ public:
 CpuNeighborList::CpuNeighborList(int blockSize) : blockSize(blockSize) {
 }
 
-void CpuNeighborList::computeNeighborList(int numAtoms, const AlignedArray<float>& atomLocations, const vector<set<int> >& exclusions,
-            const float* periodicBoxSize, bool usePeriodic, float maxDistance, ThreadPool& threads) {
-    RealVec vectors[3];
-    vectors[0] = RealVec(periodicBoxSize[0], 0, 0);
-    vectors[1] = RealVec(0, periodicBoxSize[1], 0);
-    vectors[2] = RealVec(0, 0, periodicBoxSize[2]);
-    computeNeighborList(numAtoms, atomLocations, exclusions, vectors, usePeriodic, maxDistance, threads);
-}
-
 void CpuNeighborList::computeNeighborList(int numAtoms, const AlignedArray<float>& atomLocations, const vector<set<int> >& exclusions,
             const RealVec* periodicBoxVectors, bool usePeriodic, float maxDistance, ThreadPool& threads) {
     int numBlocks = (numAtoms+blockSize-1)/blockSize;

platforms/cpu/src/CpuNonbondedForce.cpp

@@ -352,8 +352,6 @@ void CpuNonbondedForce::threadComputeDirect(ThreadPool& threads, int threadIndex
 
         // Now subtract off the exclusions, since they were implicitly included in the reciprocal space sum.
 
-        fvec4 boxSize(periodicBoxVectors[0][0], periodicBoxVectors[1][1], periodicBoxVectors[2][2], 0);
-        fvec4 invBoxSize(recipBoxSize[0], recipBoxSize[1], recipBoxSize[2], 0);
         for (int i = threadIndex; i < numberOfAtoms; i += numThreads) {
             fvec4 posI((float) atomCoordinates[i][0], (float) atomCoordinates[i][1], (float) atomCoordinates[i][2], 0.0f);
             for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter) {

platforms/cpu/tests/TestCpuCustomManyParticleForce.cpp

@@ -53,7 +53,17 @@ using namespace std;
 
 const double TOL = 1e-5;
 
-void validateAxilrodTeller(CustomManyParticleForce* force, const vector<Vec3>& positions, const vector<const int*>& expectedSets, double boxSize) {
+Vec3 computeDelta(const Vec3& pos1, const Vec3& pos2, bool periodic, const Vec3* periodicBoxVectors) {
+    Vec3 diff = pos1-pos2;
+    if (periodic) {
+        diff -= periodicBoxVectors[2]*floor(diff[2]/periodicBoxVectors[2][2]+0.5);
+        diff -= periodicBoxVectors[1]*floor(diff[1]/periodicBoxVectors[1][1]+0.5);
+        diff -= periodicBoxVectors[0]*floor(diff[0]/periodicBoxVectors[0][0]+0.5);
+    }
+    return diff;
+}
+
+void validateAxilrodTeller(CustomManyParticleForce* force, const vector<Vec3>& positions, const vector<const int*>& expectedSets, double boxSize, bool triclinic) {
     // Create a System and Context.
 
     int numParticles = force->getNumParticles();
@@ -61,7 +71,18 @@ void validateAxilrodTeller(CustomManyParticleForce* force, const vector<Vec3>& p
     System system;
     for (int i = 0; i < numParticles; i++)
         system.addParticle(1.0);
-    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    Vec3 boxVectors[3];
+    if (triclinic) {
+        boxVectors[0] = Vec3(boxSize, 0, 0);
+        boxVectors[1] = Vec3(0.2*boxSize, boxSize, 0);
+        boxVectors[2] = Vec3(-0.3*boxSize, -0.1*boxSize, boxSize);
+    }
+    else {
+        boxVectors[0] = Vec3(boxSize, 0, 0);
+        boxVectors[1] = Vec3(0, boxSize, 0);
+        boxVectors[2] = Vec3(0, 0, boxSize);
+    }
+    system.setDefaultPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
     system.addForce(force);
     VerletIntegrator integrator(0.001);
     CpuPlatform platform;
@@ -73,20 +94,14 @@ void validateAxilrodTeller(CustomManyParticleForce* force, const vector<Vec3>& p
     // See if the energy matches the expected value.
 
     double expectedEnergy = 0;
+    bool periodic = (nonbondedMethod == CustomManyParticleForce::CutoffPeriodic);
     for (int i = 0; i < (int) expectedSets.size(); i++) {
         int p1 = expectedSets[i][0];
         int p2 = expectedSets[i][1];
         int p3 = expectedSets[i][2];
-        Vec3 d12 = positions[p2]-positions[p1];
-        Vec3 d13 = positions[p3]-positions[p1];
-        Vec3 d23 = positions[p3]-positions[p2];
-        if (nonbondedMethod == CustomManyParticleForce::CutoffPeriodic) {
-            for (int j = 0; j < 3; j++) {
-                d12[j] -= floor(d12[j]/boxSize+0.5f)*boxSize;
-                d13[j] -= floor(d13[j]/boxSize+0.5f)*boxSize;
-                d23[j] -= floor(d23[j]/boxSize+0.5f)*boxSize;
-            }
-        }
+        Vec3 d12 = computeDelta(positions[p2], positions[p1], periodic, boxVectors);
+        Vec3 d13 = computeDelta(positions[p3], positions[p1], periodic, boxVectors);
+        Vec3 d23 = computeDelta(positions[p3], positions[p2], periodic, boxVectors);
         double r12 = sqrt(d12.dot(d12));
         double r13 = sqrt(d13.dot(d13));
         double r23 = sqrt(d23.dot(d23));
@@ -210,7 +225,7 @@ void testNoCutoff() {
     positions.push_back(Vec3(0.4, 0, -0.8));
     int sets[4][3] = {{0,1,2}, {1,2,3}, {2,3,0}, {3,0,1}};
     vector<const int*> expectedSets(&sets[0], &sets[4]);
-    validateAxilrodTeller(force, positions, expectedSets, 2.0);
+    validateAxilrodTeller(force, positions, expectedSets, 2.0, false);
 }
 
 void testCutoff() {
@@ -235,7 +250,7 @@ void testCutoff() {
     positions.push_back(Vec3(0.2, 0.5, -0.1));
     int sets[7][3] = {{0,1,2}, {0,1,3}, {0,1,4}, {0,2,4}, {0,3,4}, {1,2,4}, {1,3,4}};
     vector<const int*> expectedSets(&sets[0], &sets[7]);
-    validateAxilrodTeller(force, positions, expectedSets, 2.0);
+    validateAxilrodTeller(force, positions, expectedSets, 2.0, false);
 }
 
 void testPeriodic() {
@@ -261,7 +276,33 @@ void testPeriodic() {
     double boxSize = 2.1;
     int sets[5][3] = {{0,1,3}, {0,1,4}, {0,2,4}, {0,3,4}, {1,3,4}};
     vector<const int*> expectedSets(&sets[0], &sets[5]);
-    validateAxilrodTeller(force, positions, expectedSets, boxSize);
+    validateAxilrodTeller(force, positions, expectedSets, boxSize, false);
+}
+
+void testTriclinic() {
+    CustomManyParticleForce* force = new CustomManyParticleForce(3,
+        "C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
+        "theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
+        "r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
+    force->addGlobalParameter("C", 1.5);
+    force->setNonbondedMethod(CustomManyParticleForce::CutoffPeriodic);
+    force->setCutoffDistance(1.05);
+    vector<double> params;
+    force->addParticle(params);
+    force->addParticle(params);
+    force->addParticle(params);
+    force->addParticle(params);
+    force->addParticle(params);
+    vector<Vec3> positions;
+    positions.push_back(Vec3(0, 0, 0));
+    positions.push_back(Vec3(1, 0, 0));
+    positions.push_back(Vec3(0, 1.1, 0.3));
+    positions.push_back(Vec3(0.4, 0, -0.8));
+    positions.push_back(Vec3(0.2, 0.5, -0.1));
+    double boxSize = 2.1;
+    int sets[4][3] = {{0,1,3}, {0,1,4}, {0,3,4}, {1,3,4}};
+    vector<const int*> expectedSets(&sets[0], &sets[4]);
+    validateAxilrodTeller(force, positions, expectedSets, boxSize, true);
 }
 
 void testExclusions() {
@@ -286,7 +327,7 @@ void testExclusions() {
     force->addExclusion(0, 3);
     int sets[5][3] = {{0,1,4}, {1,2,3}, {1,2,4}, {1,3,4}, {2,3,4}};
     vector<const int*> expectedSets(&sets[0], &sets[5]);
-    validateAxilrodTeller(force, positions, expectedSets, 2.0);
+    validateAxilrodTeller(force, positions, expectedSets, 2.0, false);
 }
 
 void testAllTerms() {
@@ -680,6 +721,7 @@ int main() {
         testNoCutoff();
         testCutoff();
         testPeriodic();
+        testTriclinic();
         testExclusions();
         testAllTerms();
         testParameters();

platforms/cpu/tests/TestCpuCustomNonbondedForce.cpp

@@ -224,6 +224,65 @@ void testPeriodic() {
     ASSERT_EQUAL_TOL(1.9+1+0.9, state.getPotentialEnergy(), TOL);
 }
 
+void testTriclinic() {
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    Vec3 a(3.1, 0, 0);
+    Vec3 b(0.4, 3.5, 0);
+    Vec3 c(-0.1, -0.5, 4.0);
+    system.setDefaultPeriodicBoxVectors(a, b, c);
+    VerletIntegrator integrator(0.01);
+    CustomNonbondedForce* nonbonded = new CustomNonbondedForce("r");
+    nonbonded->addParticle(vector<double>());
+    nonbonded->addParticle(vector<double>());
+    nonbonded->setNonbondedMethod(CustomNonbondedForce::CutoffPeriodic);
+    const double cutoff = 1.5;
+    nonbonded->setCutoffDistance(cutoff);
+    system.addForce(nonbonded);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    for (int iteration = 0; iteration < 50; iteration++) {
+        // Generate random positions for the two particles.
+
+        positions[0] = a*genrand_real2(sfmt) + b*genrand_real2(sfmt) + c*genrand_real2(sfmt);
+        positions[1] = a*genrand_real2(sfmt) + b*genrand_real2(sfmt) + c*genrand_real2(sfmt);
+        context.setPositions(positions);
+
+        // Loop over all possible periodic copies and find the nearest one.
+
+        Vec3 delta;
+        double distance2 = 100.0;
+        for (int i = -1; i < 2; i++)
+            for (int j = -1; j < 2; j++)
+                for (int k = -1; k < 2; k++) {
+                    Vec3 d = positions[1]-positions[0]+a*i+b*j+c*k;
+                    if (d.dot(d) < distance2) {
+                        delta = d;
+                        distance2 = d.dot(d);
+                    }
+                }
+        double distance = sqrt(distance2);
+
+        // See if the force and energy are correct.
+
+        State state = context.getState(State::Forces | State::Energy);
+        if (distance >= cutoff) {
+            ASSERT_EQUAL(0.0, state.getPotentialEnergy());
+            ASSERT_EQUAL_VEC(Vec3(0, 0, 0), state.getForces()[0], 0);
+            ASSERT_EQUAL_VEC(Vec3(0, 0, 0), state.getForces()[1], 0);
+        }
+        else {
+            const Vec3 force = delta/sqrt(delta.dot(delta));
+            ASSERT_EQUAL_TOL(distance, state.getPotentialEnergy(), TOL);
+            ASSERT_EQUAL_VEC(force, state.getForces()[0], TOL);
+            ASSERT_EQUAL_VEC(-force, state.getForces()[1], TOL);
+        }
+    }
+}
+
 void testContinuous1DFunction() {
     System system;
     system.addParticle(1.0);
@@ -852,6 +911,7 @@ int main() {
         testExclusions();
         testCutoff();
         testPeriodic();
+        testTriclinic();
         testContinuous1DFunction();
         testContinuous2DFunction();
         testContinuous3DFunction();

platforms/reference/tests/TestReferenceCustomManyParticleForce.cpp

@@ -53,7 +53,17 @@ using namespace std;
 
 const double TOL = 1e-5;
 
-void validateAxilrodTeller(CustomManyParticleForce* force, const vector<Vec3>& positions, const vector<const int*>& expectedSets, double boxSize) {
+Vec3 computeDelta(const Vec3& pos1, const Vec3& pos2, bool periodic, const Vec3* periodicBoxVectors) {
+    Vec3 diff = pos1-pos2;
+    if (periodic) {
+        diff -= periodicBoxVectors[2]*floor(diff[2]/periodicBoxVectors[2][2]+0.5);
+        diff -= periodicBoxVectors[1]*floor(diff[1]/periodicBoxVectors[1][1]+0.5);
+        diff -= periodicBoxVectors[0]*floor(diff[0]/periodicBoxVectors[0][0]+0.5);
+    }
+    return diff;
+}
+
+void validateAxilrodTeller(CustomManyParticleForce* force, const vector<Vec3>& positions, const vector<const int*>& expectedSets, double boxSize, bool triclinic) {
     // Create a System and Context.
     
     int numParticles = force->getNumParticles();
@@ -61,7 +71,18 @@ void validateAxilrodTeller(CustomManyParticleForce* force, const vector<Vec3>& p
     System system;
     for (int i = 0; i < numParticles; i++)
         system.addParticle(1.0);
-    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    Vec3 boxVectors[3];
+    if (triclinic) {
+        boxVectors[0] = Vec3(boxSize, 0, 0);
+        boxVectors[1] = Vec3(0.2*boxSize, boxSize, 0);
+        boxVectors[2] = Vec3(-0.3*boxSize, -0.1*boxSize, boxSize);
+    }
+    else {
+        boxVectors[0] = Vec3(boxSize, 0, 0);
+        boxVectors[1] = Vec3(0, boxSize, 0);
+        boxVectors[2] = Vec3(0, 0, boxSize);
+    }
+    system.setDefaultPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
     system.addForce(force);
     VerletIntegrator integrator(0.001);
     ReferencePlatform platform;
@@ -73,20 +94,14 @@ void validateAxilrodTeller(CustomManyParticleForce* force, const vector<Vec3>& p
     // See if the energy matches the expected value.
     
     double expectedEnergy = 0;
+    bool periodic = (nonbondedMethod == CustomManyParticleForce::CutoffPeriodic);
     for (int i = 0; i < (int) expectedSets.size(); i++) {
         int p1 = expectedSets[i][0];
         int p2 = expectedSets[i][1];
         int p3 = expectedSets[i][2];
-        Vec3 d12 = positions[p2]-positions[p1];
-        Vec3 d13 = positions[p3]-positions[p1];
-        Vec3 d23 = positions[p3]-positions[p2];
-        if (nonbondedMethod == CustomManyParticleForce::CutoffPeriodic) {
-            for (int j = 0; j < 3; j++) {
-                d12[j] -= floor(d12[j]/boxSize+0.5f)*boxSize;
-                d13[j] -= floor(d13[j]/boxSize+0.5f)*boxSize;
-                d23[j] -= floor(d23[j]/boxSize+0.5f)*boxSize;
-            }
-        }
+        Vec3 d12 = computeDelta(positions[p2], positions[p1], periodic, boxVectors);
+        Vec3 d13 = computeDelta(positions[p3], positions[p1], periodic, boxVectors);
+        Vec3 d23 = computeDelta(positions[p3], positions[p2], periodic, boxVectors);
         double r12 = sqrt(d12.dot(d12));
         double r13 = sqrt(d13.dot(d13));
         double r23 = sqrt(d23.dot(d23));
@@ -210,7 +225,7 @@ void testNoCutoff() {
     positions.push_back(Vec3(0.4, 0, -0.8));
     int sets[4][3] = {{0,1,2}, {1,2,3}, {2,3,0}, {3,0,1}};
     vector<const int*> expectedSets(&sets[0], &sets[4]);
-    validateAxilrodTeller(force, positions, expectedSets, 2.0);
+    validateAxilrodTeller(force, positions, expectedSets, 2.0, false);
 }
 
 void testCutoff() {
@@ -235,7 +250,7 @@ void testCutoff() {
     positions.push_back(Vec3(0.2, 0.5, -0.1));
     int sets[7][3] = {{0,1,2}, {0,1,3}, {0,1,4}, {0,2,4}, {0,3,4}, {1,2,4}, {1,3,4}};
     vector<const int*> expectedSets(&sets[0], &sets[7]);
-    validateAxilrodTeller(force, positions, expectedSets, 2.0);
+    validateAxilrodTeller(force, positions, expectedSets, 2.0, false);
 }
 
 void testPeriodic() {
@@ -261,7 +276,33 @@ void testPeriodic() {
     double boxSize = 2.1;
     int sets[5][3] = {{0,1,3}, {0,1,4}, {0,2,4}, {0,3,4}, {1,3,4}};
     vector<const int*> expectedSets(&sets[0], &sets[5]);
-    validateAxilrodTeller(force, positions, expectedSets, boxSize);
+    validateAxilrodTeller(force, positions, expectedSets, boxSize, false);
+}
+
+void testTriclinic() {
+    CustomManyParticleForce* force = new CustomManyParticleForce(3,
+        "C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
+        "theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
+        "r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
+    force->addGlobalParameter("C", 1.5);
+    force->setNonbondedMethod(CustomManyParticleForce::CutoffPeriodic);
+    force->setCutoffDistance(1.05);
+    vector<double> params;
+    force->addParticle(params);
+    force->addParticle(params);
+    force->addParticle(params);
+    force->addParticle(params);
+    force->addParticle(params);
+    vector<Vec3> positions;
+    positions.push_back(Vec3(0, 0, 0));
+    positions.push_back(Vec3(1, 0, 0));
+    positions.push_back(Vec3(0, 1.1, 0.3));
+    positions.push_back(Vec3(0.4, 0, -0.8));
+    positions.push_back(Vec3(0.2, 0.5, -0.1));
+    double boxSize = 2.1;
+    int sets[4][3] = {{0,1,3}, {0,1,4}, {0,3,4}, {1,3,4}};
+    vector<const int*> expectedSets(&sets[0], &sets[4]);
+    validateAxilrodTeller(force, positions, expectedSets, boxSize, true);
 }
 
 void testExclusions() {
@@ -286,7 +327,7 @@ void testExclusions() {
     force->addExclusion(0, 3);
     int sets[5][3] = {{0,1,4}, {1,2,3}, {1,2,4}, {1,3,4}, {2,3,4}};
     vector<const int*> expectedSets(&sets[0], &sets[5]);
-    validateAxilrodTeller(force, positions, expectedSets, 2.0);
+    validateAxilrodTeller(force, positions, expectedSets, 2.0, false);
 }
 
 void testAllTerms() {
@@ -592,6 +633,7 @@ int main() {
         testNoCutoff();
         testCutoff();
         testPeriodic();
+        testTriclinic();
         testExclusions();
         testAllTerms();
         testParameters();