Skip to content

GitLab

Commit 11dc1eb6 authored by Peter Eastman's avatar Peter Eastman
Browse files

Fixed error with Drude particles

parent da28ce0e
Hide whitespace changes
Inline Side-by-side
Showing with 12 additions and 10 deletions
platforms/cuda/src/kernels/pmeExclusions.cu
View file @ 11dc1eb6
const float4 exclusionParams = PARAMS[index]; const float4 exclusionParams = PARAMS[index];
real3 delta = make_real3(pos2.x-pos1.x, pos2.y-pos1.y, pos2.z-pos1.z); real3 delta = make_real3(pos2.x-pos1.x, pos2.y-pos1.y, pos2.z-pos1.z);
const real r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z; const real r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
const real invR = RSQRT(r2); const real r = SQRT(r2);
const real r = r2*invR; const real invR = RECIP(r);
const real alphaR = EWALD_ALPHA*r; const real alphaR = EWALD_ALPHA*r;
const real expAlphaRSqr = EXP(-alphaR*alphaR); const real expAlphaRSqr = EXP(-alphaR*alphaR);
const real erfAlphaR = ERF(alphaR);
real tempForce = 0.0f; real tempForce = 0.0f;
if (erfAlphaR > 1e-6f) { if (alphaR > 1e-6f) {
const real erfAlphaR = ERF(alphaR);
const real prefactor = exclusionParams.x*invR; const real prefactor = exclusionParams.x*invR;
tempForce = -prefactor*(erfAlphaR-alphaR*expAlphaRSqr*TWO_OVER_SQRT_PI); tempForce = -prefactor*(erfAlphaR-alphaR*expAlphaRSqr*TWO_OVER_SQRT_PI);
energy -= prefactor*erfAlphaR; energy -= prefactor*erfAlphaR;
...@@ -29,6 +29,7 @@ const real dprefac = eprefac + dar6/6.0f; ...@@ -29,6 +29,7 @@ const real dprefac = eprefac + dar6/6.0f;
energy += coef*(1.0f - expDar2*eprefac); energy += coef*(1.0f - expDar2*eprefac);
tempForce += 6.0f*coef*(1.0f - expDar2*dprefac); tempForce += 6.0f*coef*(1.0f - expDar2*dprefac);
#endif #endif
delta *= tempForce*invR*invR; if (r > 0)
delta *= tempForce*invR*invR;
real3 force1 = -delta; real3 force1 = -delta;
real3 force2 = delta; real3 force2 = delta;
platforms/opencl/src/kernels/pmeExclusions.cl
View file @ 11dc1eb6
const float4 exclusionParams = PARAMS[index]; const float4 exclusionParams = PARAMS[index];
real3 delta = (real3) (pos2.x-pos1.x, pos2.y-pos1.y, pos2.z-pos1.z); real3 delta = (real3) (pos2.x-pos1.x, pos2.y-pos1.y, pos2.z-pos1.z);
const real r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z; const real r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
const real invR = RSQRT(r2); const real r = SQRT(r2);
const real r = r2*invR; const real invR = RECIP(r);
const real alphaR = EWALD_ALPHA*r; const real alphaR = EWALD_ALPHA*r;
const real expAlphaRSqr = EXP(-alphaR*alphaR); const real expAlphaRSqr = EXP(-alphaR*alphaR);
const real erfAlphaR = erf(alphaR);
real tempForce = 0.0f; real tempForce = 0.0f;
if (erfAlphaR > 1e-6f) { if (alphaR > 1e-6f) {
const real erfAlphaR = erf(alphaR);
const real prefactor = exclusionParams.x*invR; const real prefactor = exclusionParams.x*invR;
tempForce = -prefactor*(erfAlphaR-alphaR*expAlphaRSqr*TWO_OVER_SQRT_PI); tempForce = -prefactor*(erfAlphaR-alphaR*expAlphaRSqr*TWO_OVER_SQRT_PI);
energy -= prefactor*erfAlphaR; energy -= prefactor*erfAlphaR;
...@@ -29,6 +29,7 @@ const real dprefac = eprefac + dar6/6.0f; ...@@ -29,6 +29,7 @@ const real dprefac = eprefac + dar6/6.0f;
energy += coef*(1.0f - expDar2*eprefac); energy += coef*(1.0f - expDar2*eprefac);
tempForce += 6.0f*coef*(1.0f - expDar2*dprefac); tempForce += 6.0f*coef*(1.0f - expDar2*dprefac);
#endif #endif
delta *= tempForce*invR*invR; if (r > 0)
delta *= tempForce*invR*invR;
real3 force1 = -delta; real3 force1 = -delta;
real3 force2 = delta; real3 force2 = delta;
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment