Implemented the Ewald summation method in the reference platform. Preliminary...

Implemented the Ewald summation method in the reference platform. Preliminary cuda prototype is included.

olla/include/kernels.h

@@ -212,7 +212,8 @@ public:
     enum NonbondedMethod {
         NoCutoff = 0,
         CutoffNonPeriodic = 1,
-        CutoffPeriodic = 2
+        CutoffPeriodic = 2,
+        Ewald = 3
     };
     static std::string Name() {
         return "CalcNonbondedForce";

openmmapi/include/NonbondedForce.h

@@ -77,7 +77,11 @@ public:
          * each other particle.  Interactions beyond the cutoff distance are ignored.  Coulomb interactions closer than the
          * cutoff distance are modified using the reaction field method.
          */
-        CutoffPeriodic = 2
+        CutoffPeriodic = 2,
+        /**
+         * Using Ewald summation method for long-range electrostatics.
+         */
+        Ewald = 3
     };
     /**
      * Create a NonbondedForce.

platforms/cuda/src/CudaKernels.cpp

@@ -276,6 +276,14 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
             gpuSetPeriodicBoxSize(gpu, boxVectors[0][0], boxVectors[1][1], boxVectors[2][2]);
             method = PERIODIC;
         }
+
+        if (force.getNonbondedMethod() == NonbondedForce::Ewald) {
+            Vec3 boxVectors[3];
+            force.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
+            gpuSetPeriodicBoxSize(gpu, boxVectors[0][0], boxVectors[1][1], boxVectors[2][2]);
+            method = EWALD;
+
+        }
         data.nonbondedMethod = method;
         gpuSetCoulombParameters(gpu, 138.935485f, particle, c6, c12, q, symbol, exclusionList, method);
     }

platforms/cuda/src/kernels/cudatypes.h

@@ -233,7 +233,8 @@ enum CudaNonbondedMethod
 {
     NO_CUTOFF,
     CUTOFF,
-    PERIODIC
+    PERIODIC,
+    EWALD
 };
 
 struct cudaGmxSimulation {

platforms/cuda/src/kernels/kCalculateCDLJForces.cu

@@ -107,6 +107,8 @@ void GetCalculateCDLJForcesSim(gpuContext gpu)
 #define METHOD_NAME(a, b) a##PeriodicByWarp##b
 #include "kCalculateCDLJForces.h"
 
+// Include version of the kernel with Ewald method
+#include "kCalculateEwald.h"
 
 __global__ extern void kCalculateCDLJCutoffForces_12_kernel();
 
@@ -173,5 +175,30 @@ void kCalculateCDLJForces(gpuContext gpu)
                 kCalculateCDLJPeriodicForces_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
                         (sizeof(Atom)+sizeof(float3))*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pInteractingWorkUnit, numWithInteractions);
             LAUNCHERROR("kCalculateCDLJPeriodicForces");
+            break;
+        case EWALD:
+            kFindBlockBoundsPeriodic_kernel<<<(gpu->psGridBoundingBox->_length+63)/64, 64>>>();
+            LAUNCHERROR("kFindBlockBoundsPeriodic");
+            kFindBlocksWithInteractionsPeriodic_kernel<<<gpu->sim.interaction_blocks, gpu->sim.interaction_threads_per_block>>>();
+            LAUNCHERROR("kFindBlocksWithInteractionsPeriodic");
+            result = cudppCompact(gpu->cudpp, gpu->sim.pInteractingWorkUnit, gpu->sim.pInteractionCount,
+                    gpu->sim.pWorkUnit, gpu->sim.pInteractionFlag, gpu->sim.workUnits);
+            if (result != CUDPP_SUCCESS)
+            {
+                printf("Error in cudppCompact: %d\n", result);
+                exit(-1);
+            }
+            gpu->psInteractionCount->Download();
+            numWithInteractions = gpu->psInteractionCount->_pSysData[0];
+            kFindInteractionsWithinBlocksPeriodic_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
+                    sizeof(unsigned int)*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pInteractingWorkUnit, numWithInteractions);
+            if (gpu->bOutputBufferPerWarp)
+                kCalculateCDLJPeriodicByWarpForces_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
+                        (sizeof(Atom)+sizeof(float3))*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pInteractingWorkUnit, numWithInteractions);
+            else
+                kCalculateCDLJEwaldForces_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
+                        (sizeof(Atom)+sizeof(float3))*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pInteractingWorkUnit, numWithInteractions);
+            LAUNCHERROR("kCalculateCDLJPeriodicForces");
+
     }
-}
+}
\ No newline at end of file

platforms/cuda/src/kernels/kCalculateEwald.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) 2009 Stanford University and the Authors.           *
+ * Authors: Scott Le Grand, 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.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This file contains the kernel for evalauating nonbonded forces
+ * using the Ewald summation method.
+ */
+
+#include <cuComplex.h>
+
+/* Define complex multiply operations */
+
+__device__ cuComplex ComplexMul(cuComplex a, cuComplex b)
+{
+     cuComplex c;
+     c.x = a.x * b.x - a.y * b.y;
+     c.y = a.x * b.y + a.y * b.x;
+     return c;
+
+}
+
+__device__ cuComplex ComplexConjMul(cuComplex a, cuComplex b)
+{
+     cuComplex c;
+     c.x = a.x*b.x + a.y*b.y;
+     c.y = a.y*b.x - a.x*b.y;
+     return c;
+
+}
+
+__device__ cuComplex FloatComplexMul(float r, cuComplex a)
+{
+     cuComplex b;
+     b.x = r*a.x;
+     b.y = r*a.y;
+     return b;
+
+}
+
+
+/* This kernel is under development */
+
+
+__global__ void  kCalculateCDLJEwaldForces_kernel(unsigned int* workUnit, int numWorkUnits)
+{
+/*
+
+    extern __shared__ Atom sA[];
+    unsigned int totalWarps = cSim.nonbond_blocks*cSim.nonbond_threads_per_block/GRID;
+    unsigned int warp = (blockIdx.x*blockDim.x+threadIdx.x)/GRID;
+    int pos = warp*numWorkUnits/totalWarps;
+    int end = (warp+1)*numWorkUnits/totalWarps;
+
+//###########################################################################
+// EWALD RECIP SPACE
+//###########################################################################
+
+// *******************************************************************
+    float alphaEwald       =  3.123413f;
+    float factorEwald = -1.0 / (4*alphaEwald*alphaEwald);
+    float PI = 3.14159265358979323846f;
+
+
+    float SQRT_PI  = sqrt(PI);
+    float TWO_PI   = 2.0f * PI;
+    float epsilon  =  1.0f;
+
+/////##############################################################################
+
+    float recipCoeff = 4.0*PI/(cSim.periodicBoxSizeX * cSim.periodicBoxSizeY * cSim.periodicBoxSizeZ) /epsilon;
+
+// setup reciprocal box
+
+    float recipBoxSizeX = TWO_PI / cSim.periodicBoxSizeX;
+    float recipBoxSizeY = TWO_PI / cSim.periodicBoxSizeY;
+    float recipBoxSizeZ = TWO_PI / cSim.periodicBoxSizeZ;
+
+// setup K-vectors
+
+    unsigned int numRx = 60+1;
+    unsigned int numRy = 60+1;
+    unsigned int numRz = 60+1;
+    const int kmax = 61;
+
+    cuComplex eir[kmax][numberOfAtoms][3];
+    cuComplex tab_xy[numberOfAtoms];
+    cuComplex tab_qxyz[numberOfAtoms];
+
+
+  for(unsigned int i = 0; i < numberOfAtoms; i++) {
+    for(unsigned int m = 0; (m < 3); m++) {
+      eir[0][i][m].x = 1;
+      eir[0][i][m].y = 0;
+    }
+
+
+      eir[1][i][0].x = cos(apos.x*recipBoxSizeX);
+      eir[1][i][0].y = sin(apos.x*recipBoxSizeX);
+
+      eir[1][i][1].x = cos(apos.y*recipBoxSizeY);
+      eir[1][i][1].y = sin(apos.y*recipBoxSizeY);
+
+      eir[1][i][2].x = cos(apos.z*recipBoxSizeZ);
+      eir[1][i][2].y = sin(apos.z*recipBoxSizeZ);
+
+    for(unsigned int j=2; (j<kmax); j++) {
+      for(unsigned int m=0; (m<3); m++) {
+        eir[j][i][m] = ComplexMul (eir[j-1][i][m] , eir[1][i][m]);
+      }
+    }
+  }
+
+// *******************************************************************
+
+
+    int lowry = 0;
+    int lowrz = 1;
+
+    for(int rx = 0; rx < numRx; rx++) {
+
+      float kx = rx * recipBoxSizeX;
+
+      for(int ry = lowry; ry < numRy; ry++) {
+
+        float ky = ry * recipBoxSizeY;
+
+        if(ry >= 0) {
+          for(int n = 0; n < numberOfAtoms; n++) {
+            tab_xy[n] = ComplexMul (eir[rx][n][0] , eir[ry][n][1]);
+          }
+        }
+
+        else {
+          for(int n = 0; n < numberOfAtoms; n++) {
+            tab_xy[n]= ComplexConjMul (eir[rx][n][0] , (eir[-ry][n][1]));
+          }
+        }
+
+        for (int rz = lowrz; rz < numRz; rz++) {
+
+          float kz = rz * recipBoxSizeZ;
+          float k2 = kx * kx + ky * ky + kz * kz;
+          float ak = exp(k2*factorEwald) / k2;
+          float akv = 2.0 * ak * (1.0/k2 - factorEwald);
+
+          if( rz >= 0) {
+
+           for( int n = 0; n < numberOfAtoms; n++) {
+             tab_qxyz[n] = FloatComplexMul ( Charge[n] * ComplexMul (tab_xy[n] , eir[rz][n][2]));
+           }
+          }
+
+          else {
+
+            for( int n = 0; n < numberOfAtoms; n++) {
+              tab_qxyz[n] = FloatComplexMul( Charge[n] * ComplexConjMul (tab_xy[n] , conj(eir[-rz][n][2])) );
+            }
+          }
+
+          float cs = 0.0f;
+          float ss = 0.0f;
+
+          for( int n = 0; n < numberOfAtoms; n++) {
+            cs += tab_qxyz[n].x;
+            ss += tab_qxyz[n].y;
+          }
+
+          recipEnergy += ak * ( cs * cs + ss * ss);
+          float vir =  akv * ( cs * cs + ss * ss);
+
+          for(int n = 0; n < numberOfAtoms; n++) {
+            float force = ak * (cs * tab_qxyz[n].y - ss * tab_qxyz[n].x);
+            forces[n][0] += 2.0 * recipCoeff * force * kx ;
+            forces[n][1] += 2.0 * recipCoeff * force * ky ;
+            forces[n][2] += 2.0 * recipCoeff * force * kz ;
+          } 
+
+          lowrz = 1 - numRz;
+        }
+        lowry = 1 - numRy;
+      }
+    }
+
+
+//###########################################################################
+
+
+//###########################################################################
+// END EWALD RECIP SPACE
+//###########################################################################
+
+
+    while (pos < end)
+    {
+
+        // Extract cell coordinates from appropriate work unit
+        unsigned int x = workUnit[pos];
+        unsigned int y = ((x >> 2) & 0x7fff) << GRIDBITS;
+        bool bExclusionFlag = (x & 0x1);
+        x = (x >> 17) << GRIDBITS;
+        float4      apos;   // Local atom x, y, z, q
+        float3      af;     // Local atom fx, fy, fz
+        float dx;
+        float dy;
+        float dz;
+        float r2;
+        float r;
+        float invR;
+        float sig;
+        float sig2;
+        float sig6;
+        float eps;
+        float dEdR;
+        unsigned int tgx = threadIdx.x & (GRID - 1);
+        unsigned int tbx = threadIdx.x - tgx;
+        int tj = tgx;
+        Atom* psA = &sA[tbx];
+        unsigned int i      = x + tgx;
+        apos                = cSim.pPosq[i];
+        float2 a            = cSim.pAttr[i];
+        af.x                = 0.0f;
+        af.y                = 0.0f;
+        af.z                = 0.0f;
+        
+                int j                   = y + tgx;
+                float4 temp             = cSim.pPosq[j];
+                float2 temp1            = cSim.pAttr[j];
+                sA[threadIdx.x].x       = temp.x;
+                sA[threadIdx.x].y       = temp.y;
+                sA[threadIdx.x].z       = temp.z;
+                sA[threadIdx.x].q       = temp.w;
+                sA[threadIdx.x].sig     = temp1.x;
+                sA[threadIdx.x].eps     = temp1.y;
+
+            sA[threadIdx.x].fx      = 0.0f;
+            sA[threadIdx.x].fy      = 0.0f;
+            sA[threadIdx.x].fz      = 0.0f;
+            apos.w                 *= cSim.epsfac;
+
+                // Read fixed atom data into registers and GRF
+                unsigned int excl       = cSim.pExclusion[x * cSim.exclusionStride + y + tgx];
+                excl                    = (excl >> tgx) | (excl << (GRID - tgx));
+
+                for (unsigned int j = 0; j < GRID; j++)
+                {
+                    dx              = psA[tj].x - apos.x;
+                    dy              = psA[tj].y - apos.y;
+                    dz              = psA[tj].z - apos.z;
+
+                    dx -= floor(dx/cSim.periodicBoxSizeX+0.5f)*cSim.periodicBoxSizeX;
+                    dy -= floor(dy/cSim.periodicBoxSizeY+0.5f)*cSim.periodicBoxSizeY;
+                    dz -= floor(dz/cSim.periodicBoxSizeZ+0.5f)*cSim.periodicBoxSizeZ;
+
+                    r2              = dx * dx + dy * dy + dz * dz;
+                    r               = sqrt(r2);
+                    invR            = 1.0f / sqrt(r2);
+                    sig             = a.x + psA[tj].sig;
+                    sig2            = invR * sig;
+                    sig2           *= sig2;
+                    sig6            = sig2 * sig2 * sig2;
+                    eps             = a.y * psA[tj].eps;
+
+// ##### LJ
+                    dEdR            = eps * (12.0f * sig6 - 6.0f) * sig6;
+
+// ##### SHORT RANGE EWALD
+                    float alphaR    = alphaEwald * r;
+                    dEdR += apos.w * psA[tj].q * invR * invR * invR * (erfc(alphaR) + 2.0 * alphaR * exp ( - alphaR * alphaR) / SQRT_PI );
+
+
+                    if (!(excl & 0x1) || r2 > cSim.nonbondedCutoffSqr)
+                    {
+                        dEdR = 0.0f;
+                    }
+
+
+                    dx             *= dEdR;
+                    dy             *= dEdR;
+                    dz             *= dEdR;
+                    af.x           -= dx;
+                    af.y           -= dy;
+                    af.z           -= dz;
+                    psA[tj].fx     += dx;
+                    psA[tj].fy     += dy;
+                    psA[tj].fz     += dz;
+                    excl          >>= 1;
+                    tj              = (tj + 1) & (GRID - 1);
+                }
+
+            // Write results
+            float4 of;
+
+            of.x                                = af.x;
+            of.y                                = af.y;
+            of.z                                = af.z;
+            of.w                                = 0.0f;
+            int offset                          = x + tgx + (y >> GRIDBITS) * cSim.stride;
+            cSim.pForce4a[offset]               = of;
+            of.x                                = sA[threadIdx.x].fx;
+            of.y                                = sA[threadIdx.x].fy;
+            of.z                                = sA[threadIdx.x].fz;
+            offset                              = y + tgx + (x >> GRIDBITS) * cSim.stride;
+            cSim.pForce4a[offset]               = of;
+
+        pos++;
+    }
+*/
+}

