Fixed bugs in CUDA implementation of AmoebaWcaDispersionForce

plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.cpp

@@ -2231,7 +2231,7 @@ void CudaCalcAmoebaWcaDispersionForceKernel::initialize(const System& system, co
     defines["SHCTD"] = cu.doubleToString(force.getShctd());
     CUmodule module = cu.createModule(CudaKernelSources::vectorOps+CudaAmoebaKernelSources::amoebaWcaForce, defines);
     forceKernel = cu.getKernel(module, "computeWCAForce");
-    totalMaximumDispersionEnergy =  AmoebaWcaDispersionForceImpl::getTotalMaximumDispersionEnergy(force);
+    totalMaximumDispersionEnergy = AmoebaWcaDispersionForceImpl::getTotalMaximumDispersionEnergy(force);
 
     // Add an interaction to the default nonbonded kernel.  This doesn't actually do any calculations.  It's
     // just so that CudaNonbondedUtilities will keep track of the tiles.
@@ -2269,5 +2269,6 @@ void CudaCalcAmoebaWcaDispersionForceKernel::copyParametersToContext(ContextImpl
         radiusEpsilonVec[i] = make_float2((float) radius, (float) epsilon);
     }
     radiusEpsilon->upload(radiusEpsilonVec);
+    totalMaximumDispersionEnergy = AmoebaWcaDispersionForceImpl::getTotalMaximumDispersionEnergy(force);
     cu.invalidateMolecules();
 }

plugins/amoeba/platforms/cuda/src/kernels/amoebaWcaForce.cu

@@ -226,22 +226,24 @@ extern "C" __global__ void computeWCAForce(unsigned long long* __restrict__ forc
 
             // Compute forces.
 
+            unsigned int tj = tgx;
             for (unsigned int j = 0; j < TILE_SIZE; j++) {
-                int atom2 = y*TILE_SIZE+j;
+                int atom2 = y*TILE_SIZE+tj;
                 if (atom1 != atom2 && atom1 < NUM_ATOMS && atom2 < NUM_ATOMS) {
                     real3 tempForce;
                     real tempEnergy;
-                    computeOneInteraction(data, localData[tbx+j], rmixo, rmixh, emixo, emixh, tempForce, tempEnergy);
+                    computeOneInteraction(data, localData[tbx+tj], rmixo, rmixh, emixo, emixh, tempForce, tempEnergy);
                     data.force += tempForce;
-                    localData[tbx+j].force -= tempForce;
+                    localData[tbx+tj].force -= tempForce;
                     energy += (x == y ? 0.5f*tempEnergy : tempEnergy);
                     real emjxo, emjxh, rmjxo, rmjxh;
-                    initParticleParameters(localData[tbx+j].radius, localData[tbx+j].epsilon, rmjxo, rmjxh, emjxo, emjxh);
-                    computeOneInteraction(localData[tbx+j], data, rmjxo, rmjxh, emjxo, emjxh, tempForce, tempEnergy);
+                    initParticleParameters(localData[tbx+tj].radius, localData[tbx+tj].epsilon, rmjxo, rmjxh, emjxo, emjxh);
+                    computeOneInteraction(localData[tbx+tj], data, rmjxo, rmjxh, emjxo, emjxh, tempForce, tempEnergy);
                     data.force -= tempForce;
-                    localData[tbx+j].force += tempForce;
+                    localData[tbx+tj].force += tempForce;
                     energy += (x == y ? 0.5f*tempEnergy : tempEnergy);
                 }
+                tj = (tj+1) & (TILE_SIZE-1);
             }
             unsigned int offset = x*TILE_SIZE + tgx;
             atomicAdd(&forceBuffers[offset], static_cast<unsigned long long>((long long) (data.force.x*0xFFFFFFFF)));