Fixed bugs in GK force

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

@@ -1365,31 +1365,19 @@ double CudaCalcAmoebaMultipoleForceKernel::execute(ContextImpl& context, bool in
             else {
                 cu.clearBuffer(*gkKernel->getInducedField());
                 cu.clearBuffer(*gkKernel->getInducedFieldPolar());
-                vector<float> d, dp;
-                gkKernel->getInducedDipoles()->download(d);
-                gkKernel->getInducedDipolesPolar()->download(dp);
-                printf("dipole\n");
-                for (int i = 0; i < cu.getNumAtoms(); i++)
-                    printf("%d %g %g %g, %g %g %g\n", i, d[3*i], d[3*i+1], d[3*i+2], dp[3*i], dp[3*i+1], dp[3*i+2]);
                 void* computeInducedFieldArgs[] = {&inducedField->getDevicePointer(), &inducedFieldPolar->getDevicePointer(), &cu.getPosq().getDevicePointer(),
                     &inducedDipole->getDevicePointer(), &inducedDipolePolar->getDevicePointer(), &startTileIndex, &numTileIndices,
                     &gkKernel->getInducedField()->getDevicePointer(), &gkKernel->getInducedFieldPolar()->getDevicePointer(),
                     &gkKernel->getInducedDipoles()->getDevicePointer(), &gkKernel->getInducedDipolesPolar()->getDevicePointer(),
                     &gkKernel->getBornRadii()->getDevicePointer(), &dampingAndThole->getDevicePointer()};
                 cu.executeKernel(computeInducedFieldKernel, computeInducedFieldArgs, numForceThreadBlocks*forceThreadBlockSize, forceThreadBlockSize);
-                vector<long long> f, fp;
-                gkKernel->getInducedField()->download(f);
-                gkKernel->getInducedFieldPolar()->download(fp);
-                printf("field\n");
-                for (int i = 0; i < cu.getNumAtoms(); i++)
-                    printf("%d %g %g %g, %g %g %g\n", i, f[i]/(double) 0xFFFFFFFF, f[i+cu.getPaddedNumAtoms()]/(double) 0xFFFFFFFF, f[i+2*cu.getPaddedNumAtoms()]/(double) 0xFFFFFFFF, fp[i]/(double) 0xFFFFFFFF, fp[i+cu.getPaddedNumAtoms()]/(double) 0xFFFFFFFF, fp[i+2*cu.getPaddedNumAtoms()]/(double) 0xFFFFFFFF);
                 void* updateInducedGkFieldArgs[] = {&field->getDevicePointer(), &fieldPolar->getDevicePointer(),
-                    &gkKernel->getField()->getDevicePointer(), &gkKernel->getFieldPolar()->getDevicePointer(), &gkKernel->getInducedField()->getDevicePointer(),
-                    &gkKernel->getInducedFieldPolar()->getDevicePointer(), &gkKernel->getInducedDipoles()->getDevicePointer(), &gkKernel->getInducedDipolesPolar()->getDevicePointer(),
-                    &polarizability->getDevicePointer(), &inducedDipoleErrors->getDevicePointer()};
+                    &gkKernel->getField()->getDevicePointer(), &gkKernel->getInducedField()->getDevicePointer(),
+                    &gkKernel->getInducedFieldPolar()->getDevicePointer(), &gkKernel->getInducedDipoles()->getDevicePointer(),
+                    &gkKernel->getInducedDipolesPolar()->getDevicePointer(), &polarizability->getDevicePointer(), &inducedDipoleErrors->getDevicePointer()};
                 cu.executeKernel(updateInducedFieldKernel, updateInducedGkFieldArgs, cu.getNumThreadBlocks()*cu.ThreadBlockSize, cu.ThreadBlockSize, cu.ThreadBlockSize*elementSize*2);
             }