platforms/cuda/tests/TestCudaEwald.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) 2008 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.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests all the different force terms in the reference implementation of NonbondedForce.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "OpenMMContext.h"
+#include "CudaPlatform.h"
+#include "ReferencePlatform.h"
+#include "HarmonicBondForce.h"
+#include "NonbondedForce.h"
+#include "System.h"
+#include "LangevinIntegrator.h"
+#include "VerletIntegrator.h"
+#include "internal/OpenMMContextImpl.h"
+#include "kernels/gputypes.h"
+#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
+#include "../src/sfmt/SFMT.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+const double TOL = 1e-5;
+
+void testEwald() {
+    CudaPlatform platform;
+    System system(2, 0);
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* nonbonded = new NonbondedForce(2, 0);
+    nonbonded->setParticleParameters(0, 1.0, 1, 0);
+    nonbonded->setParticleParameters(1, -1.0, 1, 0);
+    nonbonded->setNonbondedMethod(NonbondedForce::Ewald);
+    const double cutoff = 2.0;
+    nonbonded->setCutoffDistance(cutoff);
+    nonbonded->setPeriodicBoxVectors(Vec3(6, 0, 0), Vec3(0, 6, 0), Vec3(0, 0, 6));
+    system.addForce(nonbonded);
+    OpenMMContext context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(3.048000,2.764000,3.156000);
+    positions[1] = Vec3(2.809000,2.888000,2.571000);
+    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 << "    energyPoten: " << state.getPotentialEnergy() << endl;
+
+    ASSERT_EQUAL_VEC(Vec3(-123.711, 64.1877, -302.716), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(123.711, -64.1877, 302.716), forces[1], TOL);
+
+    //ASSERT_EQUAL_TOL(2*138.935485*(1.0)*(1.0+krf*1.0-crf), state.getPotentialEnergy(), TOL);
+}
+
+void testPeriodic() {
+    CudaPlatform platform;
+    System system(3, 0);
+    LangevinIntegrator integrator(0.0, 0.1, 0.01);
+    HarmonicBondForce* bonds = new HarmonicBondForce(1);
+    bonds->setBondParameters(0, 0, 1, 1, 0);
+    system.addForce(bonds);
+    NonbondedForce* nonbonded = new NonbondedForce(3, 0);
+    nonbonded->setParticleParameters(0, 1.0, 1, 0);
+    nonbonded->setParticleParameters(1, 1.0, 1, 0);
+    nonbonded->setParticleParameters(2, 1.0, 1, 0);
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    const double cutoff = 2.0;
+    nonbonded->setCutoffDistance(cutoff);
+    nonbonded->setPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 4, 0), Vec3(0, 0, 4));
+    system.addForce(nonbonded);
+    OpenMMContext context(system, integrator, platform);
+    vector<Vec3> positions(3);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(2, 0, 0);
+    positions[2] = Vec3(3, 0, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Forces | State::Energy);
+    const vector<Vec3>& forces = state.getForces();
+    cout << ""  << endl;
+    cout << "    force 0: " << forces[0] << endl;
+    cout << "    force 1: " << forces[1] << endl;
+    cout << "    force 2: " << forces[2] << endl;
+    cout << "energyPoten: " << state.getPotentialEnergy() << endl;
+    cout << ""  << endl;
+    cout << "    no cutoff force: 138.935485"  << endl;
+    cout << ""  << endl;
+    const double eps = 78.3;
+    const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0);
+    const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0);
+    const double force = 138.935485*(1.0)*(1.0-2.0*krf*1.0);
+    ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[1], TOL);
+    ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[2], TOL);
+    ASSERT_EQUAL_TOL(2*138.935485*(1.0)*(1.0+krf*1.0-crf), state.getPotentialEnergy(), TOL);
+}
+
+
+
+int main() {
+    try {
+    cout << ""  << endl;
+    cout << "Executing Periodic..." << endl;
+        testPeriodic();
+    cout << ""  << endl;
+    cout << "Executing Ewald..." << endl;
+        testEwald();
+    cout << ""  << endl;
+    cout << "Done" << endl;
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}

platforms/reference/src/ReferenceKernels.cpp

@@ -357,13 +357,20 @@ void ReferenceCalcNonbondedForceKernel::executeForces(OpenMMContextImpl& context
     RealOpenMM** forceData = ((ReferenceFloatStreamImpl&) context.getForces().getImpl()).getData();
     ReferenceLJCoulombIxn clj;
     bool periodic = (nonbondedMethod == CutoffPeriodic);
+    bool ewald  = (nonbondedMethod == Ewald);
     if (nonbondedMethod != NoCutoff) {
-        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, periodic ? periodicBoxSize : NULL, nonbondedCutoff, 0.0);
+        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, (periodic || ewald) ? periodicBoxSize : NULL, nonbondedCutoff, 0.0);
         clj.setUseCutoff(nonbondedCutoff, *neighborList, 78.3f);
     }
