Merge pull request #1174 from peastman/amoebabug

Fixed a bug in AmoebaMultipoleForce

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

@@ -1040,9 +1040,9 @@ extern "C" __global__ void computeInducedDipoleForceAndEnergy(real4* __restrict_
         inducedDipolePolar[2] = cinducedDipolePolar[0]*fracToCart[0][2] + cinducedDipolePolar[1]*fracToCart[1][2] + cinducedDipolePolar[2]*fracToCart[2][2];
         real4 f = make_real4(0, 0, 0, 0);
 
-        energy += inducedDipole[0]*phi[i+NUM_ATOMS];
-        energy += inducedDipole[1]*phi[i+NUM_ATOMS*2];
-        energy += inducedDipole[2]*phi[i+NUM_ATOMS*3];
+        energy += (inducedDipole[0]+inducedDipolePolar[0])*phi[i+NUM_ATOMS];
+        energy += (inducedDipole[1]+inducedDipolePolar[1])*phi[i+NUM_ATOMS*2];
+        energy += (inducedDipole[2]+inducedDipolePolar[2])*phi[i+NUM_ATOMS*3];
 
         for (int k = 0; k < 3; k++) {
             int j1 = deriv1[k+1];
@@ -1070,7 +1070,7 @@ extern "C" __global__ void computeInducedDipoleForceAndEnergy(real4* __restrict_
         forceBuffers[i+PADDED_NUM_ATOMS] -= static_cast<unsigned long long>((long long) (f.y*0x100000000));
         forceBuffers[i+PADDED_NUM_ATOMS*2] -= static_cast<unsigned long long>((long long) (f.z*0x100000000));
     }
-    energyBuffer[blockIdx.x*blockDim.x+threadIdx.x] += 0.5f*EPSILON_FACTOR*energy;
+    energyBuffer[blockIdx.x*blockDim.x+threadIdx.x] += 0.25f*EPSILON_FACTOR*energy;
 }
 
 extern "C" __global__ void recordInducedFieldDipoles(const real* __restrict__ phid, real* const __restrict__ phip,

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

@@ -414,7 +414,7 @@ __device__ void computeOneInteraction(AtomData& atom1, AtomData& atom2, bool has
 __device__ void computeSelfEnergyAndTorque(AtomData& atom1, mixed& energy) {
     real cii = atom1.q*atom1.q;
     real3 dipole = make_real3(atom1.sphericalDipole.y, atom1.sphericalDipole.z, atom1.sphericalDipole.x);
-    real dii = dot(dipole, dipole+atom1.inducedDipole);
+    real dii = dot(dipole, dipole+(atom1.inducedDipole+atom1.inducedDipolePolar)*0.5);
 #ifdef INCLUDE_QUADRUPOLES
     real qii = (atom1.sphericalQuadrupole[0]*atom1.sphericalQuadrupole[0] +
                 atom1.sphericalQuadrupole[1]*atom1.sphericalQuadrupole[1] +

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceMultipoleForce.cpp

@@ -5805,9 +5805,9 @@ RealOpenMM AmoebaReferencePmeMultipoleForce::computeReciprocalSpaceInducedDipole
         inducedDipolePolar[1] = _inducedDipolePolar[i][0]*cartToFrac[1][0] + _inducedDipolePolar[i][1]*cartToFrac[1][1] + _inducedDipolePolar[i][2]*cartToFrac[1][2];
         inducedDipolePolar[2] = _inducedDipolePolar[i][0]*cartToFrac[2][0] + _inducedDipolePolar[i][1]*cartToFrac[2][1] + _inducedDipolePolar[i][2]*cartToFrac[2][2];
 
-        energy += inducedDipole[0]*_phi[20*i+1];
-        energy += inducedDipole[1]*_phi[20*i+2];
-        energy += inducedDipole[2]*_phi[20*i+3];
+        energy += (inducedDipole[0]+inducedDipolePolar[0])*_phi[20*i+1];
+        energy += (inducedDipole[1]+inducedDipolePolar[1])*_phi[20*i+2];
+        energy += (inducedDipole[2]+inducedDipolePolar[2])*_phi[20*i+3];
 
         RealVec f        = RealVec(0.0, 0.0, 0.0);
 
@@ -5841,7 +5841,7 @@ RealOpenMM AmoebaReferencePmeMultipoleForce::computeReciprocalSpaceInducedDipole
                                   f[0]*fracToCart[1][0] + f[1]*fracToCart[1][1] + f[2]*fracToCart[1][2],
                                   f[0]*fracToCart[2][0] + f[1]*fracToCart[2][1] + f[2]*fracToCart[2][2]);
     }
-    return (0.5*_electric*energy);
+    return (0.25*_electric*energy);
 }
 
 void AmoebaReferencePmeMultipoleForce::recordFixedMultipoleField()
@@ -6028,7 +6028,7 @@ RealOpenMM AmoebaReferencePmeMultipoleForce::calculatePmeSelfEnergy(const vector
         cii += particleI.charge*particleI.charge;
 
         RealVec dipole(particleI.sphericalDipole[1], particleI.sphericalDipole[2], particleI.sphericalDipole[0]);
-        dii += dipole.dot(dipole + _inducedDipole[ii]);
+        dii += dipole.dot(dipole + (_inducedDipole[ii]+_inducedDipolePolar[ii])*0.5);
 
         qii += (particleI.sphericalQuadrupole[0]*particleI.sphericalQuadrupole[0]
                +particleI.sphericalQuadrupole[1]*particleI.sphericalQuadrupole[1]