Cleaned up Cuda-Ewald

platforms/cuda/src/CudaKernels.cpp

@@ -291,6 +291,7 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
             Vec3 boxVectors[3];
             force.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
             gpuSetPeriodicBoxSize(gpu, (float)boxVectors[0][0], (float)boxVectors[1][1], (float)boxVectors[2][2]);
+            gpuSetEwaldParameters(gpu);//, (float)alphaEwald, (int)kmax);
             method = EWALD;
 
         }

platforms/cuda/src/kernels/cudatypes.h

@@ -292,6 +292,13 @@ struct cudaGmxSimulation {
     float           periodicBoxSizeX;               // The X dimension of the periodic box
     float           periodicBoxSizeY;               // The Y dimension of the periodic box
     float           periodicBoxSizeZ;               // The Z dimension of the periodic box
+    float           recipBoxSizeX;                  // The X dimension of the reciprocal box for Ewald summation
+    float           recipBoxSizeY;                  // The Y dimension of the reciprocal box for Ewald summation
+    float           recipBoxSizeZ;                  // The Z dimension of the reciprocal box for Ewald summation
+    float           cellVolume;                     // Ewald parameter alpha (a.k.a. kappa)
+    float           alphaEwald;                     // Ewald parameter alpha (a.k.a. kappa)
+    float           factorEwald;                    // - 1 ( 4 * alphaEwald * alphaEwald)
+    float           kmax;                           // Maximum number of reciprocal vectors
     float           reactionFieldK;                 // Constant for reaction field correction
     float           probeRadius;                    // SASA probe radius
     float           surfaceAreaFactor;              // ACE approximation surface area factor

platforms/cuda/src/kernels/gpu.cpp

@@ -413,6 +413,26 @@ void gpuSetNonbondedCutoff(gpuContext gpu, float cutoffDistance, float solventDi
     gpu->sim.reactionFieldK = pow(cutoffDistance, -3.0f)*(solventDielectric-1.0f)/(2.0f*solventDielectric+1.0f);
 }
 
+extern "C"
+void gpuSetEwaldParameters(gpuContext gpu)//, float alphaEwald, int kmax )
+{
+
+    // hard coded alphaEwald and kmax, no interface yet
+    float alpha            = 3.123413;
+    float PI               = 3.14159265358979323846f;
+    float TWO_PI           = 2.0 * PI;
+
+    gpu->sim.recipBoxSizeX = TWO_PI / gpu->sim.periodicBoxSizeX ;
+    gpu->sim.recipBoxSizeY = TWO_PI / gpu->sim.periodicBoxSizeY ;
+    gpu->sim.recipBoxSizeZ = TWO_PI / gpu->sim.periodicBoxSizeZ ;
+
+    gpu->sim.cellVolume        = gpu->sim.periodicBoxSizeX * gpu->sim.periodicBoxSizeY * gpu->sim.periodicBoxSizeZ;
+
+    gpu->sim.alphaEwald        = alpha;
+    gpu->sim.factorEwald       = -1 / (4*alpha*alpha);
+    gpu->sim.kmax              = 20+1;
+}
+
 extern "C"
 void gpuSetPeriodicBoxSize(gpuContext gpu, float xsize, float ysize, float zsize)
 {

platforms/cuda/src/kernels/gputypes.h

@@ -176,6 +176,9 @@ void gpuSetCoulombParameters(gpuContext gpu, float epsfac, const std::vector<int
 extern "C"
 void gpuSetNonbondedCutoff(gpuContext gpu, float cutoffDistance, float solventDielectric);
 
+extern "C"
+void gpuSetEwaldParameters(gpuContext gpu);//, float alphaEwald, int kmax);
+
 extern "C"
 void gpuSetPeriodicBoxSize(gpuContext gpu, float xsize, float ysize, float zsize);
 

platforms/cuda/src/kernels/kCalculateCDLJEwaldReciprocal.h

@@ -32,10 +32,7 @@
 
 __global__ void kCalculateCDLJEwaldReciprocalForces_kernel()
 {
-    float alphaEwald       =  3.123413;
-    float PI               = 3.14159265358979323846f;
-    float TWO_PI           = 2.0 * PI;
-    float SQRT_PI          = sqrt(PI);
+    float alphaEwald       = cSim.alphaEwald;
     float eps0             = 5.72765E-4;
 
     int numRx              = 20+1;
@@ -47,18 +44,13 @@ __global__ void kCalculateCDLJEwaldReciprocalForces_kernel()
     float kx, ky, kz, k2, ek;
     float Qi, Qj, SinI, SinJ, CosI, CosJ;
 
-    float factorEwald   = -1 / (4*alphaEwald*alphaEwald);
-
-    float recipBoxSizeX = TWO_PI / cSim.periodicBoxSizeX;
-    float recipBoxSizeY = TWO_PI / cSim.periodicBoxSizeY;
-    float recipBoxSizeZ = TWO_PI / cSim.periodicBoxSizeZ;
-    float V = cSim.periodicBoxSizeX * cSim.periodicBoxSizeY * cSim.periodicBoxSizeZ;
+    float V = cSim.cellVolume;
 
     float4 apos1, apos2 ;
 
     unsigned int atomID1    = threadIdx.x + blockIdx.x * blockDim.x;
 
-    while (atomID1 < cSim.atoms)
+    while (atomID1 < cSim.stride * cSim.outputBuffers)
     {
         apos1             = cSim.pPosq[atomID1];
 
@@ -72,18 +64,18 @@ __global__ void kCalculateCDLJEwaldReciprocalForces_kernel()
 
                 for(int rx = 0; rx < numRx; rx++) {
 
-                  kx = rx * recipBoxSizeX;
+                  kx = rx * cSim.recipBoxSizeX;
 
                   for(int ry = lowry; ry < numRy; ry++) {
 
-                    ky = ry * recipBoxSizeY;
+                    ky = ry * cSim.recipBoxSizeY;
 
                     for (int rz = lowrz; rz < numRz; rz++) {
 
-                      kz = rz * recipBoxSizeZ;
+                      kz = rz * cSim.recipBoxSizeZ;
 
                       k2 = kx * kx + ky * ky + kz * kz;
-                      ek = exp (  k2 * factorEwald);
+                      ek = exp (  k2 * cSim.factorEwald);
 
                       Qi = apos1.w ;
                       Qj = apos2.w ;