-            void* updateInducedFieldArgs[] = {&field->getDevicePointer(), &fieldPolar->getDevicePointer(), &npt, &npt, &inducedField->getDevicePointer(),
+            void* updateInducedFieldArgs[] = {&field->getDevicePointer(), &fieldPolar->getDevicePointer(), &npt, &inducedField->getDevicePointer(),
                 &inducedFieldPolar->getDevicePointer(), &inducedDipole->getDevicePointer(), &inducedDipolePolar->getDevicePointer(),
                 &polarizability->getDevicePointer(), &inducedDipoleErrors->getDevicePointer()};
             cu.executeKernel(updateInducedFieldKernel, updateInducedFieldArgs, cu.getNumThreadBlocks()*cu.ThreadBlockSize, cu.ThreadBlockSize, cu.ThreadBlockSize*elementSize*2);
@@ -1517,7 +1505,7 @@ double CudaCalcAmoebaMultipoleForceKernel::execute(ContextImpl& context, bool in
             void* pmeRecordInducedFieldDipolesArgs[] = {&pmePhid->getDevicePointer(), &pmePhip->getDevicePointer(),
                 &inducedField->getDevicePointer(), &inducedFieldPolar->getDevicePointer(), cu.getInvPeriodicBoxSizePointer()};
             cu.executeKernel(pmeRecordInducedFieldDipolesKernel, pmeRecordInducedFieldDipolesArgs, cu.getNumAtoms());
-            void* updateInducedFieldArgs[] = {&field->getDevicePointer(), &fieldPolar->getDevicePointer(), &npt, &npt, &inducedField->getDevicePointer(),
+            void* updateInducedFieldArgs[] = {&field->getDevicePointer(), &fieldPolar->getDevicePointer(), &npt, &inducedField->getDevicePointer(),
                 &inducedFieldPolar->getDevicePointer(), &inducedDipole->getDevicePointer(), &inducedDipolePolar->getDevicePointer(),
                 &polarizability->getDevicePointer(), &inducedDipoleErrors->getDevicePointer()};
             cu.executeKernel(updateInducedFieldKernel, updateInducedFieldArgs, cu.getNumThreadBlocks()*cu.ThreadBlockSize, cu.ThreadBlockSize, cu.ThreadBlockSize*elementSize*2);
@@ -1734,7 +1722,7 @@ private:
 };
 
 CudaCalcAmoebaGeneralizedKirkwoodForceKernel::CudaCalcAmoebaGeneralizedKirkwoodForceKernel(std::string name, const Platform& platform, CudaContext& cu, System& system) : 
-           CalcAmoebaGeneralizedKirkwoodForceKernel(name, platform), cu(cu), system(system), params(NULL), bornRadii(NULL), field(NULL), fieldPolar(NULL),
+           CalcAmoebaGeneralizedKirkwoodForceKernel(name, platform), cu(cu), system(system), params(NULL), bornRadii(NULL), field(NULL),
            inducedField(NULL), inducedFieldPolar(NULL), inducedDipoleS(NULL), inducedDipolePolarS(NULL), bornSum(NULL), bornForce(NULL) {
 }
 
@@ -1746,8 +1734,6 @@ CudaCalcAmoebaGeneralizedKirkwoodForceKernel::~CudaCalcAmoebaGeneralizedKirkwood
         delete bornRadii;
     if (field != NULL)
         delete field;
-    if (fieldPolar != NULL)
-        delete fieldPolar;
     if (inducedField != NULL)
         delete inducedField;
     if (inducedFieldPolar != NULL)
@@ -1777,7 +1763,6 @@ void CudaCalcAmoebaGeneralizedKirkwoodForceKernel::initialize(const System& syst
     params = CudaArray::create<float2>(cu, paddedNumAtoms, "amoebaGkParams");
     bornRadii = new CudaArray(cu, paddedNumAtoms, elementSize, "bornRadii");
     field = new CudaArray(cu, 3*paddedNumAtoms, sizeof(long long), "gkField");
-    fieldPolar = new CudaArray(cu, 3*paddedNumAtoms, sizeof(long long), "gkFieldPolar");
     bornSum = CudaArray::create<long long>(cu, paddedNumAtoms, "bornSum");
     bornForce = CudaArray::create<long long>(cu, paddedNumAtoms, "bornForce");
     inducedDipoleS = new CudaArray(cu, 3*paddedNumAtoms, elementSize, "inducedDipoleS");
@@ -1787,7 +1772,6 @@ void CudaCalcAmoebaGeneralizedKirkwoodForceKernel::initialize(const System& syst
         inducedFieldPolar = new CudaArray(cu, 3*paddedNumAtoms, sizeof(long long), "gkInducedFieldPolar");
     }
     cu.addAutoclearBuffer(*field);
-    cu.addAutoclearBuffer(*fieldPolar);
     cu.addAutoclearBuffer(*bornSum);
     cu.addAutoclearBuffer(*bornForce);
     vector<float2> paramsVector(paddedNumAtoms);

plugins/amoeba/platforms/cuda2/src/AmoebaCudaKernels.h

@@ -470,9 +470,6 @@ public:
     CudaArray* getField() {
         return field;
     }
-    CudaArray* getFieldPolar() {
-        return fieldPolar;
-    }
     CudaArray* getInducedField() {
         return inducedField;
     }
@@ -494,7 +491,6 @@ private:
     CudaArray* bornRadii;
     CudaArray* bornForce;
     CudaArray* field;
-    CudaArray* fieldPolar;
     CudaArray* inducedField;
     CudaArray* inducedFieldPolar;
     CudaArray* inducedDipoleS;

plugins/amoeba/platforms/cuda2/src/kernels/multipoleInducedField.cu

@@ -149,32 +149,16 @@ __device__ void computeOneInteraction(AtomData& atom1, AtomData& atom2, real3 de
         atom2.fieldPolarS += rr3*atom1.inducedDipolePolarS + dDotDelta*deltaR;
     }
 
-
-
-
-
-
     real rb2 = atom1.bornRadius*atom2.bornRadius;
-
     real expterm = EXP(-r2/(GK_C*rb2));
     real expc = expterm/GK_C; 
-
     real gf2 = RECIP(r2+rb2*expterm);
     real gf = SQRT(gf2);
     real gf3 = gf2*gf;
     real gf5 = gf3*gf2;
-
-    // reaction potential auxiliary terms
- 
     real a10 = -gf3;
-
-    // reaction potential gradient auxiliary terms
-
     real expc1 = 1 - expc;
     real a11 = expc1 * 3 * gf5;
- 
-    // unweighted dipole reaction potential gradient tensor
-
     real3 gux = GK_FD*make_real3(a10+deltaR.x*deltaR.x*a11, deltaR.x*deltaR.y*a11, deltaR.x*deltaR.z*a11);
     real3 guy = make_real3(gux.y, GK_FD*(a10+deltaR.y*deltaR.y*a11), GK_FD*deltaR.y*deltaR.z*a11);
     real3 guz = make_real3(gux.z, guy.z, GK_FD*(a10+deltaR.z*deltaR.z*a11));
@@ -280,6 +264,11 @@ extern "C" __global__ void computeInducedField(
                 localData[threadIdx.x].inducedDipolePolar = data.inducedDipolePolar;
                 localData[threadIdx.x].thole = data.thole;
                 localData[threadIdx.x].damp = data.damp;
+#ifdef USE_GK
+                localData[threadIdx.x].inducedDipoleS = data.inducedDipoleS;
+                localData[threadIdx.x].inducedDipolePolarS = data.inducedDipolePolarS;
+                localData[threadIdx.x].bornRadius = data.bornRadius;
+#endif
                 for (unsigned int j = 0; j < TILE_SIZE; j++) {
                     real3 delta = localData[tbx+j].pos-data.pos;
 #ifdef USE_PERIODIC
@@ -435,9 +424,8 @@ extern "C" __global__ void computeInducedField(
 }
 
 extern "C" __global__ void updateInducedFieldBySOR(const long long* __restrict__ fixedField, const long long* __restrict__ fixedFieldPolar,
-        const long long* __restrict__ fixedFieldS, const long long* __restrict__ fixedFieldPolarS,
-        const long long* __restrict__ inducedField, const long long* __restrict__ inducedFieldPolar, real* __restrict__ inducedDipole,
-        real* __restrict__ inducedDipolePolar, const float* __restrict__ polarizability, float2* __restrict__ errors) {
+        const long long* __restrict__ fixedFieldS, const long long* __restrict__ inducedField, const long long* __restrict__ inducedFieldPolar,
+        real* __restrict__ inducedDipole, real* __restrict__ inducedDipolePolar, const float* __restrict__ polarizability, float2* __restrict__ errors) {
     extern __shared__ real2 buffer[];
     const float polarSOR = 0.55f;
 #ifdef USE_EWALD
@@ -455,10 +443,9 @@ extern "C" __global__ void updateInducedFieldBySOR(const long long* __restrict__
             int fieldIndex = atom+component*PADDED_NUM_ATOMS;
             real previousDipole = inducedDipole[dipoleIndex];
             real previousDipolePolar = inducedDipolePolar[dipoleIndex];
-            long long fixed = fixedField[fieldIndex] + (fixedFieldS == NULL ? (long long) 0 : fixedFieldS[fieldIndex]);
-            long long fixedPolar = fixedFieldPolar[fieldIndex] + (fixedFieldPolarS == NULL ? (long long) 0 : fixedFieldPolarS[fieldIndex]);
-            real newDipole = scale*((fixed+inducedField[fieldIndex])*fieldScale+ewaldScale*previousDipole);
-            real newDipolePolar = scale*((fixedPolar+inducedFieldPolar[fieldIndex])*fieldScale+ewaldScale*previousDipolePolar);
+            long long fixedS = (fixedFieldS == NULL ? (long long) 0 : fixedFieldS[fieldIndex]);
+            real newDipole = scale*((fixedField[fieldIndex]+fixedS+inducedField[fieldIndex])*fieldScale+ewaldScale*previousDipole);
+            real newDipolePolar = scale*((fixedFieldPolar[fieldIndex]+fixedS+inducedFieldPolar[fieldIndex])*fieldScale+ewaldScale*previousDipolePolar);
             newDipole = previousDipole + polarSOR*(newDipole-previousDipole);
             newDipolePolar = previousDipolePolar + polarSOR*(newDipolePolar-previousDipolePolar);
             inducedDipole[dipoleIndex] = newDipole;