Fixed an error in the computation of Drude screening force

plugins/drude/platforms/cuda/src/kernels/drudePairForce.cu

@@ -13,7 +13,7 @@ real u = drudeParams.x*r;
 real screening = 1-(1+0.5f*u)*EXP(-u);
 real pairEnergy = drudeParams.y*screening*rInv;
 energy += pairEnergy;
-real3 f = delta*(pairEnergy*rInv*rInv);
+real3 f = delta*(drudeParams.y*rInv*rInv*(screening*rInv-0.5f*(1+u)*EXP(-u)*drudeParams.x));
 force1 += f;
 force3 -= f;
 
@@ -26,7 +26,7 @@ u = drudeParams.x*r;
 screening = 1-(1+0.5f*u)*EXP(-u);
 pairEnergy = -drudeParams.y*screening*rInv;
 energy += pairEnergy;
-f = delta*(pairEnergy*rInv*rInv);
+f = delta*(-drudeParams.y*rInv*rInv*(screening*rInv-0.5f*(1+u)*EXP(-u)*drudeParams.x));
 force1 += f;
 force4 -= f;
 
@@ -39,7 +39,7 @@ u = drudeParams.x*r;
 screening = 1-(1+0.5f*u)*EXP(-u);
 pairEnergy = -drudeParams.y*screening*rInv;
 energy += pairEnergy;
-f = delta*(pairEnergy*rInv*rInv);
+f = delta*(-drudeParams.y*rInv*rInv*(screening*rInv-0.5f*(1+u)*EXP(-u)*drudeParams.x));
 force2 += f;
 force3 -= f;
 
@@ -52,6 +52,6 @@ u = drudeParams.x*r;
 screening = 1-(1+0.5f*u)*EXP(-u);
 pairEnergy = drudeParams.y*screening*rInv;
 energy += pairEnergy;
-f = delta*(pairEnergy*rInv*rInv);
+f = delta*(drudeParams.y*rInv*rInv*(screening*rInv-0.5f*(1+u)*EXP(-u)*drudeParams.x));
 force2 += f;
 force4 -= f;

plugins/drude/platforms/opencl/src/kernels/drudePairForce.cl

@@ -13,7 +13,7 @@ real u = drudeParams.x*r;
 real screening = 1-(1+0.5f*u)*EXP(-u);
 real pairEnergy = drudeParams.y*screening*rInv;
 energy += pairEnergy;
-real4 f = delta*(pairEnergy*rInv*rInv);
+real4 f = delta*(drudeParams.y*rInv*rInv)*(screening*rInv-0.5f*(1+u)*EXP(-u)*drudeParams.x);
 force1 += f;
 force3 -= f;
 
@@ -26,7 +26,7 @@ u = drudeParams.x*r;
 screening = 1-(1+0.5f*u)*EXP(-u);
 pairEnergy = -drudeParams.y*screening*rInv;
 energy += pairEnergy;
-f = delta*(pairEnergy*rInv*rInv);
+f = delta*(-drudeParams.y*rInv*rInv)*(screening*rInv-0.5f*(1+u)*EXP(-u)*drudeParams.x);
 force1 += f;
 force4 -= f;
 
@@ -39,7 +39,7 @@ u = drudeParams.x*r;
 screening = 1-(1+0.5f*u)*EXP(-u);
 pairEnergy = -drudeParams.y*screening*rInv;
 energy += pairEnergy;
-f = delta*(pairEnergy*rInv*rInv);
+f = delta*(-drudeParams.y*rInv*rInv)*(screening*rInv-0.5f*(1+u)*EXP(-u)*drudeParams.x);
 force2 += f;
 force3 -= f;
 
@@ -52,6 +52,6 @@ u = drudeParams.x*r;
 screening = 1-(1+0.5f*u)*EXP(-u);
 pairEnergy = drudeParams.y*screening*rInv;
 energy += pairEnergy;
-f = delta*(pairEnergy*rInv*rInv);
+f = delta*(drudeParams.y*rInv*rInv)*(screening*rInv-0.5f*(1+u)*EXP(-u)*drudeParams.x);
 force2 += f;
 force4 -= f;

plugins/drude/platforms/reference/src/ReferenceDrudeKernels.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2011-2013 Stanford University and the Authors.      *
+ * Portions copyright (c) 2011-2014 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -188,6 +188,7 @@ double ReferenceCalcDrudeForceKernel::execute(ContextImpl& context, bool include
         int dipole2 = pair2[i];
         int dipole1Particles[] = {particle[dipole1], particle1[dipole1]};
         int dipole2Particles[] = {particle[dipole2], particle1[dipole2]};
+        RealOpenMM uscale = pairThole[i]/pow(polarizability[dipole1]*polarizability[dipole2], 1.0/6.0);
         for (int j = 0; j < 2; j++)
             for (int k = 0; k < 2; k++) {
                 int p1 = dipole1Particles[j];
@@ -195,10 +196,10 @@ double ReferenceCalcDrudeForceKernel::execute(ContextImpl& context, bool include
                 RealOpenMM chargeProduct = charge[dipole1]*charge[dipole2]*(j == k ? 1 : -1);
                 RealVec delta = pos[p1]-pos[p2];
                 RealOpenMM r = sqrt(delta.dot(delta));
-                RealOpenMM u = r*pairThole[i]/pow(polarizability[dipole1]*polarizability[dipole2], 1.0/6.0);
+                RealOpenMM u = r*uscale;
                 RealOpenMM screening = 1.0 - (1.0+0.5*u)*exp(-u);
                 energy += ONE_4PI_EPS0*chargeProduct*screening/r;
-                RealVec f = delta*(ONE_4PI_EPS0*chargeProduct*screening/(r*r*r));
+                RealVec f = delta*(ONE_4PI_EPS0*chargeProduct/(r*r))*(screening/r-0.5*(1+u)*exp(-u)*uscale);
                 force[p1] += f;
                 force[p2] -= f;
             }
@@ -461,4 +462,4 @@ void ReferenceIntegrateDrudeSCFStepKernel::minimize(ContextImpl& context, double
     lbfgsfloatval_t fx;
     MinimizerData data(context, drudeParticles);
     lbfgs(numDrudeParticles*3, minimizerPos, &fx, evaluate, NULL, &data, &minimizerParams);
-}
\ No newline at end of file
+}

plugins/drude/platforms/reference/tests/TestReferenceDrudeForce.cpp

@@ -71,7 +71,7 @@ void validateForce(System& system, vector<Vec3>& positions, double expectedEnerg
             offsetPos[i][j] = positions[i][j]+offset;
             context.setPositions(offsetPos);
             double e2 = context.getState(State::Energy | State::Forces).getPotentialEnergy();
-            ASSERT_EQUAL_TOL(state.getForces()[i][j], (e1-e2)/(2*offset), 1e-3);
+            ASSERT_EQUAL_TOL(state.getForces()[i][j], (e1-e2)/(2*offset), 1e-5);
         }
 }