-    if (periodic)
+    if (periodic||ewald)
         clj.setPeriodic(periodicBoxSize);
-    clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, 0);
+    if (ewald) {
+        clj.setRecipVectors();
+        clj.calculateEwaldIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, 0);
+    }
+    else {
+        clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, 0);
+    }
     ReferenceBondForce refBondForce;
     ReferenceLJCoulomb14 nonbonded14;
     if (nonbondedMethod != NoCutoff)
@@ -377,13 +384,20 @@ double ReferenceCalcNonbondedForceKernel::executeEnergy(OpenMMContextImpl& conte
     RealOpenMM energy = 0;
     ReferenceLJCoulombIxn clj;
     bool periodic = (nonbondedMethod == CutoffPeriodic);
+    bool ewald  = (nonbondedMethod == Ewald);
     if (nonbondedMethod != NoCutoff) {
-        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, periodic ? periodicBoxSize : NULL, nonbondedCutoff, 0.0);
+        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, (periodic || ewald) ? periodicBoxSize : NULL, nonbondedCutoff, 0.0);
         clj.setUseCutoff(nonbondedCutoff, *neighborList, 78.3f);
     }
-    if (periodic)
+    if (periodic || ewald)
         clj.setPeriodic(periodicBoxSize);
-    clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, &energy);
+    if (ewald) {
+        clj.setRecipVectors();
+        clj.calculateEwaldIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, &energy);
+    }
+    else {
+        clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, &energy);
+    }
     ReferenceBondForce refBondForce;
     ReferenceLJCoulomb14 nonbonded14;
     if (nonbondedMethod != NoCutoff)

platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.cpp

@@ -24,6 +24,7 @@
 
 #include <string.h>
 #include <sstream>
+#include <complex>
 
 #include "../SimTKUtilities/SimTKOpenMMCommon.h"
 #include "../SimTKUtilities/SimTKOpenMMLog.h"
@@ -164,6 +165,295 @@ int ReferenceLJCoulombIxn::getDerivedParameters( RealOpenMM c6, RealOpenMM c12,
    return ReferenceForce::DefaultReturn;
 }
 
+  /**---------------------------------------------------------------------------------------
+
+     Set the reciprocal space vectors to use with Ewald
+
+// Currently a dumb routine, vectors are set in calculateEwaldIxn
+
+     @return ReferenceForce::DefaultReturn
+
+     --------------------------------------------------------------------------------------- */
+
+ int ReferenceLJCoulombIxn::setRecipVectors() {
+    
+    return ReferenceForce::DefaultReturn;
+  }
+
+  /**---------------------------------------------------------------------------------------
+
+     Calculate the Ewald parameter based on the cutoff and the desired tolerance using
+
+     erfc( alpha*cutoff )/cutoff < tolerance
+
+     @return ReferenceForce::DefaultReturn
+
+     --------------------------------------------------------------------------------------- */
+
+/* int ReferenceLJCoulombIxn::setAlphaEwald(RealOpenMM cutoff, RealOpenMM tolerance,
+                                          RealOpenMM alphaEwald) {
+ 
+    alphaEwald = 1.0f;
+
+    while ( erfc(alphaEwald*cutoff)  >=  tolerance*cutoff) {
+      alphaEwald= 1.2 * alphaEwald;
+    }
+
+    return ReferenceForce::DefaultReturn;
+  }
+*/
+
+/**---------------------------------------------------------------------------------------
+
+   Calculate Ewald 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::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** atomCoordinates,
+                                             RealOpenMM** atomParameters, int** exclusions,
+                                             RealOpenMM* fixedParameters, RealOpenMM** forces,
+                                             RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const {
+
+    #include "../SimTKUtilities/RealTypeSimTk.h"
+    typedef std::complex<RealOpenMM> d_complex;
+    typedef std::complex<int> int_complex;
+
+// Number of R-vectors (real space vectors)
+// to be calculated automatically eventually from alphaEwald and desired precision
+
+    int numRx = 60+1;
+    int numRy = 60+1;
+    int numRz = 60+1;
+    int kmax = std::max(numRx, std::max(numRy,numRz));
+
+    static const RealOpenMM alphaEwald       =  3.123413;
+    RealOpenMM  factorEwald = -1.0 / (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;
+
+    RealOpenMM recipCoeff = 4.0*M_PI/(periodicBoxSize[0] * periodicBoxSize[1] * periodicBoxSize[2]) /epsilon;
+
+    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
+// **************************************************************************************
+
+// setup reciprocal box
+
+           RealOpenMM recipBoxSize[3] = { TWO_PI / periodicBoxSize[0], TWO_PI / periodicBoxSize[1], TWO_PI / periodicBoxSize[2]};
+
+
+// setup K-vectors
+
+  d_complex eir[kmax][numberOfAtoms][3];
+  d_complex tab_xy[numberOfAtoms];
+  d_complex tab_qxyz[numberOfAtoms];
+  d_complex a,b,c;
+
+  int  i,j,m;
+
+  if (kmax < 1) {
+      std::stringstream message;
+      message << " kmax < 1 , Aborting" << std::endl;
+      SimTKOpenMMLog::printError( message );
+  }
+
+  for(i = 0; (i < numberOfAtoms); i++) {
+    for(m = 0; (m < 3); m++) {
+      eir[0][i][m].real() = 1;
+      eir[0][i][m].imag() = 0;
+    }
+
+
+    for(m=0; (m<3); m++) {
+      eir[1][i][m].real() = cos(atomCoordinates[i][m]*recipBoxSize[m]);
+      eir[1][i][m].imag() = sin(atomCoordinates[i][m]*recipBoxSize[m]);
+    }
+
+    for(j=2; (j<kmax); j++) {
+      for(m=0; (m<3); m++) {
+        eir[j][i][m] = eir[j-1][i][m] * eir[1][i][m];
+      }
+    }
+  }
+
+// calculate reciprocal space energy and forces
+
+    int lowry = 0;
+    int lowrz = 1;
+
+    for(int rx = 0; rx < numRx; rx++) {
+
+      RealOpenMM kx = rx * recipBoxSize[0];
+
+      for(int ry = lowry; ry < numRy; ry++) {
+
+        RealOpenMM ky = ry * recipBoxSize[1];
+
+        if(ry >= 0) {
+          for(int n = 0; n < numberOfAtoms; n++) {
+            tab_xy[n] = eir[rx][n][0] * eir[ry][n][1];
+          }
+        }
+
+        else {
+          for(int n = 0; n < numberOfAtoms; n++) {
+            tab_xy[n]= eir[rx][n][0] * conj (eir[-ry][n][1]);
+          }
+        }
+
+        for (int rz = lowrz; rz < numRz; rz++) {
+
+          if( rz >= 0) {
+
+           for( int n = 0; n < numberOfAtoms; n++) {
+             tab_qxyz[n].real() = atomParameters[n][QIndex] * (tab_xy[n] * eir[rz][n][2]).real();
+             tab_qxyz[n].imag() = atomParameters[n][QIndex] * (tab_xy[n] * eir[rz][n][2]).imag();
+           }
+          }
+
+          else {
+
+            for( int n = 0; n < numberOfAtoms; n++) {
+              tab_qxyz[n].real() = atomParameters[n][QIndex] * (tab_xy[n] * conj(eir[-rz][n][2])).real() ;
+              tab_qxyz[n].imag() = atomParameters[n][QIndex] * (tab_xy[n] * conj(eir[-rz][n][2])).imag() ;
+            }
+          }
+
+          RealOpenMM cs = 0.0f;
+          RealOpenMM ss = 0.0f;
+
+          for( int n = 0; n < numberOfAtoms; n++) {
+            cs += tab_qxyz[n].real();
+            ss += tab_qxyz[n].imag();
+          }
+
+
+          RealOpenMM kz = rz * recipBoxSize[2];
+          RealOpenMM k2 = kx * kx + ky * ky + kz * kz;
+          RealOpenMM ak = exp(k2*factorEwald) / k2;
+          RealOpenMM akv = 2.0 * ak * (1.0/k2 - factorEwald);
+
+          recipEnergy += ak * ( cs * cs + ss * ss);
+          RealOpenMM vir =  akv * ( 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.0 * recipCoeff * force * kx ;
+            forces[n][1] += 2.0 * recipCoeff * force * ky ;
+            forces[n][2] += 2.0 * recipCoeff * force * kz ;
+          }
+          lowrz = 1 - numRz;
+        }
+        lowry = 1 - numRy;
+      }
+    }
+
+  recipEnergy *= recipCoeff;
+
+
+
+// **************************************************************************************
+// 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]);
+
+          realSpaceEwaldEnergy = realSpaceEwaldEnergy + atomParameters[atomID1][QIndex]*atomParameters[atomID2][QIndex]*inverseR*erfc(alphaEwald*r); 
+        }
+       }
+
+       // allocate and initialize exclusion array
+
+       int* exclusionIndices = new int[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 r         = deltaR[0][ReferenceForce::RIndex];
+       RealOpenMM r2        = deltaR[0][ReferenceForce::R2Index];
+       RealOpenMM inverseR  = one/(deltaR[0][ReferenceForce::RIndex]);
+       RealOpenMM alphaR    = alphaEwald * r;
+ 
+       realSpaceEwaldEnergy = realSpaceEwaldEnergy + atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR*erfc(alphaR); 
+       RealOpenMM dEdR = atomParameters[ii][QIndex] * atomParameters[jj][QIndex] * inverseR * inverseR * inverseR;
+       dEdR = dEdR * (erfc(alphaR) + 2.0 * 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;
+                } 
+             }
+          }
+       }
+
+       delete[] exclusionIndices;
+
+// ***********************************************************************
+
+        if( totalEnergy ) {
+              *totalEnergy += recipEnergy + realSpaceEwaldEnergy - selfEwaldEnergy;
+        }
+
+
+   return ReferenceForce::DefaultReturn;
+}
+
+
 /**---------------------------------------------------------------------------------------
 
    Calculate LJ Coulomb pair ixn

platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.h

@@ -113,6 +113,19 @@ class ReferenceLJCoulombIxn : public ReferencePairIxn {
          --------------------------------------------------------------------------------------- */
       
       int setPeriodic( RealOpenMM* boxSize );      
