CPU PME works with triclinic boxes

olla/include/openmm/kernels.h

@@ -1223,12 +1223,12 @@ public:
     virtual void initialize(int gridx, int gridy, int gridz, int numParticles, double alpha) = 0;
     /**
      * Begin computing the force and energy.
-     * 
-     * @param io               an object that coordinates data transfer
-     * @param periodicBoxSize  the size of the periodic box (measured in nm)
-     * @param includeEnergy    true if potential energy should be computed
+     *
+     * @param io                  an object that coordinates data transfer
+     * @param periodicBoxVectors  the vectors defining the periodic box (measured in nm)
+     * @param includeEnergy       true if potential energy should be computed
      */
-    virtual void beginComputation(IO& io, Vec3 periodicBoxSize, bool includeEnergy) = 0;
+    virtual void beginComputation(IO& io, const Vec3* periodicBoxVectors, bool includeEnergy) = 0;
     /**
      * Finish computing the force and energy.
      * 

platforms/cpu/src/CpuKernels.cpp

@@ -582,8 +582,8 @@ double CpuCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeFo
     if (includeReciprocal) {
         if (useOptimizedPme) {
             PmeIO io(&posq[0], &data.threadForce[0][0], numParticles);
-            Vec3 periodicBoxSize(boxVectors[0][0], boxVectors[1][1], boxVectors[2][2]);
-            optimizedPme.getAs<CalcPmeReciprocalForceKernel>().beginComputation(io, periodicBoxSize, includeEnergy);
+            Vec3 periodicBoxVectors[3] = {boxVectors[0], boxVectors[1], boxVectors[2]};
+            optimizedPme.getAs<CalcPmeReciprocalForceKernel>().beginComputation(io, periodicBoxVectors, includeEnergy);
             nonbondedEnergy += optimizedPme.getAs<CalcPmeReciprocalForceKernel>().finishComputation(io);
         }
         else

platforms/cpu/tests/TestCpuEwald.cpp

@@ -201,6 +201,57 @@ void testEwald2Ions() {
     ASSERT_EQUAL_TOL(-217.276, state.getPotentialEnergy(), 0.01/*10*TOL*/);
 }
 
