Combined the Ewald and PME routines into a single one and added calculation of vdW interactions.

platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.cpp

@@ -31,8 +31,7 @@
 #include "../SimTKUtilities/SimTKOpenMMUtilities.h"
 #include "ReferenceLJCoulombIxn.h"
 #include "ReferenceForce.h"
-
-#include "pme.h"
+#include "PME.h"
 
 // In case we're using some primitive version of Visual Studio this will
 // make sure that erf() and erfc() are defined.
@@ -85,13 +84,13 @@ ReferenceLJCoulombIxn::~ReferenceLJCoulombIxn( ){
      --------------------------------------------------------------------------------------- */
 
   int ReferenceLJCoulombIxn::setUseCutoff( RealOpenMM distance, const OpenMM::NeighborList& neighbors, RealOpenMM solventDielectric ) {
-    
+
     cutoff = true;
     cutoffDistance = distance;
     neighborList = &neighbors;
     krf = pow(cutoffDistance, -3.0f)*(solventDielectric-1.0f)/(2.0f*solventDielectric+1.0f);
     crf = (1.0f/cutoffDistance)*(3.0f*solventDielectric)/(2.0f*solventDielectric+1.0f);
-            
+
     return ReferenceForce::DefaultReturn;
   }
 
@@ -193,9 +192,9 @@ int ReferenceLJCoulombIxn::getDerivedParameters( RealOpenMM c6, RealOpenMM c12,
       parameters[SigIndex] = half;
 
    } else {
- 
+
       parameters[EpsIndex]    = c6*SQRT( one/c12 );
-   
+
       parameters[SigIndex]    = POW( (c12/c6), oneSixth );
       parameters[SigIndex]   *= half;
    }
@@ -222,9 +221,9 @@ int ReferenceLJCoulombIxn::getDerivedParameters( RealOpenMM c6, RealOpenMM c12,
    @param totalEnergy      total energy
 
    @return ReferenceForce::DefaultReturn
-      
+
    --------------------------------------------------------------------------------------- */
- 
+
 int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** atomCoordinates,
                                              RealOpenMM** atomParameters, int** exclusions,
                                              RealOpenMM* fixedParameters, RealOpenMM** forces,
@@ -233,40 +232,80 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
     #include "../SimTKUtilities/RealTypeSimTk.h"
     typedef std::complex<RealOpenMM> d_complex;
 
-    int kmax = std::max(numRx, std::max(numRy,numRz));
-
-    RealOpenMM  factorEwald = -1 / (4*alphaEwald*alphaEwald);
-
-    RealOpenMM SQRT_PI = sqrt(PI);
-    RealOpenMM TWO_PI = 2.0 * PI;
     static const RealOpenMM epsilon     =  1.0;
     static const RealOpenMM one         =  1.0;
+    static const RealOpenMM six         =  6.0;
+    static const RealOpenMM twelve      = 12.0;
 
-    RealOpenMM recipCoeff = (RealOpenMM)(4*PI/(periodicBoxSize[0] * periodicBoxSize[1] * periodicBoxSize[2]) /epsilon);
+    int kmax                            = std::max(numRx, std::max(numRy,numRz));
+    RealOpenMM  factorEwald             = -1 / (4*alphaEwald*alphaEwald);
+    RealOpenMM SQRT_PI                  = sqrt(PI);
+    RealOpenMM TWO_PI                   = 2.0 * PI;
+    RealOpenMM recipCoeff               = (RealOpenMM)(4*PI/(periodicBoxSize[0] * periodicBoxSize[1] * periodicBoxSize[2]) /epsilon);
 
-    RealOpenMM selfEwaldEnergy = 0.0;
-    RealOpenMM realSpaceEwaldEnergy = 0.0;
-    RealOpenMM recipEnergy = 0.0;
+    RealOpenMM selfEwaldEnergy          = 0.0;
+    RealOpenMM realSpaceEwaldEnergy     = 0.0;
+    RealOpenMM recipEnergy              = 0.0;
+    RealOpenMM vdwEnergy                = 0.0;
 
 // **************************************************************************************
 // SELF ENERGY
 // **************************************************************************************
 
     for( int atomID = 0; atomID < numberOfAtoms; atomID++ ){
-        selfEwaldEnergy = selfEwaldEnergy + atomParameters[atomID][QIndex]*atomParameters[atomID][QIndex];
+        selfEwaldEnergy = atomParameters[atomID][QIndex]*atomParameters[atomID][QIndex] * alphaEwald/SQRT_PI;
+
+        if( totalEnergy )
+              *totalEnergy -= selfEwaldEnergy;
+
+        if( energyByAtom ){
+           energyByAtom[atomID] -= selfEwaldEnergy;
+        }
     }
-       selfEwaldEnergy = selfEwaldEnergy * alphaEwald/SQRT_PI ;
 
 // **************************************************************************************
 // RECIPROCAL SPACE EWALD ENERGY AND FORCES
 // **************************************************************************************
 
-// setup reciprocal box
+    // PME
+
+  if (pme) {
+	pme_t          pmedata; /* abstract handle for PME data */
+	int            ngrid[3];
+	RealOpenMM virial[3][3];
+
+	/* PME grid dimensions.
+	 * We typically want to set this as the spacing rather than absolute dimensions, but
+	 * to be able to reproduce results from other programs (e.g. Gromacs) we need to be
+	 * able to set exact grid dimenisions occasionally.
+	 */
+	ngrid[0] = 16;
+	ngrid[1] = 16;
+	ngrid[2] = 16;
+
+	pme_init(&pmedata,alphaEwald,numberOfAtoms,ngrid,4,1);
+
+	pme_exec(pmedata,atomCoordinates,forces,atomParameters,periodicBoxSize,&recipEnergy,virial);
+
+	if( totalEnergy )
+       *totalEnergy += recipEnergy;
+
+    if( energyByAtom )
+        for(int n = 0; n < numberOfAtoms; n++)
+            energyByAtom[n] += recipEnergy;
+
+  }
+
+    // Ewald method
+
+  else if (ewald) {
+
+    // setup reciprocal box
 
            RealOpenMM recipBoxSize[3] = { TWO_PI / periodicBoxSize[0], TWO_PI / periodicBoxSize[1], TWO_PI / periodicBoxSize[2]};
 
 
-// setup K-vectors
+    // setup K-vectors
 
   #define EIR(x, y, z) eir[(x)*numberOfAtoms*3+(y)*3+z]
   vector<d_complex> eir(kmax*numberOfAtoms*3);
@@ -292,7 +331,8 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
         EIR(j, i, m) = EIR(j-1, i, m) * EIR(1, i, m);
   }
 
-// calculate reciprocal space energy and forces
+    // calculate reciprocal space energy and forces
+
 
     int lowry = 0;
     int lowrz = 1;
@@ -335,231 +375,92 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
             ss += tab_qxyz[n].imag();
           }
 
-
           RealOpenMM kz = rz * recipBoxSize[2];
           RealOpenMM k2 = kx * kx + ky * ky + kz * kz;
           RealOpenMM ak = exp(k2*factorEwald) / k2;
 
-          recipEnergy += ak * ( cs * cs + ss * ss);
-
           for(int n = 0; n < numberOfAtoms; n++) {
             RealOpenMM force = ak * (cs * tab_qxyz[n].imag() - ss * tab_qxyz[n].real());
             forces[n][0] += 2 * recipCoeff * force * kx ;
             forces[n][1] += 2 * recipCoeff * force * ky ;
             forces[n][2] += 2 * recipCoeff * force * kz ;
           }
+
+          recipEnergy = recipCoeff * ak * ( cs * cs + ss * ss);
+
+          if( totalEnergy )
+             *totalEnergy += recipEnergy;
+
+          if( energyByAtom )
+             for(int n = 0; n < numberOfAtoms; n++)
+               energyByAtom[n] += recipEnergy;
+
           lowrz = 1 - numRz;
         }
         lowry = 1 - numRy;
       }
     }