+      
+      /**---------------------------------------------------------------------------------------
+      
+         Set the reciprocal vectors to use with Ewald
+      
+         @param
+      
+         @return ReferenceForce::DefaultReturn
+      
+         --------------------------------------------------------------------------------------- */
+      
+      int setRecipVectors();      
+
 
       /**---------------------------------------------------------------------------------------
       
@@ -159,6 +172,31 @@ class ReferenceLJCoulombIxn : public ReferencePairIxn {
                             RealOpenMM** atomParameters, int** exclusions,
                             RealOpenMM* fixedParameters, RealOpenMM** forces,
                             RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const;
+       
+      /**---------------------------------------------------------------------------------------
+      
+         Calculate Ewald ixn
+      
+         @param numberOfAtoms    number of atoms
+         @param atomCoordinates  atom coordinates
+         @param atomParameters   atom parameters (charges, c6, c12, ...)     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 calculateEwaldIxn( int numberOfAtoms, RealOpenMM** atomCoordinates,
+                            RealOpenMM** atomParameters, int** exclusions,
+                            RealOpenMM* fixedParameters, RealOpenMM** forces,
+                            RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const;
       
 };
 

platforms/reference/src/SimTKReference/ReferencePairIxn.h

@@ -72,7 +72,30 @@ class ReferencePairIxn {
                             RealOpenMM** atomParameters, int** exclusions,
                             RealOpenMM* fixedParameters, RealOpenMM** forces,
                             RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const = 0;
+ 
+      /**---------------------------------------------------------------------------------------
+      
+         Calculate Ewald ixn
+      
+         @param numberOfAtoms    number of atoms
+         @param atomCoordinates  atom coordinates
+         @param atomParameters   atom parameters (charges, c6, c12, ...)     atomParameters[atomIndex][paramterIndex]
+         @param exclusions       atom exclusion indices                      exclusions[atomIndex][atomToExcludeIndex]
+         @param fixedParameters  non-atom parameters
+         @param forces           force array (forces added)
+         @param energyByAtom     atom energy
+         @param totalEnergy      total energy
+      
+         @return ReferenceForce::DefaultReturn
+            
+         --------------------------------------------------------------------------------------- */
+          
+      virtual int calculateEwaldIxn( int numberOfAtoms, RealOpenMM** atomCoordinates,
+                            RealOpenMM** atomParameters, int** exclusions,
+                            RealOpenMM* fixedParameters, RealOpenMM** forces,
+                            RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const = 0;
       
+     
 };
 
 // ---------------------------------------------------------------------------------------