+void testTriclinic() {
+    // Create a triclinic box containing eight particles.
+
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(2.5, 0, 0), Vec3(0.5, 3.0, 0), Vec3(0.7, 0.9, 3.5));
+    for (int i = 0; i < 8; i++)
+        system.addParticle(1.0);
+    NonbondedForce* force = new NonbondedForce();
+    system.addForce(force);
+    force->setNonbondedMethod(NonbondedForce::PME);
+    force->setCutoffDistance(1.0);
+    force->setPMEParameters(3.45891, 32, 40, 48);
+    for (int i = 0; i < 4; i++)
+        force->addParticle(-1, 0.440104, 0.4184); // Cl parameters
+    for (int i = 0; i < 4; i++)
+        force->addParticle(1, 0.332840, 0.0115897); // Na parameters
+    vector<Vec3> positions(8);
+    positions[0] = Vec3(1.744, 2.788, 3.162);
+    positions[1] = Vec3(1.048, 0.762, 2.340);
+    positions[2] = Vec3(2.489, 1.570, 2.817);
+    positions[3] = Vec3(1.027, 1.893, 3.271);
+    positions[4] = Vec3(0.937, 0.825, 0.009);
+    positions[5] = Vec3(2.290, 1.887, 3.352);
+    positions[6] = Vec3(1.266, 1.111, 2.894);
+    positions[7] = Vec3(0.933, 1.862, 3.490);
+
+    // Compute the forces and energy.
+
+    VerletIntegrator integ(0.001);
+    Context context(system, integ, platform);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+
+    // Compare them to values computed by Gromacs.
+
+    double expectedEnergy = -963.370;
+    vector<Vec3> expectedForce(8);
+    expectedForce[0] = Vec3(4.25253e+01, -1.23503e+02, 1.22139e+02);
+    expectedForce[1] = Vec3(9.74752e+01, 1.68213e+02, 1.93169e+02);
+    expectedForce[2] = Vec3(-1.50348e+02, 1.29165e+02, 3.70435e+02);
+    expectedForce[3] = Vec3(9.18644e+02, -3.52571e+00, -1.34772e+03);
+    expectedForce[4] = Vec3(-1.61193e+02, 9.01528e+01, -7.12904e+01);
+    expectedForce[5] = Vec3(2.82630e+02, 2.78029e+01, -3.72864e+02);
+    expectedForce[6] = Vec3(-1.47454e+02, -2.14448e+02, -3.55789e+02);
+    expectedForce[7] = Vec3(-8.82195e+02, -7.39132e+01, 1.46202e+03);
+    for (int i = 0; i < 8; i++) {
+        ASSERT_EQUAL_VEC(expectedForce[i], state.getForces()[i], 1e-4);
+    }
+    ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), 1e-4);
+}
+
 void testErrorTolerance(NonbondedForce::NonbondedMethod method) {
     // Create a cloud of random point charges.
 
@@ -261,6 +312,7 @@ int main(int argc, char* argv[]) {
         testEwaldPME(false);
         testEwaldPME(true);
 //        testEwald2Ions();
+        testTriclinic();
         testErrorTolerance(NonbondedForce::Ewald);
         testErrorTolerance(NonbondedForce::PME);
     }

platforms/cuda/src/CudaKernels.cpp

@@ -1377,8 +1377,8 @@ public:
     PmePreComputation(CudaContext& cu, Kernel& pme, CalcPmeReciprocalForceKernel::IO& io) : cu(cu), pme(pme), io(io) {
     }
     void computeForceAndEnergy(bool includeForces, bool includeEnergy, int groups) {
-        Vec3 boxSize(cu.getPeriodicBoxSize().x, cu.getPeriodicBoxSize().y, cu.getPeriodicBoxSize().z);
-        pme.getAs<CalcPmeReciprocalForceKernel>().beginComputation(io, boxSize, includeEnergy);
+        Vec3 boxVectors[3] = {Vec3(cu.getPeriodicBoxSize().x, 0, 0), Vec3(0, cu.getPeriodicBoxSize().y, 0), Vec3(0, 0, cu.getPeriodicBoxSize().z)};
+        pme.getAs<CalcPmeReciprocalForceKernel>().beginComputation(io, boxVectors, includeEnergy);
     }
 private:
     CudaContext& cu;

platforms/opencl/src/OpenCLKernels.cpp

@@ -1363,8 +1363,8 @@ public:
     PmePreComputation(OpenCLContext& cl, Kernel& pme, CalcPmeReciprocalForceKernel::IO& io) : cl(cl), pme(pme), io(io) {
     }
     void computeForceAndEnergy(bool includeForces, bool includeEnergy, int groups) {
-        Vec3 boxSize(cl.getPeriodicBoxSize().x, cl.getPeriodicBoxSize().y, cl.getPeriodicBoxSize().z);
-        pme.getAs<CalcPmeReciprocalForceKernel>().beginComputation(io, boxSize, includeEnergy);
+        Vec3 boxVectors[3] = {Vec3(cl.getPeriodicBoxSize().x, 0, 0), Vec3(0, cl.getPeriodicBoxSize().y, 0), Vec3(0, 0, cl.getPeriodicBoxSize().z)};
+        pme.getAs<CalcPmeReciprocalForceKernel>().beginComputation(io, boxVectors, includeEnergy);
     }
 private:
     OpenCLContext& cl;

platforms/reference/src/SimTKReference/ReferencePME.cpp

@@ -256,7 +256,8 @@ pme_update_grid_index_and_fraction(pme_t    pme,
             coord[d] -= floor(coord[d]*recipBoxVectors[d][d])*periodicBoxVectors[d][d];
         for(d=0;d<3;d++)
         {
-            t  = (coord[0]*recipBoxVectors[0][d]+coord[1]*recipBoxVectors[1][d]+coord[2]*recipBoxVectors[2][d])*pme->ngrid[d];
+            t = coord[0]*recipBoxVectors[0][d]+coord[1]*recipBoxVectors[1][d]+coord[2]*recipBoxVectors[2][d];
+            t = (t-floor(t))*pme->ngrid[d];
             ti = (int) t;
 
             pme->particlefraction[i][d] = t - ti;

plugins/cpupme/src/CpuPmeKernels.cpp

@@ -47,10 +47,13 @@ static const int PME_ORDER = 5;
 bool CpuCalcPmeReciprocalForceKernel::hasInitializedThreads = false;
 int CpuCalcPmeReciprocalForceKernel::numThreads = 0;
 
-static void spreadCharge(int start, int end, float* posq, float* grid, int gridx, int gridy, int gridz, int numParticles, Vec3 periodicBoxSize) {
+static void spreadCharge(int start, int end, float* posq, float* grid, int gridx, int gridy, int gridz, int numParticles, Vec3* periodicBoxVectors, Vec3* recipBoxVectors) {
     float temp[4];
-    fvec4 boxSize((float) periodicBoxSize[0], (float) periodicBoxSize[1], (float) periodicBoxSize[2], 0);
-    fvec4 invBoxSize((float) (1/periodicBoxSize[0]), (float) (1/periodicBoxSize[1]), (float) (1/periodicBoxSize[2]), 0);
+    fvec4 boxSize((float) periodicBoxVectors[0][0], (float) periodicBoxVectors[1][1], (float) periodicBoxVectors[2][2], 0);
+    fvec4 invBoxSize((float) recipBoxVectors[0][0], (float) recipBoxVectors[1][1], (float) recipBoxVectors[2][2], 0);
+    fvec4 recipBoxVec0((float) recipBoxVectors[0][0], (float) recipBoxVectors[0][1], (float) recipBoxVectors[0][2], 0);
+    fvec4 recipBoxVec1((float) recipBoxVectors[1][0], (float) recipBoxVectors[1][1], (float) recipBoxVectors[1][2], 0);
+    fvec4 recipBoxVec2((float) recipBoxVectors[2][0], (float) recipBoxVectors[2][1], (float) recipBoxVectors[2][2], 0);
     fvec4 gridSize(gridx, gridy, gridz, 0);
     ivec4 gridSizeInt(gridx, gridy, gridz, 0);
     fvec4 one(1);
@@ -61,9 +64,10 @@ static void spreadCharge(int start, int end, float* posq, float* grid, int gridx
         // Find the position relative to the nearest grid point.
         
         fvec4 pos(&posq[4*i]);
-        fvec4 posFloor = floor(pos*invBoxSize);
-        fvec4 posInBox = pos-boxSize*posFloor;
-        fvec4 t = posInBox*invBoxSize*gridSize;
+        float posInBox[4];
+        (pos-boxSize*floor(pos*invBoxSize)).store(posInBox);
+        fvec4 t = posInBox[0]*recipBoxVec0 + posInBox[1]*recipBoxVec1 + posInBox[2]*recipBoxVec2;
+        t = (t-floor(t))*gridSize;
         ivec4 ti = t;
         fvec4 dr = t-ti;
         ivec4 gridIndex = ti-(gridSizeInt&ti==gridSizeInt);
@@ -142,29 +146,26 @@ static void spreadCharge(int start, int end, float* posq, float* grid, int gridx
     }
 }
 
-static void computeReciprocalEterm(int start, int end, int gridx, int gridy, int gridz, vector<float>& recipEterm, double alpha, vector<float>* bsplineModuli, Vec3 periodicBoxSize) {
+static void computeReciprocalEterm(int start, int end, int gridx, int gridy, int gridz, vector<float>& recipEterm, double alpha, vector<float>* bsplineModuli, Vec3* periodicBoxVectors, Vec3* recipBoxVectors) {
     const unsigned int zsize = gridz/2+1;
     const unsigned int yzsize = gridy*zsize;
-    const float scaleFactor = (float) (M_PI*periodicBoxSize[0]*periodicBoxSize[1]*periodicBoxSize[2]);
+    const float scaleFactor = (float) (M_PI*periodicBoxVectors[0][0]*periodicBoxVectors[1][1]*periodicBoxVectors[2][2]);
     const float recipExpFactor = (float) (M_PI*M_PI/(alpha*alpha));
-    const float invPeriodicBoxSizeX = (float) (1.0/periodicBoxSize[0]);
-    const float invPeriodicBoxSizeY = (float) (1.0/periodicBoxSize[1]);
-    const float invPeriodicBoxSizeZ = (float) (1.0/periodicBoxSize[2]);
 
     int firstz = (start == 0 ? 1 : 0);
     for (int kx = start; kx < end; kx++) {
         int mx = (kx < (gridx+1)/2) ? kx : kx-gridx;
-        float mhx = mx*invPeriodicBoxSizeX;
+        float mhx = mx*(float)recipBoxVectors[0][0];
         float bx = scaleFactor*bsplineModuli[0][kx];
         for (int ky = 0; ky < gridy; ky++) {
             int my = (ky < (gridy+1)/2) ? ky : ky-gridy;
-            float mhy = my*invPeriodicBoxSizeY;
+            float mhy = mx*(float)recipBoxVectors[1][0] + my*(float)recipBoxVectors[1][1];
             float mhx2y2 = mhx*mhx + mhy*mhy;
             float bxby = bx*bsplineModuli[1][ky];
             for (int kz = firstz; kz < zsize; kz++) {
                 int index = kx*yzsize + ky*zsize + kz;
                 int mz = (kz < (gridz+1)/2) ? kz : kz-gridz;
-                float mhz = mz*invPeriodicBoxSizeZ;
+                float mhz = mx*(float)recipBoxVectors[2][0] + my*(float)recipBoxVectors[2][1] + mz*(float)recipBoxVectors[2][2];
                 float bz = bsplineModuli[2][kz];
                 float m2 = mhx2y2 + mhz*mhz;
                 float denom = m2*bxby*bz;
@@ -175,29 +176,26 @@ static void computeReciprocalEterm(int start, int end, int gridx, int gridy, int
     }
 }
 
-static float reciprocalEnergy(int start, int end, fftwf_complex* grid, int gridx, int gridy, int gridz, double alpha, vector<float>* bsplineModuli, Vec3 periodicBoxSize) {
+static float reciprocalEnergy(int start, int end, fftwf_complex* grid, int gridx, int gridy, int gridz, double alpha, vector<float>* bsplineModuli, Vec3* periodicBoxVectors, Vec3* recipBoxVectors) {
     const unsigned int zsizeHalf = gridz/2+1;
     const unsigned int yzsizeHalf = gridy*zsizeHalf;
-    const float scaleFactor = (float) (M_PI*periodicBoxSize[0]*periodicBoxSize[1]*periodicBoxSize[2]);
+    const float scaleFactor = (float) (M_PI*periodicBoxVectors[0][0]*periodicBoxVectors[1][1]*periodicBoxVectors[2][2]);
     const float recipExpFactor = (float) (M_PI*M_PI/(alpha*alpha));
-    const float invPeriodicBoxSizeX = (float) (1.0/periodicBoxSize[0]);
-    const float invPeriodicBoxSizeY = (float) (1.0/periodicBoxSize[1]);
-    const float invPeriodicBoxSizeZ = (float) (1.0/periodicBoxSize[2]);
     float energy = 0.0f;
 
     int firstz = (start == 0 ? 1 : 0);
     for (int kx = start; kx < end; kx++) {
         int mx = (kx < (gridx+1)/2) ? kx : kx-gridx;
-        float mhx = mx*invPeriodicBoxSizeX;
+        float mhx = mx*(float)recipBoxVectors[0][0];
         float bx = scaleFactor*bsplineModuli[0][kx];
         for (int ky = 0; ky < gridy; ky++) {
             int my = (ky < (gridy+1)/2) ? ky : ky-gridy;
-            float mhy = my*invPeriodicBoxSizeY;
+            float mhy = mx*(float)recipBoxVectors[1][0] + my*(float)recipBoxVectors[1][1];
             float mhx2y2 = mhx*mhx + mhy*mhy;
             float bxby = bx*bsplineModuli[1][ky];
             for (int kz = firstz; kz < gridz; kz++) {
                 int mz = (kz < (gridz+1)/2) ? kz : kz-gridz;
-                float mhz = mz*invPeriodicBoxSizeZ;
+                float mhz = mx*(float)recipBoxVectors[2][0] + my*(float)recipBoxVectors[2][1] + mz*(float)recipBoxVectors[2][2];
                 float bz = bsplineModuli[2][kz];
                 float m2 = mhx2y2 + mhz*mhz;
                 float denom = m2*bxby*bz;
@@ -242,9 +240,12 @@ static void reciprocalConvolution(int start, int end, fftwf_complex* grid, int g
     }
 }
 
-static void interpolateForces(int start, int end, float* posq, float* force, float* grid, int gridx, int gridy, int gridz, int numParticles, Vec3 periodicBoxSize) {
-    fvec4 boxSize((float) periodicBoxSize[0], (float) periodicBoxSize[1], (float) periodicBoxSize[2], 0);
-    fvec4 invBoxSize((float) (1/periodicBoxSize[0]), (float) (1/periodicBoxSize[1]), (float) (1/periodicBoxSize[2]), 0);
+static void interpolateForces(int start, int end, float* posq, float* force, float* grid, int gridx, int gridy, int gridz, int numParticles, Vec3* periodicBoxVectors, Vec3* recipBoxVectors) {
+    fvec4 boxSize((float) periodicBoxVectors[0][0], (float) periodicBoxVectors[1][1], (float) periodicBoxVectors[2][2], 0);
+    fvec4 invBoxSize((float) recipBoxVectors[0][0], (float) recipBoxVectors[1][1], (float) recipBoxVectors[2][2], 0);
+    fvec4 recipBoxVec0((float) recipBoxVectors[0][0], (float) recipBoxVectors[0][1], (float) recipBoxVectors[0][2], 0);
+    fvec4 recipBoxVec1((float) recipBoxVectors[1][0], (float) recipBoxVectors[1][1], (float) recipBoxVectors[1][2], 0);
+    fvec4 recipBoxVec2((float) recipBoxVectors[2][0], (float) recipBoxVectors[2][1], (float) recipBoxVectors[2][2], 0);
     fvec4 gridSize(gridx, gridy, gridz, 0);
     ivec4 gridSizeInt(gridx, gridy, gridz, 0);
     fvec4 one(1);
@@ -254,9 +255,10 @@ static void interpolateForces(int start, int end, float* posq, float* force, flo
         // Find the position relative to the nearest grid point.
         
         fvec4 pos(&posq[4*i]);
-        fvec4 posFloor = floor(pos*invBoxSize);
-        fvec4 posInBox = pos-boxSize*posFloor;
-        fvec4 t = posInBox*invBoxSize*gridSize;
+        float posInBox[4];
+        (pos-boxSize*floor(pos*invBoxSize)).store(posInBox);
+        fvec4 t = posInBox[0]*recipBoxVec0 + posInBox[1]*recipBoxVec1 + posInBox[2]*recipBoxVec2;
+        t = (t-floor(t))*gridSize;
         ivec4 ti = t;
         fvec4 dr = t-ti;
         ivec4 gridIndex = ti-(gridSizeInt&ti==gridSizeInt);
@@ -321,8 +323,12 @@ static void interpolateForces(int start, int end, float* posq, float* force, flo
                 }
             }
         }
-        f = invBoxSize*gridSize*f*(-epsilonFactor*posq[4*i+3]);
-        f.store(&force[4*i]);
+        f *= -epsilonFactor*posq[4*i+3];
+        float fc[4];
+        f.store(fc);
+        force[4*i+0] = fc[0]*gridx*(float)recipBoxVectors[0][0];
+        force[4*i+1] = fc[0]*gridx*(float)recipBoxVectors[1][0]+fc[1]*gridy*(float)recipBoxVectors[1][1];
+        force[4*i+2] = fc[0]*gridx*(float)recipBoxVectors[2][0]+fc[1]*gridy*(float)recipBoxVectors[2][1]+fc[2]*gridz*(float)recipBoxVectors[2][2];
     }
 }
 
@@ -476,7 +482,7 @@ void CpuCalcPmeReciprocalForceKernel::runThread(int index) {
             advanceThreads(); // Signal threads to perform charge spreading.
             advanceThreads(); // Signal threads to sum the charge grids.
             fftwf_execute_dft_r2c(forwardFFT, realGrid, complexGrid);
-            if (lastBoxSize != periodicBoxSize)
+            if (lastBoxVectors[0] != periodicBoxVectors[0] || lastBoxVectors[1] != periodicBoxVectors[1] || lastBoxVectors[2] != periodicBoxVectors[2])
                 advanceThreads(); // Signal threads to compute the reciprocal scale factors.
             if (includeEnergy)
                 advanceThreads(); // Signal threads to compute energy.
@@ -484,7 +490,9 @@ void CpuCalcPmeReciprocalForceKernel::runThread(int index) {
             fftwf_execute_dft_c2r(backwardFFT, complexGrid, realGrid);
             advanceThreads(); // Signal threads to interpolate forces.
             isFinished = true;
-            lastBoxSize = periodicBoxSize;
+            lastBoxVectors[0] = periodicBoxVectors[0];
+            lastBoxVectors[1] = periodicBoxVectors[1];
+            lastBoxVectors[2] = periodicBoxVectors[2];
             pthread_cond_signal(&mainThreadEndCondition);
         }
         pthread_mutex_unlock(&lock);
@@ -503,7 +511,7 @@ void CpuCalcPmeReciprocalForceKernel::runThread(int index) {
             threadWait();
             if (isDeleted)
                 break;
-            spreadCharge(particleStart, particleEnd, posq, threadData[index]->tempGrid, gridx, gridy, gridz, numParticles, periodicBoxSize);
+            spreadCharge(particleStart, particleEnd, posq, threadData[index]->tempGrid, gridx, gridy, gridz, numParticles, periodicBoxVectors, recipBoxVectors);
             threadWait();
             int numGrids = threadData.size();
             for (int i = gridStart; i < gridEnd; i += 4) {
@@ -513,12 +521,12 @@ void CpuCalcPmeReciprocalForceKernel::runThread(int index) {
                 sum.store(&realGrid[i]);
             }
             threadWait();
-            if (lastBoxSize != periodicBoxSize) {
-                computeReciprocalEterm(gridxStart, gridxEnd, gridx, gridy, gridz, recipEterm, alpha, bsplineModuli, periodicBoxSize);
+            if (lastBoxVectors[0] != periodicBoxVectors[0] || lastBoxVectors[1] != periodicBoxVectors[1] || lastBoxVectors[2] != periodicBoxVectors[2]) {
+                computeReciprocalEterm(gridxStart, gridxEnd, gridx, gridy, gridz, recipEterm, alpha, bsplineModuli, periodicBoxVectors, recipBoxVectors);
                 threadWait();
             }
             if (includeEnergy) {
-                double threadEnergy = reciprocalEnergy(gridxStart, gridxEnd, complexGrid, gridx, gridy, gridz, alpha, bsplineModuli, periodicBoxSize);
+                double threadEnergy = reciprocalEnergy(gridxStart, gridxEnd, complexGrid, gridx, gridy, gridz, alpha, bsplineModuli, periodicBoxVectors, recipBoxVectors);
                 pthread_mutex_lock(&lock);
                 energy += threadEnergy;
                 pthread_mutex_unlock(&lock);
@@ -526,7 +534,7 @@ void CpuCalcPmeReciprocalForceKernel::runThread(int index) {
             }
             reciprocalConvolution(gridxStart, gridxEnd, complexGrid, gridx, gridy, gridz, recipEterm);
             threadWait();
-            interpolateForces(particleStart, particleEnd, posq, &force[0], realGrid, gridx, gridy, gridz, numParticles, periodicBoxSize);
+            interpolateForces(particleStart, particleEnd, posq, &force[0], realGrid, gridx, gridy, gridz, numParticles, periodicBoxVectors, recipBoxVectors);
         }
     }
 }
@@ -547,11 +555,24 @@ void CpuCalcPmeReciprocalForceKernel::advanceThreads() {
     }
 }
 
-void CpuCalcPmeReciprocalForceKernel::beginComputation(IO& io, Vec3 periodicBoxSize, bool includeEnergy) {
+void CpuCalcPmeReciprocalForceKernel::beginComputation(IO& io, const Vec3* periodicBoxVectors, bool includeEnergy) {
     this->io = &io;
-    this->periodicBoxSize = periodicBoxSize;
+    this->periodicBoxVectors[0] = periodicBoxVectors[0];
+    this->periodicBoxVectors[1] = periodicBoxVectors[1];
+    this->periodicBoxVectors[2] = periodicBoxVectors[2];
     this->includeEnergy = includeEnergy;
     energy = 0.0;
+
+    // Invert the box vectors.
+
+    double determinant = periodicBoxVectors[0][0]*periodicBoxVectors[1][1]*periodicBoxVectors[2][2];
+    double scale = 1.0/determinant;
+    recipBoxVectors[0] = Vec3(periodicBoxVectors[1][1]*periodicBoxVectors[2][2], 0, 0)*scale;
+    recipBoxVectors[1] = Vec3(-periodicBoxVectors[1][0]*periodicBoxVectors[2][2], periodicBoxVectors[0][0]*periodicBoxVectors[2][2], 0)*scale;
+    recipBoxVectors[2] = Vec3(periodicBoxVectors[1][0]*periodicBoxVectors[2][1]-periodicBoxVectors[1][1]*periodicBoxVectors[2][0], -periodicBoxVectors[0][0]*periodicBoxVectors[2][1], periodicBoxVectors[0][0]*periodicBoxVectors[1][1])*scale;
+
+    // Do the calculation.
+
     pthread_mutex_lock(&lock);
     isFinished = false;
     pthread_cond_signal(&mainThreadStartCondition);

plugins/cpupme/src/CpuPmeKernels.h

@@ -66,11 +66,11 @@ public:
     /**
      * Begin computing the force and energy.
      * 
-     * @param io               an object that coordinates data transfer
-     * @param periodicBoxSize  the size of the periodic box (measured in nm)
-     * @param includeEnergy    true if potential energy should be computed
+     * @param io                  an object that coordinates data transfer
+     * @param periodicBoxVectors  the vectors defining the periodic box (measured in nm)
+     * @param includeEnergy       true if potential energy should be computed
      */
-    void beginComputation(IO& io, Vec3 periodicBoxSize, bool includeEnergy);
+    void beginComputation(IO& io, const Vec3* periodicBoxVectors, bool includeEnergy);
     /**
      * Finish computing the force and energy.
      * 
@@ -107,7 +107,7 @@ private:
     std::vector<float> force;
     std::vector<float> bsplineModuli[3];
     std::vector<float> recipEterm;
-    Vec3 lastBoxSize;
+    Vec3 lastBoxVectors[3];
     float* realGrid;
     fftwf_complex* complexGrid;
     fftwf_plan forwardFFT, backwardFFT;
@@ -122,7 +122,7 @@ private:
     IO* io;
     float energy;
     float* posq;
-    Vec3 periodicBoxSize;
+    Vec3 periodicBoxVectors[3], recipBoxVectors[3];
     bool includeEnergy;
 };
 

plugins/cpupme/tests/TestCpuPme.cpp

@@ -61,14 +61,25 @@ public:
     }
 };
 
-void testPME() {
+void testPME(bool triclinic) {
     // Create a cloud of random point charges.
 
     const int numParticles = 51;
     const double boxWidth = 5.0;
     const double cutoff = 1.0;
+    Vec3 boxVectors[3];
+    if (triclinic) {
+        boxVectors[0] = Vec3(boxWidth, 0, 0);
+        boxVectors[1] = Vec3(0.2*boxWidth, boxWidth, 0);
+        boxVectors[2] = Vec3(-0.3*boxWidth, -0.1*boxWidth, boxWidth);
+    }
+    else {
+        boxVectors[0] = Vec3(boxWidth, 0, 0);
+        boxVectors[1] = Vec3(0, boxWidth, 0);
+        boxVectors[2] = Vec3(0, 0, boxWidth);
+    }
     System system;
-    system.setDefaultPeriodicBoxVectors(Vec3(boxWidth, 0, 0), Vec3(0, boxWidth, 0), Vec3(0, 0, boxWidth));
+    system.setDefaultPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
     NonbondedForce* force = new NonbondedForce();
     system.addForce(force);
     vector<Vec3> positions(numParticles);
@@ -112,7 +123,7 @@ void testPME() {
     }
     double ewaldSelfEnergy = -ONE_4PI_EPS0*alpha*sumSquaredCharges/sqrt(M_PI);
     pme.initialize(gridx, gridy, gridz, numParticles, alpha);
-    pme.beginComputation(io, Vec3(boxWidth, boxWidth, boxWidth), true);
+    pme.beginComputation(io, boxVectors, true);
     double energy = pme.finishComputation(io);
     
     // See if they match.
@@ -128,7 +139,8 @@ int main(int argc, char* argv[]) {
             cout << "CPU is not supported.  Exiting." << endl;
             return 0;
         }
-        testPME();
+        testPME(false);
+        testPME(true);
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;