Fixes to dispersion PME

openmmapi/include/openmm/NonbondedForce.h

@@ -101,18 +101,18 @@ public:
          */
         CutoffPeriodic = 2,
         /**
-         * Periodic boundary conditions are used, and Ewald summation is used to compute the interaction of each particle
+         * Periodic boundary conditions are used, and Ewald summation is used to compute the Coulomb interaction of each particle
          * with all periodic copies of every other particle.
          */
         Ewald = 3,
         /**
-         * Periodic boundary conditions are used, and Particle-Mesh Ewald (PME) summation is used to compute the interaction of each particle
+         * Periodic boundary conditions are used, and Particle-Mesh Ewald (PME) summation is used to compute the Coulomb interaction of each particle
          * with all periodic copies of every other particle.
          */
         PME = 4,
         /**
          * Periodic boundary conditions are used, and Particle-Mesh Ewald (PME) summation is used to compute the interaction of each particle
-         * with all periodic copies of every other particle for both electrostatics and dispersion.  No switching is used for either interaction.
+         * with all periodic copies of every other particle for both Coulomb and Lennard-Jones.  No switching is used for either interaction.
          */
         LJPME = 5
     };

openmmapi/src/NonbondedForce.cpp

@@ -48,7 +48,8 @@ using std::stringstream;
 using std::vector;
 
 NonbondedForce::NonbondedForce() : nonbondedMethod(NoCutoff), cutoffDistance(1.0), switchingDistance(-1.0), rfDielectric(78.3),