+  }
 
-  recipEnergy *= recipCoeff;
-
-
+  else {
+      std::stringstream message;
+      message << " Wrong method for Ewald summation, Aborting" << std::endl;
+      SimTKOpenMMLog::printError( message );
+  }
 
 // **************************************************************************************
 // SHORT-RANGE ENERGY AND FORCES
 // **************************************************************************************
 
-       RealOpenMM deltaR[2][ReferenceForce::LastDeltaRIndex];
-
-       for( int atomID1 = 0; atomID1 < numberOfAtoms; atomID1++ ){
-        for( int atomID2 = atomID1 + 1; atomID2 < numberOfAtoms; atomID2++ ){
-
-          ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomID2], atomCoordinates[atomID1], periodicBoxSize, deltaR[0] );  
-          RealOpenMM r         = deltaR[0][ReferenceForce::RIndex];
-          RealOpenMM r2        = deltaR[0][ReferenceForce::R2Index];
-          RealOpenMM inverseR  = one/(deltaR[0][ReferenceForce::RIndex]);
+    for (int i = 0; i < (int) neighborList->size(); i++) {
+       OpenMM::AtomPair pair = (*neighborList)[i];
+       int ii = pair.first;
+       int jj = pair.second;
 
-          realSpaceEwaldEnergy = 
-              (RealOpenMM)(realSpaceEwaldEnergy + atomParameters[atomID1][QIndex]*atomParameters[atomID2][QIndex]*inverseR*erfc(alphaEwald*r)); 
-        }
-       }
-
-       // allocate and initialize exclusion array
-
-       vector<int> exclusionIndices(numberOfAtoms);
-       for( int ii = 0; ii < numberOfAtoms; ii++ ){
-          exclusionIndices[ii] = -1;
-       }
-
-       for( int ii = 0; ii < numberOfAtoms; ii++ ){
-
-          // set exclusions
-
-          for( int jj = 1; jj <= exclusions[ii][0]; jj++ ){
-             exclusionIndices[exclusions[ii][jj]] = ii;
-          }
-
-          // loop over atom pairs
-
-          for( int jj = ii+1; jj < numberOfAtoms; jj++ ){
-
-             if( exclusionIndices[jj] != ii ){
-
-       ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxSize, deltaR[0] );  
+       RealOpenMM deltaR[2][ReferenceForce::LastDeltaRIndex];
+       ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxSize, deltaR[0] );
        RealOpenMM r         = deltaR[0][ReferenceForce::RIndex];
        RealOpenMM r2        = deltaR[0][ReferenceForce::R2Index];
        RealOpenMM inverseR  = one/(deltaR[0][ReferenceForce::RIndex]);
        RealOpenMM alphaR    = alphaEwald * r;
- 
-       realSpaceEwaldEnergy = 
-           (RealOpenMM)(realSpaceEwaldEnergy + atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR*erfc(alphaR)); 
-       RealOpenMM dEdR = atomParameters[ii][QIndex] * atomParameters[jj][QIndex] * inverseR * inverseR * inverseR;
-       dEdR = (RealOpenMM)(dEdR * (erfc(alphaR) + 2 * alphaR * exp ( - alphaR * alphaR) / SQRT_PI ));
- 
-                for( int kk = 0; kk < 3; kk++ ){
-                   RealOpenMM force  = dEdR*deltaR[0][kk];
-                   forces[ii][kk]   += force;
-                   forces[jj][kk]   -= force;
-                } 
-             }
-          }
-       }
-
-// ***********************************************************************
-
-        if( totalEnergy ) {
-              *totalEnergy += recipEnergy + realSpaceEwaldEnergy - selfEwaldEnergy;
-        }
-
-
-   return ReferenceForce::DefaultReturn;
-}
-
-
-/**---------------------------------------------------------------------------------------
-
-   Calculate PME ixn
-
-   @param numberOfAtoms    number of atoms
-   @param atomCoordinates  atom coordinates
-   @param atomParameters   atom parameters                             atomParameters[atomIndex][paramterIndex]
-   @param exclusions       atom exclusion indices                      exclusions[atomIndex][atomToExcludeIndex]
-                           exclusions[atomIndex][0] = number of exclusions
-                           exclusions[atomIndex][1-no.] = atom indices of atoms to excluded from
-                           interacting w/ atom atomIndex
-   @param fixedParameters  non atom parameters (not currently used)
-   @param forces           force array (forces added)
-   @param energyByAtom     atom energy
-   @param totalEnergy      total energy
-
-   @return ReferenceForce::DefaultReturn
-      
-   --------------------------------------------------------------------------------------- */
- 
-int ReferenceLJCoulombIxn::calculatePMEIxn( int numberOfAtoms, RealOpenMM** atomCoordinates,
-                                             RealOpenMM** atomParameters, int** exclusions,
-                                             RealOpenMM* fixedParameters, RealOpenMM** forces,
-                                             RealOpenMM* energyByAtom, RealOpenMM* totalEnergy) const {
-
-
-
-    RealOpenMM SQRT_PI = sqrt(PI);
-    static const RealOpenMM one         =  1.0;
-
-    RealOpenMM selfEwaldEnergy = 0.0;
-    RealOpenMM realSpaceEwaldEnergy = 0.0;
-    RealOpenMM recipEnergy = 0.0;
-
-
-// **************************************************************************************
-// SELF ENERGY
-// **************************************************************************************
-
-    for( int atomID = 0; atomID < numberOfAtoms; atomID++ ){
-        selfEwaldEnergy = selfEwaldEnergy + atomParameters[atomID][QIndex]*atomParameters[atomID][QIndex];
-    }
-       selfEwaldEnergy = selfEwaldEnergy * alphaEwald/SQRT_PI ;
-
-// **************************************************************************************
-// RECIPROCAL SPACE EWALD ENERGY AND FORCES
-// **************************************************************************************
-
-	pme_t          pmedata; /* abstract handle for PME data */
-	int            ngrid[3];
-	RealOpenMM virial[3][3];
-	
-	/* PME grid dimensions.
-	 * We typically want to set this as the spacing rather than absolute dimensions, but
-	 * to be able to reproduce results from other programs (e.g. Gromacs) we need to be
-	 * able to set exact grid dimenisions occasionally.
-	 */
-	ngrid[0] = 16;
-	ngrid[1] = 16;
-	ngrid[2] = 16;
-	
-	pme_init(&pmedata,alphaEwald,numberOfAtoms,ngrid,4,1);
-	
-	pme_exec(pmedata,atomCoordinates,forces,atomParameters,periodicBoxSize,&recipEnergy,virial);
-	
-// **************************************************************************************
-// SHORT-RANGE ENERGY AND FORCES
-// **************************************************************************************
 
-       RealOpenMM deltaR[2][ReferenceForce::LastDeltaRIndex];
 
-       for( int atomID1 = 0; atomID1 < numberOfAtoms; atomID1++ ){
-        for( int atomID2 = atomID1 + 1; atomID2 < numberOfAtoms; atomID2++ ){
+       RealOpenMM dEdR      = atomParameters[ii][QIndex] * atomParameters[jj][QIndex] * inverseR * inverseR * inverseR;
+                  dEdR      = (RealOpenMM)(dEdR * (erfc(alphaR) + 2 * alphaR * exp ( - alphaR * alphaR) / SQRT_PI ));
 
-          ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomID2], atomCoordinates[atomID1], periodicBoxSize, deltaR[0] );  
-          RealOpenMM r         = deltaR[0][ReferenceForce::RIndex];
-          RealOpenMM r2        = deltaR[0][ReferenceForce::R2Index];
-          RealOpenMM inverseR  = one/(deltaR[0][ReferenceForce::RIndex]);
+       RealOpenMM sig       = atomParameters[ii][SigIndex] +  atomParameters[jj][SigIndex];
+       RealOpenMM sig2      = inverseR*sig;
+                  sig2     *= sig2;
+       RealOpenMM sig6      = sig2*sig2*sig2;
+       RealOpenMM eps       = atomParameters[ii][EpsIndex]*atomParameters[jj][EpsIndex];
 
-          realSpaceEwaldEnergy = 
-              (RealOpenMM)(realSpaceEwaldEnergy + atomParameters[atomID1][QIndex]*atomParameters[atomID2][QIndex]*inverseR*erfc(alphaEwald*r)); 
-        }
-       }
+                  dEdR     += eps*( twelve*sig6 - six )*sig6;
 
