Skip to content

GitLab

Commit 2e3d7cf7 authored by peastman's avatar peastman
Browse files

Improved accuracy of correction for exclusions with PME

parent fa80fbf7
Hide whitespace changes
Inline Side-by-side
Showing with 5 additions and 5 deletions
platforms/cpu/src/CpuNonbondedForce.cpp
View file @ 2e3d7cf7
/* Portions copyright (c) 2006-2013 Stanford University and Simbios. /* Portions copyright (c) 2006-2015 Stanford University and Simbios.
* Contributors: Pande Group * Contributors: Pande Group
* *
* Permission is hereby granted, free of charge, to any person obtaining * Permission is hereby granted, free of charge, to any person obtaining
...@@ -364,15 +364,15 @@ void CpuNonbondedForce::threadComputeDirect(ThreadPool& threads, int threadIndex ...@@ -364,15 +364,15 @@ void CpuNonbondedForce::threadComputeDirect(ThreadPool& threads, int threadIndex
float inverseR = 1/r; float inverseR = 1/r;
float chargeProd = ONE_4PI_EPS0*posq[4*i+3]*posq[4*j+3]; float chargeProd = ONE_4PI_EPS0*posq[4*i+3]*posq[4*j+3];
float alphaR = alphaEwald*r; float alphaR = alphaEwald*r;
float erfcAlphaR = erfcApprox(alphaR); float erfAlphaR = erf(alphaR);
if (1-erfcAlphaR > 1e-6f) { if (erfAlphaR > 1e-6f) {
float dEdR = (float) (chargeProd * inverseR * inverseR * inverseR); float dEdR = (float) (chargeProd * inverseR * inverseR * inverseR);
dEdR = (float) (dEdR * (1.0f-erfcAlphaR-TWO_OVER_SQRT_PI*alphaR*exp(-alphaR*alphaR))); dEdR = (float) (dEdR * (erfAlphaR-TWO_OVER_SQRT_PI*alphaR*exp(-alphaR*alphaR)));
fvec4 result = deltaR*dEdR; fvec4 result = deltaR*dEdR;
(fvec4(forces+4*i)-result).store(forces+4*i); (fvec4(forces+4*i)-result).store(forces+4*i);
(fvec4(forces+4*j)+result).store(forces+4*j); (fvec4(forces+4*j)+result).store(forces+4*j);
if (includeEnergy) if (includeEnergy)
threadEnergy[threadIndex] -= chargeProd*inverseR*(1.0f-erfcAlphaR); threadEnergy[threadIndex] -= chargeProd*inverseR*erfAlphaR;
} }
} }
} }
......
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