-        ewaldErrorTol(5e-4), alpha(0.0), useSwitchingFunction(false), useDispersionCorrection(true), recipForceGroup(-1), nx(0), ny(0), nz(0) {
+        ewaldErrorTol(5e-4), alpha(0.0), dalpha(0.0), useSwitchingFunction(false), useDispersionCorrection(true), recipForceGroup(-1),
+        nx(0), ny(0), nz(0), dnx(0), dny(0), dnz(0) {
 }
 
 NonbondedForce::NonbondedMethod NonbondedForce::getNonbondedMethod() const {

openmmapi/src/NonbondedForceImpl.cpp

@@ -90,9 +90,7 @@ void NonbondedForceImpl::initialize(ContextImpl& context) {
         exceptions[particle1].insert(particle2);
         exceptions[particle2].insert(particle1);
     }
-    if (owner.getNonbondedMethod() == NonbondedForce::CutoffPeriodic ||
-            owner.getNonbondedMethod() == NonbondedForce::Ewald ||
-            owner.getNonbondedMethod() == NonbondedForce::PME) {
+    if (owner.getNonbondedMethod() != NonbondedForce::NoCutoff) {
         Vec3 boxVectors[3];
         system.getDefaultPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
         double cutoff = owner.getCutoffDistance();
@@ -161,12 +159,19 @@ void NonbondedForceImpl::calcPMEParameters(const System& system, const Nonbonded
         system.getDefaultPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
         double tol = force.getEwaldErrorTolerance();
         alpha = (1.0/force.getCutoffDistance())*std::sqrt(-log(2.0*tol));
-        xsize = (int) ceil(2*alpha*boxVectors[0][0]/(3*pow(tol, 0.2)));
-        ysize = (int) ceil(2*alpha*boxVectors[1][1]/(3*pow(tol, 0.2)));
-        zsize = (int) ceil(2*alpha*boxVectors[2][2]/(3*pow(tol, 0.2)));
-        xsize = max(xsize, 5);
-        ysize = max(ysize, 5);
-        zsize = max(zsize, 5);
+        if (lj) {
+            xsize = (int) ceil(alpha*boxVectors[0][0]/(3*pow(tol, 0.05)));
+            ysize = (int) ceil(alpha*boxVectors[1][1]/(3*pow(tol, 0.05)));
+            zsize = (int) ceil(alpha*boxVectors[2][2]/(3*pow(tol, 0.05)));
+        }
+        else {
+            xsize = (int) ceil(2*alpha*boxVectors[0][0]/(3*pow(tol, 0.2)));
+            ysize = (int) ceil(2*alpha*boxVectors[1][1]/(3*pow(tol, 0.2)));
+            zsize = (int) ceil(2*alpha*boxVectors[2][2]/(3*pow(tol, 0.2)));
+        }
+        xsize = max(xsize, 6);
+        ysize = max(ysize, 6);
+        zsize = max(zsize, 6);
     }
 }
 

platforms/cpu/tests/TestCpuDispersionPME.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2017 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "CpuTests.h"
+#include "TestDispersionPME.h"
+
+void runPlatformTests() {
+}
+

platforms/cuda/src/kernels/coulombLennardJones.cu

@@ -30,7 +30,7 @@
             const real dar6 = dar4*dar2;
             const real invR2 = invR*invR;
             const real expDar2 = EXP(-dar2);
-            const real2 sigExpProd = sigmaEpsilon1*sigmaEpsilon2;
+            const float2 sigExpProd = sigmaEpsilon1*sigmaEpsilon2;
             const real c6 = 64*sigExpProd.x*sigExpProd.x*sigExpProd.x*sigExpProd.y;
             const real coef = invR2*invR2*invR2*c6;
             const real eprefac = 1.0f + dar2 + 0.5f*dar4;

platforms/cuda/src/kernels/pme.cu

@@ -23,7 +23,7 @@ extern "C" __global__ void gridSpreadCharge(const real4* __restrict__ posq, real
         real4 periodicBoxSize, real4 invPeriodicBoxSize, real4 periodicBoxVecX, real4 periodicBoxVecY, real4 periodicBoxVecZ,
         real3 recipBoxVecX, real3 recipBoxVecY, real3 recipBoxVecZ, const int2* __restrict__ pmeAtomGridIndex
 #ifdef USE_LJPME
-        , const real2* __restrict__ sigmaEpsilon
+        , const float2* __restrict__ sigmaEpsilon
 #endif
         ) {
     real3 data[PME_ORDER];
@@ -68,7 +68,7 @@ extern "C" __global__ void gridSpreadCharge(const real4* __restrict__ posq, real
         // Spread the charge from this atom onto each grid point.
         
 #ifdef USE_LJPME
-        const real2 sigEps = sigmaEpsilon[atom];
+        const float2 sigEps = sigmaEpsilon[atom];
         const real charge = 8*sigEps.x*sigEps.x*sigEps.x*sigEps.y;
 #else
         const real charge = pos.w;
@@ -250,7 +250,7 @@ void gridInterpolateForce(const real4* __restrict__ posq, unsigned long long* __
         real4 periodicBoxSize, real4 invPeriodicBoxSize, real4 periodicBoxVecX, real4 periodicBoxVecY, real4 periodicBoxVecZ,
         real3 recipBoxVecX, real3 recipBoxVecY, real3 recipBoxVecZ, const int2* __restrict__ pmeAtomGridIndex
 #ifdef USE_LJPME
-        , const real2* __restrict__ sigmaEpsilon
+        , const float2* __restrict__ sigmaEpsilon
 #endif
         ) {
     real3 data[PME_ORDER];
@@ -325,7 +325,7 @@ void gridInterpolateForce(const real4* __restrict__ posq, unsigned long long* __
             }
         }
 #ifdef USE_LJPME
-        const real2 sigEps = sigmaEpsilon[atom];
+        const float2 sigEps = sigmaEpsilon[atom];
         real q = 8*sigEps.x*sigEps.x*sigEps.x*sigEps.y;
 #else
         real q = pos.w*EPSILON_FACTOR;

platforms/cuda/tests/TestCudaDispersionPME.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2017 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "CudaTests.h"
+#include "TestDispersionPME.h"
+
+void runPlatformTests() {
+}

platforms/opencl/tests/TestOpenCLDispersionPME.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2017 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "OpenCLTests.h"
+#include "TestDispersionPME.h"
+
+void runPlatformTests() {
+}

platforms/reference/tests/TestReferenceDispersionPME.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2017 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "ReferenceTests.h"
+#include "TestDispersionPME.h"
+
+void runPlatformTests() {
+}

tests/TestDispersionPME.h

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008-2017 Stanford University and the Authors.      *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "openmm/NonbondedForce.h"
+#include "openmm/System.h"
+#include "openmm/Units.h"
+#include "openmm/VerletIntegrator.h"
+#include "openmm/internal/NonbondedForceImpl.h"
+#include "sfmt/SFMT.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+void testErrorTolerance() {
+    // Create a cloud of random point charges.
+
+    const int numParticles = 200;
+    const double boxWidth = 5.0;
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(boxWidth, 0, 0), Vec3(0, boxWidth, 0), Vec3(0, 0, boxWidth));
+    NonbondedForce* force = new NonbondedForce();
+    system.addForce(force);
+    vector<Vec3> positions(numParticles);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+
+    for (int i = 0; i < numParticles; i++) {
+        system.addParticle(1.0);
+        force->addParticle(0.0, 0.1+0.2*genrand_real2(sfmt), genrand_real2(sfmt));
+        while (true) {
+            Vec3 pos = Vec3(boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt));
+            double minDist = boxWidth;
+            for (int j = 0; j < i; j++) {
+                Vec3 delta = pos-positions[j];
+                minDist = min(minDist, sqrt(delta.dot(delta)));
+            }
+            if (minDist > 0.1) {
+                positions[i] = pos;
+                break;
+            }
+        }
+    }
+    force->setNonbondedMethod(NonbondedForce::LJPME);
+
+    // For various values of the cutoff and error tolerance, see if the actual error is reasonable.
+
+    for (double cutoff = 0.8; cutoff < boxWidth/2; cutoff *= 1.3) {
+        force->setCutoffDistance(cutoff);
+        vector<Vec3> refForces;
+        double norm = 0.0;
+        for (double tol = 5e-5; tol < 1e-3; tol *= 2.0) {
+            force->setEwaldErrorTolerance(tol);
+            VerletIntegrator integrator(0.01);
+            Context context(system, integrator, platform);
+            context.setPositions(positions);
+            State state = context.getState(State::Forces);
+            if (refForces.size() == 0) {
+                refForces = state.getForces();
+                for (int i = 0; i < numParticles; i++)
+                    norm += refForces[i].dot(refForces[i]);
+                norm = sqrt(norm);
+            }
+            else {
+                double diff = 0.0;
+                for (int i = 0; i < numParticles; i++) {
+                    Vec3 delta = refForces[i]-state.getForces()[i];
+                    diff += delta.dot(delta);
+                }
+                diff = sqrt(diff)/norm;
+                ASSERT(diff < 2*tol);
+            }
+
+            // See if the PME parameters were calculated correctly.
+
+            double expectedAlpha, actualAlpha;
+            int expectedSize[3], actualSize[3];
+            NonbondedForceImpl::calcPMEParameters(system, *force, expectedAlpha, expectedSize[0], expectedSize[1], expectedSize[2], true);
+            force->getLJPMEParametersInContext(context, actualAlpha, actualSize[0], actualSize[1], actualSize[2]);
+            ASSERT_EQUAL_TOL(expectedAlpha, actualAlpha, 1e-5);
+            for (int i = 0; i < 3; i++) {
+                ASSERT(actualSize[i] >= expectedSize[i]);
+                ASSERT(actualSize[i] < expectedSize[i]+10);
+            }
+        }
+    }
+}
+
+void testPMEParameters() {
+    // Create a cloud of random point charges.
+
+    const int numParticles = 51;
+    const double boxWidth = 4.7;
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(boxWidth, 0, 0), Vec3(0, boxWidth, 0), Vec3(0, 0, boxWidth));
+    NonbondedForce* force = new NonbondedForce();
+    system.addForce(force);
+    vector<Vec3> positions(numParticles);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+
+    for (int i = 0; i < numParticles; i++) {
+        system.addParticle(1.0);
+        force->addParticle(0.0, 0.1+0.2*genrand_real2(sfmt), genrand_real2(sfmt));
+        positions[i] = Vec3(boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt));
+    }
+    force->setNonbondedMethod(NonbondedForce::LJPME);
+    
+    // Compute the energy with an error tolerance of 1e-3.
+
+    force->setEwaldErrorTolerance(1e-3);
+    VerletIntegrator integrator1(0.01);
+    Context context1(system, integrator1, platform);
+    context1.setPositions(positions);
+    double energy1 = context1.getState(State::Energy).getPotentialEnergy();
+    double alpha;
+    int gridx, gridy, gridz;
+    force->getLJPMEParametersInContext(context1, alpha, gridx, gridy, gridz);
+    
+    // Try again with an error tolerance of 1e-4.
+
+    force->setEwaldErrorTolerance(1e-4);
+    VerletIntegrator integrator2(0.01);
+    Context context2(system, integrator2, platform);
+    context2.setPositions(positions);
+    double energy2 = context2.getState(State::Energy).getPotentialEnergy();
+    
+    // Now explicitly set the parameters.  These should match the values that were
+    // used for tolerance 1e-3.
+
+    force->setLJPMEParameters(alpha, gridx, gridy, gridz);
+    VerletIntegrator integrator3(0.01);
+    Context context3(system, integrator3, platform);
+    context3.setPositions(positions);
+    double energy3 = context3.getState(State::Energy).getPotentialEnergy();
+    ASSERT_EQUAL_TOL(energy1, energy3, 1e-5);
+    ASSERT(fabs((energy1-energy2)/energy1) > 1e-5);
+    force->getLJPMEParametersInContext(context2, alpha, gridx, gridy, gridz);
+}
+
+void make_waterbox(int natoms, double boxEdgeLength, NonbondedForce *forceField,  vector<Vec3> &positions, vector<double>& eps, vector<double>& sig,
+                   vector<pair<int, int> >& bonds, System &system, bool do_electrostatics) {
+    const int RESSIZE = 3;
+    const double masses[RESSIZE]    = {     8.0,    1.0,    1.0 };
+    const double charges[RESSIZE]   = {  -0.834,  0.417,  0.417 };
+    // Values from the CHARMM force field, in AKMA units
+    const double epsilons[RESSIZE]  = { -0.1521, -0.046, -0.046 };
+    const double halfrmins[RESSIZE] = {  1.7682, 0.2245, 0.2245 };
+    positions.clear();
+    if (natoms == 6) {
+        const double coords[6][3] = {
+            {  2.000000, 2.000000, 2.000000},
+            {  2.500000, 2.000000, 3.000000},
+            {  1.500000, 2.000000, 3.000000},
+            {  0.000000, 0.000000, 0.000000},
+            {  0.500000, 0.000000, 1.000000},
+            { -0.500000, 0.000000, 1.000000}
+        };
+        for (int atom = 0; atom < natoms; ++atom)
+            positions.push_back(Vec3(coords[atom][0], coords[atom][1], coords[atom][2])*OpenMM::NmPerAngstrom);
+    }
+    else if (natoms == 375) {
+        const double coords[375][3] = {
+            { -6.22, -6.25, -6.24 },
+            { -5.32, -6.03, -6.00 },
+            { -6.75, -5.56, -5.84 },
+            { -3.04, -6.23, -6.19 },
+            { -3.52, -5.55, -5.71 },
+            { -3.59, -6.43, -6.94 },
+            {  0.02, -6.23, -6.14 },
+            { -0.87, -5.97, -6.37 },
+            {  0.53, -6.03, -6.93 },
+            {  3.10, -6.20, -6.27 },
+            {  3.87, -6.35, -5.72 },
+            {  2.37, -6.11, -5.64 },
+            {  6.18, -6.14, -6.20 },
+            {  6.46, -6.66, -5.44 },
+            {  6.26, -6.74, -6.94 },
+            { -6.21, -3.15, -6.24 },
+            { -6.23, -3.07, -5.28 },
+            { -6.02, -2.26, -6.55 },
+            { -3.14, -3.07, -6.16 },
+            { -3.38, -3.63, -6.90 },
+            { -2.18, -3.05, -6.17 },
+            { -0.00, -3.16, -6.23 },
+            { -0.03, -2.30, -6.67 },
+            {  0.05, -2.95, -5.29 },
+            {  3.08, -3.11, -6.14 },
+            {  2.65, -2.55, -6.79 },
+            {  3.80, -3.53, -6.62 },
+            {  6.16, -3.14, -6.16 },
+            {  7.04, -3.32, -6.51 },
+            {  5.95, -2.27, -6.51 },
+            { -6.20, -0.04, -6.15 },
+            { -5.43,  0.32, -6.59 },
+            { -6.95,  0.33, -6.62 },
+            { -3.10, -0.06, -6.19 },
+            { -3.75,  0.42, -6.69 },
+            { -2.46,  0.60, -5.93 },
+            {  0.05, -0.01, -6.17 },
+            { -0.10,  0.02, -7.12 },
+            { -0.79,  0.16, -5.77 },
+            {  3.03,  0.00, -6.19 },
+            {  3.54,  0.08, -7.01 },
+            {  3.69, -0.22, -5.53 },
+            {  6.17,  0.05, -6.19 },
+            {  5.78, -0.73, -6.57 },
+            {  7.09, -0.17, -6.04 },
+            { -6.20,  3.15, -6.25 },
+            { -6.59,  3.18, -5.37 },
+            { -5.87,  2.25, -6.33 },
+            { -3.09,  3.04, -6.17 },
+            { -3.88,  3.58, -6.26 },
+            { -2.41,  3.54, -6.63 },
+            {  0.00,  3.06, -6.26 },
+            { -0.71,  3.64, -6.00 },
+            {  0.65,  3.15, -5.55 },
+            {  3.14,  3.06, -6.23 },
+            {  3.11,  3.31, -5.30 },
+            {  2.38,  3.49, -6.63 },
+            {  6.19,  3.14, -6.25 },
+            {  6.82,  3.25, -5.54 },
+            {  5.76,  2.30, -6.07 },
+            { -6.22,  6.26, -6.19 },
+            { -6.22,  5.74, -7.00 },
+            { -5.89,  5.67, -5.52 },
+            { -3.04,  6.24, -6.20 },
+            { -3.08,  5.28, -6.17 },
+            { -3.96,  6.52, -6.25 },
+            { -0.05,  6.21, -6.16 },
+            {  0.82,  6.58, -6.06 },
+            {  0.01,  5.64, -6.93 },
+            {  3.10,  6.25, -6.15 },
+            {  3.64,  5.47, -6.31 },
+            {  2.46,  6.24, -6.87 },
+            {  6.22,  6.20, -6.27 },
+            {  5.37,  6.42, -5.88 },
+            {  6.80,  6.07, -5.51 },
+            { -6.19, -6.15, -3.13 },
+            { -6.37, -7.01, -3.51 },
+            { -6.25, -6.29, -2.18 },
+            { -3.10, -6.27, -3.11 },
+            { -2.29, -5.77, -2.99 },
+            { -3.80, -5.62, -2.98 },
+            { -0.03, -6.18, -3.15 },
+            { -0.07, -7.05, -2.75 },
+            {  0.68, -5.74, -2.70 },
+            {  3.10, -6.14, -3.07 },
+            {  2.35, -6.72, -3.23 },
+            {  3.86, -6.65, -3.37 },
+            {  6.22, -6.20, -3.16 },
+            {  6.82, -6.36, -2.43 },
+            {  5.35, -6.13, -2.75 },
+            { -6.26, -3.13, -3.12 },
+            { -6.16, -2.27, -2.70 },
+            { -5.36, -3.47, -3.18 },
+            { -3.11, -3.05, -3.14 },
+            { -3.31, -3.96, -3.34 },
+            { -2.77, -3.06, -2.24 },
+            {  0.00, -3.13, -3.16 },
+            {  0.48, -2.37, -2.81 },
+            { -0.57, -3.40, -2.44 },
+            {  3.09, -3.09, -3.16 },
+            {  2.41, -3.19, -2.49 },
+            {  3.91, -3.07, -2.67 },
+            {  6.19, -3.04, -3.08 },
+            {  5.64, -3.61, -3.61 },
+            {  6.93, -3.58, -2.82 },
+            { -6.18, -0.00, -3.04 },
+            { -6.00, -0.59, -3.78 },
+            { -6.79,  0.64, -3.39 },
+            { -3.05, -0.03, -3.07 },
+            { -2.95,  0.80, -3.52 },
+            { -4.00, -0.20, -3.07 },
+            { -0.03,  0.03, -3.06 },
+            { -0.33, -0.37, -3.87 },
+            {  0.89, -0.21, -2.99 },
+            {  3.13, -0.05, -3.10 },
+            {  3.44,  0.81, -3.34 },
+            {  2.21,  0.07, -2.86 },
+            {  6.20, -0.05, -3.13 },
+            {  6.89,  0.60, -3.20 },
+            {  5.58,  0.30, -2.49 },
+            { -6.23,  3.09, -3.16 },
+            { -5.62,  3.79, -2.94 },
+            { -6.33,  2.60, -2.33 },
+            { -3.10,  3.08, -3.04 },
+            { -3.84,  3.47, -3.51 },
+            { -2.40,  3.01, -3.69 },
+            {  0.01,  3.04, -3.11 },
+            { -0.56,  3.59, -3.64 },
+            {  0.28,  3.60, -2.38 },
+            {  3.04,  3.11, -3.09 },
+            {  3.49,  2.30, -2.87 },
+            {  3.70,  3.66, -3.51 },
+            {  6.15,  3.14, -3.11 },
+            {  6.52,  2.52, -3.74 },
+            {  6.72,  3.06, -2.34 },
+            { -6.22,  6.15, -3.13 },
+            { -5.49,  6.21, -2.51 },
+            { -6.56,  7.04, -3.18 },
+            { -3.11,  6.24, -3.05 },
+            { -3.76,  5.83, -3.62 },
+            { -2.26,  5.92, -3.37 },
+            {  0.03,  6.25, -3.07 },
+            {  0.34,  5.63, -3.73 },
+            { -0.87,  6.00, -2.91 },
+            {  3.07,  6.15, -3.08 },
+            {  3.29,  6.92, -2.56 },
+            {  3.39,  6.35, -3.96 },
+            {  6.22,  6.14, -3.12 },
+            {  5.79,  6.38, -2.29 },
+            {  6.25,  6.96, -3.62 },
+            { -6.21, -6.20, -0.06 },
+            { -5.79, -6.87,  0.48 },
+            { -6.43, -5.50,  0.54 },
+            { -3.16, -6.21, -0.02 },
+            { -2.50, -6.87,  0.20 },
+            { -2.77, -5.37,  0.23 },
+            { -0.00, -6.14, -0.00 },
+            {  0.68, -6.72, -0.33 },
+            { -0.64, -6.73,  0.38 },
+            {  3.03, -6.20, -0.01 },
+            {  3.77, -6.56, -0.51 },
+            {  3.43, -5.85,  0.78 },
+            {  6.25, -6.16, -0.00 },
+            {  5.36, -6.09, -0.36 },
+            {  6.24, -6.97,  0.49 },
+            { -6.24, -3.05, -0.01 },
+            { -6.35, -3.64,  0.73 },
+            { -5.42, -3.33, -0.42 },
+            { -3.09, -3.06,  0.05 },
+            { -2.44, -3.62, -0.38 },
+            { -3.90, -3.21, -0.43 },
+            {  0.05, -3.10,  0.02 },
+            { -0.31, -2.35, -0.43 },
+            { -0.63, -3.77,  0.01 },
+            {  3.05, -3.09, -0.04 },
+            {  3.28, -3.90,  0.41 },
+            {  3.65, -2.43,  0.30 },
+            {  6.20, -3.04, -0.03 },
+            {  5.66, -3.31,  0.71 },
+            {  6.78, -3.79, -0.19 },
+            { -6.18,  0.04, -0.04 },
+            { -6.73, -0.73, -0.15 },
+            { -5.98,  0.06,  0.89 },
+            { -3.11, -0.04, -0.04 },
+            { -3.36, -0.08,  0.87 },
+            { -2.70,  0.81, -0.14 },
+            { -0.02, -0.02, -0.05 },
+            { -0.45,  0.28,  0.75 },
+            {  0.90,  0.15,  0.07 },
+            {  3.04,  0.02, -0.01 },
+            {  3.26, -0.82,  0.38 },
+            {  3.89,  0.45, -0.13 },
+            {  6.19,  0.05, -0.03 },
+            {  5.52, -0.56,  0.25 },
+            {  7.01, -0.29,  0.32 },
+            { -6.14,  3.08,  0.00 },
+            { -6.83,  2.82,  0.61 },
+            { -6.59,  3.64, -0.64 },
+            { -3.05,  3.09, -0.04 },
+            { -3.79,  2.50,  0.09 },
+            { -3.18,  3.80,  0.59 },
+            {  0.02,  3.14,  0.04 },
+            { -0.89,  3.04, -0.19 },
+            {  0.49,  2.57, -0.57 },
+            {  3.14,  3.15,  0.00 },
+            {  3.28,  2.28,  0.37 },
+            {  2.30,  3.08, -0.45 },
+            {  6.27,  3.08, -0.00 },
+            {  5.55,  2.54, -0.33 },
+            {  5.83,  3.87,  0.34 },
+            { -6.18,  6.15, -0.03 },
+            { -6.45,  6.21,  0.88 },
+            { -6.26,  7.05, -0.36 },
+            { -3.06,  6.19, -0.05 },
+            { -2.84,  6.64,  0.76 },
+            { -3.99,  5.96,  0.03 },
+            { -0.00,  6.20,  0.06 },
+            { -0.67,  5.99, -0.59 },
+            {  0.76,  6.46, -0.44 },
+            {  3.10,  6.26, -0.03 },
+            {  3.57,  6.09,  0.78 },
+            {  2.57,  5.47, -0.18 },
+            {  6.26,  6.18,  0.02 },
+            {  5.53,  5.64, -0.29 },
+            {  5.95,  7.08, -0.06 },
+            { -6.26, -6.21,  3.07 },
+            { -5.98, -6.38,  3.97 },
+            { -5.46, -5.94,  2.62 },
+            { -3.10, -6.24,  3.04 },
+            { -2.69, -6.51,  3.87 },
+            { -3.43, -5.35,  3.21 },
+            { -0.03, -6.16,  3.06 },
+            {  0.83, -6.00,  3.42 },
+            { -0.30, -6.99,  3.45 },
+            {  3.15, -6.25,  3.11 },
+            {  2.77, -5.60,  3.72 },
+            {  2.68, -6.10,  2.28 },
+            {  6.20, -6.21,  3.16 },
+            {  5.75, -6.73,  2.50 },
+            {  6.69, -5.56,  2.66 },
+            { -6.17, -3.10,  3.04 },
+            { -6.82, -2.44,  3.28 },
+            { -6.12, -3.69,  3.80 },
+            { -3.08, -3.04,  3.11 },
+            { -3.59, -3.56,  3.72 },
+            { -2.97, -3.61,  2.34 },
+            {  0.01, -3.04,  3.11 },
+            { -0.86, -3.41,  3.20 },
+            {  0.56, -3.78,  2.86 },
+            {  3.07, -3.07,  3.15 },
+            {  3.81, -3.68,  3.13 },
+            {  2.80, -2.98,  2.23 },
+            {  6.20, -3.04,  3.13 },
+            {  5.48, -3.64,  2.92 },
+            {  6.98, -3.49,  2.81 },
+            { -6.18, -0.05,  3.12 },
+            { -6.41,  0.66,  3.69 },
+            { -6.33,  0.28,  2.23 },
+            { -3.05,  0.03,  3.10 },
+            { -3.46, -0.42,  3.83 },
+            { -3.57, -0.19,  2.33 },
+            {  0.03, -0.02,  3.15 },
+            {  0.23, -0.08,  2.21 },
+            { -0.81,  0.41,  3.18 },
+            {  3.09,  0.00,  3.03 },
+            {  2.48, -0.29,  3.71 },
+            {  3.91,  0.16,  3.51 },
+            {  6.19, -0.06,  3.11 },
+            {  6.05,  0.47,  2.33 },
+            {  6.59,  0.52,  3.74 },
+            { -6.20,  3.05,  3.05 },
+            { -6.87,  3.73,  3.17 },
+            { -5.55,  3.24,  3.73 },
+            { -3.11,  3.06,  3.15 },
+            { -3.64,  3.74,  2.71 },
+            { -2.32,  3.00,  2.62 },
+            {  0.02,  3.05,  3.06 },
+            { -0.87,  3.14,  3.38 },
+            {  0.48,  3.82,  3.42 },
+            {  3.07,  3.10,  3.16 },
+            {  3.95,  3.44,  2.97 },
+            {  2.76,  2.73,  2.32 },
+            {  6.19,  3.07,  3.16 },
+            {  7.02,  3.30,  2.72 },
+            {  5.52,  3.27,  2.51 },
+            { -6.19,  6.24,  3.15 },
+            { -5.56,  5.88,  2.52 },
+            { -7.05,  5.96,  2.83 },
+            { -3.10,  6.14,  3.08 },
+            { -2.34,  6.69,  3.27 },
+            { -3.86,  6.69,  3.29 },
+            { -0.04,  6.24,  3.13 },
+            {  0.63,  6.54,  2.53 },
+            {  0.08,  5.29,  3.18 },
+            {  3.12,  6.24,  3.14 },
+            {  3.57,  5.82,  2.40 },
+            {  2.23,  5.90,  3.12 },
+            {  6.25,  6.19,  3.06 },
+            {  5.55,  5.59,  3.32 },
+            {  6.08,  6.99,  3.55 },
+            { -6.20, -6.16,  6.15 },
+            { -6.29, -5.99,  7.09 },
+            { -6.09, -7.11,  6.09 },
+            { -3.09, -6.19,  6.27 },
+            { -2.56, -5.90,  5.52 },
+            { -3.80, -6.69,  5.87 },
+            {  0.02, -6.25,  6.24 },
+            { -0.70, -5.70,  6.51 },
+            {  0.25, -5.93,  5.36 },
+            {  3.11, -6.18,  6.14 },
+            {  3.76, -6.54,  6.74 },
+            {  2.29, -6.20,  6.64 },
+            {  6.22, -6.17,  6.15 },
+            {  6.61, -6.98,  6.47 },
+            {  5.56, -5.94,  6.81 },
+            { -6.21, -3.10,  6.14 },
+            { -6.76, -2.66,  6.78 },
+            { -5.51, -3.50,  6.65 },
+            { -3.13, -3.05,  6.18 },
+            { -2.19, -3.14,  6.34 },
+            { -3.50, -3.89,  6.43 },
+            {  0.01, -3.06,  6.15 },
+            { -0.06, -2.81,  7.07 },
+            { -0.25, -3.98,  6.13 },
+            {  3.04, -3.09,  6.17 },
+            {  3.84, -3.51,  5.84 },
+            {  3.25, -2.85,  7.08 },
+            {  6.26, -3.13,  6.19 },
+            {  6.01, -2.20,  6.09 },
+            {  5.47, -3.55,  6.54 },
+            { -6.20,  0.01,  6.27 },
+            { -5.79, -0.70,  5.78 },
+            { -6.67,  0.51,  5.60 },
+            { -3.13,  0.01,  6.14 },
+            { -3.53, -0.35,  6.94 },
+            { -2.21,  0.17,  6.39 },
+            { -0.04, -0.04,  6.20 },
+            {  0.26,  0.47,  5.46 },
+            {  0.51,  0.22,  6.93 },
+            {  3.10, -0.05,  6.23 },
+            {  2.33,  0.44,  5.95 },
+            {  3.85,  0.45,  5.92 },
+            {  6.19, -0.01,  6.26 },
+            {  7.05,  0.16,  5.88 },
+            {  5.58,  0.02,  5.52 },
+            { -6.22,  3.04,  6.17 },
+            { -5.45,  3.57,  5.95 },
+            { -6.62,  3.50,  6.92 },
+            { -3.09,  3.16,  6.21 },
+            { -3.71,  2.75,  5.61 },
+            { -2.60,  2.43,  6.59 },
+            { -0.02,  3.10,  6.26 },
+            {  0.89,  3.27,  6.05 },
+            { -0.44,  2.94,  5.41 },
+            {  3.12,  3.04,  6.23 },
+            {  2.31,  3.53,  6.43 },
+            {  3.59,  3.60,  5.60 },
+            {  6.23,  3.05,  6.24 },
+            {  5.92,  3.91,  6.54 },
+            {  6.02,  3.03,  5.30 },
+            { -6.15,  6.21,  6.24 },
+            { -6.27,  6.46,  5.32 },
+            { -7.00,  5.85,  6.51 },
+            { -3.07,  6.15,  6.22 },
+            { -3.98,  6.27,  5.94 },
+            { -2.66,  7.01,  6.10 },
+            {  0.04,  6.20,  6.25 },
+            { -0.38,  5.50,  5.75 },
+            { -0.36,  7.00,  5.93 },
+            {  3.12,  6.15,  6.24 },
+            {  3.66,  6.88,  5.93 },
+            {  2.25,  6.33,  5.86 },
+            {  6.20,  6.27,  6.19 },
+            {  5.46,  5.65,  6.19 },
+            {  6.97,  5.73,  6.39 }
+        };
+        for (int atom = 0; atom < natoms; ++atom)
+            positions.push_back(Vec3(coords[atom][0], coords[atom][1], coords[atom][2])*OpenMM::NmPerAngstrom);
+    }
+    else
+        throw exception();
+
+    system.setDefaultPeriodicBoxVectors(Vec3(boxEdgeLength, 0, 0),
+                                        Vec3(0, boxEdgeLength, 0),
+                                        Vec3(0, 0, boxEdgeLength));
+
+    sig.clear();
+    eps.clear();
+    bonds.clear();
+    for (int atom = 0; atom < natoms; ++atom) {
+        system.addParticle(masses[atom%RESSIZE]);
+        double sigma = 2.0*pow(2.0, -1.0/6.0)*halfrmins[atom%RESSIZE]*OpenMM::NmPerAngstrom;
+        double epsilon = fabs(epsilons[atom%RESSIZE])*OpenMM::KJPerKcal;
+        sig.push_back(0.5*sigma);
+        eps.push_back(2.0*sqrt(epsilon));
+        if (atom%RESSIZE == 0) {
+            bonds.push_back(pair<int, int>(atom, atom+1));
+            bonds.push_back(pair<int, int>(atom, atom+2));
+        }
+        double charge = do_electrostatics ? charges[atom] : 0;
+        forceField->addParticle(charge, sigma, epsilon);
+    }
+}
+
+void testWater2DpmeEnergiesForcesNoExclusions() {
+    const double cutoff = 7.0*OpenMM::NmPerAngstrom;
+    const double alpha = 4.0124063605;
+    const double dalpha = 4.0124063605;
+    const int grid = 32;
+    NonbondedForce* forceField = new NonbondedForce();
+
+    vector<Vec3> positions;
+    vector<double> epsvals;
+    vector<double> sigvals;
+    vector<pair<int, int> > bonds;
+    System system;
+
+    const int numAtoms = 6;
+    double boxEdgeLength = 25*OpenMM::NmPerAngstrom;
+
+    make_waterbox(numAtoms, boxEdgeLength, forceField,  positions, epsvals, sigvals, bonds, system, false);
+    forceField->setNonbondedMethod(OpenMM::NonbondedForce::LJPME);
+    forceField->setPMEParameters(alpha, grid, grid, grid);
+    forceField->setLJPMEParameters(dalpha, grid, grid, grid);
+    forceField->setCutoffDistance(cutoff);
+    forceField->setReactionFieldDielectric(1.0);
+    system.addForce(forceField);
+
+    VerletIntegrator integrator(0.01);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+
+
+//Gromacs reference values from the following inputs
+//
+//------------------------------------------------
+// water2.gro
+//------------------------------------------------
+//MD of 2 waters, t= 0.0
+//    6
+//    1SOL    OW1    1   0.200   0.200   0.200
+//    1SOL    HW1    2   0.250   0.200   0.300
+//    1SOL    HW2    3   0.150   0.200   0.300
+//    2SOL    OW1    4   0.000   0.000   0.000
+//    2SOL    HW1    5   0.050   0.000   0.100
+//    2SOL    HW2    6  -0.050   0.000   0.100
+//    2.50000   2.50000   2.50000
+//------------------------------------------------
+//
+//------------------------------------------------
+// water2.gro
+//------------------------------------------------
+//[ defaults ]
+//; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
+//  1             2               no              1.0     1.0
+//
+//[ atomtypes ]
+//; full atom descriptions are available in ffoplsaa.atp
+//; name  bond_type    mass    charge   ptype          sigma      epsilon
+//OW   OW 8     15.99940    -0.834       A    3.15057e-01  6.36386e-01
+//HW   HW 1      1.00800     0.417       A    4.00014e-02  1.92464e-01
+//
+//
+//#ifdef NOEXCLUDE
+//
+//[ moleculetype ]
+//; molname   nrexcl
+//SOL     0
+//
+//#else
+//
+//[ moleculetype ]
+//; molname   nrexcl
+//SOL     3
+//
+//#endif
+//
+//#ifdef NOCHARGE
+//[ atoms ]
+//; id    at type res nr  residu name at name     cg nr   charge
+//1      OW  1       SOL             OW1             1        0.000
+//2      HW  1       SOL             HW1             1        0.000
+//3      HW  1       SOL             HW2             1        0.000
+//#else
+//[ atoms ]
+//; id    at type res nr  residu name at name     cg nr   charge
+//1      OW  1       SOL             OW1             1       -0.834
+//2      HW  1       SOL             HW1             1        0.417
+//3      HW  1       SOL             HW2             1        0.417
+//#endif
+//
+//[ settles ]
+//; i j   funct   length
+//1   1   0.09572 0.15139
+//
+//#ifndef NOEXCLUDE
+//[ exclusions ]
+//1   2   3
+//2   1   3
+//3   1   2
+//#endif
+//------------------------------------------------
+//
+//------------------------------------------------
+// water2.top
+//------------------------------------------------
+//#include "tip3p.tpi"
+//
+//[ system ]
+//Water dimer
+//
+//[ molecules ]
+//SOL 2
+//------------------------------------------------
+//
+//------------------------------------------------
+// water2.mdp
+//------------------------------------------------
+//define = -DNOCHARGE -DNOEXCLUDE
+//;define = -DNOCHARGE
+//rvdw             = 0.7
+//vdw-modifier     = Potential-Shift
+//;vdw-modifier     = None
+//rlist            = 0.9
+//rcoulomb         = 0.7
+//nstfout          = 1
+//fourierspacing   = 0.08
+//pme-order        = 5
+//ewald-rtol       = 1.5E-2
+//ewald-rtol-lj    = 1.5E-2
+//vdwtype          = PME
+//lj-pme-comb-rule = geometric
+//continuation     = yes
+//------------------------------------------------
+//
+//run commands:-
+//    gmx grompp -c water2.gro -p water2.top -f water2.mdp -o run.tpr
+//    gmx mdrun -s run.tpr -o traj.trr -pforce 1E-16
+//    gmx dump -f traj.trr
+
+    double refenergy = 1.34804e+03;
+    vector<Vec3> refforces(6);
+    refforces[0] = Vec3( 3.15301e-01, -3.91114e-03, -6.68965e+04);
+    refforces[1] = Vec3( 1.67241e+04, -3.79846e-02,  3.34487e+04);
+    refforces[2] = Vec3(-1.67243e+04, -9.55931e-02,  3.34486e+04);
+    refforces[3] = Vec3(-1.71643e+00, -1.76696e+00, -6.68993e+04);
+    refforces[4] = Vec3( 1.67254e+04,  1.59080e+00,  3.34495e+04);
+    refforces[5] = Vec3(-1.67239e+04,  3.13897e-01,  3.34491e+04);
+
+    State state = context.getState(State::Forces | State::Energy);
+    double energy = state.getPotentialEnergy();
+    const vector<Vec3>& forces = state.getForces();
+
+    ASSERT_EQUAL_TOL(refenergy, energy, 1E-4);
+    for (int n = 0; n < numAtoms; ++n)
+        ASSERT_EQUAL_VEC(refforces[n], forces[n], 1E-4);
+}
+
+
+void testWater2DpmeEnergiesForcesWithExclusions() {
+    const double cutoff = 7.0*OpenMM::NmPerAngstrom;
+    const double alpha = 4.0124063605;
+    const double dalpha = 4.0124063605;
+    const int grid = 32;
+    NonbondedForce* forceField = new NonbondedForce();
+
+    vector<Vec3> positions;
+    vector<double> epsvals;
+    vector<double> sigvals;
+    vector<pair<int, int> > bonds;
+    System system;
+
+    const int numAtoms = 6;
+    double boxEdgeLength = 25*OpenMM::NmPerAngstrom;
+
+    make_waterbox(numAtoms, boxEdgeLength, forceField,  positions, epsvals, sigvals, bonds, system, false);
+
+    forceField->setNonbondedMethod(OpenMM::NonbondedForce::LJPME);
+    forceField->createExceptionsFromBonds(bonds, 1.0, 1.0);
+    forceField->setPMEParameters(alpha, grid, grid, grid);
+    forceField->setLJPMEParameters(dalpha, grid, grid, grid);
+    forceField->setCutoffDistance(cutoff);
+    forceField->setUseSwitchingFunction(false);
+    forceField->setUseDispersionCorrection(false);
+
+    forceField->setReactionFieldDielectric(1.0);
+    system.addForce(forceField);
+
+    VerletIntegrator integrator(0.01);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+
+    double refenergy = -7.78060e-01;
+    vector<Vec3> refforces(6);
+
+    // Gromacs reference values.  See comments in testWater2DpmeEnergiesForcesNoExclusions() for details.
+    refforces[0] = Vec3( 3.15301e-01, -3.91114e-03,  9.48373e-01);
+    refforces[1] = Vec3(-4.74745e-02, -3.79846e-02, -5.69616e-02);
+    refforces[2] = Vec3(-7.16866e-02, -9.55931e-02, -1.43348e-01);
+    refforces[3] = Vec3(-1.78136e+00, -1.76696e+00, -1.70022e+00);
+    refforces[4] = Vec3( 1.19309e+00,  1.59080e+00,  7.95443e-01);
+    refforces[5] = Vec3( 3.92366e-01,  3.13897e-01,  1.56962e-01);
+
+    State state = context.getState(State::Forces | State::Energy);
+    double energy = state.getPotentialEnergy();
+    const vector<Vec3>& forces = state.getForces();
+
+
+    ASSERT_EQUAL_TOL(refenergy, energy, 1E-4);
+    for (int n = 0; n < numAtoms; ++n)
+        ASSERT_EQUAL_VEC(refforces[n], forces[n], 5E-4);
+}
+
+
+void testWater125DpmeVsLongCutoffNoExclusions() {
+    const double cutoff = 8.5*OpenMM::NmPerAngstrom;
+    const double alpha = 0.45*OpenMM::AngstromsPerNm;
+    const double dalpha = 0.45*OpenMM::AngstromsPerNm;
+    const int grid = 32;
+    NonbondedForce* forceField = new NonbondedForce();
+
+    vector<Vec3> positions;
+    vector<double> epsvals;
+    vector<double> sigvals;
+    vector<pair<int, int> > bonds;
+    System system;
+
+    const int numAtoms = 375;
+    double boxEdgeLength = 17.01*OpenMM::NmPerAngstrom;
+
+    make_waterbox(numAtoms, boxEdgeLength, forceField,  positions, epsvals, sigvals, bonds, system, false);
+
+    forceField->setNonbondedMethod(OpenMM::NonbondedForce::LJPME);
+    forceField->setPMEParameters(alpha, grid, grid, grid);
+    forceField->setLJPMEParameters(dalpha, grid, grid, grid);
+    forceField->setCutoffDistance(cutoff);
+    forceField->setReactionFieldDielectric(1.0);
+    system.addForce(forceField);
+
+    VerletIntegrator integrator(0.01);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+
+    State state = context.getState(State::Forces | State::Energy);
+    double energy = state.getPotentialEnergy();
+
+//Gromacs reference values.  See comments in testWater2DpmeEnergiesForcesNoExclusions() for details.
+//Coordinates are from make_waterbox, and the .gro file looks like
+//
+//;define = -DNOCHARGE -DNOEXCLUDE
+//define = -DNOCHARGE
+//vdw-modifier     = Potential-Shift
+//rvdw             = 1.15
+//rlist            = 1.15
+//rcoulomb         = 1.15
+//fourierspacing   = 0.04
+//nstfout          = 1
+//pme-order        = 5
+//ewald-rtol       = 1E-8
+//ewald-rtol-lj    = 1E-8
+//vdwtype          = PME
+//lj-pme-comb-rule = geometric
+//continuation     = yes
+
+    double gromacs_energy = 5.63157e+05;
+    const vector<Vec3>& forces = state.getForces();
+    vector<Vec3> refforces(numAtoms);
+    refforces[  0] = Vec3(-1.12446e+05, -2.71952e+05, -1.91471e+05);
+    refforces[  1] = Vec3( 2.70479e+05,  6.61172e+04,  7.21277e+04);
+    refforces[  2] = Vec3(-1.58058e+05,  2.05779e+05,  1.19292e+05);
+    refforces[  3] = Vec3( 3.16438e+05, -1.27994e+05,  1.08811e+05);
+    refforces[  4] = Vec3(-1.36520e+05,  1.93405e+05,  1.36521e+05);
+    refforces[  5] = Vec3(-1.79957e+05, -6.54374e+04, -2.45393e+05);
+    refforces[  6] = Vec3( 1.30626e+05, -1.36728e+05,  2.93833e+05);
+    refforces[  7] = Vec3(-2.74564e+05,  8.02095e+04, -7.09524e+04);
+    refforces[  8] = Vec3( 1.43952e+05,  5.64523e+04, -2.22980e+05);
+    refforces[  9] = Vec3(-4.22852e+04,  2.14660e+04, -3.23254e+05);
+    refforces[ 10] = Vec3( 2.28059e+05, -4.44251e+04,  1.62898e+05);
+    refforces[ 11] = Vec3(-1.85776e+05,  2.29044e+04,  1.60326e+05);
+    refforces[ 12] = Vec3(-1.02075e+05,  3.27346e+05,  1.47993e+04);
+    refforces[ 13] = Vec3( 7.77195e+04, -1.44336e+05,  2.10959e+05);
+    refforces[ 14] = Vec3( 2.44142e+04, -1.83111e+05, -2.25837e+05);
+    refforces[ 15] = Vec3(-4.81843e+04, -2.72889e+05, -1.75068e+05);
+    refforces[ 16] = Vec3(-5.46666e+03,  2.18711e+04,  2.62456e+05);
+    refforces[ 17] = Vec3( 5.35890e+04,  2.51021e+05, -8.74311e+04);
+    refforces[ 18] = Vec3(-2.04721e+05,  1.58918e+05,  2.20448e+05);
+    refforces[ 19] = Vec3(-7.05932e+04, -1.64717e+05, -2.17659e+05);
+    refforces[ 20] = Vec3( 2.75339e+05,  5.73611e+03, -2.86718e+03);
+    refforces[ 21] = Vec3(-5.65480e+03, -2.81809e+05, -1.38358e+05);
+    refforces[ 22] = Vec3(-7.85576e+03,  2.25204e+05, -1.15217e+05);
+    refforces[ 23] = Vec3( 1.34846e+04,  5.66345e+04,  2.53512e+05);
+    refforces[ 24] = Vec3(-7.73167e+04, -4.43114e+04,  3.22354e+05);
+    refforces[ 25] = Vec3(-1.24372e+05,  1.61973e+05, -1.88001e+05);
+    refforces[ 26] = Vec3( 2.01689e+05, -1.17651e+05, -1.34455e+05);
+    refforces[ 27] = Vec3(-1.79300e+05, -1.97922e+05,  1.94294e+05);
+    refforces[ 28] = Vec3( 2.38945e+05, -4.88761e+04, -9.50350e+04);
+    refforces[ 29] = Vec3(-5.95911e+04,  2.46878e+05, -9.93159e+04);
+    refforces[ 30] = Vec3(-1.31892e+04, -2.15675e+05,  2.68757e+05);
+    refforces[ 31] = Vec3( 2.31232e+05,  1.08107e+05, -1.32129e+05);
+    refforces[ 32] = Vec3(-2.18089e+05,  1.07592e+05, -1.36669e+05);
+    refforces[ 33] = Vec3( 1.90632e+04, -3.62889e+05,  8.61608e+04);
+    refforces[ 34] = Vec3(-2.16178e+05,  1.59639e+05, -1.66288e+05);
+    refforces[ 35] = Vec3( 1.97134e+05,  2.03295e+05,  8.00862e+04);
+    refforces[ 36] = Vec3( 3.40071e+05, -6.87734e+04,  1.22584e+05);
+    refforces[ 37] = Vec3(-4.17943e+04,  8.35901e+03, -2.64693e+05);
+    refforces[ 38] = Vec3(-2.98358e+05,  6.03813e+04,  1.42075e+05);
+    refforces[ 39] = Vec3(-3.21693e+05,  4.40885e+04,  1.41154e+04);
+    refforces[ 40] = Vec3( 1.28817e+05,  2.02067e+04, -2.07114e+05);
+    refforces[ 41] = Vec3( 1.92948e+05, -6.43160e+04,  1.92949e+05);
+    refforces[ 42] = Vec3(-1.46001e+05,  3.20840e+05,  7.97305e+04);
+    refforces[ 43] = Vec3(-1.27705e+05, -2.55411e+05, -1.24428e+05);
+    refforces[ 44] = Vec3( 2.73737e+05, -6.54594e+04,  4.46323e+04);
+    refforces[ 45] = Vec3( 1.50212e+04,  2.43629e+05, -2.20084e+05);
+    refforces[ 46] = Vec3(-1.07432e+05,  8.26408e+03,  2.42419e+05);
+    refforces[ 47] = Vec3( 9.23685e+04, -2.51915e+05, -2.23909e+04);
+    refforces[ 48] = Vec3( 3.14327e+04, -2.94199e+05,  1.55484e+05);
+    refforces[ 49] = Vec3(-2.23665e+05,  1.52883e+05, -2.54793e+04);
+    refforces[ 50] = Vec3( 1.92227e+05,  1.41342e+05, -1.30033e+05);
+    refforces[ 51] = Vec3( 5.73542e+04, -2.08689e+05, -2.68170e+05);
+    refforces[ 52] = Vec3(-2.26784e+05,  1.85259e+05,  8.30474e+04);
+    refforces[ 53] = Vec3( 1.69446e+05,  2.34613e+04,  1.85087e+05);
+    refforces[ 54] = Vec3( 2.25490e+05, -1.91311e+05, -1.40103e+05);
+    refforces[ 55] = Vec3(-8.20802e+03,  6.83988e+04,  2.54448e+05);
+    refforces[ 56] = Vec3(-2.17315e+05,  1.22953e+05, -1.14373e+05);
+    refforces[ 57] = Vec3(-7.09510e+04,  2.05639e+05, -2.69540e+05);
+    refforces[ 58] = Vec3( 1.93603e+05,  3.38041e+04,  2.18192e+05);
+    refforces[ 59] = Vec3(-1.22579e+05, -2.39458e+05,  5.13124e+04);
+    refforces[ 60] = Vec3(-1.07250e+05,  3.35982e+05,  6.89600e+03);
+    refforces[ 61] = Vec3( 6.90672e-01, -1.44228e+05, -2.24659e+05);
+    refforces[ 62] = Vec3( 1.07222e+05, -1.91701e+05,  2.17695e+05);
+    refforces[ 63] = Vec3( 2.64846e+05,  1.94002e+05,  5.27262e+03);
+    refforces[ 64] = Vec3(-1.12985e+04, -2.71176e+05,  8.47510e+03);
+    refforces[ 65] = Vec3(-2.53630e+05,  7.71894e+04, -1.37831e+04);
+    refforces[ 66] = Vec3(-3.04551e+05,  4.21918e+04,  1.88948e+05);
+    refforces[ 67] = Vec3( 2.87326e+05,  1.22193e+05,  3.30263e+04);
+    refforces[ 68] = Vec3( 1.73006e+04, -1.64358e+05, -2.22023e+05);
+    refforces[ 69] = Vec3( 2.45556e+04,  2.20596e+05,  2.41913e+05);
+    refforces[ 70] = Vec3( 1.50804e+05, -2.17829e+05, -4.46813e+04);
+    refforces[ 71] = Vec3(-1.75369e+05, -2.74087e+03, -1.97287e+05);
+    refforces[ 72] = Vec3( 8.65787e+04, -2.77188e+04, -3.15014e+05);
+    refforces[ 73] = Vec3(-2.42127e+05,  6.26659e+04,  1.11092e+05);
+    refforces[ 74] = Vec3( 1.55591e+05, -3.48752e+04,  2.03883e+05);
+    refforces[ 75] = Vec3( 7.06224e+04,  2.96642e+05, -1.51267e+05);
+    refforces[ 76] = Vec3(-5.39280e+04, -2.57657e+05, -1.13850e+05);
+    refforces[ 77] = Vec3(-1.67422e+04, -3.90659e+04,  2.65102e+05);
+    refforces[ 78] = Vec3(-4.52577e+04, -3.21552e+05, -7.01075e+04);
+    refforces[ 79] = Vec3( 2.35182e+05,  1.45173e+05,  3.48411e+04);
+    refforces[ 80] = Vec3(-1.89931e+05,  1.76363e+05,  3.52722e+04);
+    refforces[ 81] = Vec3(-2.29347e+05,  1.06977e+05, -2.70703e+05);
+    refforces[ 82] = Vec3(-1.17929e+04, -2.56493e+05,  1.17929e+05);
+    refforces[ 83] = Vec3( 2.41161e+05,  1.49452e+05,  1.52847e+05);
+    refforces[ 84] = Vec3( 2.32633e+03,  3.03443e+05,  1.27473e+05);
+    refforces[ 85] = Vec3(-2.11215e+05, -1.63335e+05, -4.50583e+04);
+    refforces[ 86] = Vec3( 2.08887e+05, -1.40170e+05, -8.24544e+04);
+    refforces[ 87] = Vec3( 5.83172e+04,  2.82133e+04, -3.26001e+05);
+    refforces[ 88] = Vec3( 1.76890e+05, -4.71698e+04,  2.15220e+05);
+    refforces[ 89] = Vec3(-2.35168e+05,  1.89225e+04,  1.10825e+05);
+    refforces[ 90] = Vec3(-2.72444e+05, -1.46998e+05, -1.00637e+05);
+    refforces[ 91] = Vec3( 2.78426e+04,  2.39442e+05,  1.16936e+05);
+    refforces[ 92] = Vec3( 2.44570e+05, -9.23920e+04, -1.63042e+04);
+    refforces[ 93] = Vec3(-3.10100e+04,  2.93390e+05, -1.87196e+05);
+    refforces[ 94] = Vec3(-6.38830e+04, -2.90671e+05, -6.38833e+04);
+    refforces[ 95] = Vec3( 9.48775e+04, -2.79002e+03,  2.51149e+05);
+    refforces[ 96] = Vec3( 4.18753e+04, -1.23438e+05, -3.10188e+05);
+    refforces[ 97] = Vec3( 1.29157e+05,  2.04500e+05,  9.41770e+04);
+    refforces[ 98] = Vec3(-1.71038e+05, -8.10180e+04,  2.16049e+05);
+    refforces[ 99] = Vec3(-5.61490e+04,  2.26993e+04, -3.44090e+05);
+    refforces[100] = Vec3(-1.96230e+05, -2.88572e+04,  1.93344e+05);
+    refforces[101] = Vec3( 2.52383e+05,  6.15570e+03,  1.50813e+05);
+    refforces[102] = Vec3(-6.33126e+04,  3.55941e+05,  8.51300e+04);
+    refforces[103] = Vec3(-1.75405e+05, -1.81783e+05, -1.69027e+05);
+    refforces[104] = Vec3( 2.38759e+05, -1.74232e+05,  8.38891e+04);
+    refforces[105] = Vec3( 1.51480e+05, -4.90606e+04,  3.17956e+05);
+    refforces[106] = Vec3( 4.93229e+04, -1.61668e+05, -2.02771e+05);
+    refforces[107] = Vec3(-2.00831e+05,  2.10712e+05, -1.15233e+05);
+    refforces[108] = Vec3( 2.20204e+05, -2.33820e+05,  1.51386e+05);
+    refforces[109] = Vec3( 3.36498e+04,  2.79296e+05, -1.51424e+05);
+    refforces[110] = Vec3(-2.53896e+05, -4.54334e+04,  2.10242e-01);
+    refforces[111] = Vec3(-1.94665e+05,  2.05890e+05,  2.40510e+05);
+    refforces[112] = Vec3(-9.73233e+04, -1.29764e+05, -2.62775e+05);
+    refforces[113] = Vec3( 2.92049e+05, -7.61856e+04,  2.22208e+04);
+    refforces[114] = Vec3( 1.60261e+05, -3.43716e+05,  1.52456e+04);
+    refforces[115] = Vec3( 1.11151e+05,  3.08358e+05, -8.60527e+04);
+    refforces[116] = Vec3(-2.71443e+05,  3.54056e+04,  7.08103e+04);
+    refforces[117] = Vec3(-4.27331e+04, -3.19871e+05, -1.68415e+05);
+    refforces[118] = Vec3( 2.28408e+05,  2.15171e+05, -2.31720e+04);
+    refforces[119] = Vec3(-1.85616e+05,  1.04782e+05,  1.91602e+05);
+    refforces[120] = Vec3(-1.66057e+05, -9.58191e+04, -2.78448e+05);
+    refforces[121] = Vec3( 1.91258e+05,  2.19476e+05,  6.89779e+04);
+    refforces[122] = Vec3(-2.52379e+04, -1.23674e+05,  2.09489e+05);
+    refforces[123] = Vec3( 6.55172e+03, -9.23173e+04,  3.29554e+05);
+    refforces[124] = Vec3(-2.14685e+05,  1.13143e+05, -1.36353e+05);
+    refforces[125] = Vec3( 2.08124e+05, -2.08124e+04, -1.93257e+05);
+    refforces[126] = Vec3( 1.02696e+05, -3.39299e+05, -4.47114e+04);
+    refforces[127] = Vec3(-1.81786e+05,  1.75407e+05, -1.69029e+05);
+    refforces[128] = Vec3( 7.90534e+04,  1.63962e+05,  2.13738e+05);
+    refforces[129] = Vec3(-3.45588e+05,  9.36838e+04,  5.67939e+04);
+    refforces[130] = Vec3( 1.44967e+05, -2.60943e+05,  7.08730e+04);
+    refforces[131] = Vec3( 2.00649e+05,  1.67207e+05, -1.27685e+05);
+    refforces[132] = Vec3(-2.70181e+05,  2.05707e+05, -3.11852e+04);
+    refforces[133] = Vec3( 1.09331e+05, -1.83207e+05, -1.86162e+05);
+    refforces[134] = Vec3( 1.60883e+05, -2.25804e+04,  2.17340e+05);
+    refforces[135] = Vec3(-1.04496e+05, -2.97000e+05, -1.63733e+05);
+    refforces[136] = Vec3( 2.11303e+05,  1.73661e+04,  1.79462e+05);
+    refforces[137] = Vec3(-1.06849e+05,  2.79694e+05, -1.57134e+04);
+    refforces[138] = Vec3(-3.82271e+04,  2.11942e+05,  2.60118e+05);
+    refforces[139] = Vec3(-1.96098e+05, -1.23692e+05, -1.71962e+05);
+    refforces[140] = Vec3( 2.34334e+05, -8.82195e+04, -8.82188e+04);
+    refforces[141] = Vec3( 2.17634e+05,  2.72527e+05,  1.42967e+05);
+    refforces[142] = Vec3( 9.30924e+04, -1.86185e+05, -1.98197e+05);
+    refforces[143] = Vec3(-3.10770e+05, -8.63246e+04,  5.52480e+04);
+    refforces[144] = Vec3(-1.63880e+05, -2.98869e+05,  1.01299e+05);
+    refforces[145] = Vec3( 6.83424e+04,  2.39196e+05,  1.61538e+05);
+    refforces[146] = Vec3( 9.55787e+04,  5.97354e+04, -2.62846e+05);
+    refforces[147] = Vec3( 1.06430e+05, -2.97083e+05, -7.97261e+04);
+    refforces[148] = Vec3(-1.14920e+05,  6.41394e+04,  2.21823e+05);
+    refforces[149] = Vec3( 8.52531e+03,  2.33043e+05, -1.42101e+05);
+    refforces[150] = Vec3(-4.96390e+04, -4.12076e+04, -3.67834e+05);
+    refforces[151] = Vec3( 1.25352e+05, -1.99962e+05,  1.61166e+05);
+    refforces[152] = Vec3(-7.57841e+04,  2.41138e+05,  2.06689e+05);
+    refforces[153] = Vec3(-3.06397e+05, -5.15222e+04, -1.37080e+05);
+    refforces[154] = Vec3( 1.92949e+05, -1.92945e+05,  6.43163e+04);
+    refforces[155] = Vec3( 1.13490e+05,  2.44443e+05,  7.27504e+04);
+    refforces[156] = Vec3(-3.74471e+03,  3.83223e+05, -2.14754e+04);
+    refforces[157] = Vec3( 2.17881e+05, -1.85835e+05, -1.05735e+05);
+    refforces[158] = Vec3(-2.14191e+05, -1.97454e+05,  1.27176e+05);
+    refforces[159] = Vec3(-3.33357e+05, -1.38307e+04, -1.17320e+05);
+    refforces[160] = Vec3( 2.04157e+05, -9.93173e+04, -1.37945e+05);
+    refforces[161] = Vec3( 1.29262e+05,  1.13105e+05,  2.55294e+05);
+    refforces[162] = Vec3( 2.50185e+05,  2.64217e+05, -7.18219e+04);
+    refforces[163] = Vec3(-2.46668e+05,  1.94013e+04, -9.97745e+04);
+    refforces[164] = Vec3(-3.50266e+03, -2.83658e+05,  1.71598e+05);
+    refforces[165] = Vec3(-2.05824e+05,  2.71177e+05, -1.16442e+05);
+    refforces[166] = Vec3(-3.52163e+04, -1.88897e+05,  2.36923e+05);
+    refforces[167] = Vec3( 2.41012e+05, -8.22954e+04, -1.20505e+05);
+    refforces[168] = Vec3( 6.89329e+04,  2.09489e+05,  2.76583e+05);
+    refforces[169] = Vec3( 1.87999e+05, -1.61968e+05, -1.24368e+05);
+    refforces[170] = Vec3(-2.56932e+05, -4.75793e+04, -1.52255e+05);
+    refforces[171] = Vec3( 3.39431e+05, -5.75509e+04,  1.62795e+05);
+    refforces[172] = Vec3(-1.27767e+05,  2.66184e+05, -1.59709e+05);
+    refforces[173] = Vec3(-2.11726e+05, -2.08613e+05, -3.11349e+03);
+    refforces[174] = Vec3(-2.58601e+05,  4.61975e+04, -2.46016e+05);
+    refforces[175] = Vec3( 7.15587e+04, -2.52015e+05,  1.40007e+05);
+    refforces[176] = Vec3( 1.87109e+05,  2.05820e+05,  1.06028e+05);
+    refforces[177] = Vec3( 3.92176e+03,  2.94819e+05, -1.84065e+05);
+    refforces[178] = Vec3(-1.67230e+05, -8.36141e+04,  2.29167e+05);
+    refforces[179] = Vec3( 1.63334e+05, -2.11213e+05, -4.50586e+04);
+    refforces[180] = Vec3( 1.11151e+05,  2.40551e+05, -2.68210e+05);
+    refforces[181] = Vec3(-1.76495e+05, -2.47099e+05, -3.53002e+04);
+    refforces[182] = Vec3( 6.52872e+04,  6.52855e+03,  3.03586e+05);
+    refforces[183] = Vec3(-4.84049e+04, -2.73305e+05, -2.95224e+05);
+    refforces[184] = Vec3(-9.04206e+04, -1.44672e+04,  3.29135e+05);
+    refforces[185] = Vec3( 1.38830e+05,  2.87820e+05, -3.38610e+04);
+    refforces[186] = Vec3(-2.09072e+05, -1.53610e+05, -2.86654e+05);
+    refforces[187] = Vec3(-1.30322e+05,  9.09223e+04,  2.42461e+05);
+    refforces[188] = Vec3( 3.39380e+05,  6.27113e+04,  4.42669e+04);
+    refforces[189] = Vec3(-3.15692e+05,  1.48900e+05, -9.20420e+04);
+    refforces[190] = Vec3( 7.13687e+04, -2.72503e+05,  1.26518e+05);
+    refforces[191] = Vec3( 2.44351e+05,  1.23612e+05, -3.44964e+04);
+    refforces[192] = Vec3(-2.80727e+04,  3.15077e+05, -2.05427e+05);
+    refforces[193] = Vec3(-2.29004e+05, -2.08496e+05,  9.57028e+04);
+    refforces[194] = Vec3( 2.57096e+05, -1.06603e+05,  1.09738e+05);
+    refforces[195] = Vec3( 3.33213e+05, -7.80040e+04, -5.04749e+03);
+    refforces[196] = Vec3(-2.07680e+05, -7.82574e+04,  1.83605e+05);
+    refforces[197] = Vec3(-1.25574e+05,  1.56274e+05, -1.78599e+05);
+    refforces[198] = Vec3( 2.66795e+05, -2.82879e+04, -2.26620e+05);
+    refforces[199] = Vec3(-2.28291e+05, -1.82015e+05,  4.01046e+04);
+    refforces[200] = Vec3(-3.85027e+04,  2.10288e+05,  1.86592e+05);
+    refforces[201] = Vec3( 1.93076e+05,  2.05300e+05,  2.64591e+05);
+    refforces[202] = Vec3(-3.32259e+05, -3.65117e+04, -8.39768e+04);
+    refforces[203] = Vec3( 1.39204e+05, -1.68822e+05, -1.80669e+05);
+    refforces[204] = Vec3( 2.15422e+05,  2.88264e+05,  2.49786e+04);
+    refforces[205] = Vec3( 4.29235e+04, -2.66744e+05,  1.13441e+05);
+    refforces[206] = Vec3(-2.58339e+05, -2.15280e+04, -1.38395e+05);
+    refforces[207] = Vec3( 3.27580e+05, -4.93781e+04,  7.43625e+03);
+    refforces[208] = Vec3(-2.11622e+05, -1.58716e+05, -9.69925e+04);
+    refforces[209] = Vec3(-1.15917e+05,  2.08125e+05,  8.95715e+04);
+    refforces[210] = Vec3( 1.10967e+05, -2.71544e+05, -2.06279e+05);
+    refforces[211] = Vec3(-8.86155e+04,  1.96918e+04,  2.98676e+05);
+    refforces[212] = Vec3(-2.23913e+04,  2.51912e+05, -9.23690e+04);
+    refforces[213] = Vec3( 1.91611e+05, -7.99998e+04, -2.83463e+05);
+    refforces[214] = Vec3( 7.08725e+04,  1.44963e+05,  2.60941e+05);
+    refforces[215] = Vec3(-2.62503e+05, -6.49202e+04,  2.25808e+04);
+    refforces[216] = Vec3(-6.63168e+04, -2.84265e+04,  3.72700e+05);
+    refforces[217] = Vec3(-2.02131e+05, -6.33547e+04, -1.96097e+05);
+    refforces[218] = Vec3( 2.68466e+05,  9.18414e+04, -1.76621e+05);
+    refforces[219] = Vec3(-6.80052e+03,  2.72744e+05, -2.21849e+05);
+    refforces[220] = Vec3( 1.52709e+05, -5.52356e+04,  2.63181e+05);
+    refforces[221] = Vec3(-1.45890e+05, -2.17461e+05, -4.12893e+04);
+    refforces[222] = Vec3( 3.07522e+05, -1.21146e+05,  1.14593e+05);
+    refforces[223] = Vec3(-2.11787e+05, -1.56664e+05, -8.99363e+04);
+    refforces[224] = Vec3(-9.57076e+04,  2.77857e+05, -2.46987e+04);
+    refforces[225] = Vec3(-3.24882e+05, -3.09840e+04, -1.31947e+05);
+    refforces[226] = Vec3( 8.33124e+04, -5.05809e+04,  2.67793e+05);
+    refforces[227] = Vec3( 2.41536e+05,  8.15185e+04, -1.35863e+05);
+    refforces[228] = Vec3(-2.16561e+04, -1.67607e+05, -2.70252e+05);
+    refforces[229] = Vec3( 1.10825e+05, -7.29802e+04,  2.24354e+05);
+    refforces[230] = Vec3(-8.91995e+04,  2.40572e+05,  4.59512e+04);
+    refforces[231] = Vec3(-2.22246e+05,  1.96766e+05, -2.46638e+05);
+    refforces[232] = Vec3( 3.04748e+05,  5.66971e+04,  1.27568e+05);
+    refforces[233] = Vec3(-8.24632e+04, -2.53495e+05,  1.19113e+05);
+    refforces[234] = Vec3( 2.20621e+05, -2.03896e+05,  6.63128e+04);
+    refforces[235] = Vec3(-9.59082e+04,  1.64055e+05,  1.53960e+05);
+    refforces[236] = Vec3(-1.24758e+05,  3.98170e+04, -2.20320e+05);
+    refforces[237] = Vec3(-7.79583e+03, -3.49688e+04,  3.64330e+05);
+    refforces[238] = Vec3(-1.43293e+05, -1.65580e+05, -2.10163e+05);
+    refforces[239] = Vec3( 1.51159e+05,  2.00522e+05, -1.54247e+05);
+    refforces[240] = Vec3( 1.82058e+05, -3.72909e+04, -2.80373e+05);
+    refforces[241] = Vec3(-1.95792e+05,  1.98809e+05,  7.22936e+04);
+    refforces[242] = Vec3( 1.36910e+04, -1.61546e+05,  2.08094e+05);
+    refforces[243] = Vec3( 1.40288e+05,  3.27379e+05,  4.78569e+03);
+    refforces[244] = Vec3(-1.70016e+05, -1.73349e+05,  2.03353e+05);
+    refforces[245] = Vec3( 2.97331e+04, -1.54072e+05, -2.08134e+05);
+    refforces[246] = Vec3( 1.21932e+05,  3.52267e+05,  4.69291e+04);
+    refforces[247] = Vec3(-2.91621e+05, -1.24021e+05,  3.01674e+04);
+    refforces[248] = Vec3( 1.69673e+05, -2.28288e+05, -7.71238e+04);
+    refforces[249] = Vec3(-1.41105e+05,  1.52818e+05,  2.59197e+05);
+    refforces[250] = Vec3( 2.15511e+05, -1.77650e+05, -5.82442e+03);
+    refforces[251] = Vec3(-7.43779e+04,  2.47926e+04, -2.53438e+05);
+    refforces[252] = Vec3(-3.34214e+04,  3.09554e+05,  1.58194e+05);
+    refforces[253] = Vec3(-2.05877e+05, -1.71563e+05, -6.00467e+04);
+    refforces[254] = Vec3( 2.39329e+05, -1.38076e+05, -9.81875e+04);
+    refforces[255] = Vec3( 1.32793e+05, -3.72271e+05,  2.88491e+04);
+    refforces[256] = Vec3(-9.02633e+04,  2.78646e+05,  2.23700e+05);
+    refforces[257] = Vec3(-4.25632e+04,  9.36427e+04, -2.52553e+05);
+    refforces[258] = Vec3( 2.97252e+05,  2.17294e+05, -2.50945e+03);
+    refforces[259] = Vec3(-1.35724e+05, -1.48966e+05,  2.41658e+05);
+    refforces[260] = Vec3(-1.61536e+05, -6.83423e+04, -2.39199e+05);
+    refforces[261] = Vec3( 2.50550e+05, -1.39926e+05,  2.48375e+05);
+    refforces[262] = Vec3( 5.51783e+04, -1.65530e+04, -2.59341e+05);
+    refforces[263] = Vec3(-3.05734e+05,  1.56505e+05,  1.09189e+04);
+    refforces[264] = Vec3(-4.44616e+04,  4.17040e+04, -3.31507e+05);
+    refforces[265] = Vec3(-1.79702e+05, -8.54319e+04,  2.00323e+05);
+    refforces[266] = Vec3( 2.24182e+05,  4.37421e+04,  1.31227e+05);
+    refforces[267] = Vec3(-9.89367e+04, -3.76136e+05,  2.17157e+04);
+    refforces[268] = Vec3(-4.44424e+04,  1.68244e+05, -2.47606e+05);
+    refforces[269] = Vec3( 1.43421e+05,  2.07965e+05,  2.25893e+05);
+    refforces[270] = Vec3(-1.09030e+03, -2.44686e+05, -2.30145e+05);
+    refforces[271] = Vec3(-1.86963e+05,  1.89757e+05,  3.34863e+04);
+    refforces[272] = Vec3( 1.88003e+05,  5.49541e+04,  1.96680e+05);
+    refforces[273] = Vec3(-1.15452e+05, -1.55243e+05,  2.81412e+05);
+    refforces[274] = Vec3(-1.35893e+05,  1.74355e+05, -1.12817e+05);
+    refforces[275] = Vec3( 2.51361e+05, -1.90905e+04, -1.68633e+05);
+    refforces[276] = Vec3( 1.76202e+05, -2.31601e+05, -2.00886e+05);
+    refforces[277] = Vec3(-2.96697e+05,  3.00025e+04,  1.06676e+05);
+    refforces[278] = Vec3( 1.20457e+05,  2.01635e+05,  9.42700e+04);
+    refforces[279] = Vec3(-1.66320e+05, -9.08355e+02,  2.65478e+05);
+    refforces[280] = Vec3( 2.44823e+05,  9.45895e+04, -5.28589e+04);
+    refforces[281] = Vec3(-7.84765e+04, -9.36654e+04, -2.12648e+05);
+    refforces[282] = Vec3(-6.58002e+03, -1.21918e+05,  3.16459e+05);
+    refforces[283] = Vec3( 2.15225e+05,  5.96416e+04, -1.14097e+05);
+    refforces[284] = Vec3(-2.08615e+05,  6.22725e+04, -2.02387e+05);
+    refforces[285] = Vec3( 7.09186e+04,  1.83619e+05,  2.71862e+05);
+    refforces[286] = Vec3( 1.78911e+05, -1.02234e+05, -1.78910e+05);
+    refforces[287] = Vec3(-2.49882e+05, -8.13582e+04, -9.29803e+04);
+    refforces[288] = Vec3(-1.35529e+04, -3.20226e+05, -1.16293e+05);
+    refforces[289] = Vec3( 2.28053e+05,  1.65034e+05,  5.70125e+04);
+    refforces[290] = Vec3(-2.14515e+05,  1.55237e+05,  5.92732e+04);
+    refforces[291] = Vec3(-2.65006e+05,  1.75233e+05,  1.91268e+05);
+    refforces[292] = Vec3( 2.29905e+05,  1.02940e+05, -2.05886e+05);
+    refforces[293] = Vec3( 3.51345e+04, -2.78156e+05,  1.46393e+04);
+    refforces[294] = Vec3( 1.59458e+05,  2.25128e+05,  2.11611e+05);
+    refforces[295] = Vec3( 1.24813e+05, -1.16492e+05, -2.05249e+05);
+    refforces[296] = Vec3(-2.84281e+05, -1.08601e+05, -6.38824e+03);
+    refforces[297] = Vec3( 2.61767e+05, -7.55993e+04, -2.32576e+05);
+    refforces[298] = Vec3(-2.07803e+05, -1.78116e+05,  7.71831e+04);
+    refforces[299] = Vec3(-5.39242e+04,  2.53755e+05,  1.55427e+05);
+    refforces[300] = Vec3(-6.44117e+03,  2.31328e+05, -2.54922e+05);
+    refforces[301] = Vec3(-2.61097e+04,  4.93208e+04,  2.72709e+05);
+    refforces[302] = Vec3( 3.25053e+04, -2.80718e+05, -1.77304e+04);
+    refforces[303] = Vec3( 7.06280e+04,  7.26140e+04,  3.28446e+05);
+    refforces[304] = Vec3( 1.45890e+05,  7.98271e+04, -2.06449e+05);
+    refforces[305] = Vec3(-2.16516e+05, -1.52473e+05, -1.21980e+05);
+    refforces[306] = Vec3( 1.94879e+05, -2.83084e+05,  1.41824e+05);
+    refforces[307] = Vec3(-2.57149e+05,  1.96432e+05,  9.64284e+04);
+    refforces[308] = Vec3( 6.22633e+04,  8.66280e+04, -2.38228e+05);
+    refforces[309] = Vec3( 3.26456e+04,  1.17141e+05, -3.28374e+05);
+    refforces[310] = Vec3( 2.01298e+05, -1.11486e+05,  1.85810e+05);
+    refforces[311] = Vec3(-2.33937e+05, -5.70482e+03,  1.42641e+05);
+    refforces[312] = Vec3( 6.42964e+04,  1.88719e+05, -2.86572e+05);
+    refforces[313] = Vec3( 1.22184e+05, -2.53768e+05,  1.00254e+05);
+    refforces[314] = Vec3(-1.86433e+05,  6.49694e+04,  1.86429e+05);
+    refforces[315] = Vec3(-4.12317e+04, -1.73607e+04, -3.68620e+05);
+    refforces[316] = Vec3(-1.78981e+05,  1.43186e+05,  2.08269e+05);
+    refforces[317] = Vec3( 2.20155e+05, -1.25802e+05,  1.60395e+05);
+    refforces[318] = Vec3(-1.58613e+05,  3.01594e+05, -1.29345e+05);
+    refforces[319] = Vec3( 2.79698e+05, -2.67790e+04,  4.76063e+04);
+    refforces[320] = Vec3(-1.21063e+05, -2.74847e+05,  8.17976e+04);
+    refforces[321] = Vec3( 1.00451e+05,  2.03427e+05, -2.75044e+05);
+    refforces[322] = Vec3(-2.13956e+04,  7.64144e+04,  2.81199e+05);
+    refforces[323] = Vec3(-7.91039e+04, -2.79910e+05, -6.08590e+03);
+    refforces[324] = Vec3(-2.80664e+05,  5.26142e+04, -1.53706e+05);
+    refforces[325] = Vec3( 2.23923e+05, -1.17559e+05, -9.23685e+04);
+    refforces[326] = Vec3( 5.68062e+04,  6.49218e+04,  2.46158e+05);
+    refforces[327] = Vec3( 2.88911e+05, -1.17893e+05, -7.40828e+04);
+    refforces[328] = Vec3(-6.38612e+04,  2.37565e+05, -2.55450e+04);
+    refforces[329] = Vec3(-2.25033e+05, -1.19637e+05,  9.96952e+04);
+    refforces[330] = Vec3( 1.03510e+04,  7.35559e+04,  3.47129e+05);
+    refforces[331] = Vec3( 1.26778e+05, -2.19542e+05, -1.51515e+05);
+    refforces[332] = Vec3(-1.37190e+05,  1.45950e+05, -1.95573e+05);
+    refforces[333] = Vec3(-1.31816e+05,  5.57951e+04, -2.81935e+05);
+    refforces[334] = Vec3(-1.08367e+05, -9.75306e+04,  2.16731e+05);
+    refforces[335] = Vec3( 2.40191e+05,  4.17715e+04,  6.52667e+04);
+    refforces[336] = Vec3(-2.86680e+05, -2.63007e+05,  1.37964e+04);
+    refforces[337] = Vec3( 1.03914e+05,  1.76653e+05, -2.56324e+05);
+    refforces[338] = Vec3( 1.82779e+05,  8.64040e+04,  2.42594e+05);
+    refforces[339] = Vec3( 1.09800e+03, -3.11636e+05,  1.85833e+05);
+    refforces[340] = Vec3(-2.39750e+05,  1.52568e+05, -8.71817e+04);
+    refforces[341] = Vec3( 2.38635e+05,  1.59089e+05, -9.86351e+04);
+    refforces[342] = Vec3(-8.76706e+04, -6.10520e+04,  3.31680e+05);
+    refforces[343] = Vec3( 2.64692e+05,  5.23233e+04, -1.16959e+05);
+    refforces[344] = Vec3(-1.76973e+05,  8.70353e+03, -2.14689e+05);
+    refforces[345] = Vec3(-1.15976e+05, -2.72303e+05, -1.33307e+05);
+    refforces[346] = Vec3( 2.20682e+05,  1.51896e+05, -6.30516e+04);
+    refforces[347] = Vec3(-1.04744e+05,  1.20457e+05,  1.96395e+05);
+    refforces[348] = Vec3( 4.57064e+04,  3.43550e+05,  7.23245e+04);
+    refforces[349] = Vec3(-1.91414e+05, -1.26580e+05, -1.85241e+05);
+    refforces[350] = Vec3( 1.45684e+05, -2.17040e+05,  1.12977e+05);
+    refforces[351] = Vec3(-1.88623e+05, -1.22936e+04,  3.10223e+05);
+    refforces[352] = Vec3( 3.06931e+05,  5.73380e+04, -7.08300e+04);
+    refforces[353] = Vec3(-1.18279e+05, -4.50586e+04, -2.39376e+05);
+    refforces[354] = Vec3( 8.31603e+04, -2.75958e+05,  1.16877e+05);
+    refforces[355] = Vec3(-2.08784e+05,  1.26300e+05,  5.15494e+04);
+    refforces[356] = Vec3( 1.25608e+05,  1.49661e+05, -1.68370e+05);
+    refforces[357] = Vec3( 1.44100e+05, -2.34453e+05,  1.73916e+05);
+    refforces[358] = Vec3(-8.65072e+04,  2.39988e+05,  8.37141e+04);
+    refforces[359] = Vec3(-5.75432e+04, -5.48037e+03, -2.57574e+05);
+    refforces[360] = Vec3( 2.72504e+05,  2.99781e+04,  1.85778e+05);
+    refforces[361] = Vec3(-3.41026e+04,  7.10483e+04, -2.61466e+05);
+    refforces[362] = Vec3(-2.38458e+05, -1.00996e+05,  7.57451e+04);
+    refforces[363] = Vec3( 1.45826e+05, -2.81284e+05,  1.15502e+05);
+    refforces[364] = Vec3(-2.63407e+05,  3.47332e+04, -8.10477e+04);
+    refforces[365] = Vec3( 1.17594e+05,  2.46657e+05, -3.44181e+04);
+    refforces[366] = Vec3( 2.59381e+05, -5.73297e+04,  2.56690e+05);
+    refforces[367] = Vec3(-1.25834e+05, -2.09725e+05, -1.49803e+05);
+    refforces[368] = Vec3(-1.33554e+05,  2.67105e+05, -1.06844e+05);
+    refforces[369] = Vec3( 7.18473e+04, -2.59021e+05,  1.89797e+05);
+    refforces[370] = Vec3( 1.56664e+05,  2.11783e+05, -8.99369e+04);
+    refforces[371] = Vec3(-2.28511e+05,  4.72761e+04, -9.98083e+04);
+    refforces[372] = Vec3(-1.99083e+04,  3.17042e+05, -5.62792e+04);
+    refforces[373] = Vec3(-1.96875e+05, -1.64948e+05, -7.96758e-01);
+    refforces[374] = Vec3( 2.16840e+05, -1.52072e+05,  5.63215e+04);
+
+
+    const double longcutoff = 30.0*OpenMM::NmPerAngstrom;
+    const double longcutoff2 = longcutoff*longcutoff;
+    const double cutoff2 = cutoff*cutoff;
+    const double cutoff6inv = 1.0 / (cutoff2*cutoff2*cutoff2);
+    const int nboxes = ceil(longcutoff/boxEdgeLength);
+
+    double refenergy = 0.0;
+    // Loop over home box first...
+    for (int i = 0; i < numAtoms; ++ i) {
+        for (int j = i+1; j < numAtoms; ++j) {
+            Vec3 dR = positions[i] - positions[j];
+            double R2 = dR[0]*dR[0] + dR[1]*dR[1] + dR[2]*dR[2];
+            double sig2 = (sigvals[i] + sigvals[j])*(sigvals[i] + sigvals[j]) / R2;
+            double sig6 = sig2*sig2*sig2;
+            double eps = epsvals[i]*epsvals[j];
+            refenergy += 2.0*eps*(sig6-1.0)*sig6;
+            if (R2 < cutoff2) {
+                // Add a shift term for direct space parts withing t
+                refenergy += 2.0*eps*(pow(sigvals[i]+sigvals[j], 6) - 64.0*pow(sigvals[i]*sigvals[j], 3))*cutoff6inv;
+            }
+        }
+    }
+
+    // ... and now add in the image box terms
+    for (int bx = -nboxes; bx <= nboxes; ++bx) {
+        for (int by = -nboxes; by <= nboxes; ++by) {
+            for (int bz = -nboxes; bz <= nboxes; ++bz) {
+                if (bx==0 && by==0 && bz==0) continue;
+                Vec3 offset(bx*boxEdgeLength, by*boxEdgeLength, bz*boxEdgeLength);
+                for (int i = 0; i < numAtoms; ++ i) {
+                    for (int j = 0; j < numAtoms; ++j) {
+                        Vec3 dR = positions[i] - positions[j] + offset;
+                        double R2 = dR[0]*dR[0] + dR[1]*dR[1] + dR[2]*dR[2];
+                        if (R2 > longcutoff2) continue;
+                        double sig2 = (sigvals[i] + sigvals[j])*(sigvals[i] + sigvals[j]) / R2;
+                        double sig6 = sig2*sig2*sig2;
+                        double eps = epsvals[i]*epsvals[j];
+                        refenergy += eps*(sig6-1.0)*sig6;
+                        if (R2 < cutoff2) {
+                            // Add a shift term for direct space parts withing teh
+                            refenergy += eps*(pow(sigvals[i]+sigvals[j], 6) - 64.0*pow(sigvals[i]*sigvals[j], 3))*cutoff6inv;
+                        }
+                    }
+                }
+            }
+        }
+    }
+    refenergy *= 0.5;
+
+    // For this test the reference energy is 545636 kJ/mol, while the difference between DPME and 30A cutoffs
+    // is just 1.062 kJ/mol.  The difference is due to the fact that arithmetic mean combination rules are used
+    // up to the cutoff, while the reciprocal space uses the geometric mean.  See DOI: 10.1021/acs.jctc.5b00726
+
+    ASSERT_EQUAL_TOL(refenergy, energy, 5E-6);
+    ASSERT_EQUAL_TOL(gromacs_energy, energy, 5E-4);
+
+    for (int n = 0; n < numAtoms; ++n)
+        ASSERT_EQUAL_VEC(refforces[n], forces[n], 5E-4);
+}
+
+
+
+void testWater125DpmeVsLongCutoffWithExclusions() {
+    const double cutoff = 11.5*OpenMM::NmPerAngstrom;
+    const double alpha = 0.45*OpenMM::AngstromsPerNm;
+    const double dalpha = 4.2760443169;
+    const int grid = 60;
+    NonbondedForce* forceField = new NonbondedForce();
+
+    vector<Vec3> positions;
+    vector<double> epsvals;
+    vector<double> sigvals;
+    vector<pair<int, int> > bonds;
+    System system;
+
+    const int numAtoms = 375;
+    double boxEdgeLength = 23.01*OpenMM::NmPerAngstrom;
+
+    make_waterbox(numAtoms, boxEdgeLength, forceField,  positions, epsvals, sigvals, bonds, system, false);
+
+    forceField->setNonbondedMethod(OpenMM::NonbondedForce::LJPME);
+    forceField->createExceptionsFromBonds(bonds, 1.0, 1.0);
+    forceField->setPMEParameters(alpha, grid, grid, grid);
+    forceField->setLJPMEParameters(dalpha, grid, grid, grid);
+    forceField->setCutoffDistance(cutoff);
+    forceField->setReactionFieldDielectric(1.0);
+    system.addForce(forceField);
+
+    VerletIntegrator integrator(0.01);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+
+    State state = context.getState(State::Forces | State::Energy);
+    double energy = state.getPotentialEnergy();
+    const vector<Vec3>& forces = state.getForces();
+    // Gromacs reference values.  See comments in testWater2DpmeEnergiesForcesNoExclusions() for details.
+    double gromacs_energy = -2.17481e+02;
+    vector<Vec3> refforces(numAtoms);
+    refforces[  0] = Vec3(-3.10130e+01, -6.01219e+01, -5.31302e+01);
+    refforces[  1] = Vec3( 4.21376e+00,  9.13811e-01,  9.61319e-01);
+    refforces[  2] = Vec3( 1.42363e+00,  2.97219e+00,  1.27185e+00);
+    refforces[  3] = Vec3(-3.70172e+01, -2.95228e+01, -6.27987e+01);
+    refforces[  4] = Vec3(-6.57386e-01,  2.39279e+00,  1.53243e+00);
+    refforces[  5] = Vec3(-1.61426e+00,  7.11728e-01,  9.99752e-01);
+    refforces[  6] = Vec3( 1.60813e+01, -6.77979e+01, -1.01302e+02);
+    refforces[  7] = Vec3(-4.23118e+00,  9.25550e-01,  1.28003e+00);
+    refforces[  8] = Vec3( 1.81589e+00,  1.10083e+00,  1.28615e+00);
+    refforces[  9] = Vec3(-2.24427e+00, -5.71071e+01, -3.18064e+01);
+    refforces[ 10] = Vec3( 3.41715e+00,  1.32918e+00,  1.27908e+00);
+    refforces[ 11] = Vec3(-3.09930e+00,  9.91868e-01,  1.69245e+00);
+    refforces[ 12] = Vec3( 6.08033e+01, -1.02990e+02, -8.26733e+01);
+    refforces[ 13] = Vec3(-1.01346e+00,  1.35475e+00,  4.13725e+00);
+    refforces[ 14] = Vec3(-6.71645e-01,  4.69536e-01,  4.64774e-01);
+    refforces[ 15] = Vec3(-6.39056e+01, -9.85756e-01, -4.97837e+01);
+    refforces[ 16] = Vec3( 9.97113e-01,  1.13195e-01,  4.49387e+00);
+    refforces[ 17] = Vec3( 9.48864e-01,  4.22217e+00,  1.56364e+00);
+    refforces[ 18] = Vec3( 2.18124e+01, -6.14078e+01, -8.10129e+01);
+    refforces[ 19] = Vec3(-5.59701e-01, -1.47528e+00,  1.40157e+00);
+    refforces[ 20] = Vec3( 4.29539e+00,  4.27738e-02,  1.03009e+00);
+    refforces[ 21] = Vec3(-2.62377e+01,  2.55532e+01, -6.89576e+01);
+    refforces[ 22] = Vec3( 1.17222e-01,  3.81618e+00,  1.74502e+00);
+    refforces[ 23] = Vec3( 1.46448e-01,  1.06214e-01,  4.27946e+00);
+    refforces[ 24] = Vec3(-1.90819e+00,  9.86513e+00, -1.05364e+02);
+    refforces[ 25] = Vec3(-1.18537e+00,  1.75021e+00,  1.54448e+00);
+    refforces[ 26] = Vec3( 2.96182e+00, -9.48323e-01,  1.61299e+00);
+    refforces[ 27] = Vec3( 5.61257e+01,  7.69427e+01, -5.84534e+01);
+    refforces[ 28] = Vec3(-8.11821e-01, -7.15455e-01,  6.60049e-01);
+    refforces[ 29] = Vec3(-9.20901e-01,  3.76217e+00,  1.35685e+00);
+    refforces[ 30] = Vec3(-5.03791e+01,  2.24324e+01, -4.31606e+01);
+    refforces[ 31] = Vec3( 3.72530e+00,  3.48275e-01,  1.40060e+00);
+    refforces[ 32] = Vec3( 7.61265e-01,  8.56416e-01,  6.23359e-01);
+    refforces[ 33] = Vec3( 1.95683e+01,  4.16950e+01, -4.33266e+01);
+    refforces[ 34] = Vec3(-2.19677e+00,  1.06287e+00,  1.58260e+00);
+    refforces[ 35] = Vec3( 1.70923e+00,  2.09511e+00,  9.84492e-01);
+    refforces[ 36] = Vec3(-7.84273e+01, -3.34396e+01, -3.66486e+01);
+    refforces[ 37] = Vec3(-1.53567e-01,  1.74934e-01,  1.15266e+00);
+    refforces[ 38] = Vec3(-3.45548e+00,  1.68633e-01,  1.15923e+00);
+    refforces[ 39] = Vec3( 6.90725e+01, -2.09426e+01, -5.31532e+01);
+    refforces[ 40] = Vec3( 1.33274e+00,  3.09504e-01,  1.40065e+00);
+    refforces[ 41] = Vec3( 1.65828e+00, -1.96378e-01,  2.68350e+00);
+    refforces[ 42] = Vec3( 3.39954e+01, -2.65418e+01, -6.81184e+01);
+    refforces[ 43] = Vec3(-1.17491e+00, -2.88793e+00,  1.40541e+00);
+    refforces[ 44] = Vec3(-9.63132e-01, -3.78178e-01,  1.06756e+00);
+    refforces[ 45] = Vec3(-4.44610e+01, -1.81736e+01, -6.19409e+01);
+    refforces[ 46] = Vec3( 1.48529e+00,  8.95548e-06,  4.60732e+00);
+    refforces[ 47] = Vec3( 1.00451e+00, -3.88830e+00,  1.27726e+00);
+    refforces[ 48] = Vec3(-4.23022e+00,  2.52963e+01, -3.04756e+01);
+    refforces[ 49] = Vec3(-3.45228e+00,  7.61191e-01,  9.40731e-01);
+    refforces[ 50] = Vec3( 2.78622e+00,  8.32081e-01,  1.37588e+00);
+    refforces[ 51] = Vec3( 1.63360e+01,  3.02230e+01, -3.91446e+01);
+    refforces[ 52] = Vec3(-2.64761e+00,  1.18509e+00,  9.81393e-01);
+    refforces[ 53] = Vec3( 1.74883e+00, -1.51958e-01,  2.56504e+00);
+    refforces[ 54] = Vec3(-3.01111e+01,  4.08021e+01, -3.45058e+01);
+    refforces[ 55] = Vec3(-1.17129e-01,  9.74282e-02,  4.56402e+00);
+    refforces[ 56] = Vec3(-3.03696e+00,  4.44012e-01,  1.44269e+00);
+    refforces[ 57] = Vec3( 7.57956e+01, -1.12043e+01, -3.92084e+01);
+    refforces[ 58] = Vec3(-1.55617e+00,  2.26444e-01,  2.88234e+00);
+    refforces[ 59] = Vec3(-1.37287e+00, -3.79045e+00,  8.55326e-01);
+    refforces[ 60] = Vec3(-2.88831e+01,  5.72484e+01, -7.17030e+01);
+    refforces[ 61] = Vec3( 6.88458e-01, -2.06343e+00,  1.03552e+00);
+    refforces[ 62] = Vec3( 6.94194e-01, -1.83673e+00,  2.90979e+00);
+    refforces[ 63] = Vec3(-7.87489e+01,  2.17729e+01, -3.72371e+01);
+    refforces[ 64] = Vec3( 1.39267e-01, -4.55150e+00,  9.53893e-01);
+    refforces[ 65] = Vec3(-4.27928e+00, -1.14139e+00,  9.25405e-01);
+    refforces[ 66] = Vec3( 7.13947e+01,  3.02722e+01, -5.16551e+01);
+    refforces[ 67] = Vec3( 4.07933e+00, -1.28199e+00,  9.34707e-01);
+    refforces[ 68] = Vec3(-7.53668e-02, -2.14888e+00,  1.27713e+00);
+    refforces[ 69] = Vec3(-9.36762e+00,  2.98321e+01, -5.70761e+01);
+    refforces[ 70] = Vec3( 1.35989e+00, -3.09988e+00,  9.89945e-01);
+    refforces[ 71] = Vec3(-2.20314e+00, -7.81892e-01,  1.21781e+00);
+    refforces[ 72] = Vec3( 4.89194e+01,  7.40749e+01, -4.22017e+01);
+    refforces[ 73] = Vec3(-4.25696e+00, -1.18954e+00,  1.04383e+00);
+    refforces[ 74] = Vec3(-1.40510e+00, -9.87719e-01,  3.30559e+00);
+    refforces[ 75] = Vec3(-4.95912e+01, -8.24838e+01, -1.82510e+01);
+    refforces[ 76] = Vec3( 7.39933e-01,  9.29483e-01, -1.09883e+00);
+    refforces[ 77] = Vec3( 9.70359e-01,  1.03012e+00,  4.01759e+00);
+    refforces[ 78] = Vec3(-7.87499e+00, -1.87981e+01,  5.77017e+00);
+    refforces[ 79] = Vec3( 3.82847e+00,  1.01191e+00, -7.30715e-02);
+    refforces[ 80] = Vec3(-2.64831e+00,  1.68402e+00,  1.06474e-01);
+    refforces[ 81] = Vec3( 1.81335e+01, -6.72766e+01,  7.21397e+01);
+    refforces[ 82] = Vec3(-1.89188e-01,  1.27803e+00,  8.45103e-01);
+    refforces[ 83] = Vec3( 2.47584e+00,  1.54044e+00,  6.55440e-01);
+    refforces[ 84] = Vec3(-1.05824e+00, -6.50833e+01, -3.88413e+01);
+    refforces[ 85] = Vec3(-3.21762e+00,  1.58809e+00, -1.32879e-01);
+    refforces[ 86] = Vec3( 3.30285e+00,  1.52154e+00, -3.56471e-01);
+    refforces[ 87] = Vec3( 4.41972e+01, -3.58815e+01,  4.01227e+01);
+    refforces[ 88] = Vec3(-1.22923e+00,  8.91826e-01,  2.84505e+00);
+    refforces[ 89] = Vec3(-4.31649e+00,  1.00895e+00,  4.43366e-01);
+    refforces[ 90] = Vec3(-3.65176e+01,  4.86616e+01, -6.44163e+00);
+    refforces[ 91] = Vec3( 1.03989e+00,  3.91910e+00,  7.93373e-01);
+    refforces[ 92] = Vec3( 4.11089e+00, -4.60733e-01, -7.78461e-02);
+    refforces[ 93] = Vec3(-1.58332e+01, -6.82045e+01,  6.62934e+01);
+    refforces[ 94] = Vec3(-7.75645e-02, -3.66897e+00, -3.55044e-01);
+    refforces[ 95] = Vec3( 4.08399e-01,  2.84950e-01,  3.80641e+00);
+    refforces[ 96] = Vec3(-5.05975e+00,  4.20346e+01,  3.49518e+01);
+    refforces[ 97] = Vec3( 9.82673e-01,  2.64451e+00,  4.02971e-01);
+    refforces[ 98] = Vec3(-1.24992e+00, -3.34479e-01,  1.95012e+00);
+    refforces[ 99] = Vec3( 2.31815e+00, -2.30056e+00,  6.49918e+01);
+    refforces[100] = Vec3(-2.16986e+00, -1.05184e-01,  1.89395e+00);
+    refforces[101] = Vec3( 3.51302e+00,  1.83706e-01,  7.08534e-01);
+    refforces[102] = Vec3( 4.46929e+01, -7.11604e+01, -6.62201e+00);
+    refforces[103] = Vec3(-1.87730e+00, -1.10206e+00, -1.47504e+00);
+    refforces[104] = Vec3(-1.41639e+00, -1.36347e+00,  6.10811e-01);
+    refforces[105] = Vec3(-3.12256e+01, -1.73923e+01, -4.47428e+01);
+    refforces[106] = Vec3( 7.92004e-01, -1.39066e+00, -2.68309e+00);
+    refforces[107] = Vec3( 1.54891e+00,  2.50018e+00, -6.21395e-01);
+    refforces[108] = Vec3(-3.83585e+01,  3.90231e+01, -3.68876e+01);
+    refforces[109] = Vec3( 1.86446e-01,  3.49120e+00, -4.00101e-01);
+    refforces[110] = Vec3(-4.05592e+00, -2.21609e-01,  2.13257e-01);
+    refforces[111] = Vec3( 5.70991e+01, -5.95502e+01, -4.08285e+01);
+    refforces[112] = Vec3(-4.25039e-01, -3.96206e-01, -3.23550e+00);
+    refforces[113] = Vec3( 4.14811e+00, -1.47859e-01,  1.11071e-01);
+    refforces[114] = Vec3(-2.74979e+01,  4.40260e+01,  3.06766e+00);
+    refforces[115] = Vec3( 2.22804e-01,  3.48934e+00, -1.96162e-01);
+    refforces[116] = Vec3(-4.21754e+00,  1.18655e-01,  3.36974e-01);
+    refforces[117] = Vec3( 6.32783e+01,  8.02719e+01,  1.38219e+01);
+    refforces[118] = Vec3(-1.49058e+00,  1.87045e+00, -1.32164e-01);
+    refforces[119] = Vec3(-2.30776e+00,  1.73240e-01,  1.92445e+00);
+    refforces[120] = Vec3(-4.04272e+01, -1.77143e+01,  1.57461e+01);
+    refforces[121] = Vec3( 1.82378e+00,  2.60376e+00,  2.33866e-01);
+    refforces[122] = Vec3( 1.33379e+00, -1.16675e+00,  3.34265e+00);
+    refforces[123] = Vec3(-8.78170e+00,  1.35313e+01, -5.40700e+01);
+    refforces[124] = Vec3(-2.72817e+00,  3.27413e-01, -6.95012e-01);
+    refforces[125] = Vec3( 2.41935e+00, -2.66049e-01, -1.62940e+00);
+    refforces[126] = Vec3(-3.63273e+01,  6.91610e+01, -3.72663e+00);
+    refforces[127] = Vec3(-1.40546e+00,  7.48503e-01, -9.73230e-01);
+    refforces[128] = Vec3( 5.39530e-01,  8.56302e-01,  2.57617e+00);
+    refforces[129] = Vec3( 2.54426e+01, -5.04169e+01, -1.76229e+01);
+    refforces[130] = Vec3( 6.82595e-01, -3.17362e+00,  2.20610e-01);
+    refforces[131] = Vec3( 1.94539e+00,  1.53325e+00, -4.24624e-01);
+    refforces[132] = Vec3( 3.61678e+01, -7.89030e+01, -9.61808e+00);
+    refforces[133] = Vec3(-1.29295e+00, -1.55223e+00, -1.88241e+00);
+    refforces[134] = Vec3(-1.51438e+00,  3.14594e-02,  3.60043e+00);
+    refforces[135] = Vec3(-4.60860e+01,  6.18971e+01,  1.38188e+01);
+    refforces[136] = Vec3( 3.07506e+00, -9.60659e-01,  1.81485e+00);
+    refforces[137] = Vec3( 6.46515e-01, -7.30203e-01, -2.17471e-01);
+    refforces[138] = Vec3( 7.76539e+00,  3.18751e+01, -6.18739e+01);
+    refforces[139] = Vec3(-2.12790e+00, -1.06577e+00, -1.03098e+00);
+    refforces[140] = Vec3( 3.70469e+00, -7.70591e-01, -2.00619e-01);
+    refforces[141] = Vec3(-4.02748e+01,  1.93804e+01,  1.92911e+01);
+    refforces[142] = Vec3( 5.74429e-01, -1.58840e+00, -2.07083e+00);
+    refforces[143] = Vec3(-4.21689e+00, -7.25280e-01,  1.61839e-01);
+    refforces[144] = Vec3( 4.10890e+01,  6.50772e+01, -6.64661e+00);
+    refforces[145] = Vec3( 2.02022e-01, -1.46124e+00,  1.99785e+00);
+    refforces[146] = Vec3( 1.91733e-01, -1.14906e+00, -4.18291e+00);
+    refforces[147] = Vec3( 3.77785e+01,  1.01191e+02, -6.07436e+00);
+    refforces[148] = Vec3(-1.03030e+00, -1.26650e+00,  3.57188e+00);
+    refforces[149] = Vec3(-7.42295e-01, -1.05437e+00, -1.47491e+00);
+    refforces[150] = Vec3(-7.40150e+01, -3.53477e+01,  1.81946e+01);
+    refforces[151] = Vec3( 1.65267e+00,  1.39873e+00,  1.46780e+00);
+    refforces[152] = Vec3( 1.18632e+00,  2.66165e+00,  1.61236e+00);
+    refforces[153] = Vec3( 3.89536e+01, -2.96071e+01, -1.38610e+01);
+    refforces[154] = Vec3( 2.05061e+00,  1.54050e+00,  5.35476e-01);
+    refforces[155] = Vec3( 4.98336e-01,  3.85036e+00,  4.21301e-01);
+    refforces[156] = Vec3(-5.53456e+01, -6.87157e+01, -3.55693e+01);
+    refforces[157] = Vec3( 3.17070e+00,  1.65800e+00, -3.99476e-01);
+    refforces[158] = Vec3(-2.02395e+00,  1.48143e+00,  1.00879e+00);
+    refforces[159] = Vec3( 5.94087e+01, -4.57694e+01,  2.77585e+01);
+    refforces[160] = Vec3( 2.05570e+00,  1.46552e+00, -1.41447e+00);
+    refforces[161] = Vec3( 2.97741e-01,  9.12876e-01,  3.27902e+00);
+    refforces[162] = Vec3( 1.93442e+01, -4.16391e+01,  3.34803e-01);
+    refforces[163] = Vec3(-3.93865e+00,  8.46322e-01, -4.47285e-01);
+    refforces[164] = Vec3(-7.34898e-01,  9.82679e-01,  1.44017e+00);
+    refforces[165] = Vec3(-3.34824e+01, -1.38834e+01, -2.62952e+01);
+    refforces[166] = Vec3( 1.26693e+00, -1.10171e+00,  3.22806e+00);
+    refforces[167] = Vec3( 3.49758e+00, -1.84534e-01, -4.62412e-01);
+    refforces[168] = Vec3( 1.82122e+00, -5.69776e+01, -3.86605e+01);
+    refforces[169] = Vec3( 1.92404e+00, -1.05303e+00, -4.52186e-01);
+    refforces[170] = Vec3(-3.25172e+00,  9.44050e-02, -5.87318e-01);
+    refforces[171] = Vec3(-5.98475e+01,  1.85708e+01, -2.67060e+01);
+    refforces[172] = Vec3(-3.77230e-01,  3.29609e+00, -4.81590e-01);
+    refforces[173] = Vec3(-1.66002e+00, -2.41098e+00,  1.03425e-01);
+    refforces[174] = Vec3( 6.47852e+01,  4.38661e+00,  1.77482e+01);
+    refforces[175] = Vec3(-1.29162e-01, -3.40031e+00,  3.91521e-01);
+    refforces[176] = Vec3( 1.28253e+00,  2.14633e+00,  2.64792e-01);
+    refforces[177] = Vec3( 2.95010e+01, -5.60376e+00,  4.16733e+01);
+    refforces[178] = Vec3(-1.43378e+00, -1.96178e-01,  2.40594e+00);
+    refforces[179] = Vec3(-1.56255e+00, -3.08967e+00, -3.14335e-01);
+    refforces[180] = Vec3(-6.02384e+01, -1.59585e+01,  7.16798e+01);
+    refforces[181] = Vec3( 1.68669e+00, -3.54228e+00, -2.98122e-01);
+    refforces[182] = Vec3( 9.46412e-01, -8.03852e-02,  4.17305e+00);
+    refforces[183] = Vec3( 3.87313e+00,  4.46632e+01,  4.59130e+01);
+    refforces[184] = Vec3(-9.83299e-02, -8.11669e-02,  4.04515e+00);
+    refforces[185] = Vec3( 6.64172e-01,  3.66814e+00, -2.58278e-01);
+    refforces[186] = Vec3(-1.65450e+01,  2.35336e+01,  7.09485e+01);
+    refforces[187] = Vec3(-6.98821e-01,  2.45516e-01,  2.69560e+00);
+    refforces[188] = Vec3( 4.04145e+00,  5.46258e-02, -1.82210e-02);
+    refforces[189] = Vec3( 2.39311e+01,  1.22785e+01, -2.13260e+01);
+    refforces[190] = Vec3(-3.33802e-02, -3.60205e+00,  6.21747e-01);
+    refforces[191] = Vec3( 3.50616e+00,  6.09623e-01, -7.40284e-02);
+    refforces[192] = Vec3( 2.32693e+01, -1.88836e+01,  1.24193e+01);
+    refforces[193] = Vec3(-2.25948e+00, -1.74718e+00,  2.65705e-01);
+    refforces[194] = Vec3(-1.32949e+00, -8.67446e-01,  7.75386e-01);
+    refforces[195] = Vec3(-4.40626e+01,  1.17032e+01, -4.18284e+01);
+    refforces[196] = Vec3( 1.61723e+00, -6.18309e-01,  2.45312e+00);
+    refforces[197] = Vec3( 1.42535e+00,  1.71801e+00, -1.86477e+00);
+    refforces[198] = Vec3( 3.37157e+00, -1.65818e+01,  7.59773e+01);
+    refforces[199] = Vec3(-2.81989e+00, -1.37877e+00, -9.60015e-02);
+    refforces[200] = Vec3(-8.77602e-02,  2.44299e+00,  1.17362e+00);
+    refforces[201] = Vec3( 2.44644e+01, -3.26831e+01, -5.34506e+01);
+    refforces[202] = Vec3(-4.14853e+00, -9.48873e-02, -2.18126e-01);
+    refforces[203] = Vec3( 8.13810e-01, -1.30429e+00, -1.62539e+00);
+    refforces[204] = Vec3( 1.03414e+01, -4.14898e+00,  2.56841e+01);
+    refforces[205] = Vec3(-1.50990e-01, -3.90514e+00,  2.98649e-01);
+    refforces[206] = Vec3(-3.60035e+00, -7.83546e-02, -7.95300e-01);
+    refforces[207] = Vec3( 4.49059e+01,  2.88959e+01,  1.53626e+01);
+    refforces[208] = Vec3(-2.55855e+00, -1.55231e+00, -3.96734e-01);
+    refforces[209] = Vec3(-1.14072e+00,  3.25255e+00,  2.39860e-01);
+    refforces[210] = Vec3(-4.12850e+01,  6.11654e+01,  2.43179e+01);
+    refforces[211] = Vec3( 1.20781e+00, -8.28196e-01,  4.26212e+00);
+    refforces[212] = Vec3( 7.30387e-01, -9.26965e-01, -9.32605e-01);
+    refforces[213] = Vec3(-1.64972e+01,  4.58667e+01,  5.58154e+01);
+    refforces[214] = Vec3( 4.43484e-01, -1.67205e+00,  3.42164e+00);
+    refforces[215] = Vec3(-4.14170e+00, -9.68400e-01, -9.21231e-02);
+    refforces[216] = Vec3( 1.79252e+01,  6.24205e+01, -1.47777e+01);
+    refforces[217] = Vec3(-2.53437e+00, -9.16373e-01, -1.63636e+00);
+    refforces[218] = Vec3( 3.32247e+00, -1.29776e+00, -1.07427e+00);
+    refforces[219] = Vec3( 1.79752e+01,  5.16430e+01,  3.90249e+01);
+    refforces[220] = Vec3( 7.14180e-01, -1.00018e+00,  3.14118e+00);
+    refforces[221] = Vec3(-1.07965e+00, -3.41285e+00, -2.24161e-01);
+    refforces[222] = Vec3( 3.16317e+01,  4.98025e+01, -4.17518e+01);
+    refforces[223] = Vec3(-2.55899e+00, -1.63433e+00, -3.76149e-01);
+    refforces[224] = Vec3(-1.19443e+00, -1.02227e+00, -5.37349e-02);
+    refforces[225] = Vec3(-3.77249e+01, -4.85845e+01, -2.11073e+01);
+    refforces[226] = Vec3( 7.98995e-01,  1.26056e+00,  4.28442e+00);
+    refforces[227] = Vec3( 3.35648e+00,  9.11653e-01, -8.07025e-01);
+    refforces[228] = Vec3(-3.16004e+01, -2.04156e+01,  5.01542e+01);
+    refforces[229] = Vec3( 8.97766e-01,  1.34997e+00,  2.86585e+00);
+    refforces[230] = Vec3(-2.09603e-01,  4.05651e+00, -1.91043e-02);
+    refforces[231] = Vec3( 3.61582e+01, -3.35593e+01,  4.19940e+01);
+    refforces[232] = Vec3( 3.66075e+00,  7.71021e-01,  3.45598e-01);
+    refforces[233] = Vec3(-7.50361e-01,  1.22549e+00,  7.98677e-01);
+    refforces[234] = Vec3(-4.42511e+01, -2.71049e+01, -4.56576e+01);
+    refforces[235] = Vec3(-3.30181e-01,  1.98824e+00,  2.21005e+00);
+    refforces[236] = Vec3(-7.68719e-01,  9.09063e-01, -3.77633e+00);
+    refforces[237] = Vec3( 7.33817e+01, -3.03886e+01, -7.59758e+01);
+    refforces[238] = Vec3(-1.71701e+00,  1.41048e+00, -1.84182e+00);
+    refforces[239] = Vec3(-1.35491e+00,  2.27816e+00, -1.08799e+00);
+    refforces[240] = Vec3(-4.58218e+01, -2.91931e+01,  1.12369e+01);
+    refforces[241] = Vec3( 1.68546e+00,  2.81258e+00,  3.93755e-01);
+    refforces[242] = Vec3( 8.28917e-01, -1.32237e+00,  2.91182e+00);
+    refforces[243] = Vec3( 5.44446e+00, -3.99466e+01,  6.77923e+00);
+    refforces[244] = Vec3(-1.09231e+00, -5.21927e-01,  1.83581e+00);
+    refforces[245] = Vec3( 5.24849e-02, -7.21832e-01, -3.35896e+00);
+    refforces[246] = Vec3(-1.29863e+01, -3.94851e+01, -2.70347e+01);
+    refforces[247] = Vec3(-3.94725e+00, -2.01688e-01,  5.93703e-02);
+    refforces[248] = Vec3( 1.30447e+00, -2.35362e+00, -2.94278e-01);
+    refforces[249] = Vec3( 2.56542e+01, -3.79434e+01, -6.19782e+01);
+    refforces[250] = Vec3( 2.48025e+00, -1.19994e+00,  1.41640e-01);
+    refforces[251] = Vec3(-5.01379e-01,  2.48607e-01, -3.97535e+00);
+    refforces[252] = Vec3( 3.46461e+01, -8.37214e+01, -3.92225e+01);
+    refforces[253] = Vec3(-2.80216e+00, -1.13524e+00,  1.85582e-02);
+    refforces[254] = Vec3(-1.28421e+00, -9.50174e-01, -5.56828e-01);
+    refforces[255] = Vec3(-3.56755e+01,  1.44176e+01, -7.63345e-01);
+    refforces[256] = Vec3( 1.25547e+00,  2.64846e+00,  1.29226e+00);
+    refforces[257] = Vec3( 1.22691e+00,  5.72301e-01, -4.09853e+00);
+    refforces[258] = Vec3(-6.02997e+00, -1.31579e+01, -5.03304e+01);
+    refforces[259] = Vec3(-4.59947e-01, -8.09031e-01,  3.08383e+00);
+    refforces[260] = Vec3(-8.93872e-01, -3.70797e-01, -2.89773e+00);
+    refforces[261] = Vec3(-1.92458e+00,  2.50787e+01, -4.30937e+01);
+    refforces[262] = Vec3( 2.10381e-01, -1.24125e-01, -4.01554e+00);
+    refforces[263] = Vec3(-3.85853e+00,  6.72311e-01,  4.06873e-02);
+    refforces[264] = Vec3( 1.66952e+01,  1.47383e+01,  4.13458e+01);
+    refforces[265] = Vec3(-2.12261e+00, -5.40916e-01,  1.54532e+00);
+    refforces[266] = Vec3( 3.54041e+00, -4.73740e-02,  4.73489e-01);
+    refforces[267] = Vec3( 4.43789e+01,  6.96271e+01,  4.06335e+00);
+    refforces[268] = Vec3(-9.59563e-01,  9.94730e-01, -3.12464e+00);
+    refforces[269] = Vec3(-1.42363e+00,  1.62107e+00,  1.87365e+00);
+    refforces[270] = Vec3(-5.38380e+01,  2.27365e+01,  1.95666e+01);
+    refforces[271] = Vec3( 1.47405e+00,  2.11523e+00,  1.97777e-01);
+    refforces[272] = Vec3( 2.54377e+00, -4.04439e-02,  1.97927e+00);
+    refforces[273] = Vec3( 1.33568e+01,  1.88750e+01, -3.59233e+01);
+    refforces[274] = Vec3(-1.09293e+00,  2.25841e+00, -4.61222e-01);
+    refforces[275] = Vec3( 3.27762e+00, -1.71208e-01, -8.67824e-01);
+    refforces[276] = Vec3(-3.41864e+01,  3.38579e+01,  5.95219e+01);
+    refforces[277] = Vec3(-4.18249e+00, -8.81154e-02,  1.97498e-01);
+    refforces[278] = Vec3( 1.06177e+00,  2.34089e+00,  2.42037e-01);
+    refforces[279] = Vec3( 2.35073e+01,  1.61625e+01, -2.58858e+01);
+    refforces[280] = Vec3( 3.87903e+00,  1.67313e-01, -4.03473e-02);
+    refforces[281] = Vec3(-5.12341e-01, -5.62857e-01, -3.19237e+00);
+    refforces[282] = Vec3( 3.39154e+01, -5.09470e+00, -2.29146e+01);
+    refforces[283] = Vec3(-1.24992e+00,  5.24960e-01, -1.16022e+00);
+    refforces[284] = Vec3(-2.39537e+00,  1.82818e-01, -1.53421e+00);
+    refforces[285] = Vec3(-5.56195e+01,  2.89266e+01, -2.61137e+01);
+    refforces[286] = Vec3( 2.49596e+00, -1.03213e+00, -1.51258e+00);
+    refforces[287] = Vec3( 1.02471e+00, -1.05342e+00, -6.79809e-01);
+    refforces[288] = Vec3(-1.52895e+01,  4.83734e+01, -7.33233e+00);
+    refforces[289] = Vec3( 3.77837e+00, -1.59224e+00,  2.02675e-01);
+    refforces[290] = Vec3(-3.51854e+00, -1.54234e+00,  2.89348e-01);
+    refforces[291] = Vec3( 3.22262e+01,  2.32542e+01,  2.37269e+01);
+    refforces[292] = Vec3( 1.88520e+00, -1.38396e+00, -1.98750e+00);
+    refforces[293] = Vec3( 4.99036e-02, -4.53129e+00, -2.13859e-01);
+    refforces[294] = Vec3(-7.48274e+00,  3.71318e+01, -2.35667e+01);
+    refforces[295] = Vec3( 7.81671e-01, -1.07371e+00, -2.37009e+00);
+    refforces[296] = Vec3(-3.98896e+00, -9.36146e-01,  2.62843e-02);
+    refforces[297] = Vec3( 4.33081e+01,  4.24806e+01,  3.21190e+01);
+    refforces[298] = Vec3(-2.52531e+00, -1.90285e+00,  2.72091e-01);
+    refforces[299] = Vec3(-9.70927e-01, -1.13219e+00,  1.47351e+00);
+    refforces[300] = Vec3(-4.71070e+01, -7.19253e+01,  5.82986e+01);
+    refforces[301] = Vec3( 6.82929e-01,  1.00679e+00, -6.94522e-01);
+    refforces[302] = Vec3( 8.31359e-01,  7.09804e-01, -8.12652e-01);
+    refforces[303] = Vec3( 3.62309e+00, -3.42931e+01,  1.98231e+01);
+    refforces[304] = Vec3( 1.23138e+00,  8.49297e-01, -2.39283e+00);
+    refforces[305] = Vec3(-2.87428e+00,  1.41050e+00, -1.03948e+00);
+    refforces[306] = Vec3(-3.83149e+00, -2.61638e+01,  2.98550e+01);
+    refforces[307] = Vec3(-2.90010e+00,  1.38356e+00, -1.22906e+00);
+    refforces[308] = Vec3( 2.62337e-01,  9.72795e-01, -4.02660e+00);
+    refforces[309] = Vec3( 7.84655e+00, -5.09531e+01,  7.96762e+01);
+    refforces[310] = Vec3( 2.46967e+00,  9.48901e-01, -1.17317e+00);
+    refforces[311] = Vec3(-4.11596e+00,  7.95633e-01, -1.39335e+00);
+    refforces[312] = Vec3( 5.03705e+01, -8.09993e+01,  1.13172e+02);
+    refforces[313] = Vec3(-4.48925e-01,  4.29309e-01, -5.50572e-01);
+    refforces[314] = Vec3(-2.66201e+00,  1.05662e+00, -1.44247e+00);
+    refforces[315] = Vec3(-6.12660e+01,  2.31218e+01,  4.57287e+01);
+    refforces[316] = Vec3( 8.42017e-01,  1.41042e+00, -7.78795e-01);
+    refforces[317] = Vec3( 3.11298e+00, -8.04273e-01, -1.51273e+00);
+    refforces[318] = Vec3( 1.83823e+01, -2.83588e+01,  6.13597e+01);
+    refforces[319] = Vec3( 4.37781e+00,  8.12091e-02, -1.33992e+00);
+    refforces[320] = Vec3(-5.81396e-01, -3.61991e+00, -1.22489e+00);
+    refforces[321] = Vec3(-4.85270e+01, -6.50314e+01,  7.16441e+01);
+    refforces[322] = Vec3(-3.40093e-02,  9.01514e-01, -1.26573e+00);
+    refforces[323] = Vec3(-8.67193e-02, -3.94828e+00, -8.97127e-01);
+    refforces[324] = Vec3( 6.25833e+01, -2.31957e+01,  8.66089e+01);
+    refforces[325] = Vec3( 2.49809e+00, -7.64197e-01, -1.27775e+00);
+    refforces[326] = Vec3( 2.11019e-01,  7.81764e-01, -1.23664e+00);
+    refforces[327] = Vec3( 2.04801e+01,  3.14002e+01,  6.96727e+01);
+    refforces[328] = Vec3(-7.99932e-01,  4.35586e+00, -8.77957e-01);
+    refforces[329] = Vec3(-2.87129e+00, -1.02321e+00, -1.42702e+00);
+    refforces[330] = Vec3(-6.36347e+01, -3.58428e+01,  4.52682e+01);
+    refforces[331] = Vec3( 1.32210e+00, -2.58500e+00, -1.43327e+00);
+    refforces[332] = Vec3( 1.52218e+00,  1.63513e+00, -2.29784e+00);
+    refforces[333] = Vec3( 4.36105e+00,  3.70776e+01,  6.46133e+01);
+    refforces[334] = Vec3(-1.19452e+00, -1.01252e+00, -1.43172e+00);
+    refforces[335] = Vec3( 4.34983e+00, -6.52342e-02, -1.33662e+00);
+    refforces[336] = Vec3( 1.13050e+01,  4.99187e+01,  7.16119e+01);
+    refforces[337] = Vec3( 2.29221e-01,  6.35843e-01, -3.59259e+00);
+    refforces[338] = Vec3( 1.65239e+00,  4.20995e-01, -1.42878e+00);
+    refforces[339] = Vec3(-1.79210e+01,  1.92423e+01,  1.77516e+01);
+    refforces[340] = Vec3(-2.91278e+00,  9.22214e-01, -8.41281e-01);
+    refforces[341] = Vec3( 2.96465e+00,  1.05962e+00, -6.85280e-01);
+    refforces[342] = Vec3( 5.19378e+01, -2.55266e+01,  3.63518e+01);
+    refforces[343] = Vec3(-1.13088e+00,  4.54600e-01, -1.30862e+00);
+    refforces[344] = Vec3(-2.05369e+00, -9.69459e-03, -2.69351e+00);
+    refforces[345] = Vec3(-4.22214e+01,  4.88282e+01,  3.76741e+01);
+    refforces[346] = Vec3( 3.46996e+00,  9.57280e-01, -7.60128e-01);
+    refforces[347] = Vec3( 7.62585e-01,  1.14548e+00, -8.06697e-01);
+    refforces[348] = Vec3(-2.35024e+01, -6.73585e+01,  6.49050e+01);
+    refforces[349] = Vec3(-1.58984e+00, -6.44551e-01, -2.37758e+00);
+    refforces[350] = Vec3( 1.32604e+00, -2.37576e+00, -1.54939e+00);
+    refforces[351] = Vec3( 2.57776e+01, -1.40743e+01,  2.06580e+01);
+    refforces[352] = Vec3( 4.24916e+00,  6.10665e-02, -8.19360e-01);
+    refforces[353] = Vec3(-7.94423e-01, -1.75764e-01, -3.51414e+00);
+    refforces[354] = Vec3(-1.24260e+01,  6.96334e-01,  5.97705e+01);
+    refforces[355] = Vec3(-3.38899e+00,  9.44919e-01, -1.25780e+00);
+    refforces[356] = Vec3( 6.81664e-01,  1.16965e+00, -2.34256e+00);
+    refforces[357] = Vec3( 5.15292e+01,  5.13086e+01,  6.16724e+01);
+    refforces[358] = Vec3(-1.04216e+00,  3.34685e+00, -1.33787e+00);
+    refforces[359] = Vec3(-9.57017e-01, -1.59610e-01, -4.31963e+00);
+    refforces[360] = Vec3(-5.88173e+01,  3.30524e+01,  6.23031e+01);
+    refforces[361] = Vec3( 9.35413e-01, -1.09829e+00, -4.65202e+00);
+    refforces[362] = Vec3( 6.95250e-01, -1.24375e+00, -6.50832e-01);
+    refforces[363] = Vec3( 1.64454e+01,  1.08071e+02,  3.76968e+01);
+    refforces[364] = Vec3(-4.25249e+00, -1.10053e+00, -7.56357e-01);
+    refforces[365] = Vec3( 1.17898e+00, -1.10215e+00, -8.21366e-01);
+    refforces[366] = Vec3(-5.20527e+00,  5.49522e+01,  4.55303e+01);
+    refforces[367] = Vec3(-7.98794e-01, -2.78246e+00, -1.27815e+00);
+    refforces[368] = Vec3(-1.02109e+00, -1.21921e+00, -1.26706e+00);
+    refforces[369] = Vec3( 3.05185e+00,  3.99846e+01,  5.39902e+01);
+    refforces[370] = Vec3( 1.93675e+00, -1.28986e+00, -1.21655e+00);
+    refforces[371] = Vec3(-4.14952e+00, -1.06509e+00, -9.52154e-01);
+    refforces[372] = Vec3( 6.06000e+01,  2.50872e+01,  4.29623e+01);
+    refforces[373] = Vec3(-3.50284e+00, -1.47848e+00, -7.98049e-01);
+    refforces[374] = Vec3(-8.66087e-01, -1.73056e+00, -6.59554e-01);
+
+    // Find the exclusion information
+    vector<set<int> > exclusions(numAtoms);
+    for (int i = 0; i < forceField->getNumExceptions(); i++) {
+        int particle1, particle2;
+        double chargeProd, sigma, epsilon;
+        forceField->getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
+        exclusions[particle1].insert(particle2);
+        exclusions[particle2].insert(particle1);
+    }
+
+    const double longcutoff = 35.0*OpenMM::NmPerAngstrom;
+    const double longcutoff2 = longcutoff*longcutoff;
+    const double cutoff2 = cutoff*cutoff;
+    const double cutoff6inv = 1.0 / (cutoff2*cutoff2*cutoff2);
+    const int nboxes = ceil(longcutoff/boxEdgeLength);
+
+    double refenergy = 0.0;
+    // Loop over home box first...
+    for (int i = 0; i < numAtoms; ++ i) {
+        for (int j = i+1; j < numAtoms; ++j) {
+            Vec3 dR = positions[i] - positions[j];
+            double R2 = dR[0]*dR[0] + dR[1]*dR[1] + dR[2]*dR[2];
+            double sig2 = (sigvals[i] + sigvals[j])*(sigvals[i] + sigvals[j]) / R2;
+            double sig6 = sig2*sig2*sig2;
+            double eps = epsvals[i]*epsvals[j];
+            refenergy += 2.0*eps*(sig6-1.0)*sig6;
+            if (R2 < cutoff2) {
+                // Add a shift term for direct space parts withing t
+                refenergy += 2.0*eps*(pow(sigvals[i]+sigvals[j], 6) - 64.0*pow(sigvals[i]*sigvals[j], 3))*cutoff6inv;
+            }
+        }
+    }
+
+    // ... back out exclusions
+    for (int ii = 0; ii < numAtoms; ii++) {
+        for (set<int>::const_iterator iter = exclusions[ii].begin(); iter != exclusions[ii].end(); ++iter) {
+            if (*iter > ii) {
+                int i = ii;
+                int j = *iter;
+                Vec3 dR = positions[i] - positions[j];
+                double R2 = dR[0]*dR[0] + dR[1]*dR[1] + dR[2]*dR[2];
+                double sig2 = (sigvals[i] + sigvals[j])*(sigvals[i] + sigvals[j]) / R2;
+                double sig6 = sig2*sig2*sig2;
+                double eps = epsvals[i]*epsvals[j];
+                refenergy -= 2.0*eps*(sig6-1.0)*sig6;
+                if (R2 < cutoff2) {
+                    // Add a shift term for direct space parts withing t
+                    refenergy -= 2.0*eps*(pow(sigvals[i]+sigvals[j], 6) - 64.0*pow(sigvals[i]*sigvals[j], 3))*cutoff6inv;
+                }
+            }
+        }
+    }
+
+    // ... and now add in the image box terms
+    for (int bx = -nboxes; bx <= nboxes; ++bx) {
+        for (int by = -nboxes; by <= nboxes; ++by) {
+            for (int bz = -nboxes; bz <= nboxes; ++bz) {
+                if (bx==0 && by==0 && bz==0) continue;
+                Vec3 offset(bx*boxEdgeLength, by*boxEdgeLength, bz*boxEdgeLength);
+                for (int i = 0; i < numAtoms; ++ i) {
+                    for (int j = 0; j < numAtoms; ++j) {
+                        Vec3 dR = positions[i] - positions[j] + offset;
+                        double R2 = dR[0]*dR[0] + dR[1]*dR[1] + dR[2]*dR[2];
+                        if (R2 > longcutoff2) continue;
+                        double sig2 = (sigvals[i] + sigvals[j])*(sigvals[i] + sigvals[j]) / R2;
+                        double sig6 = sig2*sig2*sig2;
+                        double eps = epsvals[i]*epsvals[j];
+                        refenergy += eps*(sig6-1.0)*sig6;
+                        if (R2 < cutoff2) {
+                            // Add a shift term for direct space parts withing teh
+                            refenergy += eps*(pow(sigvals[i]+sigvals[j], 6) - 64.0*pow(sigvals[i]*sigvals[j], 3))*cutoff6inv;
+                        }
+                    }
+                }
+            }
+        }
+    }
+    refenergy *= 0.5;
+    // For this test the reference energy is -294.078 kJ/mol, while the difference between DPME and 30A cutoffs
+    // is just 0.064 kJ/mol.  The difference is due to the fact that arithmetic mean combination rules are used
+    // up to the cutoff, while the reciprocal space uses the geometric mean.  See DOI: 10.1021/acs.jctc.5b00726
+    ASSERT_EQUAL_TOL(refenergy, energy, 5E-4);
+    ASSERT_EQUAL_TOL(gromacs_energy, energy, 5E-5);
+
+    // Forces accumulated in single precision are tested to a more permissive criterion; the double
+    // precision platform can match to 5E-5.
+    for (int n = 0; n < numAtoms; ++n)
+        ASSERT_EQUAL_VEC(refforces[n], forces[n], 1E-4);
+
+}
+
+void runPlatformTests();
+
+int main(int argc, char* argv[]) {
+    try {
+        initializeTests(argc, argv);
+        testErrorTolerance();
+        testPMEParameters();
+        testWater2DpmeEnergiesForcesNoExclusions();
+        testWater2DpmeEnergiesForcesWithExclusions();
+        testWater125DpmeVsLongCutoffNoExclusions();
+        testWater125DpmeVsLongCutoffWithExclusions();
+        runPlatformTests();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}

tests/TestNonbondedForce.h

@@ -35,7 +35,6 @@
 #include "openmm/HarmonicBondForce.h"
 #include "openmm/NonbondedForce.h"
 #include "openmm/System.h"
-#include "openmm/Units.h"
 #include "openmm/VerletIntegrator.h"
 #include "SimTKOpenMMRealType.h"
 #include "sfmt/SFMT.h"
@@ -48,1661 +47,6 @@ using namespace std;
 
 const double TOL = 1e-5;
 
-void make_waterbox(int natoms, double boxEdgeLength, NonbondedForce *forceField,  vector<Vec3> &positions, vector<double>& eps, vector<double>& sig,
-                   vector<pair<int, int> >& bonds, System &system, bool do_electrostatics) {
-    const int RESSIZE = 3;
-    const double masses[RESSIZE]    = {     8.0,    1.0,    1.0 };
-    const double charges[RESSIZE]   = {  -0.834,  0.417,  0.417 };
-    // Values from the CHARMM force field, in AKMA units
-    const double epsilons[RESSIZE]  = { -0.1521, -0.046, -0.046 };
-    const double halfrmins[RESSIZE] = {  1.7682, 0.2245, 0.2245 };
-    positions.clear();
-    if(natoms == 6){
-        const double coords[6][3] = {
-            {  2.000000, 2.000000, 2.000000},
-            {  2.500000, 2.000000, 3.000000},
-            {  1.500000, 2.000000, 3.000000},
-            {  0.000000, 0.000000, 0.000000},
-            {  0.500000, 0.000000, 1.000000},
-            { -0.500000, 0.000000, 1.000000}
-        };
-        for(int atom = 0; atom < natoms; ++atom)
-            positions.push_back(Vec3(coords[atom][0], coords[atom][1], coords[atom][2])*OpenMM::NmPerAngstrom);
-    }else if(natoms == 375){
-        const double coords[375][3] = {
-            { -6.22, -6.25, -6.24 },
-            { -5.32, -6.03, -6.00 },
-            { -6.75, -5.56, -5.84 },
-            { -3.04, -6.23, -6.19 },
-            { -3.52, -5.55, -5.71 },
-            { -3.59, -6.43, -6.94 },
-            {  0.02, -6.23, -6.14 },
-            { -0.87, -5.97, -6.37 },
-            {  0.53, -6.03, -6.93 },
-            {  3.10, -6.20, -6.27 },
-            {  3.87, -6.35, -5.72 },
-            {  2.37, -6.11, -5.64 },
-            {  6.18, -6.14, -6.20 },
-            {  6.46, -6.66, -5.44 },
-            {  6.26, -6.74, -6.94 },
-            { -6.21, -3.15, -6.24 },
-            { -6.23, -3.07, -5.28 },
-            { -6.02, -2.26, -6.55 },
-            { -3.14, -3.07, -6.16 },
-            { -3.38, -3.63, -6.90 },
-            { -2.18, -3.05, -6.17 },
-            { -0.00, -3.16, -6.23 },
-            { -0.03, -2.30, -6.67 },
-            {  0.05, -2.95, -5.29 },
-            {  3.08, -3.11, -6.14 },
-            {  2.65, -2.55, -6.79 },
-            {  3.80, -3.53, -6.62 },
-            {  6.16, -3.14, -6.16 },
-            {  7.04, -3.32, -6.51 },
-            {  5.95, -2.27, -6.51 },
-            { -6.20, -0.04, -6.15 },
-            { -5.43,  0.32, -6.59 },
-            { -6.95,  0.33, -6.62 },
-            { -3.10, -0.06, -6.19 },
-            { -3.75,  0.42, -6.69 },
-            { -2.46,  0.60, -5.93 },
-            {  0.05, -0.01, -6.17 },
-            { -0.10,  0.02, -7.12 },
-            { -0.79,  0.16, -5.77 },
-            {  3.03,  0.00, -6.19 },
-            {  3.54,  0.08, -7.01 },
-            {  3.69, -0.22, -5.53 },
-            {  6.17,  0.05, -6.19 },
-            {  5.78, -0.73, -6.57 },
-            {  7.09, -0.17, -6.04 },
-            { -6.20,  3.15, -6.25 },
-            { -6.59,  3.18, -5.37 },
-            { -5.87,  2.25, -6.33 },
-            { -3.09,  3.04, -6.17 },
-            { -3.88,  3.58, -6.26 },
-            { -2.41,  3.54, -6.63 },
-            {  0.00,  3.06, -6.26 },
-            { -0.71,  3.64, -6.00 },
-            {  0.65,  3.15, -5.55 },
-            {  3.14,  3.06, -6.23 },
-            {  3.11,  3.31, -5.30 },
-            {  2.38,  3.49, -6.63 },
-            {  6.19,  3.14, -6.25 },
-            {  6.82,  3.25, -5.54 },
-            {  5.76,  2.30, -6.07 },
-            { -6.22,  6.26, -6.19 },
-            { -6.22,  5.74, -7.00 },
-            { -5.89,  5.67, -5.52 },
-            { -3.04,  6.24, -6.20 },
-            { -3.08,  5.28, -6.17 },
-            { -3.96,  6.52, -6.25 },
-            { -0.05,  6.21, -6.16 },
-            {  0.82,  6.58, -6.06 },
-            {  0.01,  5.64, -6.93 },
-            {  3.10,  6.25, -6.15 },
-            {  3.64,  5.47, -6.31 },
-            {  2.46,  6.24, -6.87 },
-            {  6.22,  6.20, -6.27 },
-            {  5.37,  6.42, -5.88 },
-            {  6.80,  6.07, -5.51 },
-            { -6.19, -6.15, -3.13 },
-            { -6.37, -7.01, -3.51 },
-            { -6.25, -6.29, -2.18 },
-            { -3.10, -6.27, -3.11 },
-            { -2.29, -5.77, -2.99 },
-            { -3.80, -5.62, -2.98 },
-            { -0.03, -6.18, -3.15 },
-            { -0.07, -7.05, -2.75 },
-            {  0.68, -5.74, -2.70 },
-            {  3.10, -6.14, -3.07 },
-            {  2.35, -6.72, -3.23 },
-            {  3.86, -6.65, -3.37 },
-            {  6.22, -6.20, -3.16 },
-            {  6.82, -6.36, -2.43 },
-            {  5.35, -6.13, -2.75 },
-            { -6.26, -3.13, -3.12 },
-            { -6.16, -2.27, -2.70 },
-            { -5.36, -3.47, -3.18 },
-            { -3.11, -3.05, -3.14 },
-            { -3.31, -3.96, -3.34 },
-            { -2.77, -3.06, -2.24 },
-            {  0.00, -3.13, -3.16 },
-            {  0.48, -2.37, -2.81 },
-            { -0.57, -3.40, -2.44 },
-            {  3.09, -3.09, -3.16 },
-            {  2.41, -3.19, -2.49 },
-            {  3.91, -3.07, -2.67 },
-            {  6.19, -3.04, -3.08 },
-            {  5.64, -3.61, -3.61 },
-            {  6.93, -3.58, -2.82 },
-            { -6.18, -0.00, -3.04 },
-            { -6.00, -0.59, -3.78 },
-            { -6.79,  0.64, -3.39 },
-            { -3.05, -0.03, -3.07 },
-            { -2.95,  0.80, -3.52 },
-            { -4.00, -0.20, -3.07 },
-            { -0.03,  0.03, -3.06 },
-            { -0.33, -0.37, -3.87 },
-            {  0.89, -0.21, -2.99 },
-            {  3.13, -0.05, -3.10 },
-            {  3.44,  0.81, -3.34 },
-            {  2.21,  0.07, -2.86 },
-            {  6.20, -0.05, -3.13 },
-            {  6.89,  0.60, -3.20 },
-            {  5.58,  0.30, -2.49 },
-            { -6.23,  3.09, -3.16 },
-            { -5.62,  3.79, -2.94 },
-            { -6.33,  2.60, -2.33 },
-            { -3.10,  3.08, -3.04 },
-            { -3.84,  3.47, -3.51 },
-            { -2.40,  3.01, -3.69 },
-            {  0.01,  3.04, -3.11 },
-            { -0.56,  3.59, -3.64 },
-            {  0.28,  3.60, -2.38 },
-            {  3.04,  3.11, -3.09 },
-            {  3.49,  2.30, -2.87 },
-            {  3.70,  3.66, -3.51 },
-            {  6.15,  3.14, -3.11 },
-            {  6.52,  2.52, -3.74 },
-            {  6.72,  3.06, -2.34 },
-            { -6.22,  6.15, -3.13 },
-            { -5.49,  6.21, -2.51 },
-            { -6.56,  7.04, -3.18 },
-            { -3.11,  6.24, -3.05 },
-            { -3.76,  5.83, -3.62 },
-            { -2.26,  5.92, -3.37 },
-            {  0.03,  6.25, -3.07 },
-            {  0.34,  5.63, -3.73 },
-            { -0.87,  6.00, -2.91 },
-            {  3.07,  6.15, -3.08 },
-            {  3.29,  6.92, -2.56 },
-            {  3.39,  6.35, -3.96 },
-            {  6.22,  6.14, -3.12 },
-            {  5.79,  6.38, -2.29 },
-            {  6.25,  6.96, -3.62 },
-            { -6.21, -6.20, -0.06 },
-            { -5.79, -6.87,  0.48 },
-            { -6.43, -5.50,  0.54 },
-            { -3.16, -6.21, -0.02 },
-            { -2.50, -6.87,  0.20 },
-            { -2.77, -5.37,  0.23 },
-            { -0.00, -6.14, -0.00 },
-            {  0.68, -6.72, -0.33 },
-            { -0.64, -6.73,  0.38 },
-            {  3.03, -6.20, -0.01 },
-            {  3.77, -6.56, -0.51 },
-            {  3.43, -5.85,  0.78 },
-            {  6.25, -6.16, -0.00 },
-            {  5.36, -6.09, -0.36 },
-            {  6.24, -6.97,  0.49 },
-            { -6.24, -3.05, -0.01 },
-            { -6.35, -3.64,  0.73 },
-            { -5.42, -3.33, -0.42 },
-            { -3.09, -3.06,  0.05 },
-            { -2.44, -3.62, -0.38 },
-            { -3.90, -3.21, -0.43 },
-            {  0.05, -3.10,  0.02 },
-            { -0.31, -2.35, -0.43 },
-            { -0.63, -3.77,  0.01 },
-            {  3.05, -3.09, -0.04 },
-            {  3.28, -3.90,  0.41 },
-            {  3.65, -2.43,  0.30 },
-            {  6.20, -3.04, -0.03 },
-            {  5.66, -3.31,  0.71 },
-            {  6.78, -3.79, -0.19 },
-            { -6.18,  0.04, -0.04 },
-            { -6.73, -0.73, -0.15 },
-            { -5.98,  0.06,  0.89 },
-            { -3.11, -0.04, -0.04 },
-            { -3.36, -0.08,  0.87 },
-            { -2.70,  0.81, -0.14 },
-            { -0.02, -0.02, -0.05 },
-            { -0.45,  0.28,  0.75 },
-            {  0.90,  0.15,  0.07 },
-            {  3.04,  0.02, -0.01 },
-            {  3.26, -0.82,  0.38 },
-            {  3.89,  0.45, -0.13 },
-            {  6.19,  0.05, -0.03 },
-            {  5.52, -0.56,  0.25 },
-            {  7.01, -0.29,  0.32 },
-            { -6.14,  3.08,  0.00 },
-            { -6.83,  2.82,  0.61 },
-            { -6.59,  3.64, -0.64 },
-            { -3.05,  3.09, -0.04 },
-            { -3.79,  2.50,  0.09 },
-            { -3.18,  3.80,  0.59 },
-            {  0.02,  3.14,  0.04 },
-            { -0.89,  3.04, -0.19 },
-            {  0.49,  2.57, -0.57 },
-            {  3.14,  3.15,  0.00 },
-            {  3.28,  2.28,  0.37 },
-            {  2.30,  3.08, -0.45 },
-            {  6.27,  3.08, -0.00 },
-            {  5.55,  2.54, -0.33 },
-            {  5.83,  3.87,  0.34 },
-            { -6.18,  6.15, -0.03 },
-            { -6.45,  6.21,  0.88 },
-            { -6.26,  7.05, -0.36 },
-            { -3.06,  6.19, -0.05 },
-            { -2.84,  6.64,  0.76 },
-            { -3.99,  5.96,  0.03 },
-            { -0.00,  6.20,  0.06 },
-            { -0.67,  5.99, -0.59 },
-            {  0.76,  6.46, -0.44 },
-            {  3.10,  6.26, -0.03 },
-            {  3.57,  6.09,  0.78 },
-            {  2.57,  5.47, -0.18 },
-            {  6.26,  6.18,  0.02 },
-            {  5.53,  5.64, -0.29 },
-            {  5.95,  7.08, -0.06 },
-            { -6.26, -6.21,  3.07 },
-            { -5.98, -6.38,  3.97 },
-            { -5.46, -5.94,  2.62 },
-            { -3.10, -6.24,  3.04 },
-            { -2.69, -6.51,  3.87 },
-            { -3.43, -5.35,  3.21 },
-            { -0.03, -6.16,  3.06 },
-            {  0.83, -6.00,  3.42 },
-            { -0.30, -6.99,  3.45 },
-            {  3.15, -6.25,  3.11 },
-            {  2.77, -5.60,  3.72 },
-            {  2.68, -6.10,  2.28 },
-            {  6.20, -6.21,  3.16 },
-            {  5.75, -6.73,  2.50 },
-            {  6.69, -5.56,  2.66 },
-            { -6.17, -3.10,  3.04 },
-            { -6.82, -2.44,  3.28 },
-            { -6.12, -3.69,  3.80 },
-            { -3.08, -3.04,  3.11 },
-            { -3.59, -3.56,  3.72 },
-            { -2.97, -3.61,  2.34 },
-            {  0.01, -3.04,  3.11 },
-            { -0.86, -3.41,  3.20 },
-            {  0.56, -3.78,  2.86 },
-            {  3.07, -3.07,  3.15 },
-            {  3.81, -3.68,  3.13 },
-            {  2.80, -2.98,  2.23 },
-            {  6.20, -3.04,  3.13 },
-            {  5.48, -3.64,  2.92 },
-            {  6.98, -3.49,  2.81 },
-            { -6.18, -0.05,  3.12 },
-            { -6.41,  0.66,  3.69 },
-            { -6.33,  0.28,  2.23 },
-            { -3.05,  0.03,  3.10 },
-            { -3.46, -0.42,  3.83 },
-            { -3.57, -0.19,  2.33 },
-            {  0.03, -0.02,  3.15 },
-            {  0.23, -0.08,  2.21 },
-            { -0.81,  0.41,  3.18 },
-            {  3.09,  0.00,  3.03 },
-            {  2.48, -0.29,  3.71 },
-            {  3.91,  0.16,  3.51 },
-            {  6.19, -0.06,  3.11 },
-            {  6.05,  0.47,  2.33 },
-            {  6.59,  0.52,  3.74 },
-            { -6.20,  3.05,  3.05 },
-            { -6.87,  3.73,  3.17 },
-            { -5.55,  3.24,  3.73 },
-            { -3.11,  3.06,  3.15 },
-            { -3.64,  3.74,  2.71 },
-            { -2.32,  3.00,  2.62 },
-            {  0.02,  3.05,  3.06 },
-            { -0.87,  3.14,  3.38 },
-            {  0.48,  3.82,  3.42 },
-            {  3.07,  3.10,  3.16 },
-            {  3.95,  3.44,  2.97 },
-            {  2.76,  2.73,  2.32 },
-            {  6.19,  3.07,  3.16 },
-            {  7.02,  3.30,  2.72 },
-            {  5.52,  3.27,  2.51 },
-            { -6.19,  6.24,  3.15 },
-            { -5.56,  5.88,  2.52 },
-            { -7.05,  5.96,  2.83 },
-            { -3.10,  6.14,  3.08 },
-            { -2.34,  6.69,  3.27 },
-            { -3.86,  6.69,  3.29 },
-            { -0.04,  6.24,  3.13 },
-            {  0.63,  6.54,  2.53 },
-            {  0.08,  5.29,  3.18 },
-            {  3.12,  6.24,  3.14 },
-            {  3.57,  5.82,  2.40 },
-            {  2.23,  5.90,  3.12 },
-            {  6.25,  6.19,  3.06 },
-            {  5.55,  5.59,  3.32 },
-            {  6.08,  6.99,  3.55 },
-            { -6.20, -6.16,  6.15 },
-            { -6.29, -5.99,  7.09 },
-            { -6.09, -7.11,  6.09 },
-            { -3.09, -6.19,  6.27 },
-            { -2.56, -5.90,  5.52 },
-            { -3.80, -6.69,  5.87 },
-            {  0.02, -6.25,  6.24 },
-            { -0.70, -5.70,  6.51 },
-            {  0.25, -5.93,  5.36 },
-            {  3.11, -6.18,  6.14 },
-            {  3.76, -6.54,  6.74 },
-            {  2.29, -6.20,  6.64 },
-            {  6.22, -6.17,  6.15 },
-            {  6.61, -6.98,  6.47 },
-            {  5.56, -5.94,  6.81 },
-            { -6.21, -3.10,  6.14 },
-            { -6.76, -2.66,  6.78 },
-            { -5.51, -3.50,  6.65 },
-            { -3.13, -3.05,  6.18 },
-            { -2.19, -3.14,  6.34 },
-            { -3.50, -3.89,  6.43 },
-            {  0.01, -3.06,  6.15 },
-            { -0.06, -2.81,  7.07 },
-            { -0.25, -3.98,  6.13 },
-            {  3.04, -3.09,  6.17 },
-            {  3.84, -3.51,  5.84 },
-            {  3.25, -2.85,  7.08 },
-            {  6.26, -3.13,  6.19 },
-            {  6.01, -2.20,  6.09 },
-            {  5.47, -3.55,  6.54 },
-            { -6.20,  0.01,  6.27 },
-            { -5.79, -0.70,  5.78 },
-            { -6.67,  0.51,  5.60 },
-            { -3.13,  0.01,  6.14 },
-            { -3.53, -0.35,  6.94 },
-            { -2.21,  0.17,  6.39 },
-            { -0.04, -0.04,  6.20 },
-            {  0.26,  0.47,  5.46 },
-            {  0.51,  0.22,  6.93 },
-            {  3.10, -0.05,  6.23 },
-            {  2.33,  0.44,  5.95 },
-            {  3.85,  0.45,  5.92 },
-            {  6.19, -0.01,  6.26 },
-            {  7.05,  0.16,  5.88 },
-            {  5.58,  0.02,  5.52 },
-            { -6.22,  3.04,  6.17 },
-            { -5.45,  3.57,  5.95 },
-            { -6.62,  3.50,  6.92 },
-            { -3.09,  3.16,  6.21 },
-            { -3.71,  2.75,  5.61 },
-            { -2.60,  2.43,  6.59 },
-            { -0.02,  3.10,  6.26 },
-            {  0.89,  3.27,  6.05 },
-            { -0.44,  2.94,  5.41 },
-            {  3.12,  3.04,  6.23 },
-            {  2.31,  3.53,  6.43 },
-            {  3.59,  3.60,  5.60 },
-            {  6.23,  3.05,  6.24 },
-            {  5.92,  3.91,  6.54 },
-            {  6.02,  3.03,  5.30 },
-            { -6.15,  6.21,  6.24 },
-            { -6.27,  6.46,  5.32 },
-            { -7.00,  5.85,  6.51 },
-            { -3.07,  6.15,  6.22 },
-            { -3.98,  6.27,  5.94 },
-            { -2.66,  7.01,  6.10 },
-            {  0.04,  6.20,  6.25 },
-            { -0.38,  5.50,  5.75 },
-            { -0.36,  7.00,  5.93 },
-            {  3.12,  6.15,  6.24 },
-            {  3.66,  6.88,  5.93 },
-            {  2.25,  6.33,  5.86 },
-            {  6.20,  6.27,  6.19 },
-            {  5.46,  5.65,  6.19 },
-            {  6.97,  5.73,  6.39 }
-        };
-        for(int atom = 0; atom < natoms; ++atom)
-            positions.push_back(Vec3(coords[atom][0], coords[atom][1], coords[atom][2])*OpenMM::NmPerAngstrom);
-    }else{
-        throw exception();
-    }
-
-    system.setDefaultPeriodicBoxVectors(Vec3(boxEdgeLength, 0, 0),
-                                        Vec3(0, boxEdgeLength, 0),
-                                        Vec3(0, 0, boxEdgeLength));
-
-    sig.clear();
-    eps.clear();
-    bonds.clear();
-    for(int atom = 0; atom < natoms; ++atom){
-        system.addParticle(masses[atom%RESSIZE]);
-        double sigma = 2.0*pow(2.0, -1.0/6.0)*halfrmins[atom%RESSIZE]*OpenMM::NmPerAngstrom;
-        double epsilon = fabs(epsilons[atom%RESSIZE])*OpenMM::KJPerKcal;
-        sig.push_back(0.5*sigma);
-        eps.push_back(2.0*sqrt(epsilon));
-        if(atom%RESSIZE == 0){
-            bonds.push_back(pair<int, int>(atom, atom+1));
-            bonds.push_back(pair<int, int>(atom, atom+2));
-        }
-        double charge = do_electrostatics ? charges[atom] : 0;
-        forceField->addParticle(charge, sigma, epsilon);
-    }
-}
-
-void print_forces(const vector<Vec3>& forces){
-    // Print the forces in AKMA units, to compare against other codes.
-    std::cout << "Forces:" << std::endl;
-    for(int n = 0; n < forces.size(); ++ n){
-        std::cout << setw(3)<< n+1
-                  << setw(16) << setprecision(10) << -forces[n][0]*OpenMM::NmPerAngstrom*OpenMM::KcalPerKJ
-                  << setw(16) << setprecision(10) << -forces[n][1]*OpenMM::NmPerAngstrom*OpenMM::KcalPerKJ
-                  << setw(16) << setprecision(10) << -forces[n][2]*OpenMM::NmPerAngstrom*OpenMM::KcalPerKJ << std::endl;
-    }
-}
-
-
-void test_water2_dpme_energies_forces_no_exclusions() {
-    const double cutoff = 7.0*OpenMM::NmPerAngstrom;
-    const double alpha = 4.0124063605;
-    const double dalpha = 4.0124063605;
-    const int grid = 32;
-    NonbondedForce* forceField = new NonbondedForce();
-
-    vector<Vec3> positions;
-    vector<double> epsvals;
-    vector<double> sigvals;
-    vector<pair<int, int> > bonds;
-    System system;
-
-    const int NATOMS = 6;
-    double boxEdgeLength = 25*OpenMM::NmPerAngstrom;
-
-    make_waterbox(NATOMS, boxEdgeLength, forceField,  positions, epsvals, sigvals, bonds, system, false);
-    forceField->setNonbondedMethod(OpenMM::NonbondedForce::LJPME);
-    forceField->setPMEParameters(alpha, grid, grid, grid);
-    forceField->setLJPMEParameters(dalpha, grid, grid, grid);
-    forceField->setCutoffDistance(cutoff);
-    forceField->setReactionFieldDielectric(1.0);
-    system.addForce(forceField);
-
-    VerletIntegrator integrator(0.01);
-    Context context(system, integrator, platform);
-    context.setPositions(positions);
-
-
-    /**
-     * Gromacs reference values from the following inputs
-     *
-------------------------------------------------
- water2.gro
-------------------------------------------------
-MD of 2 waters, t= 0.0
-    6
-    1SOL    OW1    1   0.200   0.200   0.200
-    1SOL    HW1    2   0.250   0.200   0.300
-    1SOL    HW2    3   0.150   0.200   0.300
-    2SOL    OW1    4   0.000   0.000   0.000
-    2SOL    HW1    5   0.050   0.000   0.100
-    2SOL    HW2    6  -0.050   0.000   0.100
-    2.50000   2.50000   2.50000
-------------------------------------------------
-
-------------------------------------------------
- water2.gro
-------------------------------------------------
-[ defaults ]
-; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
-  1             2               no              1.0     1.0
-
-[ atomtypes ]
-; full atom descriptions are available in ffoplsaa.atp
-; name  bond_type    mass    charge   ptype          sigma      epsilon
-OW   OW 8     15.99940    -0.834       A    3.15057e-01  6.36386e-01
-HW   HW 1      1.00800     0.417       A    4.00014e-02  1.92464e-01
-
-
-#ifdef NOEXCLUDE
-
-[ moleculetype ]
-; molname   nrexcl
-SOL     0
-
-#else
-
-[ moleculetype ]
-; molname   nrexcl
-SOL     3
-
-#endif
-
-#ifdef NOCHARGE
-[ atoms ]
-; id    at type res nr  residu name at name     cg nr   charge
-1      OW  1       SOL             OW1             1        0.000
-2      HW  1       SOL             HW1             1        0.000
-3      HW  1       SOL             HW2             1        0.000
-#else
-[ atoms ]
-; id    at type res nr  residu name at name     cg nr   charge
-1      OW  1       SOL             OW1             1       -0.834
-2      HW  1       SOL             HW1             1        0.417
-3      HW  1       SOL             HW2             1        0.417
-#endif
-
-[ settles ]
-; i j   funct   length
-1   1   0.09572 0.15139
-
-#ifndef NOEXCLUDE
-[ exclusions ]
-1   2   3
-2   1   3
-3   1   2
-#endif
-------------------------------------------------
-
-------------------------------------------------
- water2.top
-------------------------------------------------
-#include "tip3p.tpi"
-
-[ system ]
-Water dimer
-
-[ molecules ]
-SOL 2
-------------------------------------------------
-
-------------------------------------------------
- water2.mdp
-------------------------------------------------
-define = -DNOCHARGE -DNOEXCLUDE
-;define = -DNOCHARGE
-rvdw             = 0.7
-vdw-modifier     = Potential-Shift
-;vdw-modifier     = None
-rlist            = 0.9
-rcoulomb         = 0.7
-nstfout          = 1
-fourierspacing   = 0.08
-pme-order        = 5
-ewald-rtol       = 1.5E-2
-ewald-rtol-lj    = 1.5E-2
-vdwtype          = PME
-lj-pme-comb-rule = geometric
-continuation     = yes
-------------------------------------------------
-
-run commands:-
-    gmx grompp -c water2.gro -p water2.top -f water2.mdp -o run.tpr
-    gmx mdrun -s run.tpr -o traj.trr -pforce 1E-16
-    gmx dump -f traj.trr
-
-     */
-    double refenergy = 1.34804e+03;
-    vector<Vec3> refforces(6);
-    refforces[0] = Vec3( 3.15301e-01, -3.91114e-03, -6.68965e+04);
-    refforces[1] = Vec3( 1.67241e+04, -3.79846e-02,  3.34487e+04);
-    refforces[2] = Vec3(-1.67243e+04, -9.55931e-02,  3.34486e+04);
-    refforces[3] = Vec3(-1.71643e+00, -1.76696e+00, -6.68993e+04);
-    refforces[4] = Vec3( 1.67254e+04,  1.59080e+00,  3.34495e+04);
-    refforces[5] = Vec3(-1.67239e+04,  3.13897e-01,  3.34491e+04);
-
-    State state = context.getState(State::Forces | State::Energy);
-    double energy = state.getPotentialEnergy();
-//    std::cout << "Energy " << energy << std::endl;
-    const vector<Vec3>& forces = state.getForces();
-//    print_forces(forces);
-
-    ASSERT_EQUAL_TOL(refenergy, energy, 1E-4);
-    for(int n = 0; n < NATOMS; ++n)
-        ASSERT_EQUAL_VEC(refforces[n], forces[n], 1E-4);
-}
-
-
-void test_water2_dpme_energies_forces_with_exclusions() {
-    const double cutoff = 7.0*OpenMM::NmPerAngstrom;
-    const double alpha = 4.0124063605;
-    const double dalpha = 4.0124063605;
-    const int grid = 32;
-    NonbondedForce* forceField = new NonbondedForce();
-
-    vector<Vec3> positions;
-    vector<double> epsvals;
-    vector<double> sigvals;
-    vector<pair<int, int> > bonds;
-    System system;
-
-    const int NATOMS = 6;
-    double boxEdgeLength = 25*OpenMM::NmPerAngstrom;
-
-    make_waterbox(NATOMS, boxEdgeLength, forceField,  positions, epsvals, sigvals, bonds, system, false);
-
-    forceField->setNonbondedMethod(OpenMM::NonbondedForce::LJPME);
-    forceField->createExceptionsFromBonds(bonds, 1.0, 1.0);
-    forceField->setPMEParameters(alpha, grid, grid, grid);
-    forceField->setLJPMEParameters(dalpha, grid, grid, grid);
-    forceField->setCutoffDistance(cutoff);
-    forceField->setUseSwitchingFunction(false);
-    forceField->setUseDispersionCorrection(false);
-
-    forceField->setReactionFieldDielectric(1.0);
-    system.addForce(forceField);
-
-    VerletIntegrator integrator(0.01);
-    Context context(system, integrator, platform);
-    context.setPositions(positions);
-
-    double refenergy = -7.78060e-01;
-    vector<Vec3> refforces(6);
-
-    // Gromacs reference values.  See comments in test_water2_dpme_energies_forces_no_exclusions() for details.
-    refforces[0] = Vec3( 3.15301e-01, -3.91114e-03,  9.48373e-01);
-    refforces[1] = Vec3(-4.74745e-02, -3.79846e-02, -5.69616e-02);
-    refforces[2] = Vec3(-7.16866e-02, -9.55931e-02, -1.43348e-01);
-    refforces[3] = Vec3(-1.78136e+00, -1.76696e+00, -1.70022e+00);
-    refforces[4] = Vec3( 1.19309e+00,  1.59080e+00,  7.95443e-01);
-    refforces[5] = Vec3( 3.92366e-01,  3.13897e-01,  1.56962e-01);
-
-    State state = context.getState(State::Forces | State::Energy);
-    double energy = state.getPotentialEnergy();
-//    std::cout << "Energy " << energy << std::endl;
-    const vector<Vec3>& forces = state.getForces();
-//    print_forces(forces);
-
-
-    ASSERT_EQUAL_TOL(refenergy, energy, 1E-4);
-    for(int n = 0; n < NATOMS; ++n)
-        ASSERT_EQUAL_VEC(refforces[n], forces[n], 5E-4);
-}
-
-
-void test_water125_dpme_vs_long_cutoff_no_exclusions() {
-    const double cutoff = 8.5*OpenMM::NmPerAngstrom;
-    const double alpha = 0.45*OpenMM::AngstromsPerNm;
-    const double dalpha = 0.45*OpenMM::AngstromsPerNm;
-    const int grid = 32;
-    NonbondedForce* forceField = new NonbondedForce();
-
-    vector<Vec3> positions;
-    vector<double> epsvals;
-    vector<double> sigvals;
-    vector<pair<int, int> > bonds;
-    System system;
-
-    const int NATOMS = 375;
-    double boxEdgeLength = 17.01*OpenMM::NmPerAngstrom;
-
-    make_waterbox(NATOMS, boxEdgeLength, forceField,  positions, epsvals, sigvals, bonds, system, false);
-
-    forceField->setNonbondedMethod(OpenMM::NonbondedForce::LJPME);
-    forceField->setPMEParameters(alpha, grid, grid, grid);
-    forceField->setLJPMEParameters(dalpha, grid, grid, grid);
-    forceField->setCutoffDistance(cutoff);
-    forceField->setReactionFieldDielectric(1.0);
-    system.addForce(forceField);
-
-    VerletIntegrator integrator(0.01);
-    Context context(system, integrator, platform);
-    context.setPositions(positions);
-
-    State state = context.getState(State::Forces | State::Energy);
-    double energy = state.getPotentialEnergy();
-    /*
-     *  Gromacs reference values.  See comments in test_water2_dpme_energies_forces_no_exclusions() for details.
-     * Coordinates are from make_waterbox, and the .gro file looks like
-     *
-;define = -DNOCHARGE -DNOEXCLUDE
-define = -DNOCHARGE
-vdw-modifier     = Potential-Shift
-rvdw             = 1.15
-rlist            = 1.15
-rcoulomb         = 1.15
-fourierspacing   = 0.04
-nstfout          = 1
-pme-order        = 5
-ewald-rtol       = 1E-8
-ewald-rtol-lj    = 1E-8
-vdwtype          = PME
-lj-pme-comb-rule = geometric
-continuation     = yes
-
-     */
-    double gromacs_energy = 5.63157e+05;
-//    std::cout << "Energy " << energy << std::endl;
-    const vector<Vec3>& forces = state.getForces();
-//    print_forces(forces);
-    vector<Vec3> refforces(NATOMS);
-    refforces[  0] = Vec3(-1.12446e+05, -2.71952e+05, -1.91471e+05);
-    refforces[  1] = Vec3( 2.70479e+05,  6.61172e+04,  7.21277e+04);
-    refforces[  2] = Vec3(-1.58058e+05,  2.05779e+05,  1.19292e+05);
-    refforces[  3] = Vec3( 3.16438e+05, -1.27994e+05,  1.08811e+05);
-    refforces[  4] = Vec3(-1.36520e+05,  1.93405e+05,  1.36521e+05);
-    refforces[  5] = Vec3(-1.79957e+05, -6.54374e+04, -2.45393e+05);
-    refforces[  6] = Vec3( 1.30626e+05, -1.36728e+05,  2.93833e+05);
-    refforces[  7] = Vec3(-2.74564e+05,  8.02095e+04, -7.09524e+04);
-    refforces[  8] = Vec3( 1.43952e+05,  5.64523e+04, -2.22980e+05);
-    refforces[  9] = Vec3(-4.22852e+04,  2.14660e+04, -3.23254e+05);
-    refforces[ 10] = Vec3( 2.28059e+05, -4.44251e+04,  1.62898e+05);
-    refforces[ 11] = Vec3(-1.85776e+05,  2.29044e+04,  1.60326e+05);
-    refforces[ 12] = Vec3(-1.02075e+05,  3.27346e+05,  1.47993e+04);
-    refforces[ 13] = Vec3( 7.77195e+04, -1.44336e+05,  2.10959e+05);
-    refforces[ 14] = Vec3( 2.44142e+04, -1.83111e+05, -2.25837e+05);
-    refforces[ 15] = Vec3(-4.81843e+04, -2.72889e+05, -1.75068e+05);
-    refforces[ 16] = Vec3(-5.46666e+03,  2.18711e+04,  2.62456e+05);
-    refforces[ 17] = Vec3( 5.35890e+04,  2.51021e+05, -8.74311e+04);
-    refforces[ 18] = Vec3(-2.04721e+05,  1.58918e+05,  2.20448e+05);
-    refforces[ 19] = Vec3(-7.05932e+04, -1.64717e+05, -2.17659e+05);
-    refforces[ 20] = Vec3( 2.75339e+05,  5.73611e+03, -2.86718e+03);
-    refforces[ 21] = Vec3(-5.65480e+03, -2.81809e+05, -1.38358e+05);
-    refforces[ 22] = Vec3(-7.85576e+03,  2.25204e+05, -1.15217e+05);
-    refforces[ 23] = Vec3( 1.34846e+04,  5.66345e+04,  2.53512e+05);
-    refforces[ 24] = Vec3(-7.73167e+04, -4.43114e+04,  3.22354e+05);
-    refforces[ 25] = Vec3(-1.24372e+05,  1.61973e+05, -1.88001e+05);
-    refforces[ 26] = Vec3( 2.01689e+05, -1.17651e+05, -1.34455e+05);
-    refforces[ 27] = Vec3(-1.79300e+05, -1.97922e+05,  1.94294e+05);
-    refforces[ 28] = Vec3( 2.38945e+05, -4.88761e+04, -9.50350e+04);
-    refforces[ 29] = Vec3(-5.95911e+04,  2.46878e+05, -9.93159e+04);
-    refforces[ 30] = Vec3(-1.31892e+04, -2.15675e+05,  2.68757e+05);
-    refforces[ 31] = Vec3( 2.31232e+05,  1.08107e+05, -1.32129e+05);
-    refforces[ 32] = Vec3(-2.18089e+05,  1.07592e+05, -1.36669e+05);
-    refforces[ 33] = Vec3( 1.90632e+04, -3.62889e+05,  8.61608e+04);
-    refforces[ 34] = Vec3(-2.16178e+05,  1.59639e+05, -1.66288e+05);
-    refforces[ 35] = Vec3( 1.97134e+05,  2.03295e+05,  8.00862e+04);
-    refforces[ 36] = Vec3( 3.40071e+05, -6.87734e+04,  1.22584e+05);
-    refforces[ 37] = Vec3(-4.17943e+04,  8.35901e+03, -2.64693e+05);
-    refforces[ 38] = Vec3(-2.98358e+05,  6.03813e+04,  1.42075e+05);
-    refforces[ 39] = Vec3(-3.21693e+05,  4.40885e+04,  1.41154e+04);
-    refforces[ 40] = Vec3( 1.28817e+05,  2.02067e+04, -2.07114e+05);
-    refforces[ 41] = Vec3( 1.92948e+05, -6.43160e+04,  1.92949e+05);
-    refforces[ 42] = Vec3(-1.46001e+05,  3.20840e+05,  7.97305e+04);
-    refforces[ 43] = Vec3(-1.27705e+05, -2.55411e+05, -1.24428e+05);
-    refforces[ 44] = Vec3( 2.73737e+05, -6.54594e+04,  4.46323e+04);
-    refforces[ 45] = Vec3( 1.50212e+04,  2.43629e+05, -2.20084e+05);
-    refforces[ 46] = Vec3(-1.07432e+05,  8.26408e+03,  2.42419e+05);
-    refforces[ 47] = Vec3( 9.23685e+04, -2.51915e+05, -2.23909e+04);
-    refforces[ 48] = Vec3( 3.14327e+04, -2.94199e+05,  1.55484e+05);
-    refforces[ 49] = Vec3(-2.23665e+05,  1.52883e+05, -2.54793e+04);
-    refforces[ 50] = Vec3( 1.92227e+05,  1.41342e+05, -1.30033e+05);
-    refforces[ 51] = Vec3( 5.73542e+04, -2.08689e+05, -2.68170e+05);
-    refforces[ 52] = Vec3(-2.26784e+05,  1.85259e+05,  8.30474e+04);
-    refforces[ 53] = Vec3( 1.69446e+05,  2.34613e+04,  1.85087e+05);
-    refforces[ 54] = Vec3( 2.25490e+05, -1.91311e+05, -1.40103e+05);
-    refforces[ 55] = Vec3(-8.20802e+03,  6.83988e+04,  2.54448e+05);
-    refforces[ 56] = Vec3(-2.17315e+05,  1.22953e+05, -1.14373e+05);
-    refforces[ 57] = Vec3(-7.09510e+04,  2.05639e+05, -2.69540e+05);
-    refforces[ 58] = Vec3( 1.93603e+05,  3.38041e+04,  2.18192e+05);
-    refforces[ 59] = Vec3(-1.22579e+05, -2.39458e+05,  5.13124e+04);
-    refforces[ 60] = Vec3(-1.07250e+05,  3.35982e+05,  6.89600e+03);
-    refforces[ 61] = Vec3( 6.90672e-01, -1.44228e+05, -2.24659e+05);
-    refforces[ 62] = Vec3( 1.07222e+05, -1.91701e+05,  2.17695e+05);
-    refforces[ 63] = Vec3( 2.64846e+05,  1.94002e+05,  5.27262e+03);
-    refforces[ 64] = Vec3(-1.12985e+04, -2.71176e+05,  8.47510e+03);
-    refforces[ 65] = Vec3(-2.53630e+05,  7.71894e+04, -1.37831e+04);
-    refforces[ 66] = Vec3(-3.04551e+05,  4.21918e+04,  1.88948e+05);
-    refforces[ 67] = Vec3( 2.87326e+05,  1.22193e+05,  3.30263e+04);
-    refforces[ 68] = Vec3( 1.73006e+04, -1.64358e+05, -2.22023e+05);
-    refforces[ 69] = Vec3( 2.45556e+04,  2.20596e+05,  2.41913e+05);
-    refforces[ 70] = Vec3( 1.50804e+05, -2.17829e+05, -4.46813e+04);
-    refforces[ 71] = Vec3(-1.75369e+05, -2.74087e+03, -1.97287e+05);
-    refforces[ 72] = Vec3( 8.65787e+04, -2.77188e+04, -3.15014e+05);
-    refforces[ 73] = Vec3(-2.42127e+05,  6.26659e+04,  1.11092e+05);
-    refforces[ 74] = Vec3( 1.55591e+05, -3.48752e+04,  2.03883e+05);
-    refforces[ 75] = Vec3( 7.06224e+04,  2.96642e+05, -1.51267e+05);
-    refforces[ 76] = Vec3(-5.39280e+04, -2.57657e+05, -1.13850e+05);
-    refforces[ 77] = Vec3(-1.67422e+04, -3.90659e+04,  2.65102e+05);
-    refforces[ 78] = Vec3(-4.52577e+04, -3.21552e+05, -7.01075e+04);
-    refforces[ 79] = Vec3( 2.35182e+05,  1.45173e+05,  3.48411e+04);
-    refforces[ 80] = Vec3(-1.89931e+05,  1.76363e+05,  3.52722e+04);
-    refforces[ 81] = Vec3(-2.29347e+05,  1.06977e+05, -2.70703e+05);
-    refforces[ 82] = Vec3(-1.17929e+04, -2.56493e+05,  1.17929e+05);
-    refforces[ 83] = Vec3( 2.41161e+05,  1.49452e+05,  1.52847e+05);
-    refforces[ 84] = Vec3( 2.32633e+03,  3.03443e+05,  1.27473e+05);
-    refforces[ 85] = Vec3(-2.11215e+05, -1.63335e+05, -4.50583e+04);
-    refforces[ 86] = Vec3( 2.08887e+05, -1.40170e+05, -8.24544e+04);
-    refforces[ 87] = Vec3( 5.83172e+04,  2.82133e+04, -3.26001e+05);
-    refforces[ 88] = Vec3( 1.76890e+05, -4.71698e+04,  2.15220e+05);
-    refforces[ 89] = Vec3(-2.35168e+05,  1.89225e+04,  1.10825e+05);
-    refforces[ 90] = Vec3(-2.72444e+05, -1.46998e+05, -1.00637e+05);
-    refforces[ 91] = Vec3( 2.78426e+04,  2.39442e+05,  1.16936e+05);
-    refforces[ 92] = Vec3( 2.44570e+05, -9.23920e+04, -1.63042e+04);
-    refforces[ 93] = Vec3(-3.10100e+04,  2.93390e+05, -1.87196e+05);
-    refforces[ 94] = Vec3(-6.38830e+04, -2.90671e+05, -6.38833e+04);
-    refforces[ 95] = Vec3( 9.48775e+04, -2.79002e+03,  2.51149e+05);
-    refforces[ 96] = Vec3( 4.18753e+04, -1.23438e+05, -3.10188e+05);
-    refforces[ 97] = Vec3( 1.29157e+05,  2.04500e+05,  9.41770e+04);
-    refforces[ 98] = Vec3(-1.71038e+05, -8.10180e+04,  2.16049e+05);
-    refforces[ 99] = Vec3(-5.61490e+04,  2.26993e+04, -3.44090e+05);
-    refforces[100] = Vec3(-1.96230e+05, -2.88572e+04,  1.93344e+05);
-    refforces[101] = Vec3( 2.52383e+05,  6.15570e+03,  1.50813e+05);
-    refforces[102] = Vec3(-6.33126e+04,  3.55941e+05,  8.51300e+04);
-    refforces[103] = Vec3(-1.75405e+05, -1.81783e+05, -1.69027e+05);
-    refforces[104] = Vec3( 2.38759e+05, -1.74232e+05,  8.38891e+04);
-    refforces[105] = Vec3( 1.51480e+05, -4.90606e+04,  3.17956e+05);
-    refforces[106] = Vec3( 4.93229e+04, -1.61668e+05, -2.02771e+05);
-    refforces[107] = Vec3(-2.00831e+05,  2.10712e+05, -1.15233e+05);
-    refforces[108] = Vec3( 2.20204e+05, -2.33820e+05,  1.51386e+05);
-    refforces[109] = Vec3( 3.36498e+04,  2.79296e+05, -1.51424e+05);
-    refforces[110] = Vec3(-2.53896e+05, -4.54334e+04,  2.10242e-01);
-    refforces[111] = Vec3(-1.94665e+05,  2.05890e+05,  2.40510e+05);
-    refforces[112] = Vec3(-9.73233e+04, -1.29764e+05, -2.62775e+05);
-    refforces[113] = Vec3( 2.92049e+05, -7.61856e+04,  2.22208e+04);
-    refforces[114] = Vec3( 1.60261e+05, -3.43716e+05,  1.52456e+04);
-    refforces[115] = Vec3( 1.11151e+05,  3.08358e+05, -8.60527e+04);
-    refforces[116] = Vec3(-2.71443e+05,  3.54056e+04,  7.08103e+04);
-    refforces[117] = Vec3(-4.27331e+04, -3.19871e+05, -1.68415e+05);
-    refforces[118] = Vec3( 2.28408e+05,  2.15171e+05, -2.31720e+04);
-    refforces[119] = Vec3(-1.85616e+05,  1.04782e+05,  1.91602e+05);
-    refforces[120] = Vec3(-1.66057e+05, -9.58191e+04, -2.78448e+05);
-    refforces[121] = Vec3( 1.91258e+05,  2.19476e+05,  6.89779e+04);
-    refforces[122] = Vec3(-2.52379e+04, -1.23674e+05,  2.09489e+05);
-    refforces[123] = Vec3( 6.55172e+03, -9.23173e+04,  3.29554e+05);
-    refforces[124] = Vec3(-2.14685e+05,  1.13143e+05, -1.36353e+05);
-    refforces[125] = Vec3( 2.08124e+05, -2.08124e+04, -1.93257e+05);
-    refforces[126] = Vec3( 1.02696e+05, -3.39299e+05, -4.47114e+04);
-    refforces[127] = Vec3(-1.81786e+05,  1.75407e+05, -1.69029e+05);
-    refforces[128] = Vec3( 7.90534e+04,  1.63962e+05,  2.13738e+05);
-    refforces[129] = Vec3(-3.45588e+05,  9.36838e+04,  5.67939e+04);
-    refforces[130] = Vec3( 1.44967e+05, -2.60943e+05,  7.08730e+04);
-    refforces[131] = Vec3( 2.00649e+05,  1.67207e+05, -1.27685e+05);
-    refforces[132] = Vec3(-2.70181e+05,  2.05707e+05, -3.11852e+04);
-    refforces[133] = Vec3( 1.09331e+05, -1.83207e+05, -1.86162e+05);
-    refforces[134] = Vec3( 1.60883e+05, -2.25804e+04,  2.17340e+05);
-    refforces[135] = Vec3(-1.04496e+05, -2.97000e+05, -1.63733e+05);
-    refforces[136] = Vec3( 2.11303e+05,  1.73661e+04,  1.79462e+05);
-    refforces[137] = Vec3(-1.06849e+05,  2.79694e+05, -1.57134e+04);
-    refforces[138] = Vec3(-3.82271e+04,  2.11942e+05,  2.60118e+05);
-    refforces[139] = Vec3(-1.96098e+05, -1.23692e+05, -1.71962e+05);
-    refforces[140] = Vec3( 2.34334e+05, -8.82195e+04, -8.82188e+04);
-    refforces[141] = Vec3( 2.17634e+05,  2.72527e+05,  1.42967e+05);
-    refforces[142] = Vec3( 9.30924e+04, -1.86185e+05, -1.98197e+05);
-    refforces[143] = Vec3(-3.10770e+05, -8.63246e+04,  5.52480e+04);
-    refforces[144] = Vec3(-1.63880e+05, -2.98869e+05,  1.01299e+05);
-    refforces[145] = Vec3( 6.83424e+04,  2.39196e+05,  1.61538e+05);
-    refforces[146] = Vec3( 9.55787e+04,  5.97354e+04, -2.62846e+05);
-    refforces[147] = Vec3( 1.06430e+05, -2.97083e+05, -7.97261e+04);
-    refforces[148] = Vec3(-1.14920e+05,  6.41394e+04,  2.21823e+05);
-    refforces[149] = Vec3( 8.52531e+03,  2.33043e+05, -1.42101e+05);
-    refforces[150] = Vec3(-4.96390e+04, -4.12076e+04, -3.67834e+05);
-    refforces[151] = Vec3( 1.25352e+05, -1.99962e+05,  1.61166e+05);
-    refforces[152] = Vec3(-7.57841e+04,  2.41138e+05,  2.06689e+05);
-    refforces[153] = Vec3(-3.06397e+05, -5.15222e+04, -1.37080e+05);
-    refforces[154] = Vec3( 1.92949e+05, -1.92945e+05,  6.43163e+04);
-    refforces[155] = Vec3( 1.13490e+05,  2.44443e+05,  7.27504e+04);
-    refforces[156] = Vec3(-3.74471e+03,  3.83223e+05, -2.14754e+04);
-    refforces[157] = Vec3( 2.17881e+05, -1.85835e+05, -1.05735e+05);
-    refforces[158] = Vec3(-2.14191e+05, -1.97454e+05,  1.27176e+05);
-    refforces[159] = Vec3(-3.33357e+05, -1.38307e+04, -1.17320e+05);
-    refforces[160] = Vec3( 2.04157e+05, -9.93173e+04, -1.37945e+05);
-    refforces[161] = Vec3( 1.29262e+05,  1.13105e+05,  2.55294e+05);
-    refforces[162] = Vec3( 2.50185e+05,  2.64217e+05, -7.18219e+04);
-    refforces[163] = Vec3(-2.46668e+05,  1.94013e+04, -9.97745e+04);
-    refforces[164] = Vec3(-3.50266e+03, -2.83658e+05,  1.71598e+05);
-    refforces[165] = Vec3(-2.05824e+05,  2.71177e+05, -1.16442e+05);
-    refforces[166] = Vec3(-3.52163e+04, -1.88897e+05,  2.36923e+05);
-    refforces[167] = Vec3( 2.41012e+05, -8.22954e+04, -1.20505e+05);
-    refforces[168] = Vec3( 6.89329e+04,  2.09489e+05,  2.76583e+05);
-    refforces[169] = Vec3( 1.87999e+05, -1.61968e+05, -1.24368e+05);
-    refforces[170] = Vec3(-2.56932e+05, -4.75793e+04, -1.52255e+05);
-    refforces[171] = Vec3( 3.39431e+05, -5.75509e+04,  1.62795e+05);
-    refforces[172] = Vec3(-1.27767e+05,  2.66184e+05, -1.59709e+05);
-    refforces[173] = Vec3(-2.11726e+05, -2.08613e+05, -3.11349e+03);
-    refforces[174] = Vec3(-2.58601e+05,  4.61975e+04, -2.46016e+05);
-    refforces[175] = Vec3( 7.15587e+04, -2.52015e+05,  1.40007e+05);
-    refforces[176] = Vec3( 1.87109e+05,  2.05820e+05,  1.06028e+05);
-    refforces[177] = Vec3( 3.92176e+03,  2.94819e+05, -1.84065e+05);
-    refforces[178] = Vec3(-1.67230e+05, -8.36141e+04,  2.29167e+05);
-    refforces[179] = Vec3( 1.63334e+05, -2.11213e+05, -4.50586e+04);
-    refforces[180] = Vec3( 1.11151e+05,  2.40551e+05, -2.68210e+05);
-    refforces[181] = Vec3(-1.76495e+05, -2.47099e+05, -3.53002e+04);
-    refforces[182] = Vec3( 6.52872e+04,  6.52855e+03,  3.03586e+05);
-    refforces[183] = Vec3(-4.84049e+04, -2.73305e+05, -2.95224e+05);
-    refforces[184] = Vec3(-9.04206e+04, -1.44672e+04,  3.29135e+05);
-    refforces[185] = Vec3( 1.38830e+05,  2.87820e+05, -3.38610e+04);
-    refforces[186] = Vec3(-2.09072e+05, -1.53610e+05, -2.86654e+05);
-    refforces[187] = Vec3(-1.30322e+05,  9.09223e+04,  2.42461e+05);
-    refforces[188] = Vec3( 3.39380e+05,  6.27113e+04,  4.42669e+04);
-    refforces[189] = Vec3(-3.15692e+05,  1.48900e+05, -9.20420e+04);
-    refforces[190] = Vec3( 7.13687e+04, -2.72503e+05,  1.26518e+05);
-    refforces[191] = Vec3( 2.44351e+05,  1.23612e+05, -3.44964e+04);
-    refforces[192] = Vec3(-2.80727e+04,  3.15077e+05, -2.05427e+05);
-    refforces[193] = Vec3(-2.29004e+05, -2.08496e+05,  9.57028e+04);
-    refforces[194] = Vec3( 2.57096e+05, -1.06603e+05,  1.09738e+05);
-    refforces[195] = Vec3( 3.33213e+05, -7.80040e+04, -5.04749e+03);
-    refforces[196] = Vec3(-2.07680e+05, -7.82574e+04,  1.83605e+05);
-    refforces[197] = Vec3(-1.25574e+05,  1.56274e+05, -1.78599e+05);
-    refforces[198] = Vec3( 2.66795e+05, -2.82879e+04, -2.26620e+05);
-    refforces[199] = Vec3(-2.28291e+05, -1.82015e+05,  4.01046e+04);
-    refforces[200] = Vec3(-3.85027e+04,  2.10288e+05,  1.86592e+05);
-    refforces[201] = Vec3( 1.93076e+05,  2.05300e+05,  2.64591e+05);
-    refforces[202] = Vec3(-3.32259e+05, -3.65117e+04, -8.39768e+04);
-    refforces[203] = Vec3( 1.39204e+05, -1.68822e+05, -1.80669e+05);
-    refforces[204] = Vec3( 2.15422e+05,  2.88264e+05,  2.49786e+04);
-    refforces[205] = Vec3( 4.29235e+04, -2.66744e+05,  1.13441e+05);
-    refforces[206] = Vec3(-2.58339e+05, -2.15280e+04, -1.38395e+05);
-    refforces[207] = Vec3( 3.27580e+05, -4.93781e+04,  7.43625e+03);
-    refforces[208] = Vec3(-2.11622e+05, -1.58716e+05, -9.69925e+04);
-    refforces[209] = Vec3(-1.15917e+05,  2.08125e+05,  8.95715e+04);
-    refforces[210] = Vec3( 1.10967e+05, -2.71544e+05, -2.06279e+05);
-    refforces[211] = Vec3(-8.86155e+04,  1.96918e+04,  2.98676e+05);
-    refforces[212] = Vec3(-2.23913e+04,  2.51912e+05, -9.23690e+04);
-    refforces[213] = Vec3( 1.91611e+05, -7.99998e+04, -2.83463e+05);
-    refforces[214] = Vec3( 7.08725e+04,  1.44963e+05,  2.60941e+05);
-    refforces[215] = Vec3(-2.62503e+05, -6.49202e+04,  2.25808e+04);
-    refforces[216] = Vec3(-6.63168e+04, -2.84265e+04,  3.72700e+05);
-    refforces[217] = Vec3(-2.02131e+05, -6.33547e+04, -1.96097e+05);
-    refforces[218] = Vec3( 2.68466e+05,  9.18414e+04, -1.76621e+05);
-    refforces[219] = Vec3(-6.80052e+03,  2.72744e+05, -2.21849e+05);
-    refforces[220] = Vec3( 1.52709e+05, -5.52356e+04,  2.63181e+05);
-    refforces[221] = Vec3(-1.45890e+05, -2.17461e+05, -4.12893e+04);
-    refforces[222] = Vec3( 3.07522e+05, -1.21146e+05,  1.14593e+05);
-    refforces[223] = Vec3(-2.11787e+05, -1.56664e+05, -8.99363e+04);
-    refforces[224] = Vec3(-9.57076e+04,  2.77857e+05, -2.46987e+04);
-    refforces[225] = Vec3(-3.24882e+05, -3.09840e+04, -1.31947e+05);
-    refforces[226] = Vec3( 8.33124e+04, -5.05809e+04,  2.67793e+05);
-    refforces[227] = Vec3( 2.41536e+05,  8.15185e+04, -1.35863e+05);
-    refforces[228] = Vec3(-2.16561e+04, -1.67607e+05, -2.70252e+05);
-    refforces[229] = Vec3( 1.10825e+05, -7.29802e+04,  2.24354e+05);
-    refforces[230] = Vec3(-8.91995e+04,  2.40572e+05,  4.59512e+04);
-    refforces[231] = Vec3(-2.22246e+05,  1.96766e+05, -2.46638e+05);
-    refforces[232] = Vec3( 3.04748e+05,  5.66971e+04,  1.27568e+05);
-    refforces[233] = Vec3(-8.24632e+04, -2.53495e+05,  1.19113e+05);
-    refforces[234] = Vec3( 2.20621e+05, -2.03896e+05,  6.63128e+04);
-    refforces[235] = Vec3(-9.59082e+04,  1.64055e+05,  1.53960e+05);
-    refforces[236] = Vec3(-1.24758e+05,  3.98170e+04, -2.20320e+05);
-    refforces[237] = Vec3(-7.79583e+03, -3.49688e+04,  3.64330e+05);
-    refforces[238] = Vec3(-1.43293e+05, -1.65580e+05, -2.10163e+05);
-    refforces[239] = Vec3( 1.51159e+05,  2.00522e+05, -1.54247e+05);
-    refforces[240] = Vec3( 1.82058e+05, -3.72909e+04, -2.80373e+05);
-    refforces[241] = Vec3(-1.95792e+05,  1.98809e+05,  7.22936e+04);
-    refforces[242] = Vec3( 1.36910e+04, -1.61546e+05,  2.08094e+05);
-    refforces[243] = Vec3( 1.40288e+05,  3.27379e+05,  4.78569e+03);
-    refforces[244] = Vec3(-1.70016e+05, -1.73349e+05,  2.03353e+05);
-    refforces[245] = Vec3( 2.97331e+04, -1.54072e+05, -2.08134e+05);
-    refforces[246] = Vec3( 1.21932e+05,  3.52267e+05,  4.69291e+04);
-    refforces[247] = Vec3(-2.91621e+05, -1.24021e+05,  3.01674e+04);
-    refforces[248] = Vec3( 1.69673e+05, -2.28288e+05, -7.71238e+04);
-    refforces[249] = Vec3(-1.41105e+05,  1.52818e+05,  2.59197e+05);
-    refforces[250] = Vec3( 2.15511e+05, -1.77650e+05, -5.82442e+03);
-    refforces[251] = Vec3(-7.43779e+04,  2.47926e+04, -2.53438e+05);
-    refforces[252] = Vec3(-3.34214e+04,  3.09554e+05,  1.58194e+05);
-    refforces[253] = Vec3(-2.05877e+05, -1.71563e+05, -6.00467e+04);
-    refforces[254] = Vec3( 2.39329e+05, -1.38076e+05, -9.81875e+04);
-    refforces[255] = Vec3( 1.32793e+05, -3.72271e+05,  2.88491e+04);
-    refforces[256] = Vec3(-9.02633e+04,  2.78646e+05,  2.23700e+05);
-    refforces[257] = Vec3(-4.25632e+04,  9.36427e+04, -2.52553e+05);
-    refforces[258] = Vec3( 2.97252e+05,  2.17294e+05, -2.50945e+03);
-    refforces[259] = Vec3(-1.35724e+05, -1.48966e+05,  2.41658e+05);
-    refforces[260] = Vec3(-1.61536e+05, -6.83423e+04, -2.39199e+05);
-    refforces[261] = Vec3( 2.50550e+05, -1.39926e+05,  2.48375e+05);
-    refforces[262] = Vec3( 5.51783e+04, -1.65530e+04, -2.59341e+05);
-    refforces[263] = Vec3(-3.05734e+05,  1.56505e+05,  1.09189e+04);
-    refforces[264] = Vec3(-4.44616e+04,  4.17040e+04, -3.31507e+05);
-    refforces[265] = Vec3(-1.79702e+05, -8.54319e+04,  2.00323e+05);
-    refforces[266] = Vec3( 2.24182e+05,  4.37421e+04,  1.31227e+05);
-    refforces[267] = Vec3(-9.89367e+04, -3.76136e+05,  2.17157e+04);
-    refforces[268] = Vec3(-4.44424e+04,  1.68244e+05, -2.47606e+05);
-    refforces[269] = Vec3( 1.43421e+05,  2.07965e+05,  2.25893e+05);
-    refforces[270] = Vec3(-1.09030e+03, -2.44686e+05, -2.30145e+05);
-    refforces[271] = Vec3(-1.86963e+05,  1.89757e+05,  3.34863e+04);
-    refforces[272] = Vec3( 1.88003e+05,  5.49541e+04,  1.96680e+05);
-    refforces[273] = Vec3(-1.15452e+05, -1.55243e+05,  2.81412e+05);
-    refforces[274] = Vec3(-1.35893e+05,  1.74355e+05, -1.12817e+05);
-    refforces[275] = Vec3( 2.51361e+05, -1.90905e+04, -1.68633e+05);
-    refforces[276] = Vec3( 1.76202e+05, -2.31601e+05, -2.00886e+05);
-    refforces[277] = Vec3(-2.96697e+05,  3.00025e+04,  1.06676e+05);
-    refforces[278] = Vec3( 1.20457e+05,  2.01635e+05,  9.42700e+04);
-    refforces[279] = Vec3(-1.66320e+05, -9.08355e+02,  2.65478e+05);
-    refforces[280] = Vec3( 2.44823e+05,  9.45895e+04, -5.28589e+04);
-    refforces[281] = Vec3(-7.84765e+04, -9.36654e+04, -2.12648e+05);
-    refforces[282] = Vec3(-6.58002e+03, -1.21918e+05,  3.16459e+05);
-    refforces[283] = Vec3( 2.15225e+05,  5.96416e+04, -1.14097e+05);
-    refforces[284] = Vec3(-2.08615e+05,  6.22725e+04, -2.02387e+05);
-    refforces[285] = Vec3( 7.09186e+04,  1.83619e+05,  2.71862e+05);
-    refforces[286] = Vec3( 1.78911e+05, -1.02234e+05, -1.78910e+05);
-    refforces[287] = Vec3(-2.49882e+05, -8.13582e+04, -9.29803e+04);
-    refforces[288] = Vec3(-1.35529e+04, -3.20226e+05, -1.16293e+05);
-    refforces[289] = Vec3( 2.28053e+05,  1.65034e+05,  5.70125e+04);
-    refforces[290] = Vec3(-2.14515e+05,  1.55237e+05,  5.92732e+04);
-    refforces[291] = Vec3(-2.65006e+05,  1.75233e+05,  1.91268e+05);
-    refforces[292] = Vec3( 2.29905e+05,  1.02940e+05, -2.05886e+05);
-    refforces[293] = Vec3( 3.51345e+04, -2.78156e+05,  1.46393e+04);
-    refforces[294] = Vec3( 1.59458e+05,  2.25128e+05,  2.11611e+05);
-    refforces[295] = Vec3( 1.24813e+05, -1.16492e+05, -2.05249e+05);
-    refforces[296] = Vec3(-2.84281e+05, -1.08601e+05, -6.38824e+03);
-    refforces[297] = Vec3( 2.61767e+05, -7.55993e+04, -2.32576e+05);
-    refforces[298] = Vec3(-2.07803e+05, -1.78116e+05,  7.71831e+04);
-    refforces[299] = Vec3(-5.39242e+04,  2.53755e+05,  1.55427e+05);
-    refforces[300] = Vec3(-6.44117e+03,  2.31328e+05, -2.54922e+05);
-    refforces[301] = Vec3(-2.61097e+04,  4.93208e+04,  2.72709e+05);
-    refforces[302] = Vec3( 3.25053e+04, -2.80718e+05, -1.77304e+04);
-    refforces[303] = Vec3( 7.06280e+04,  7.26140e+04,  3.28446e+05);
-    refforces[304] = Vec3( 1.45890e+05,  7.98271e+04, -2.06449e+05);
-    refforces[305] = Vec3(-2.16516e+05, -1.52473e+05, -1.21980e+05);
-    refforces[306] = Vec3( 1.94879e+05, -2.83084e+05,  1.41824e+05);
-    refforces[307] = Vec3(-2.57149e+05,  1.96432e+05,  9.64284e+04);
-    refforces[308] = Vec3( 6.22633e+04,  8.66280e+04, -2.38228e+05);
-    refforces[309] = Vec3( 3.26456e+04,  1.17141e+05, -3.28374e+05);
-    refforces[310] = Vec3( 2.01298e+05, -1.11486e+05,  1.85810e+05);
-    refforces[311] = Vec3(-2.33937e+05, -5.70482e+03,  1.42641e+05);
-    refforces[312] = Vec3( 6.42964e+04,  1.88719e+05, -2.86572e+05);
-    refforces[313] = Vec3( 1.22184e+05, -2.53768e+05,  1.00254e+05);
-    refforces[314] = Vec3(-1.86433e+05,  6.49694e+04,  1.86429e+05);
-    refforces[315] = Vec3(-4.12317e+04, -1.73607e+04, -3.68620e+05);
-    refforces[316] = Vec3(-1.78981e+05,  1.43186e+05,  2.08269e+05);
-    refforces[317] = Vec3( 2.20155e+05, -1.25802e+05,  1.60395e+05);
-    refforces[318] = Vec3(-1.58613e+05,  3.01594e+05, -1.29345e+05);
-    refforces[319] = Vec3( 2.79698e+05, -2.67790e+04,  4.76063e+04);
-    refforces[320] = Vec3(-1.21063e+05, -2.74847e+05,  8.17976e+04);
-    refforces[321] = Vec3( 1.00451e+05,  2.03427e+05, -2.75044e+05);
-    refforces[322] = Vec3(-2.13956e+04,  7.64144e+04,  2.81199e+05);
-    refforces[323] = Vec3(-7.91039e+04, -2.79910e+05, -6.08590e+03);
-    refforces[324] = Vec3(-2.80664e+05,  5.26142e+04, -1.53706e+05);
-    refforces[325] = Vec3( 2.23923e+05, -1.17559e+05, -9.23685e+04);
-    refforces[326] = Vec3( 5.68062e+04,  6.49218e+04,  2.46158e+05);
-    refforces[327] = Vec3( 2.88911e+05, -1.17893e+05, -7.40828e+04);
-    refforces[328] = Vec3(-6.38612e+04,  2.37565e+05, -2.55450e+04);
-    refforces[329] = Vec3(-2.25033e+05, -1.19637e+05,  9.96952e+04);
-    refforces[330] = Vec3( 1.03510e+04,  7.35559e+04,  3.47129e+05);
-    refforces[331] = Vec3( 1.26778e+05, -2.19542e+05, -1.51515e+05);
-    refforces[332] = Vec3(-1.37190e+05,  1.45950e+05, -1.95573e+05);
-    refforces[333] = Vec3(-1.31816e+05,  5.57951e+04, -2.81935e+05);
-    refforces[334] = Vec3(-1.08367e+05, -9.75306e+04,  2.16731e+05);
-    refforces[335] = Vec3( 2.40191e+05,  4.17715e+04,  6.52667e+04);
-    refforces[336] = Vec3(-2.86680e+05, -2.63007e+05,  1.37964e+04);
-    refforces[337] = Vec3( 1.03914e+05,  1.76653e+05, -2.56324e+05);
-    refforces[338] = Vec3( 1.82779e+05,  8.64040e+04,  2.42594e+05);
-    refforces[339] = Vec3( 1.09800e+03, -3.11636e+05,  1.85833e+05);
-    refforces[340] = Vec3(-2.39750e+05,  1.52568e+05, -8.71817e+04);
-    refforces[341] = Vec3( 2.38635e+05,  1.59089e+05, -9.86351e+04);
-    refforces[342] = Vec3(-8.76706e+04, -6.10520e+04,  3.31680e+05);
-    refforces[343] = Vec3( 2.64692e+05,  5.23233e+04, -1.16959e+05);
-    refforces[344] = Vec3(-1.76973e+05,  8.70353e+03, -2.14689e+05);
-    refforces[345] = Vec3(-1.15976e+05, -2.72303e+05, -1.33307e+05);
-    refforces[346] = Vec3( 2.20682e+05,  1.51896e+05, -6.30516e+04);
-    refforces[347] = Vec3(-1.04744e+05,  1.20457e+05,  1.96395e+05);
-    refforces[348] = Vec3( 4.57064e+04,  3.43550e+05,  7.23245e+04);
-    refforces[349] = Vec3(-1.91414e+05, -1.26580e+05, -1.85241e+05);
-    refforces[350] = Vec3( 1.45684e+05, -2.17040e+05,  1.12977e+05);
-    refforces[351] = Vec3(-1.88623e+05, -1.22936e+04,  3.10223e+05);
-    refforces[352] = Vec3( 3.06931e+05,  5.73380e+04, -7.08300e+04);
-    refforces[353] = Vec3(-1.18279e+05, -4.50586e+04, -2.39376e+05);
-    refforces[354] = Vec3( 8.31603e+04, -2.75958e+05,  1.16877e+05);
-    refforces[355] = Vec3(-2.08784e+05,  1.26300e+05,  5.15494e+04);
-    refforces[356] = Vec3( 1.25608e+05,  1.49661e+05, -1.68370e+05);
-    refforces[357] = Vec3( 1.44100e+05, -2.34453e+05,  1.73916e+05);
-    refforces[358] = Vec3(-8.65072e+04,  2.39988e+05,  8.37141e+04);
-    refforces[359] = Vec3(-5.75432e+04, -5.48037e+03, -2.57574e+05);
-    refforces[360] = Vec3( 2.72504e+05,  2.99781e+04,  1.85778e+05);
-    refforces[361] = Vec3(-3.41026e+04,  7.10483e+04, -2.61466e+05);
-    refforces[362] = Vec3(-2.38458e+05, -1.00996e+05,  7.57451e+04);
-    refforces[363] = Vec3( 1.45826e+05, -2.81284e+05,  1.15502e+05);
-    refforces[364] = Vec3(-2.63407e+05,  3.47332e+04, -8.10477e+04);
-    refforces[365] = Vec3( 1.17594e+05,  2.46657e+05, -3.44181e+04);
-    refforces[366] = Vec3( 2.59381e+05, -5.73297e+04,  2.56690e+05);
-    refforces[367] = Vec3(-1.25834e+05, -2.09725e+05, -1.49803e+05);
-    refforces[368] = Vec3(-1.33554e+05,  2.67105e+05, -1.06844e+05);
-    refforces[369] = Vec3( 7.18473e+04, -2.59021e+05,  1.89797e+05);
-    refforces[370] = Vec3( 1.56664e+05,  2.11783e+05, -8.99369e+04);
-    refforces[371] = Vec3(-2.28511e+05,  4.72761e+04, -9.98083e+04);
-    refforces[372] = Vec3(-1.99083e+04,  3.17042e+05, -5.62792e+04);
-    refforces[373] = Vec3(-1.96875e+05, -1.64948e+05, -7.96758e-01);
-    refforces[374] = Vec3( 2.16840e+05, -1.52072e+05,  5.63215e+04);
-
-
-    const double longcutoff = 30.0*OpenMM::NmPerAngstrom;
-    const double longcutoff2 = longcutoff*longcutoff;
-    const double cutoff2 = cutoff*cutoff;
-    const double cutoff6inv = 1.0 / (cutoff2*cutoff2*cutoff2);
-    const int nboxes = ceil(longcutoff/boxEdgeLength);
-
-    double refenergy = 0.0;
-    // Loop over home box first...
-    for(int i = 0; i < NATOMS; ++ i){
-        for(int j = i+1; j < NATOMS; ++j){
-            Vec3 dR = positions[i] - positions[j];
-            double R2 = dR[0]*dR[0] + dR[1]*dR[1] + dR[2]*dR[2];
-            double sig2 = (sigvals[i] + sigvals[j])*(sigvals[i] + sigvals[j]) / R2;
-            double sig6 = sig2*sig2*sig2;
-            double eps = epsvals[i]*epsvals[j];
-            refenergy += 2.0*eps*(sig6-1.0)*sig6;
-            if(R2 < cutoff2){
-                // Add a shift term for direct space parts withing t
-                refenergy += 2.0*eps*(pow(sigvals[i]+sigvals[j], 6) - 64.0*pow(sigvals[i]*sigvals[j], 3))*cutoff6inv;
-            }
-        }
-    }
-
-    // ... and now add in the image box terms
-    for(int bx = -nboxes; bx <= nboxes; ++bx){
-        for(int by = -nboxes; by <= nboxes; ++by){
-            for(int bz = -nboxes; bz <= nboxes; ++bz){
-                if(bx==0 && by==0 && bz==0) continue;
-                Vec3 offset(bx*boxEdgeLength, by*boxEdgeLength, bz*boxEdgeLength);
-                for(int i = 0; i < NATOMS; ++ i){
-                    for(int j = 0; j < NATOMS; ++j){
-                        Vec3 dR = positions[i] - positions[j] + offset;
-                        double R2 = dR[0]*dR[0] + dR[1]*dR[1] + dR[2]*dR[2];
-                        if(R2 > longcutoff2) continue;
-                        double sig2 = (sigvals[i] + sigvals[j])*(sigvals[i] + sigvals[j]) / R2;
-                        double sig6 = sig2*sig2*sig2;
-                        double eps = epsvals[i]*epsvals[j];
-                        refenergy += eps*(sig6-1.0)*sig6;
-                        if(R2 < cutoff2){
-                            // Add a shift term for direct space parts withing teh
-                            refenergy += eps*(pow(sigvals[i]+sigvals[j], 6) - 64.0*pow(sigvals[i]*sigvals[j], 3))*cutoff6inv;
-                        }
-                    }
-                }
-            }
-        }
-    }
-    refenergy *= 0.5;
-
-    // For this test the reference energy is 545636 kJ/mol, while the difference between DPME and 30A cutoffs
-    // is just 1.062 kJ/mol.  The difference is due to the fact that arithmetic mean combination rules are used
-    // up to the cutoff, while the reciprocal space uses the geometric mean.  See DOI: 10.1021/acs.jctc.5b00726
-
-    ASSERT_EQUAL_TOL(refenergy, energy, 5E-6);
-    ASSERT_EQUAL_TOL(gromacs_energy, energy, 5E-4);
-
-    for(int n = 0; n < NATOMS; ++n)
-        ASSERT_EQUAL_VEC(refforces[n], forces[n], 5E-4);
-}
-
-
-
-void test_water125_dpme_vs_long_cutoff_with_exclusions() {
-    const double cutoff = 11.5*OpenMM::NmPerAngstrom;
-    const double alpha = 0.45*OpenMM::AngstromsPerNm;
-    const double dalpha = 4.2760443169;
-    const int grid = 60;
-    NonbondedForce* forceField = new NonbondedForce();
-
-    vector<Vec3> positions;
-    vector<double> epsvals;
-    vector<double> sigvals;
-    vector<pair<int, int> > bonds;
-    System system;
-
-    const int NATOMS = 375;
-    double boxEdgeLength = 23.01*OpenMM::NmPerAngstrom;
-
-    make_waterbox(NATOMS, boxEdgeLength, forceField,  positions, epsvals, sigvals, bonds, system, false);
-
-    forceField->setNonbondedMethod(OpenMM::NonbondedForce::LJPME);
-    forceField->createExceptionsFromBonds(bonds, 1.0, 1.0);
-    forceField->setPMEParameters(alpha, grid, grid, grid);
-    forceField->setLJPMEParameters(dalpha, grid, grid, grid);
-    forceField->setCutoffDistance(cutoff);
-    forceField->setReactionFieldDielectric(1.0);
-    system.addForce(forceField);
-
-    VerletIntegrator integrator(0.01);
-    Context context(system, integrator, platform);
-    context.setPositions(positions);
-
-    State state = context.getState(State::Forces | State::Energy);
-    double energy = state.getPotentialEnergy();
-//    std::cout << "Energy " << energy << std::endl;
-    const vector<Vec3>& forces = state.getForces();
-//    print_forces(forces);
-    // Gromacs reference values.  See comments in test_water2_dpme_energies_forces_no_exclusions() for details.
-    double gromacs_energy = -2.17481e+02;
-    vector<Vec3> refforces(NATOMS);
-    refforces[  0] = Vec3(-3.10130e+01, -6.01219e+01, -5.31302e+01);
-    refforces[  1] = Vec3( 4.21376e+00,  9.13811e-01,  9.61319e-01);
-    refforces[  2] = Vec3( 1.42363e+00,  2.97219e+00,  1.27185e+00);
-    refforces[  3] = Vec3(-3.70172e+01, -2.95228e+01, -6.27987e+01);
-    refforces[  4] = Vec3(-6.57386e-01,  2.39279e+00,  1.53243e+00);
-    refforces[  5] = Vec3(-1.61426e+00,  7.11728e-01,  9.99752e-01);
-    refforces[  6] = Vec3( 1.60813e+01, -6.77979e+01, -1.01302e+02);
-    refforces[  7] = Vec3(-4.23118e+00,  9.25550e-01,  1.28003e+00);
-    refforces[  8] = Vec3( 1.81589e+00,  1.10083e+00,  1.28615e+00);
-    refforces[  9] = Vec3(-2.24427e+00, -5.71071e+01, -3.18064e+01);
-    refforces[ 10] = Vec3( 3.41715e+00,  1.32918e+00,  1.27908e+00);
-    refforces[ 11] = Vec3(-3.09930e+00,  9.91868e-01,  1.69245e+00);
-    refforces[ 12] = Vec3( 6.08033e+01, -1.02990e+02, -8.26733e+01);
-    refforces[ 13] = Vec3(-1.01346e+00,  1.35475e+00,  4.13725e+00);
-    refforces[ 14] = Vec3(-6.71645e-01,  4.69536e-01,  4.64774e-01);
-    refforces[ 15] = Vec3(-6.39056e+01, -9.85756e-01, -4.97837e+01);
-    refforces[ 16] = Vec3( 9.97113e-01,  1.13195e-01,  4.49387e+00);
-    refforces[ 17] = Vec3( 9.48864e-01,  4.22217e+00,  1.56364e+00);
-    refforces[ 18] = Vec3( 2.18124e+01, -6.14078e+01, -8.10129e+01);
-    refforces[ 19] = Vec3(-5.59701e-01, -1.47528e+00,  1.40157e+00);
-    refforces[ 20] = Vec3( 4.29539e+00,  4.27738e-02,  1.03009e+00);
-    refforces[ 21] = Vec3(-2.62377e+01,  2.55532e+01, -6.89576e+01);
-    refforces[ 22] = Vec3( 1.17222e-01,  3.81618e+00,  1.74502e+00);
-    refforces[ 23] = Vec3( 1.46448e-01,  1.06214e-01,  4.27946e+00);
-    refforces[ 24] = Vec3(-1.90819e+00,  9.86513e+00, -1.05364e+02);
-    refforces[ 25] = Vec3(-1.18537e+00,  1.75021e+00,  1.54448e+00);
-    refforces[ 26] = Vec3( 2.96182e+00, -9.48323e-01,  1.61299e+00);
-    refforces[ 27] = Vec3( 5.61257e+01,  7.69427e+01, -5.84534e+01);
-    refforces[ 28] = Vec3(-8.11821e-01, -7.15455e-01,  6.60049e-01);
-    refforces[ 29] = Vec3(-9.20901e-01,  3.76217e+00,  1.35685e+00);
-    refforces[ 30] = Vec3(-5.03791e+01,  2.24324e+01, -4.31606e+01);
-    refforces[ 31] = Vec3( 3.72530e+00,  3.48275e-01,  1.40060e+00);
-    refforces[ 32] = Vec3( 7.61265e-01,  8.56416e-01,  6.23359e-01);
-    refforces[ 33] = Vec3( 1.95683e+01,  4.16950e+01, -4.33266e+01);
-    refforces[ 34] = Vec3(-2.19677e+00,  1.06287e+00,  1.58260e+00);
-    refforces[ 35] = Vec3( 1.70923e+00,  2.09511e+00,  9.84492e-01);
-    refforces[ 36] = Vec3(-7.84273e+01, -3.34396e+01, -3.66486e+01);
-    refforces[ 37] = Vec3(-1.53567e-01,  1.74934e-01,  1.15266e+00);
-    refforces[ 38] = Vec3(-3.45548e+00,  1.68633e-01,  1.15923e+00);
-    refforces[ 39] = Vec3( 6.90725e+01, -2.09426e+01, -5.31532e+01);
-    refforces[ 40] = Vec3( 1.33274e+00,  3.09504e-01,  1.40065e+00);
-    refforces[ 41] = Vec3( 1.65828e+00, -1.96378e-01,  2.68350e+00);
-    refforces[ 42] = Vec3( 3.39954e+01, -2.65418e+01, -6.81184e+01);
-    refforces[ 43] = Vec3(-1.17491e+00, -2.88793e+00,  1.40541e+00);
-    refforces[ 44] = Vec3(-9.63132e-01, -3.78178e-01,  1.06756e+00);
-    refforces[ 45] = Vec3(-4.44610e+01, -1.81736e+01, -6.19409e+01);
-    refforces[ 46] = Vec3( 1.48529e+00,  8.95548e-06,  4.60732e+00);
-    refforces[ 47] = Vec3( 1.00451e+00, -3.88830e+00,  1.27726e+00);
-    refforces[ 48] = Vec3(-4.23022e+00,  2.52963e+01, -3.04756e+01);
-    refforces[ 49] = Vec3(-3.45228e+00,  7.61191e-01,  9.40731e-01);
-    refforces[ 50] = Vec3( 2.78622e+00,  8.32081e-01,  1.37588e+00);
-    refforces[ 51] = Vec3( 1.63360e+01,  3.02230e+01, -3.91446e+01);
-    refforces[ 52] = Vec3(-2.64761e+00,  1.18509e+00,  9.81393e-01);
-    refforces[ 53] = Vec3( 1.74883e+00, -1.51958e-01,  2.56504e+00);
-    refforces[ 54] = Vec3(-3.01111e+01,  4.08021e+01, -3.45058e+01);
-    refforces[ 55] = Vec3(-1.17129e-01,  9.74282e-02,  4.56402e+00);
-    refforces[ 56] = Vec3(-3.03696e+00,  4.44012e-01,  1.44269e+00);
-    refforces[ 57] = Vec3( 7.57956e+01, -1.12043e+01, -3.92084e+01);
-    refforces[ 58] = Vec3(-1.55617e+00,  2.26444e-01,  2.88234e+00);
-    refforces[ 59] = Vec3(-1.37287e+00, -3.79045e+00,  8.55326e-01);
-    refforces[ 60] = Vec3(-2.88831e+01,  5.72484e+01, -7.17030e+01);
-    refforces[ 61] = Vec3( 6.88458e-01, -2.06343e+00,  1.03552e+00);
-    refforces[ 62] = Vec3( 6.94194e-01, -1.83673e+00,  2.90979e+00);
-    refforces[ 63] = Vec3(-7.87489e+01,  2.17729e+01, -3.72371e+01);
-    refforces[ 64] = Vec3( 1.39267e-01, -4.55150e+00,  9.53893e-01);
-    refforces[ 65] = Vec3(-4.27928e+00, -1.14139e+00,  9.25405e-01);
-    refforces[ 66] = Vec3( 7.13947e+01,  3.02722e+01, -5.16551e+01);
-    refforces[ 67] = Vec3( 4.07933e+00, -1.28199e+00,  9.34707e-01);
-    refforces[ 68] = Vec3(-7.53668e-02, -2.14888e+00,  1.27713e+00);
-    refforces[ 69] = Vec3(-9.36762e+00,  2.98321e+01, -5.70761e+01);
-    refforces[ 70] = Vec3( 1.35989e+00, -3.09988e+00,  9.89945e-01);
-    refforces[ 71] = Vec3(-2.20314e+00, -7.81892e-01,  1.21781e+00);
-    refforces[ 72] = Vec3( 4.89194e+01,  7.40749e+01, -4.22017e+01);
-    refforces[ 73] = Vec3(-4.25696e+00, -1.18954e+00,  1.04383e+00);
-    refforces[ 74] = Vec3(-1.40510e+00, -9.87719e-01,  3.30559e+00);
-    refforces[ 75] = Vec3(-4.95912e+01, -8.24838e+01, -1.82510e+01);
-    refforces[ 76] = Vec3( 7.39933e-01,  9.29483e-01, -1.09883e+00);
-    refforces[ 77] = Vec3( 9.70359e-01,  1.03012e+00,  4.01759e+00);
-    refforces[ 78] = Vec3(-7.87499e+00, -1.87981e+01,  5.77017e+00);
-    refforces[ 79] = Vec3( 3.82847e+00,  1.01191e+00, -7.30715e-02);
-    refforces[ 80] = Vec3(-2.64831e+00,  1.68402e+00,  1.06474e-01);
-    refforces[ 81] = Vec3( 1.81335e+01, -6.72766e+01,  7.21397e+01);
-    refforces[ 82] = Vec3(-1.89188e-01,  1.27803e+00,  8.45103e-01);
-    refforces[ 83] = Vec3( 2.47584e+00,  1.54044e+00,  6.55440e-01);
-    refforces[ 84] = Vec3(-1.05824e+00, -6.50833e+01, -3.88413e+01);
-    refforces[ 85] = Vec3(-3.21762e+00,  1.58809e+00, -1.32879e-01);
-    refforces[ 86] = Vec3( 3.30285e+00,  1.52154e+00, -3.56471e-01);
-    refforces[ 87] = Vec3( 4.41972e+01, -3.58815e+01,  4.01227e+01);
-    refforces[ 88] = Vec3(-1.22923e+00,  8.91826e-01,  2.84505e+00);
-    refforces[ 89] = Vec3(-4.31649e+00,  1.00895e+00,  4.43366e-01);
-    refforces[ 90] = Vec3(-3.65176e+01,  4.86616e+01, -6.44163e+00);
-    refforces[ 91] = Vec3( 1.03989e+00,  3.91910e+00,  7.93373e-01);
-    refforces[ 92] = Vec3( 4.11089e+00, -4.60733e-01, -7.78461e-02);
-    refforces[ 93] = Vec3(-1.58332e+01, -6.82045e+01,  6.62934e+01);
-    refforces[ 94] = Vec3(-7.75645e-02, -3.66897e+00, -3.55044e-01);
-    refforces[ 95] = Vec3( 4.08399e-01,  2.84950e-01,  3.80641e+00);
-    refforces[ 96] = Vec3(-5.05975e+00,  4.20346e+01,  3.49518e+01);
-    refforces[ 97] = Vec3( 9.82673e-01,  2.64451e+00,  4.02971e-01);
-    refforces[ 98] = Vec3(-1.24992e+00, -3.34479e-01,  1.95012e+00);
-    refforces[ 99] = Vec3( 2.31815e+00, -2.30056e+00,  6.49918e+01);
-    refforces[100] = Vec3(-2.16986e+00, -1.05184e-01,  1.89395e+00);
-    refforces[101] = Vec3( 3.51302e+00,  1.83706e-01,  7.08534e-01);
-    refforces[102] = Vec3( 4.46929e+01, -7.11604e+01, -6.62201e+00);
-    refforces[103] = Vec3(-1.87730e+00, -1.10206e+00, -1.47504e+00);
-    refforces[104] = Vec3(-1.41639e+00, -1.36347e+00,  6.10811e-01);
-    refforces[105] = Vec3(-3.12256e+01, -1.73923e+01, -4.47428e+01);
-    refforces[106] = Vec3( 7.92004e-01, -1.39066e+00, -2.68309e+00);
-    refforces[107] = Vec3( 1.54891e+00,  2.50018e+00, -6.21395e-01);
-    refforces[108] = Vec3(-3.83585e+01,  3.90231e+01, -3.68876e+01);
-    refforces[109] = Vec3( 1.86446e-01,  3.49120e+00, -4.00101e-01);
-    refforces[110] = Vec3(-4.05592e+00, -2.21609e-01,  2.13257e-01);
-    refforces[111] = Vec3( 5.70991e+01, -5.95502e+01, -4.08285e+01);
-    refforces[112] = Vec3(-4.25039e-01, -3.96206e-01, -3.23550e+00);
-    refforces[113] = Vec3( 4.14811e+00, -1.47859e-01,  1.11071e-01);
-    refforces[114] = Vec3(-2.74979e+01,  4.40260e+01,  3.06766e+00);
-    refforces[115] = Vec3( 2.22804e-01,  3.48934e+00, -1.96162e-01);
-    refforces[116] = Vec3(-4.21754e+00,  1.18655e-01,  3.36974e-01);
-    refforces[117] = Vec3( 6.32783e+01,  8.02719e+01,  1.38219e+01);
-    refforces[118] = Vec3(-1.49058e+00,  1.87045e+00, -1.32164e-01);
-    refforces[119] = Vec3(-2.30776e+00,  1.73240e-01,  1.92445e+00);
-    refforces[120] = Vec3(-4.04272e+01, -1.77143e+01,  1.57461e+01);
-    refforces[121] = Vec3( 1.82378e+00,  2.60376e+00,  2.33866e-01);
-    refforces[122] = Vec3( 1.33379e+00, -1.16675e+00,  3.34265e+00);
-    refforces[123] = Vec3(-8.78170e+00,  1.35313e+01, -5.40700e+01);
-    refforces[124] = Vec3(-2.72817e+00,  3.27413e-01, -6.95012e-01);
-    refforces[125] = Vec3( 2.41935e+00, -2.66049e-01, -1.62940e+00);
-    refforces[126] = Vec3(-3.63273e+01,  6.91610e+01, -3.72663e+00);
-    refforces[127] = Vec3(-1.40546e+00,  7.48503e-01, -9.73230e-01);
-    refforces[128] = Vec3( 5.39530e-01,  8.56302e-01,  2.57617e+00);
-    refforces[129] = Vec3( 2.54426e+01, -5.04169e+01, -1.76229e+01);
-    refforces[130] = Vec3( 6.82595e-01, -3.17362e+00,  2.20610e-01);
-    refforces[131] = Vec3( 1.94539e+00,  1.53325e+00, -4.24624e-01);
-    refforces[132] = Vec3( 3.61678e+01, -7.89030e+01, -9.61808e+00);
-    refforces[133] = Vec3(-1.29295e+00, -1.55223e+00, -1.88241e+00);
-    refforces[134] = Vec3(-1.51438e+00,  3.14594e-02,  3.60043e+00);
-    refforces[135] = Vec3(-4.60860e+01,  6.18971e+01,  1.38188e+01);
-    refforces[136] = Vec3( 3.07506e+00, -9.60659e-01,  1.81485e+00);
-    refforces[137] = Vec3( 6.46515e-01, -7.30203e-01, -2.17471e-01);
-    refforces[138] = Vec3( 7.76539e+00,  3.18751e+01, -6.18739e+01);
-    refforces[139] = Vec3(-2.12790e+00, -1.06577e+00, -1.03098e+00);
-    refforces[140] = Vec3( 3.70469e+00, -7.70591e-01, -2.00619e-01);
-    refforces[141] = Vec3(-4.02748e+01,  1.93804e+01,  1.92911e+01);
-    refforces[142] = Vec3( 5.74429e-01, -1.58840e+00, -2.07083e+00);
-    refforces[143] = Vec3(-4.21689e+00, -7.25280e-01,  1.61839e-01);
-    refforces[144] = Vec3( 4.10890e+01,  6.50772e+01, -6.64661e+00);
-    refforces[145] = Vec3( 2.02022e-01, -1.46124e+00,  1.99785e+00);
-    refforces[146] = Vec3( 1.91733e-01, -1.14906e+00, -4.18291e+00);
-    refforces[147] = Vec3( 3.77785e+01,  1.01191e+02, -6.07436e+00);
-    refforces[148] = Vec3(-1.03030e+00, -1.26650e+00,  3.57188e+00);
-    refforces[149] = Vec3(-7.42295e-01, -1.05437e+00, -1.47491e+00);
-    refforces[150] = Vec3(-7.40150e+01, -3.53477e+01,  1.81946e+01);
-    refforces[151] = Vec3( 1.65267e+00,  1.39873e+00,  1.46780e+00);
-    refforces[152] = Vec3( 1.18632e+00,  2.66165e+00,  1.61236e+00);
-    refforces[153] = Vec3( 3.89536e+01, -2.96071e+01, -1.38610e+01);
-    refforces[154] = Vec3( 2.05061e+00,  1.54050e+00,  5.35476e-01);
-    refforces[155] = Vec3( 4.98336e-01,  3.85036e+00,  4.21301e-01);
-    refforces[156] = Vec3(-5.53456e+01, -6.87157e+01, -3.55693e+01);
-    refforces[157] = Vec3( 3.17070e+00,  1.65800e+00, -3.99476e-01);
-    refforces[158] = Vec3(-2.02395e+00,  1.48143e+00,  1.00879e+00);
-    refforces[159] = Vec3( 5.94087e+01, -4.57694e+01,  2.77585e+01);
-    refforces[160] = Vec3( 2.05570e+00,  1.46552e+00, -1.41447e+00);
-    refforces[161] = Vec3( 2.97741e-01,  9.12876e-01,  3.27902e+00);
-    refforces[162] = Vec3( 1.93442e+01, -4.16391e+01,  3.34803e-01);
-    refforces[163] = Vec3(-3.93865e+00,  8.46322e-01, -4.47285e-01);
-    refforces[164] = Vec3(-7.34898e-01,  9.82679e-01,  1.44017e+00);
-    refforces[165] = Vec3(-3.34824e+01, -1.38834e+01, -2.62952e+01);
-    refforces[166] = Vec3( 1.26693e+00, -1.10171e+00,  3.22806e+00);
-    refforces[167] = Vec3( 3.49758e+00, -1.84534e-01, -4.62412e-01);
-    refforces[168] = Vec3( 1.82122e+00, -5.69776e+01, -3.86605e+01);
-    refforces[169] = Vec3( 1.92404e+00, -1.05303e+00, -4.52186e-01);
-    refforces[170] = Vec3(-3.25172e+00,  9.44050e-02, -5.87318e-01);
-    refforces[171] = Vec3(-5.98475e+01,  1.85708e+01, -2.67060e+01);
-    refforces[172] = Vec3(-3.77230e-01,  3.29609e+00, -4.81590e-01);
-    refforces[173] = Vec3(-1.66002e+00, -2.41098e+00,  1.03425e-01);
-    refforces[174] = Vec3( 6.47852e+01,  4.38661e+00,  1.77482e+01);
-    refforces[175] = Vec3(-1.29162e-01, -3.40031e+00,  3.91521e-01);
-    refforces[176] = Vec3( 1.28253e+00,  2.14633e+00,  2.64792e-01);
-    refforces[177] = Vec3( 2.95010e+01, -5.60376e+00,  4.16733e+01);
-    refforces[178] = Vec3(-1.43378e+00, -1.96178e-01,  2.40594e+00);
-    refforces[179] = Vec3(-1.56255e+00, -3.08967e+00, -3.14335e-01);
-    refforces[180] = Vec3(-6.02384e+01, -1.59585e+01,  7.16798e+01);
-    refforces[181] = Vec3( 1.68669e+00, -3.54228e+00, -2.98122e-01);
-    refforces[182] = Vec3( 9.46412e-01, -8.03852e-02,  4.17305e+00);
-    refforces[183] = Vec3( 3.87313e+00,  4.46632e+01,  4.59130e+01);
-    refforces[184] = Vec3(-9.83299e-02, -8.11669e-02,  4.04515e+00);
-    refforces[185] = Vec3( 6.64172e-01,  3.66814e+00, -2.58278e-01);
-    refforces[186] = Vec3(-1.65450e+01,  2.35336e+01,  7.09485e+01);
-    refforces[187] = Vec3(-6.98821e-01,  2.45516e-01,  2.69560e+00);
-    refforces[188] = Vec3( 4.04145e+00,  5.46258e-02, -1.82210e-02);
-    refforces[189] = Vec3( 2.39311e+01,  1.22785e+01, -2.13260e+01);
-    refforces[190] = Vec3(-3.33802e-02, -3.60205e+00,  6.21747e-01);
-    refforces[191] = Vec3( 3.50616e+00,  6.09623e-01, -7.40284e-02);
-    refforces[192] = Vec3( 2.32693e+01, -1.88836e+01,  1.24193e+01);
-    refforces[193] = Vec3(-2.25948e+00, -1.74718e+00,  2.65705e-01);
-    refforces[194] = Vec3(-1.32949e+00, -8.67446e-01,  7.75386e-01);
-    refforces[195] = Vec3(-4.40626e+01,  1.17032e+01, -4.18284e+01);
-    refforces[196] = Vec3( 1.61723e+00, -6.18309e-01,  2.45312e+00);
-    refforces[197] = Vec3( 1.42535e+00,  1.71801e+00, -1.86477e+00);
-    refforces[198] = Vec3( 3.37157e+00, -1.65818e+01,  7.59773e+01);
-    refforces[199] = Vec3(-2.81989e+00, -1.37877e+00, -9.60015e-02);
-    refforces[200] = Vec3(-8.77602e-02,  2.44299e+00,  1.17362e+00);
-    refforces[201] = Vec3( 2.44644e+01, -3.26831e+01, -5.34506e+01);
-    refforces[202] = Vec3(-4.14853e+00, -9.48873e-02, -2.18126e-01);
-    refforces[203] = Vec3( 8.13810e-01, -1.30429e+00, -1.62539e+00);
-    refforces[204] = Vec3( 1.03414e+01, -4.14898e+00,  2.56841e+01);
-    refforces[205] = Vec3(-1.50990e-01, -3.90514e+00,  2.98649e-01);
-    refforces[206] = Vec3(-3.60035e+00, -7.83546e-02, -7.95300e-01);
-    refforces[207] = Vec3( 4.49059e+01,  2.88959e+01,  1.53626e+01);
-    refforces[208] = Vec3(-2.55855e+00, -1.55231e+00, -3.96734e-01);
-    refforces[209] = Vec3(-1.14072e+00,  3.25255e+00,  2.39860e-01);
-    refforces[210] = Vec3(-4.12850e+01,  6.11654e+01,  2.43179e+01);
-    refforces[211] = Vec3( 1.20781e+00, -8.28196e-01,  4.26212e+00);
-    refforces[212] = Vec3( 7.30387e-01, -9.26965e-01, -9.32605e-01);
-    refforces[213] = Vec3(-1.64972e+01,  4.58667e+01,  5.58154e+01);
-    refforces[214] = Vec3( 4.43484e-01, -1.67205e+00,  3.42164e+00);
-    refforces[215] = Vec3(-4.14170e+00, -9.68400e-01, -9.21231e-02);
-    refforces[216] = Vec3( 1.79252e+01,  6.24205e+01, -1.47777e+01);
-    refforces[217] = Vec3(-2.53437e+00, -9.16373e-01, -1.63636e+00);
-    refforces[218] = Vec3( 3.32247e+00, -1.29776e+00, -1.07427e+00);
-    refforces[219] = Vec3( 1.79752e+01,  5.16430e+01,  3.90249e+01);
-    refforces[220] = Vec3( 7.14180e-01, -1.00018e+00,  3.14118e+00);
-    refforces[221] = Vec3(-1.07965e+00, -3.41285e+00, -2.24161e-01);
-    refforces[222] = Vec3( 3.16317e+01,  4.98025e+01, -4.17518e+01);
-    refforces[223] = Vec3(-2.55899e+00, -1.63433e+00, -3.76149e-01);
-    refforces[224] = Vec3(-1.19443e+00, -1.02227e+00, -5.37349e-02);
-    refforces[225] = Vec3(-3.77249e+01, -4.85845e+01, -2.11073e+01);
-    refforces[226] = Vec3( 7.98995e-01,  1.26056e+00,  4.28442e+00);
-    refforces[227] = Vec3( 3.35648e+00,  9.11653e-01, -8.07025e-01);
-    refforces[228] = Vec3(-3.16004e+01, -2.04156e+01,  5.01542e+01);
-    refforces[229] = Vec3( 8.97766e-01,  1.34997e+00,  2.86585e+00);
-    refforces[230] = Vec3(-2.09603e-01,  4.05651e+00, -1.91043e-02);
-    refforces[231] = Vec3( 3.61582e+01, -3.35593e+01,  4.19940e+01);
-    refforces[232] = Vec3( 3.66075e+00,  7.71021e-01,  3.45598e-01);
-    refforces[233] = Vec3(-7.50361e-01,  1.22549e+00,  7.98677e-01);
-    refforces[234] = Vec3(-4.42511e+01, -2.71049e+01, -4.56576e+01);
-    refforces[235] = Vec3(-3.30181e-01,  1.98824e+00,  2.21005e+00);
-    refforces[236] = Vec3(-7.68719e-01,  9.09063e-01, -3.77633e+00);
-    refforces[237] = Vec3( 7.33817e+01, -3.03886e+01, -7.59758e+01);
-    refforces[238] = Vec3(-1.71701e+00,  1.41048e+00, -1.84182e+00);
-    refforces[239] = Vec3(-1.35491e+00,  2.27816e+00, -1.08799e+00);
-    refforces[240] = Vec3(-4.58218e+01, -2.91931e+01,  1.12369e+01);
-    refforces[241] = Vec3( 1.68546e+00,  2.81258e+00,  3.93755e-01);
-    refforces[242] = Vec3( 8.28917e-01, -1.32237e+00,  2.91182e+00);
-    refforces[243] = Vec3( 5.44446e+00, -3.99466e+01,  6.77923e+00);
-    refforces[244] = Vec3(-1.09231e+00, -5.21927e-01,  1.83581e+00);
-    refforces[245] = Vec3( 5.24849e-02, -7.21832e-01, -3.35896e+00);
-    refforces[246] = Vec3(-1.29863e+01, -3.94851e+01, -2.70347e+01);
-    refforces[247] = Vec3(-3.94725e+00, -2.01688e-01,  5.93703e-02);
-    refforces[248] = Vec3( 1.30447e+00, -2.35362e+00, -2.94278e-01);
-    refforces[249] = Vec3( 2.56542e+01, -3.79434e+01, -6.19782e+01);
-    refforces[250] = Vec3( 2.48025e+00, -1.19994e+00,  1.41640e-01);
-    refforces[251] = Vec3(-5.01379e-01,  2.48607e-01, -3.97535e+00);
-    refforces[252] = Vec3( 3.46461e+01, -8.37214e+01, -3.92225e+01);
-    refforces[253] = Vec3(-2.80216e+00, -1.13524e+00,  1.85582e-02);
-    refforces[254] = Vec3(-1.28421e+00, -9.50174e-01, -5.56828e-01);
-    refforces[255] = Vec3(-3.56755e+01,  1.44176e+01, -7.63345e-01);
-    refforces[256] = Vec3( 1.25547e+00,  2.64846e+00,  1.29226e+00);
-    refforces[257] = Vec3( 1.22691e+00,  5.72301e-01, -4.09853e+00);
-    refforces[258] = Vec3(-6.02997e+00, -1.31579e+01, -5.03304e+01);
-    refforces[259] = Vec3(-4.59947e-01, -8.09031e-01,  3.08383e+00);
-    refforces[260] = Vec3(-8.93872e-01, -3.70797e-01, -2.89773e+00);
-    refforces[261] = Vec3(-1.92458e+00,  2.50787e+01, -4.30937e+01);
-    refforces[262] = Vec3( 2.10381e-01, -1.24125e-01, -4.01554e+00);
-    refforces[263] = Vec3(-3.85853e+00,  6.72311e-01,  4.06873e-02);
-    refforces[264] = Vec3( 1.66952e+01,  1.47383e+01,  4.13458e+01);
-    refforces[265] = Vec3(-2.12261e+00, -5.40916e-01,  1.54532e+00);
-    refforces[266] = Vec3( 3.54041e+00, -4.73740e-02,  4.73489e-01);
-    refforces[267] = Vec3( 4.43789e+01,  6.96271e+01,  4.06335e+00);
-    refforces[268] = Vec3(-9.59563e-01,  9.94730e-01, -3.12464e+00);
-    refforces[269] = Vec3(-1.42363e+00,  1.62107e+00,  1.87365e+00);
-    refforces[270] = Vec3(-5.38380e+01,  2.27365e+01,  1.95666e+01);
-    refforces[271] = Vec3( 1.47405e+00,  2.11523e+00,  1.97777e-01);
-    refforces[272] = Vec3( 2.54377e+00, -4.04439e-02,  1.97927e+00);
-    refforces[273] = Vec3( 1.33568e+01,  1.88750e+01, -3.59233e+01);
-    refforces[274] = Vec3(-1.09293e+00,  2.25841e+00, -4.61222e-01);
-    refforces[275] = Vec3( 3.27762e+00, -1.71208e-01, -8.67824e-01);
-    refforces[276] = Vec3(-3.41864e+01,  3.38579e+01,  5.95219e+01);
-    refforces[277] = Vec3(-4.18249e+00, -8.81154e-02,  1.97498e-01);
-    refforces[278] = Vec3( 1.06177e+00,  2.34089e+00,  2.42037e-01);
-    refforces[279] = Vec3( 2.35073e+01,  1.61625e+01, -2.58858e+01);
-    refforces[280] = Vec3( 3.87903e+00,  1.67313e-01, -4.03473e-02);
-    refforces[281] = Vec3(-5.12341e-01, -5.62857e-01, -3.19237e+00);
-    refforces[282] = Vec3( 3.39154e+01, -5.09470e+00, -2.29146e+01);
-    refforces[283] = Vec3(-1.24992e+00,  5.24960e-01, -1.16022e+00);
-    refforces[284] = Vec3(-2.39537e+00,  1.82818e-01, -1.53421e+00);
-    refforces[285] = Vec3(-5.56195e+01,  2.89266e+01, -2.61137e+01);
-    refforces[286] = Vec3( 2.49596e+00, -1.03213e+00, -1.51258e+00);
-    refforces[287] = Vec3( 1.02471e+00, -1.05342e+00, -6.79809e-01);
-    refforces[288] = Vec3(-1.52895e+01,  4.83734e+01, -7.33233e+00);
-    refforces[289] = Vec3( 3.77837e+00, -1.59224e+00,  2.02675e-01);
-    refforces[290] = Vec3(-3.51854e+00, -1.54234e+00,  2.89348e-01);
-    refforces[291] = Vec3( 3.22262e+01,  2.32542e+01,  2.37269e+01);
-    refforces[292] = Vec3( 1.88520e+00, -1.38396e+00, -1.98750e+00);
-    refforces[293] = Vec3( 4.99036e-02, -4.53129e+00, -2.13859e-01);
-    refforces[294] = Vec3(-7.48274e+00,  3.71318e+01, -2.35667e+01);
-    refforces[295] = Vec3( 7.81671e-01, -1.07371e+00, -2.37009e+00);
-    refforces[296] = Vec3(-3.98896e+00, -9.36146e-01,  2.62843e-02);
-    refforces[297] = Vec3( 4.33081e+01,  4.24806e+01,  3.21190e+01);
-    refforces[298] = Vec3(-2.52531e+00, -1.90285e+00,  2.72091e-01);
-    refforces[299] = Vec3(-9.70927e-01, -1.13219e+00,  1.47351e+00);
-    refforces[300] = Vec3(-4.71070e+01, -7.19253e+01,  5.82986e+01);
-    refforces[301] = Vec3( 6.82929e-01,  1.00679e+00, -6.94522e-01);
-    refforces[302] = Vec3( 8.31359e-01,  7.09804e-01, -8.12652e-01);
-    refforces[303] = Vec3( 3.62309e+00, -3.42931e+01,  1.98231e+01);
-    refforces[304] = Vec3( 1.23138e+00,  8.49297e-01, -2.39283e+00);
-    refforces[305] = Vec3(-2.87428e+00,  1.41050e+00, -1.03948e+00);
-    refforces[306] = Vec3(-3.83149e+00, -2.61638e+01,  2.98550e+01);
-    refforces[307] = Vec3(-2.90010e+00,  1.38356e+00, -1.22906e+00);
-    refforces[308] = Vec3( 2.62337e-01,  9.72795e-01, -4.02660e+00);
-    refforces[309] = Vec3( 7.84655e+00, -5.09531e+01,  7.96762e+01);
-    refforces[310] = Vec3( 2.46967e+00,  9.48901e-01, -1.17317e+00);
-    refforces[311] = Vec3(-4.11596e+00,  7.95633e-01, -1.39335e+00);
-    refforces[312] = Vec3( 5.03705e+01, -8.09993e+01,  1.13172e+02);
-    refforces[313] = Vec3(-4.48925e-01,  4.29309e-01, -5.50572e-01);
-    refforces[314] = Vec3(-2.66201e+00,  1.05662e+00, -1.44247e+00);
-    refforces[315] = Vec3(-6.12660e+01,  2.31218e+01,  4.57287e+01);
-    refforces[316] = Vec3( 8.42017e-01,  1.41042e+00, -7.78795e-01);
-    refforces[317] = Vec3( 3.11298e+00, -8.04273e-01, -1.51273e+00);
-    refforces[318] = Vec3( 1.83823e+01, -2.83588e+01,  6.13597e+01);
-    refforces[319] = Vec3( 4.37781e+00,  8.12091e-02, -1.33992e+00);
-    refforces[320] = Vec3(-5.81396e-01, -3.61991e+00, -1.22489e+00);
-    refforces[321] = Vec3(-4.85270e+01, -6.50314e+01,  7.16441e+01);
-    refforces[322] = Vec3(-3.40093e-02,  9.01514e-01, -1.26573e+00);
-    refforces[323] = Vec3(-8.67193e-02, -3.94828e+00, -8.97127e-01);
-    refforces[324] = Vec3( 6.25833e+01, -2.31957e+01,  8.66089e+01);
-    refforces[325] = Vec3( 2.49809e+00, -7.64197e-01, -1.27775e+00);
-    refforces[326] = Vec3( 2.11019e-01,  7.81764e-01, -1.23664e+00);
-    refforces[327] = Vec3( 2.04801e+01,  3.14002e+01,  6.96727e+01);
-    refforces[328] = Vec3(-7.99932e-01,  4.35586e+00, -8.77957e-01);
-    refforces[329] = Vec3(-2.87129e+00, -1.02321e+00, -1.42702e+00);
-    refforces[330] = Vec3(-6.36347e+01, -3.58428e+01,  4.52682e+01);
-    refforces[331] = Vec3( 1.32210e+00, -2.58500e+00, -1.43327e+00);
-    refforces[332] = Vec3( 1.52218e+00,  1.63513e+00, -2.29784e+00);
-    refforces[333] = Vec3( 4.36105e+00,  3.70776e+01,  6.46133e+01);
-    refforces[334] = Vec3(-1.19452e+00, -1.01252e+00, -1.43172e+00);
-    refforces[335] = Vec3( 4.34983e+00, -6.52342e-02, -1.33662e+00);
-    refforces[336] = Vec3( 1.13050e+01,  4.99187e+01,  7.16119e+01);
-    refforces[337] = Vec3( 2.29221e-01,  6.35843e-01, -3.59259e+00);
-    refforces[338] = Vec3( 1.65239e+00,  4.20995e-01, -1.42878e+00);
-    refforces[339] = Vec3(-1.79210e+01,  1.92423e+01,  1.77516e+01);
-    refforces[340] = Vec3(-2.91278e+00,  9.22214e-01, -8.41281e-01);
-    refforces[341] = Vec3( 2.96465e+00,  1.05962e+00, -6.85280e-01);
-    refforces[342] = Vec3( 5.19378e+01, -2.55266e+01,  3.63518e+01);
-    refforces[343] = Vec3(-1.13088e+00,  4.54600e-01, -1.30862e+00);
-    refforces[344] = Vec3(-2.05369e+00, -9.69459e-03, -2.69351e+00);
-    refforces[345] = Vec3(-4.22214e+01,  4.88282e+01,  3.76741e+01);
-    refforces[346] = Vec3( 3.46996e+00,  9.57280e-01, -7.60128e-01);
-    refforces[347] = Vec3( 7.62585e-01,  1.14548e+00, -8.06697e-01);
-    refforces[348] = Vec3(-2.35024e+01, -6.73585e+01,  6.49050e+01);
-    refforces[349] = Vec3(-1.58984e+00, -6.44551e-01, -2.37758e+00);
-    refforces[350] = Vec3( 1.32604e+00, -2.37576e+00, -1.54939e+00);
-    refforces[351] = Vec3( 2.57776e+01, -1.40743e+01,  2.06580e+01);
-    refforces[352] = Vec3( 4.24916e+00,  6.10665e-02, -8.19360e-01);
-    refforces[353] = Vec3(-7.94423e-01, -1.75764e-01, -3.51414e+00);
-    refforces[354] = Vec3(-1.24260e+01,  6.96334e-01,  5.97705e+01);
-    refforces[355] = Vec3(-3.38899e+00,  9.44919e-01, -1.25780e+00);
-    refforces[356] = Vec3( 6.81664e-01,  1.16965e+00, -2.34256e+00);
-    refforces[357] = Vec3( 5.15292e+01,  5.13086e+01,  6.16724e+01);
-    refforces[358] = Vec3(-1.04216e+00,  3.34685e+00, -1.33787e+00);
-    refforces[359] = Vec3(-9.57017e-01, -1.59610e-01, -4.31963e+00);
-    refforces[360] = Vec3(-5.88173e+01,  3.30524e+01,  6.23031e+01);
-    refforces[361] = Vec3( 9.35413e-01, -1.09829e+00, -4.65202e+00);
-    refforces[362] = Vec3( 6.95250e-01, -1.24375e+00, -6.50832e-01);
-    refforces[363] = Vec3( 1.64454e+01,  1.08071e+02,  3.76968e+01);
-    refforces[364] = Vec3(-4.25249e+00, -1.10053e+00, -7.56357e-01);
-    refforces[365] = Vec3( 1.17898e+00, -1.10215e+00, -8.21366e-01);
-    refforces[366] = Vec3(-5.20527e+00,  5.49522e+01,  4.55303e+01);
-    refforces[367] = Vec3(-7.98794e-01, -2.78246e+00, -1.27815e+00);
-    refforces[368] = Vec3(-1.02109e+00, -1.21921e+00, -1.26706e+00);
-    refforces[369] = Vec3( 3.05185e+00,  3.99846e+01,  5.39902e+01);
-    refforces[370] = Vec3( 1.93675e+00, -1.28986e+00, -1.21655e+00);
-    refforces[371] = Vec3(-4.14952e+00, -1.06509e+00, -9.52154e-01);
-    refforces[372] = Vec3( 6.06000e+01,  2.50872e+01,  4.29623e+01);
-    refforces[373] = Vec3(-3.50284e+00, -1.47848e+00, -7.98049e-01);
-    refforces[374] = Vec3(-8.66087e-01, -1.73056e+00, -6.59554e-01);
-
-    // Find the exclusion information
-    vector<set<int> > exclusions(NATOMS);
-    for (int i = 0; i < forceField->getNumExceptions(); i++) {
-        int particle1, particle2;
-        double chargeProd, sigma, epsilon;
-        forceField->getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
-        exclusions[particle1].insert(particle2);
-        exclusions[particle2].insert(particle1);
-    }
-
-    const double longcutoff = 35.0*OpenMM::NmPerAngstrom;
-    const double longcutoff2 = longcutoff*longcutoff;
-    const double cutoff2 = cutoff*cutoff;
-    const double cutoff6inv = 1.0 / (cutoff2*cutoff2*cutoff2);
-    const int nboxes = ceil(longcutoff/boxEdgeLength);
-
-    double refenergy = 0.0;
-    // Loop over home box first...
-    for(int i = 0; i < NATOMS; ++ i){
-        for(int j = i+1; j < NATOMS; ++j){
-            Vec3 dR = positions[i] - positions[j];
-            double R2 = dR[0]*dR[0] + dR[1]*dR[1] + dR[2]*dR[2];
-            double sig2 = (sigvals[i] + sigvals[j])*(sigvals[i] + sigvals[j]) / R2;
-            double sig6 = sig2*sig2*sig2;
-            double eps = epsvals[i]*epsvals[j];
-            refenergy += 2.0*eps*(sig6-1.0)*sig6;
-            if(R2 < cutoff2){
-                // Add a shift term for direct space parts withing t
-                refenergy += 2.0*eps*(pow(sigvals[i]+sigvals[j], 6) - 64.0*pow(sigvals[i]*sigvals[j], 3))*cutoff6inv;
-            }
-        }
-    }
-
-    // ... back out exclusions
-    for (int ii = 0; ii < NATOMS; ii++){
-        for (set<int>::const_iterator iter = exclusions[ii].begin(); iter != exclusions[ii].end(); ++iter) {
-            if (*iter > ii) {
-                int i = ii;
-                int j = *iter;
-                Vec3 dR = positions[i] - positions[j];
-                double R2 = dR[0]*dR[0] + dR[1]*dR[1] + dR[2]*dR[2];
-                double sig2 = (sigvals[i] + sigvals[j])*(sigvals[i] + sigvals[j]) / R2;
-                double sig6 = sig2*sig2*sig2;
-                double eps = epsvals[i]*epsvals[j];
-                refenergy -= 2.0*eps*(sig6-1.0)*sig6;
-                if(R2 < cutoff2){
-                    // Add a shift term for direct space parts withing t
-                    refenergy -= 2.0*eps*(pow(sigvals[i]+sigvals[j], 6) - 64.0*pow(sigvals[i]*sigvals[j], 3))*cutoff6inv;
-                }
-            }
-        }
-    }
-
-    // ... and now add in the image box terms
-    for(int bx = -nboxes; bx <= nboxes; ++bx){
-        for(int by = -nboxes; by <= nboxes; ++by){
-            for(int bz = -nboxes; bz <= nboxes; ++bz){
-                if(bx==0 && by==0 && bz==0) continue;
-                Vec3 offset(bx*boxEdgeLength, by*boxEdgeLength, bz*boxEdgeLength);
-                for(int i = 0; i < NATOMS; ++ i){
-                    for(int j = 0; j < NATOMS; ++j){
-                        Vec3 dR = positions[i] - positions[j] + offset;
-                        double R2 = dR[0]*dR[0] + dR[1]*dR[1] + dR[2]*dR[2];
-                        if(R2 > longcutoff2) continue;
-                        double sig2 = (sigvals[i] + sigvals[j])*(sigvals[i] + sigvals[j]) / R2;
-                        double sig6 = sig2*sig2*sig2;
-                        double eps = epsvals[i]*epsvals[j];
-                        refenergy += eps*(sig6-1.0)*sig6;
-                        if(R2 < cutoff2){
-                            // Add a shift term for direct space parts withing teh
-                            refenergy += eps*(pow(sigvals[i]+sigvals[j], 6) - 64.0*pow(sigvals[i]*sigvals[j], 3))*cutoff6inv;
-                        }
-                    }
-                }
-            }
-        }
-    }
-    refenergy *= 0.5;
-    // For this test the reference energy is -294.078 kJ/mol, while the difference between DPME and 30A cutoffs
-    // is just 0.064 kJ/mol.  The difference is due to the fact that arithmetic mean combination rules are used
-    // up to the cutoff, while the reciprocal space uses the geometric mean.  See DOI: 10.1021/acs.jctc.5b00726
-    ASSERT_EQUAL_TOL(refenergy, energy, 5E-4);
-    ASSERT_EQUAL_TOL(gromacs_energy, energy, 5E-5);
-
-    // Forces accumulated in single precision are tested to a more permissive criterion; the double
-    // precision platform can match to 5E-5.
-    for(int n = 0; n < NATOMS; ++n)
-        ASSERT_EQUAL_VEC(refforces[n], forces[n], 1E-4);
-
-}
-
-
 void testCoulomb() {
     System system;
     system.addParticle(1.0);
@@ -2346,10 +690,6 @@ void runPlatformTests();
 int main(int argc, char* argv[]) {
     try {
         initializeTests(argc, argv);
-        test_water2_dpme_energies_forces_no_exclusions();
-        test_water2_dpme_energies_forces_with_exclusions();
-        test_water125_dpme_vs_long_cutoff_no_exclusions();
-        test_water125_dpme_vs_long_cutoff_with_exclusions();
         testCoulomb();
         testLJ();
         testExclusionsAnd14();