platforms/reference/tests/TestReferenceEwald.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) 2008 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.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the Eewald summation method reference implementation of NonbondedForce.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "OpenMMContext.h"
+#include "ReferencePlatform.h"
+#include "NonbondedForce.h"
+#include "System.h"
+#include "VerletIntegrator.h"
+#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
+#include "HarmonicBondForce.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+const double TOL = 1e-5;
+
+void testEwald() {
+    ReferencePlatform platform;
+    System system(2, 0);
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* nonbonded = new NonbondedForce(2, 0);
+    nonbonded->setParticleParameters(0, 1.0, 1, 0);
+    nonbonded->setParticleParameters(1, -1.0, 1, 0);
+    nonbonded->setNonbondedMethod(NonbondedForce::Ewald);
+    const double cutoff = 2.0;
+    nonbonded->setCutoffDistance(cutoff);
+    nonbonded->setPeriodicBoxVectors(Vec3(6, 0, 0), Vec3(0, 6, 0), Vec3(0, 0, 6));
+    system.addForce(nonbonded);
+    OpenMMContext context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(3.048000,2.764000,3.156000);
+    positions[1] = Vec3(2.809000,2.888000,2.571000);
+    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_VEC(Vec3(-123.711, 64.1877, -302.716), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(123.711, -64.1877, 302.716), forces[1], TOL);
+//    const double eps = 78.3;
+//    const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0);
+//    const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0);
+//    const double force = 138.935485*(1.0)*(1.0-2.0*krf*1.0);
+//    ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[0], TOL);
+//    ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[1], TOL);
+//    ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[2], TOL);
+//    ASSERT_EQUAL_TOL(2*138.935485*(1.0)*(1.0+krf*1.0-crf), state.getPotentialEnergy(), TOL);
+}
+
+/*
+void testEwald4() {
+    ReferencePlatform platform;
+    System system(4, 0);
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* nonbonded = new NonbondedForce(4, 0);
+    nonbonded->setParticleParameters(0, 1.0, 1, 0);
+    nonbonded->setParticleParameters(1, 1.0, 1, 0);
+    nonbonded->setParticleParameters(2, -1.0, 1, 0);
+    nonbonded->setParticleParameters(3, -1.0, 1, 0);
+    nonbonded->setNonbondedMethod(NonbondedForce::Ewald);
+    const double cutoff = 2.0;
+    nonbonded->setCutoffDistance(cutoff);
+    nonbonded->setPeriodicBoxVectors(Vec3(6, 0, 0), Vec3(0, 6, 0), Vec3(0, 0, 6));
+    system.addForce(nonbonded);
+    OpenMMContext context(system, integrator, platform);
+    vector<Vec3> positions(4);
+    positions[0] = Vec3(3.048000,2.764000,3.156000);
+    positions[1] = Vec3(3.348000,2.764000,3.156000);
+    positions[2] = Vec3(2.809000,2.888000,2.571000);
+    positions[3] = Vec3(2.509000,2.888000,2.571000);
+    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 << "energyPoten: " << state.getPotentialEnergy() << endl;
+ //   const double eps = 78.3;
+ //   const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0);
+ //   const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0);
+//    const double force = 138.935485*(1.0)*(1.0-2.0*krf*1.0);
+//    ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[0], TOL);
+//    ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[1], TOL);
+//    ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[2], TOL);
+//    ASSERT_EQUAL_TOL(2*138.935485*(1.0)*(1.0+krf*1.0-crf), state.getPotentialEnergy(), TOL);
+}
+*/
+
+/*
+void testPeriodic() {
+    ReferencePlatform platform;
+    System system(2, 0);
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* nonbonded = new NonbondedForce(2, 0);
+    nonbonded->setParticleParameters(0, 1.0, 1, 0);
+    nonbonded->setParticleParameters(1, -1.0, 1, 0);
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodci);
+    const double cutoff = 1.0;
+    nonbonded->setCutoffDistance(cutoff);
+    nonbonded->setPeriodicBoxVectors(Vec3(6, 0, 0), Vec3(0, 6, 0), Vec3(0, 0, 6));
+    system.addForce(nonbonded);
+    OpenMMContext context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(3.044293,2.765923,3.146914);
+    positions[1] = Vec3(2.812707,2.886077,2.580086);
+    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 << "energyPoten: " << state.getPotentialEnergy() << endl;
+}
+*/
+
+int main() {
+    try {
+//        testPeriodic();
+        testEwald();
+//        testEwald4();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}