Began implementing triclinic boxes for CPU platform

platforms/cpu/include/CpuNonbondedForce.h

@@ -83,11 +83,11 @@ class CpuNonbondedForce {
          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(float* periodicBoxSize);
+      void setPeriodic(RealVec* periodicBoxVectors);
        
       /**---------------------------------------------------------------------------------------
       
@@ -161,11 +161,14 @@ protected:
         bool cutoff;
         bool useSwitch;
         bool periodic;
+        bool triclinic;
         bool ewald;
         bool pme;
         bool tableIsValid;
         const CpuNeighborList* neighborList;
-        float periodicBoxSize[3];
+        float recipBoxSize[3];
+        RealVec periodicBoxVectors[3];
+        AlignedArray<fvec4> periodicBoxVec4;
         float cutoffDistance, switchingDistance;
         float krf, crf;
         float alphaEwald;

platforms/cpu/src/CpuKernels.cpp

@@ -73,6 +73,11 @@ static RealVec& extractBoxSize(ContextImpl& context) {
     return *(RealVec*) data->periodicBoxSize;
 }
 
+static RealVec* extractBoxVectors(ContextImpl& context) {
+    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
+    return (RealVec*) data->periodicBoxVectors;
+}
+
 static ReferenceConstraints& extractConstraints(ContextImpl& context) {
     ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
     return *(ReferenceConstraints*) data->constraints;
@@ -147,19 +152,36 @@ public:
 
         AlignedArray<float>& posq = data.posq;
         vector<RealVec>& posData = extractPositions(context);
-        RealVec boxSize = extractBoxSize(context);
-        double invBoxSize[3] = {1/boxSize[0], 1/boxSize[1], 1/boxSize[2]};
+        RealVec* boxVectors = extractBoxVectors(context);
+        double boxSize[3] = {boxVectors[0][0], boxVectors[1][1], boxVectors[2][2]};
+        double invBoxSize[3] = {1/boxVectors[0][0], 1/boxVectors[1][1], 1/boxVectors[2][2]};
+        bool triclinic = (boxVectors[0][1] != 0 || boxVectors[0][2] != 0 || boxVectors[1][0] != 0 || boxVectors[1][2] != 0 || boxVectors[2][0] != 0 || boxVectors[2][1] != 0);
         int numParticles = context.getSystem().getNumParticles();
         int numThreads = threads.getNumThreads();
         int start = threadIndex*numParticles/numThreads;
         int end = (threadIndex+1)*numParticles/numThreads;
-        if (data.isPeriodic)
-            for (int i = start; i < end; i++)
-                for (int j = 0; j < 3; j++) {
-                    RealOpenMM x = posData[i][j];
-                    double base = floor(x*invBoxSize[j])*boxSize[j];
-                    posq[4*i+j] = (float) (x-base);
+        if (data.isPeriodic) {
+            if (triclinic) {
+                for (int i = start; i < end; i++) {
+                    RealVec pos = posData[i];
+                    pos -= boxVectors[2]*floor(pos[2]*invBoxSize[2]);
+                    pos -= boxVectors[1]*floor(pos[1]*invBoxSize[1]);
+                    pos -= boxVectors[0]*floor(pos[0]*invBoxSize[0]);
+                    posq[4*i] = (float) pos[0];
+                    posq[4*i+1] = (float) pos[1];
+                    posq[4*i+2] = (float) pos[2];
+                }
+            }
+            else {
+                for (int i = start; i < end; i++) {
+                    for (int j = 0; j < 3; j++) {
+                        RealOpenMM x = posData[i][j];
+                        double base = floor(x*invBoxSize[j])*boxSize[j];
+                        posq[4*i+j] = (float) (x-base);
+                    }
                 }
+            }
+        }
         else
             for (int i = start; i < end; i++) {
                 posq[4*i] = (float) posData[i][0];
@@ -543,10 +565,11 @@ double CpuCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeFo
         nonbonded->setUseCutoff(nonbondedCutoff, *neighborList, rfDielectric);
     }
     if (data.isPeriodic) {
+        RealVec* boxVectors = extractBoxVectors(context);
         double minAllowedSize = 1.999999*nonbondedCutoff;
-        if (boxSize[0] < minAllowedSize || boxSize[1] < minAllowedSize || boxSize[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(floatBoxSize);
+        nonbonded->setPeriodic(boxVectors);
     }
     if (ewald)
         nonbonded->setUseEwald(ewaldAlpha, kmax[0], kmax[1], kmax[2]);

platforms/cpu/src/CpuNonbondedForce.cpp

@@ -103,20 +103,30 @@ void CpuNonbondedForce::setUseSwitchingFunction(float distance) {
      also 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 CpuNonbondedForce::setPeriodic(float* periodicBoxSize) {
+  void CpuNonbondedForce::setPeriodic(RealVec* periodicBoxVectors) {
 
     assert(cutoff);
-    assert(periodicBoxSize[0] >= 2*cutoffDistance);
-    assert(periodicBoxSize[1] >= 2*cutoffDistance);
-    assert(periodicBoxSize[2] >= 2*cutoffDistance);
+    assert(periodicBoxVectors[0][0] >= 2.0*cutoffDistance);
+    assert(periodicBoxVectors[1][1] >= 2.0*cutoffDistance);
+    assert(periodicBoxVectors[2][2] >= 2.0*cutoffDistance);
     periodic = true;
-    this->periodicBoxSize[0] = periodicBoxSize[0];
-    this->periodicBoxSize[1] = periodicBoxSize[1];
-    this->periodicBoxSize[2] = periodicBoxSize[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);
   }
 
   /**---------------------------------------------------------------------------------------
@@ -186,7 +196,7 @@ void CpuNonbondedForce::calculateReciprocalIxn(int numberOfAtoms, float* posq, c
     int kmax                       = (ewald ? max(numRx, max(numRy,numRz)) : 0);
     float factorEwald              = -1 / (4*alphaEwald*alphaEwald);
     float TWO_PI                   = 2.0 * PI_M;
-    float recipCoeff               = (float)(ONE_4PI_EPS0*4*PI_M/(periodicBoxSize[0] * periodicBoxSize[1] * periodicBoxSize[2]) /epsilon);
+    float recipCoeff               = (float)(ONE_4PI_EPS0*4*PI_M/(periodicBoxVectors[0][0] * periodicBoxVectors[1][1] * periodicBoxVectors[2][2]) /epsilon);
 
     if (pme) {
         pme_t pmedata;
@@ -194,9 +204,8 @@ void CpuNonbondedForce::calculateReciprocalIxn(int numberOfAtoms, float* posq, c
         vector<RealOpenMM> charges(numberOfAtoms);
         for (int i = 0; i < numberOfAtoms; i++)
             charges[i] = posq[4*i+3];
-        RealVec boxVectors[3] = {Vec3(periodicBoxSize[0], 0, 0), Vec3(0, periodicBoxSize[1], 0), Vec3(0, 0, periodicBoxSize[2])};
         RealOpenMM recipEnergy = 0.0;
-        pme_exec(pmedata, atomCoordinates, forces, charges, boxVectors, &recipEnergy);
+        pme_exec(pmedata, atomCoordinates, forces, charges, periodicBoxVectors, &recipEnergy);
         if (totalEnergy)
             *totalEnergy += recipEnergy;
         pme_destroy(pmedata);
@@ -208,7 +217,7 @@ void CpuNonbondedForce::calculateReciprocalIxn(int numberOfAtoms, float* posq, c
 
         // setup reciprocal box
 
-        float recipBoxSize[3] = { TWO_PI / periodicBoxSize[0], TWO_PI / periodicBoxSize[1], TWO_PI / periodicBoxSize[2]};
+        float recipBoxSize[3] = { TWO_PI / periodicBoxVectors[0][0], TWO_PI / periodicBoxVectors[1][1], TWO_PI / periodicBoxVectors[2][2]};
 
 
         // setup K-vectors
@@ -329,8 +338,8 @@ void CpuNonbondedForce::threadComputeDirect(ThreadPool& threads, int threadIndex
     threadEnergy[threadIndex] = 0;
     double* energyPtr = (includeEnergy ? &threadEnergy[threadIndex] : NULL);
     float* forces = &(*threadForce)[threadIndex][0];
-    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 (ewald || pme) {
         // Compute the interactions from the neighbor list.
 
@@ -343,8 +352,8 @@ void CpuNonbondedForce::threadComputeDirect(ThreadPool& threads, int threadIndex
 
         // Now subtract off the exclusions, since they were implicitly included in the reciprocal space sum.
 
-        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);
         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) {
@@ -453,8 +462,15 @@ void CpuNonbondedForce::calculateOneIxn(int ii, int jj, float* forces, double* t
 void CpuNonbondedForce::getDeltaR(const fvec4& posI, const fvec4& posJ, fvec4& deltaR, float& r2, bool periodic, const fvec4& boxSize, const fvec4& invBoxSize) const {
     deltaR = posJ-posI;
     if (periodic) {
-        fvec4 base = round(deltaR*invBoxSize)*boxSize;
-        deltaR = deltaR-base;
+        if (triclinic) {
+            deltaR -= periodicBoxVec4[2]*floor(deltaR[2]*recipBoxSize[2]+0.5f);
+            deltaR -= periodicBoxVec4[1]*floor(deltaR[1]*recipBoxSize[1]+0.5f);
+            deltaR -= periodicBoxVec4[0]*floor(deltaR[0]*recipBoxSize[0]+0.5f);
+        }
+        else {
+            fvec4 base = round(deltaR*invBoxSize)*boxSize;
+            deltaR = deltaR-base;
+        }
     }
     r2 = dot3(deltaR, deltaR);
 }

platforms/cpu/src/CpuNonbondedForceVec4.cpp

@@ -262,9 +262,22 @@ void CpuNonbondedForceVec4::getDeltaR(const float* posI, const fvec4& x, const f
     dy = y-posI[1];
     dz = z-posI[2];
     if (periodic) {
-        dx -= round(dx*invBoxSize[0])*boxSize[0];
-        dy -= round(dy*invBoxSize[1])*boxSize[1];
-        dz -= round(dz*invBoxSize[2])*boxSize[2];
+        if (triclinic) {
+            fvec4 scale3 = floor(dz*recipBoxSize[2]+0.5f);
+            dx -= scale3*periodicBoxVectors[2][0];
+            dy -= scale3*periodicBoxVectors[2][1];
+            dz -= scale3*periodicBoxVectors[2][2];
+            fvec4 scale2 = floor(dy*recipBoxSize[1]+0.5f);
+            dx -= scale2*periodicBoxVectors[1][0];
+            dy -= scale2*periodicBoxVectors[1][1];
+            fvec4 scale1 = floor(dx*recipBoxSize[0]+0.5f);
+            dx -= scale1*periodicBoxVectors[0][0];
+        }
+        else {
+            dx -= round(dx*invBoxSize[0])*boxSize[0];
+            dy -= round(dy*invBoxSize[1])*boxSize[1];
+            dz -= round(dz*invBoxSize[2])*boxSize[2];
+        }
     }
     r2 = dx*dx + dy*dy + dz*dz;
 }

platforms/cpu/src/CpuNonbondedForceVec8.cpp

@@ -292,9 +292,22 @@ void CpuNonbondedForceVec8::getDeltaR(const float* posI, const fvec8& x, const f
     dy = y-posI[1];
     dz = z-posI[2];
     if (periodic) {
-        dx -= round(dx*invBoxSize[0])*boxSize[0];
-        dy -= round(dy*invBoxSize[1])*boxSize[1];
-        dz -= round(dz*invBoxSize[2])*boxSize[2];
+        if (triclinic) {
+            fvec8 scale3 = floor(dz*recipBoxSize[2]+0.5f);
+            dx -= scale3*periodicBoxVectors[2][0];
+            dy -= scale3*periodicBoxVectors[2][1];
+            dz -= scale3*periodicBoxVectors[2][2];
+            fvec8 scale2 = floor(dy*recipBoxSize[1]+0.5f);
+            dx -= scale2*periodicBoxVectors[1][0];
+            dy -= scale2*periodicBoxVectors[1][1];
+            fvec8 scale1 = floor(dx*recipBoxSize[0]+0.5f);
+            dx -= scale1*periodicBoxVectors[0][0];
+        }
+        else {
+            dx -= round(dx*invBoxSize[0])*boxSize[0];
+            dy -= round(dy*invBoxSize[1])*boxSize[1];
+            dz -= round(dz*invBoxSize[2])*boxSize[2];
+        }
     }
     r2 = dx*dx + dy*dy + dz*dz;
 }

platforms/cpu/tests/TestCpuNonbondedForce.cpp

@@ -353,6 +353,67 @@ void testPeriodic() {
     ASSERT_EQUAL_TOL(2*ONE_4PI_EPS0*(1.0)*(1.0+krf*1.0-crf), 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);
+    NonbondedForce* nonbonded = new NonbondedForce();
+    nonbonded->addParticle(1.0, 1, 0);
+    nonbonded->addParticle(1.0, 1, 0);
+    nonbonded->setNonbondedMethod(NonbondedForce::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);
+    const double eps = 78.3;
+    const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0);
+    const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0);
+    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*ONE_4PI_EPS0*(-1.0/(distance*distance*distance)+2.0*krf);
+            ASSERT_EQUAL_TOL(ONE_4PI_EPS0*(1.0/distance+krf*distance*distance-crf), state.getPotentialEnergy(), 1e-4);
+            ASSERT_EQUAL_VEC(force, state.getForces()[0], TOL);
+            ASSERT_EQUAL_VEC(-force, state.getForces()[1], TOL);
+        }
+    }
+}
 
 void testLargeSystem() {
     const int numMolecules = 600;
@@ -635,6 +696,7 @@ int main(int argc, char* argv[]) {
         testCutoff();
         testCutoff14();
         testPeriodic();
+        testTriclinic();
         testLargeSystem();
         testDispersionCorrection();
         testChangingParameters();