Skip to content

GitLab

Commit 8480944d authored by Peter Eastman's avatar Peter Eastman
Browse files

Bug fixes to Ewald summation

parent e5c065df
Hide whitespace changes
Inline Side-by-side
Showing with 88 additions and 3 deletions
platforms/cuda/src/kernels/kCalculateCDLJForces.h
View file @ 8480944d
...@@ -179,6 +179,14 @@ __global__ void METHOD_NAME(kCalculateCDLJ, Forces_kernel)(unsigned int* workUni ...@@ -179,6 +179,14 @@ __global__ void METHOD_NAME(kCalculateCDLJ, Forces_kernel)(unsigned int* workUni
dEdR += apos.w * psA[j].q * invR * (erfc(alphaR) + 2.0f * alphaR * exp ( - alphaR * alphaR) / SQRT_PI ); dEdR += apos.w * psA[j].q * invR * (erfc(alphaR) + 2.0f * alphaR * exp ( - alphaR * alphaR) / SQRT_PI );
/* E */ /* E */
CDLJ_energy += apos.w * psA[j].q * invR * erfc(alphaR); CDLJ_energy += apos.w * psA[j].q * invR * erfc(alphaR);
bool needCorrection = !(excl & 0x1) && x+tgx != y+j && x+tgx < cSim.atoms && y+j < cSim.atoms;
if (needCorrection)
{
// Subtract off the part of this interaction that was included in the reciprocal space contribution.
dEdR = -apos.w * psA[j].q * invR * (erf(alphaR) - 2.0f * alphaR * exp ( - alphaR * alphaR) / SQRT_PI );
CDLJ_energy = -apos.w * psA[j].q * invR * erf(alphaR);
}
#else #else
dEdR += apos.w * psA[j].q * (invR - 2.0f * cSim.reactionFieldK * r2); dEdR += apos.w * psA[j].q * (invR - 2.0f * cSim.reactionFieldK * r2);
/* E */ /* E */
...@@ -191,7 +199,11 @@ __global__ void METHOD_NAME(kCalculateCDLJ, Forces_kernel)(unsigned int* workUni ...@@ -191,7 +199,11 @@ __global__ void METHOD_NAME(kCalculateCDLJ, Forces_kernel)(unsigned int* workUni
#endif #endif
dEdR *= invR * invR; dEdR *= invR * invR;
#ifdef USE_CUTOFF #ifdef USE_CUTOFF
#ifdef USE_EWALD
if ((!(excl & 0x1) && !needCorrection) || r2 > cSim.nonbondedCutoffSqr)
#else
if (!(excl & 0x1) || r2 > cSim.nonbondedCutoffSqr) if (!(excl & 0x1) || r2 > cSim.nonbondedCutoffSqr)
#endif
#else #else
if (!(excl & 0x1)) if (!(excl & 0x1))
#endif #endif
...@@ -458,6 +470,14 @@ __global__ void METHOD_NAME(kCalculateCDLJ, Forces_kernel)(unsigned int* workUni ...@@ -458,6 +470,14 @@ __global__ void METHOD_NAME(kCalculateCDLJ, Forces_kernel)(unsigned int* workUni
dEdR += apos.w * psA[tj].q * invR * (erfc(alphaR) + 2.0f * alphaR * exp ( - alphaR * alphaR) / SQRT_PI ); dEdR += apos.w * psA[tj].q * invR * (erfc(alphaR) + 2.0f * alphaR * exp ( - alphaR * alphaR) / SQRT_PI );
/* E */ /* E */
CDLJ_energy += apos.w * psA[tj].q * invR * erfc(alphaR); CDLJ_energy += apos.w * psA[tj].q * invR * erfc(alphaR);
bool needCorrection = !(excl & 0x1) && x+tgx != y+tj && x+tgx < cSim.atoms && y+tj < cSim.atoms;
if (needCorrection)
{
// Subtract off the part of this interaction that was included in the reciprocal space contribution.
dEdR = -apos.w * psA[tj].q * invR * (erf(alphaR) - 2.0f * alphaR * exp ( - alphaR * alphaR) / SQRT_PI );
CDLJ_energy = -apos.w * psA[tj].q * invR * erf(alphaR);
}
#else #else
dEdR += apos.w * psA[tj].q * (invR - 2.0f * cSim.reactionFieldK * r2); dEdR += apos.w * psA[tj].q * (invR - 2.0f * cSim.reactionFieldK * r2);
/* E */ /* E */
...@@ -470,7 +490,11 @@ __global__ void METHOD_NAME(kCalculateCDLJ, Forces_kernel)(unsigned int* workUni ...@@ -470,7 +490,11 @@ __global__ void METHOD_NAME(kCalculateCDLJ, Forces_kernel)(unsigned int* workUni
#endif #endif
dEdR *= invR * invR; dEdR *= invR * invR;
#ifdef USE_CUTOFF #ifdef USE_CUTOFF
#ifdef USE_EWALD
if ((!(excl & 0x1) && !needCorrection) || r2 > cSim.nonbondedCutoffSqr)
#else
if (!(excl & 0x1) || r2 > cSim.nonbondedCutoffSqr) if (!(excl & 0x1) || r2 > cSim.nonbondedCutoffSqr)
#endif
#else #else
if (!(excl & 0x1)) if (!(excl & 0x1))
#endif #endif
......
platforms/cuda/src/kernels/kCalculateCDLJObcGbsaForces1.h
View file @ 8480944d
...@@ -194,6 +194,14 @@ __global__ void METHOD_NAME(kCalculateCDLJObcGbsa, Forces1_kernel)(unsigned int* ...@@ -194,6 +194,14 @@ __global__ void METHOD_NAME(kCalculateCDLJObcGbsa, Forces1_kernel)(unsigned int*
dEdR += apos.w * psA[j].q * invR * (erfc(alphaR) + 2.0f * alphaR * exp ( - alphaR * alphaR) / SQRT_PI ); dEdR += apos.w * psA[j].q * invR * (erfc(alphaR) + 2.0f * alphaR * exp ( - alphaR * alphaR) / SQRT_PI );
/* E */ /* E */
CDLJObcGbsa_energy += apos.w * psA[j].q * invR * erfc(alphaR); CDLJObcGbsa_energy += apos.w * psA[j].q * invR * erfc(alphaR);
bool needCorrection = !(excl & 0x1) && x+tgx != y+j && x+tgx < cSim.atoms && y+j < cSim.atoms;
if (needCorrection)
{
// Subtract off the part of this interaction that was included in the reciprocal space contribution.
dEdR = -apos.w * psA[j].q * invR * (erf(alphaR) - 2.0f * alphaR * exp ( - alphaR * alphaR) / SQRT_PI );
CDLJObcGbsa_energy = -apos.w * psA[j].q * invR * erf(alphaR);
}
#else #else
dEdR += apos.w * psA[j].q * (invR - 2.0f * cSim.reactionFieldK * r2); dEdR += apos.w * psA[j].q * (invR - 2.0f * cSim.reactionFieldK * r2);
/* E */ /* E */
...@@ -206,7 +214,11 @@ __global__ void METHOD_NAME(kCalculateCDLJObcGbsa, Forces1_kernel)(unsigned int* ...@@ -206,7 +214,11 @@ __global__ void METHOD_NAME(kCalculateCDLJObcGbsa, Forces1_kernel)(unsigned int*
CDLJObcGbsa_energy += factorX; CDLJObcGbsa_energy += factorX;
#endif #endif
dEdR *= invR * invR; dEdR *= invR * invR;
#ifdef USE_EWALD
if (!(excl & 0x1) && !needCorrection)
#else
if (!(excl & 0x1)) if (!(excl & 0x1))
#endif
{ {
dEdR = 0.0f; dEdR = 0.0f;
/* E */ /* E */
...@@ -556,6 +568,14 @@ __global__ void METHOD_NAME(kCalculateCDLJObcGbsa, Forces1_kernel)(unsigned int* ...@@ -556,6 +568,14 @@ __global__ void METHOD_NAME(kCalculateCDLJObcGbsa, Forces1_kernel)(unsigned int*
dEdR += apos.w * psA[tj].q * invR * (erfc(alphaR) + 2.0f * alphaR * exp ( - alphaR * alphaR) / SQRT_PI ); dEdR += apos.w * psA[tj].q * invR * (erfc(alphaR) + 2.0f * alphaR * exp ( - alphaR * alphaR) / SQRT_PI );
/* E */ /* E */
CDLJObcGbsa_energy += apos.w * psA[tj].q * invR * erfc(alphaR); CDLJObcGbsa_energy += apos.w * psA[tj].q * invR * erfc(alphaR);
bool needCorrection = !(excl & 0x1) && x+tgx != y+tj && x+tgx < cSim.atoms && y+tj < cSim.atoms;
if (needCorrection)
{
// Subtract off the part of this interaction that was included in the reciprocal space contribution.
dEdR = -apos.w * psA[tj].q * invR * (erf(alphaR) - 2.0f * alphaR * exp ( - alphaR * alphaR) / SQRT_PI );
CDLJObcGbsa_energy = -apos.w * psA[tj].q * invR * erf(alphaR);
}
#else #else
dEdR += apos.w * psA[tj].q * (invR - 2.0f * cSim.reactionFieldK * r2); dEdR += apos.w * psA[tj].q * (invR - 2.0f * cSim.reactionFieldK * r2);
/* E */ /* E */
...@@ -568,7 +588,11 @@ __global__ void METHOD_NAME(kCalculateCDLJObcGbsa, Forces1_kernel)(unsigned int* ...@@ -568,7 +588,11 @@ __global__ void METHOD_NAME(kCalculateCDLJObcGbsa, Forces1_kernel)(unsigned int*
CDLJObcGbsa_energy += factorX; CDLJObcGbsa_energy += factorX;
#endif #endif
dEdR *= invR * invR; dEdR *= invR * invR;
#ifdef USE_EWALD
if (!(excl & 0x1) && !needCorrection)
#else
if (!(excl & 0x1)) if (!(excl & 0x1))
#endif
{ {
dEdR = 0.0f; dEdR = 0.0f;
/* E */ /* E */
......
platforms/cuda/src/kernels/kCalculateLocalForces.cu
View file @ 8480944d
...@@ -446,7 +446,7 @@ __global__ void kCalculateLocalForces_kernel() ...@@ -446,7 +446,7 @@ __global__ void kCalculateLocalForces_kernel()
pos += blockDim.x * gridDim.x; pos += blockDim.x * gridDim.x;
} }
if (cSim.nonbondedMethod == NO_CUTOFF) if (cSim.nonbondedMethod == NO_CUTOFF || cSim.nonbondedMethod == EWALD)
{ {
while (pos < cSim.LJ14_offset) while (pos < cSim.LJ14_offset)
{ {
......
platforms/reference/src/ReferenceKernels.cpp
View file @ 8480944d
...@@ -436,7 +436,7 @@ void ReferenceCalcNonbondedForceKernel::executeForces(ContextImpl& context) { ...@@ -436,7 +436,7 @@ void ReferenceCalcNonbondedForceKernel::executeForces(ContextImpl& context) {
clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, 0); clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, 0);
ReferenceBondForce refBondForce; ReferenceBondForce refBondForce;
ReferenceLJCoulomb14 nonbonded14; ReferenceLJCoulomb14 nonbonded14;
if (nonbondedMethod != NoCutoff) if (nonbondedMethod == CutoffNonPeriodic || nonbondedMethod == CutoffPeriodic)
nonbonded14.setUseCutoff(nonbondedCutoff, rfDielectric); nonbonded14.setUseCutoff(nonbondedCutoff, rfDielectric);
refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, 0, 0, 0, nonbonded14); refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, 0, 0, 0, nonbonded14);
} }
...@@ -462,7 +462,7 @@ double ReferenceCalcNonbondedForceKernel::executeEnergy(ContextImpl& context) { ...@@ -462,7 +462,7 @@ double ReferenceCalcNonbondedForceKernel::executeEnergy(ContextImpl& context) {
clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, &energy); clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, &energy);
ReferenceBondForce refBondForce; ReferenceBondForce refBondForce;
ReferenceLJCoulomb14 nonbonded14; ReferenceLJCoulomb14 nonbonded14;
if (nonbondedMethod != NoCutoff) if (nonbondedMethod == CutoffNonPeriodic || nonbondedMethod == CutoffPeriodic)
nonbonded14.setUseCutoff(nonbondedCutoff, rfDielectric); nonbonded14.setUseCutoff(nonbondedCutoff, rfDielectric);
RealOpenMM* energyArray = new RealOpenMM[num14]; RealOpenMM* energyArray = new RealOpenMM[num14];
for (int i = 0; i < num14; ++i) for (int i = 0; i < num14; ++i)
......
platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.cpp
View file @ 8480944d
...@@ -459,6 +459,43 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at ...@@ -459,6 +459,43 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
} }
// Now subtract off the exclusions, since they were implicitly included in the reciprocal space sum.
for (int i = 0; i < numberOfAtoms; i++)
for (int j = 1; j <= exclusions[i][0]; j++)
if (exclusions[i][j] > i) {
int ii = i;
int jj = exclusions[i][j];
RealOpenMM deltaR[2][ReferenceForce::LastDeltaRIndex];
ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxSize, deltaR[0] );
RealOpenMM r = deltaR[0][ReferenceForce::RIndex];
RealOpenMM inverseR = one/(deltaR[0][ReferenceForce::RIndex]);
RealOpenMM alphaR = alphaEwald * r;
RealOpenMM dEdR = atomParameters[ii][QIndex] * atomParameters[jj][QIndex] * inverseR * inverseR * inverseR;
dEdR = (RealOpenMM)(dEdR * (erf(alphaR) - 2 * alphaR * exp ( - alphaR * alphaR) / SQRT_PI ));
// accumulate forces
for( int kk = 0; kk < 3; kk++ ){
RealOpenMM force = dEdR*deltaR[0][kk];
forces[ii][kk] -= force;
forces[jj][kk] += force;
}
// accumulate energies
realSpaceEwaldEnergy = atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR*erf(alphaR);
if( totalEnergy )
*totalEnergy -= realSpaceEwaldEnergy;
if( energyByAtom ){
energyByAtom[ii] -= realSpaceEwaldEnergy;
energyByAtom[jj] -= realSpaceEwaldEnergy;
}
}
// *********************************************************************** // ***********************************************************************
......
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