-       // allocate and initialize exclusion array
+       // accumulate forces
 
-       vector<int> exclusionIndices(numberOfAtoms);
-       for( int ii = 0; ii < numberOfAtoms; ii++ ){
-          exclusionIndices[ii] = -1;
+       for( int kk = 0; kk < 3; kk++ ){
+          RealOpenMM force  = dEdR*deltaR[0][kk];
+          forces[ii][kk]   += force;
+          forces[jj][kk]   -= force;
        }
 
-       for( int ii = 0; ii < numberOfAtoms; ii++ ){
+       // accumulate energies
 
-          // set exclusions
+       realSpaceEwaldEnergy = atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR*erfc(alphaR);
+       vdwEnergy = eps*(sig6-one)*sig6;
 
-          for( int jj = 1; jj <= exclusions[ii][0]; jj++ ){
-             exclusionIndices[exclusions[ii][jj]] = ii;
-          }
-
-          // loop over atom pairs
-
-          for( int jj = ii+1; jj < numberOfAtoms; jj++ ){
+        if( totalEnergy )
+              *totalEnergy += realSpaceEwaldEnergy + vdwEnergy;
 
-             if( exclusionIndices[jj] != ii ){
+        if( energyByAtom ){
+           energyByAtom[ii] += realSpaceEwaldEnergy + vdwEnergy;
+           energyByAtom[jj] += realSpaceEwaldEnergy + vdwEnergy;
+        }
 
-       ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxSize, deltaR[0] );  
-       RealOpenMM r         = deltaR[0][ReferenceForce::RIndex];
-       RealOpenMM r2        = deltaR[0][ReferenceForce::R2Index];
-       RealOpenMM inverseR  = one/(deltaR[0][ReferenceForce::RIndex]);
-       RealOpenMM alphaR    = alphaEwald * r;
- 
-       realSpaceEwaldEnergy = 
-           (RealOpenMM)(realSpaceEwaldEnergy + atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR*erfc(alphaR)); 
-       RealOpenMM dEdR = atomParameters[ii][QIndex] * atomParameters[jj][QIndex] * inverseR * inverseR * inverseR;
-       dEdR = (RealOpenMM)(dEdR * (erfc(alphaR) + 2 * alphaR * exp ( - alphaR * alphaR) / SQRT_PI ));
- 
-                for( int kk = 0; kk < 3; kk++ ){
-                   RealOpenMM force  = dEdR*deltaR[0][kk];
-                   forces[ii][kk]   += force;
-                   forces[jj][kk]   -= force;
-                } 
-             }
-          }
-       }
+    }
 
 // ***********************************************************************
 
-        if( totalEnergy ) {
-              *totalEnergy += recipEnergy + realSpaceEwaldEnergy - selfEwaldEnergy;
-        }
-
 
    return ReferenceForce::DefaultReturn;
 }
@@ -582,18 +483,16 @@ int ReferenceLJCoulombIxn::calculatePMEIxn( int numberOfAtoms, RealOpenMM** atom
    @param totalEnergy      total energy
 
    @return ReferenceForce::DefaultReturn
-      
+
    --------------------------------------------------------------------------------------- */
-    
+
 int ReferenceLJCoulombIxn::calculatePairIxn( int numberOfAtoms, RealOpenMM** atomCoordinates,
                                              RealOpenMM** atomParameters, int** exclusions,
                                              RealOpenMM* fixedParameters, RealOpenMM** forces,
                                              RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const {
 
-   if (ewald)
+   if (ewald || pme)
         return calculateEwaldIxn(numberOfAtoms, atomCoordinates, atomParameters, exclusions, fixedParameters, forces, energyByAtom, totalEnergy);
-   if (pme)
-        return calculatePMEIxn(numberOfAtoms, atomCoordinates, atomParameters, exclusions, fixedParameters, forces, energyByAtom, totalEnergy);
    if (cutoff) {
        for (int i = 0; i < (int) neighborList->size(); i++) {
            OpenMM::AtomPair pair = (*neighborList)[i];
@@ -681,9 +580,9 @@ int ReferenceLJCoulombIxn::calculateOneIxn( int ii, int jj, RealOpenMM** atomCoo
     // get deltaR, R2, and R between 2 atoms
 
     if (periodic)
-        ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxSize, deltaR[0] );  
+        ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxSize, deltaR[0] );
     else
-        ReferenceForce::getDeltaR( atomCoordinates[jj], atomCoordinates[ii], deltaR[0] );  
+        ReferenceForce::getDeltaR( atomCoordinates[jj], atomCoordinates[ii], deltaR[0] );
 
     RealOpenMM r2        = deltaR[0][ReferenceForce::R2Index];
     RealOpenMM inverseR  = one/(deltaR[0][ReferenceForce::RIndex]);
@@ -699,7 +598,7 @@ int ReferenceLJCoulombIxn::calculateOneIxn( int ii, int jj, RealOpenMM** atomCoo
                else
                    dEdR += atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR;
                dEdR     *= inverseR*inverseR;
-               
+
     // accumulate forces
 
     for( int kk = 0; kk < 3; kk++ ){
@@ -726,7 +625,7 @@ int ReferenceLJCoulombIxn::calculateOneIxn( int ii, int jj, RealOpenMM** atomCoo
         }
     }
 
-    // debug 
+    // debug
 
     if( debug == ii ){
        static bool printHeader = false;
@@ -734,11 +633,11 @@ int ReferenceLJCoulombIxn::calculateOneIxn( int ii, int jj, RealOpenMM** atomCoo
        message << methodName;
        message << std::endl;
        int pairArray[2] = { ii, jj };
-       if( !printHeader  ){  
+       if( !printHeader  ){
           printHeader = true;
           message << std::endl;
           message << methodName.c_str() << " a0 k [c q p s] r1 r2  angle dt rp p[] dot cosine angle dEdR*r F[]" << std::endl;
-       }   
+       }
 
        message << std::endl;
        for( int kk = 0; kk < 2; kk++ ){

platforms/reference/tests/TestReferenceEwald.cpp

@@ -49,6 +49,41 @@ using namespace std;
 
 const double TOL = 1e-5;
 
+void testLargeSystem() {
+    ReferencePlatform platform;
+    System system;
+    for (int i = 0; i < 500; i++)
+        system.addParticle(22.99);
+    for (int i = 0; i < 500; i++)
+        system.addParticle(35.45);
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* nonbonded = new NonbondedForce();
+    for (int i = 0; i < 500; i++)
+        nonbonded->addParticle(1.0, 1.0,0.0);
+//        nonbonded->addParticle(1.0, 0.33284,0.0115897);
+    for (int i = 0; i < 500; i++)
+        nonbonded->addParticle(-1.0, 1.0,0.0);
+//        nonbonded->addParticle(-1.0, 0.440104,0.4184);
+    nonbonded->setNonbondedMethod(NonbondedForce::PME);
+    const double cutoff = 1.0;
+    nonbonded->setCutoffDistance(cutoff);
+    nonbonded->setPeriodicBoxVectors(Vec3(2.82, 0, 0), Vec3(0, 2.82, 0), Vec3(0, 0, 2.82));
+    nonbonded->setEwaldErrorTolerance(TOL);
+    system.addForce(nonbonded);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(1000);
+    #include "nacl_crystal.dat"
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+//    cout << "force 0: " << forces[0] << endl;
+//    cout << "force 1: " << forces[1] << endl;
+    cout << "PotentialEnergy: " << state.getPotentialEnergy() << endl;
+    ASSERT_EQUAL_TOL(-430355, state.getPotentialEnergy(), 100*TOL);
+//    ASSERT_EQUAL_VEC(Vec3(-123.711, 64.1877, -302.716), forces[0], 10*TOL);
+//    ASSERT_EQUAL_VEC(Vec3(123.711, -64.1877, 302.716), forces[1], 10*TOL);
+}
+
 void testEwald() {
     ReferencePlatform platform;
     System system;
@@ -58,6 +93,9 @@ void testEwald() {
     NonbondedForce* nonbonded = new NonbondedForce();
     nonbonded->addParticle(1.0, 1, 0);
     nonbonded->addParticle(-1.0, 1, 0);
+// Sodium Chloride
+//    nonbonded->addParticle(1.0, 0.33284,0.0115897);
+//    nonbonded->addParticle(-1.0, 0.440104,0.4184);
     nonbonded->setNonbondedMethod(NonbondedForce::Ewald);
     const double cutoff = 2.0;
     nonbonded->setCutoffDistance(cutoff);
@@ -71,9 +109,9 @@ void testEwald() {
     context.setPositions(positions);
     State state = context.getState(State::Forces | State::Energy);
     const vector<Vec3>& forces = state.getForces();
-    cout << "force 0: " << forces[0] << endl;
-    cout << "force 1: " << forces[1] << endl;
-    cout << "PotentialEnergy: " << state.getPotentialEnergy() << endl;
+//    cout << "force 0: " << forces[0] << endl;
+//    cout << "force 1: " << forces[1] << endl;
+//    cout << "PotentialEnergy: " << state.getPotentialEnergy() << endl;
     ASSERT_EQUAL_VEC(Vec3(-123.711, 64.1877, -302.716), forces[0], 10*TOL);
     ASSERT_EQUAL_VEC(Vec3(123.711, -64.1877, 302.716), forces[1], 10*TOL);
 }
@@ -101,48 +139,7 @@ void testPME() {
     system.addForce(nonbonded);
     Context context(system, integrator, platform);
     vector<Vec3> positions(42);
-positions[0] = Vec3( 0.23,0.628,0.113);
-positions[1] = Vec3(0.137,0.626, 0.15);
-positions[2] = Vec3(0.231,0.589,0.021);
-positions[3] = Vec3(-0.307,-0.351,0.703);
-positions[4] = Vec3(-0.364,-0.367,0.784);
-positions[5] = Vec3(-0.366,-0.341,0.623);
-positions[6] = Vec3(-0.569,-0.634,-0.439);
-positions[7] = Vec3(-0.532,-0.707,-0.497);
-positions[8] = Vec3(-0.517,-0.629,-0.354);
-positions[9] = Vec3(-0.871, 0.41,-0.62);
-positions[10] = Vec3(-0.948,0.444,-0.566);
-positions[11] = Vec3(-0.905,0.359,-0.699);
-positions[12] = Vec3(0.249,-0.077,-0.621);
-positions[13] = Vec3(0.306,-0.142,-0.571);
-positions[14] = Vec3(0.233,-0.11,-0.714);
-positions[15] = Vec3(0.561,0.222,-0.715);
-positions[16] = Vec3(0.599,0.138,-0.678);
-positions[17] = Vec3(0.473,0.241,-0.671);
-positions[18] = Vec3(-0.515,-0.803,-0.628);
-positions[19] = Vec3(-0.491,-0.866,-0.702);
-positions[20] = Vec3(-0.605,-0.763,-0.646);
-positions[21] = Vec3(-0.021,0.175,-0.899);
-positions[22] = Vec3(0.018, 0.09,-0.935);
-positions[23] = Vec3(-0.119,0.177,-0.918);
-positions[24] = Vec3(-0.422,0.856,-0.464);
-positions[25] = Vec3(-0.479,0.908,-0.527);
-positions[26] = Vec3(-0.326,0.868,-0.488);
-positions[27] = Vec3(-0.369,-0.095,-0.903);
-positions[28] = Vec3(-0.336,-0.031,-0.972);
-positions[29] = Vec3(-0.303,-0.101,-0.828);
-positions[30] = Vec3(0.594,0.745,0.652);
-positions[31] = Vec3(0.644, 0.83,0.633);
-positions[32] = Vec3(0.506,0.747,0.604);
-positions[33] = Vec3(-0.157,-0.375,-0.758);
-positions[34] = Vec3(-0.25, -0.4,-0.785);
-positions[35] = Vec3(-0.131,-0.425,-0.676);
-positions[36] = Vec3(0.618,-0.295,-0.578);
-positions[37] = Vec3(0.613,-0.213,-0.521);
-positions[38] = Vec3(0.707,-0.298,-0.623);
-positions[39] = Vec3(0.039,-0.785,  0.3);
-positions[40] = Vec3(0.138,-0.796,0.291);
-positions[41] = Vec3(-0.001,-0.871,0.332);
+    #include "water.dat"
     context.setPositions(positions);
     State state = context.getState(State::Forces | State::Energy);
     const vector<Vec3>& forces = state.getForces();
@@ -156,8 +153,9 @@ positions[41] = Vec3(-0.001,-0.871,0.332);
 
 int main() {
     try {
-        testEwald();
-        testPME();
+     testLargeSystem();
+     testEwald();
+     testPME();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;

platforms/reference/tests/water.dat

+positions[0] = Vec3( 0.23,0.628,0.113);
+positions[1] = Vec3(0.137,0.626, 0.15);
+positions[2] = Vec3(0.231,0.589,0.021);
+positions[3] = Vec3(-0.307,-0.351,0.703);
+positions[4] = Vec3(-0.364,-0.367,0.784);
+positions[5] = Vec3(-0.366,-0.341,0.623);
+positions[6] = Vec3(-0.569,-0.634,-0.439);
+positions[7] = Vec3(-0.532,-0.707,-0.497);
+positions[8] = Vec3(-0.517,-0.629,-0.354);
+positions[9] = Vec3(-0.871, 0.41,-0.62);
+positions[10] = Vec3(-0.948,0.444,-0.566);
+positions[11] = Vec3(-0.905,0.359,-0.699);
+positions[12] = Vec3(0.249,-0.077,-0.621);
+positions[13] = Vec3(0.306,-0.142,-0.571);
+positions[14] = Vec3(0.233,-0.11,-0.714);
+positions[15] = Vec3(0.561,0.222,-0.715);
+positions[16] = Vec3(0.599,0.138,-0.678);
+positions[17] = Vec3(0.473,0.241,-0.671);
+positions[18] = Vec3(-0.515,-0.803,-0.628);
+positions[19] = Vec3(-0.491,-0.866,-0.702);
+positions[20] = Vec3(-0.605,-0.763,-0.646);
+positions[21] = Vec3(-0.021,0.175,-0.899);
+positions[22] = Vec3(0.018, 0.09,-0.935);
+positions[23] = Vec3(-0.119,0.177,-0.918);
+positions[24] = Vec3(-0.422,0.856,-0.464);
+positions[25] = Vec3(-0.479,0.908,-0.527);
+positions[26] = Vec3(-0.326,0.868,-0.488);
+positions[27] = Vec3(-0.369,-0.095,-0.903);
+positions[28] = Vec3(-0.336,-0.031,-0.972);
+positions[29] = Vec3(-0.303,-0.101,-0.828);
+positions[30] = Vec3(0.594,0.745,0.652);
+positions[31] = Vec3(0.644, 0.83,0.633);
+positions[32] = Vec3(0.506,0.747,0.604);
+positions[33] = Vec3(-0.157,-0.375,-0.758);
+positions[34] = Vec3(-0.25, -0.4,-0.785);
+positions[35] = Vec3(-0.131,-0.425,-0.676);
+positions[36] = Vec3(0.618,-0.295,-0.578);
+positions[37] = Vec3(0.613,-0.213,-0.521);
+positions[38] = Vec3(0.707,-0.298,-0.623);
+positions[39] = Vec3(0.039,-0.785,  0.3);
+positions[40] = Vec3(0.138,-0.796,0.291);
+positions[41] = Vec3(-0.001,-0.871,0.332);