Commit 7ebaa816 authored by Peter Eastman's avatar Peter Eastman
Browse files

Fixes and minor improvements

parent 73714806
...@@ -709,7 +709,7 @@ private: ...@@ -709,7 +709,7 @@ private:
double ewaldSelfEnergy, dispersionCoefficient, alpha, dispersionAlpha; double ewaldSelfEnergy, dispersionCoefficient, alpha, dispersionAlpha;
int gridSizeX, gridSizeY, gridSizeZ; int gridSizeX, gridSizeY, gridSizeZ;
int dispersionGridSizeX, dispersionGridSizeY, dispersionGridSizeZ; int dispersionGridSizeX, dispersionGridSizeY, dispersionGridSizeZ;
bool hasCoulomb, hasLJ, usePmeQueue, doLJPME, usePosqCharges; bool hasCoulomb, hasLJ, usePmeQueue, doLJPME, usePosqCharges, recomputeParams, hasOffsets;
NonbondedMethod nonbondedMethod; NonbondedMethod nonbondedMethod;
static const int PmeOrder = 5; static const int PmeOrder = 5;
}; };
......
...@@ -1636,6 +1636,8 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb ...@@ -1636,6 +1636,8 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
if (cl.getContextIndex() == 0) { if (cl.getContextIndex() == 0) {
paramsDefines["INCLUDE_EWALD"] = "1"; paramsDefines["INCLUDE_EWALD"] = "1";
paramsDefines["EWALD_SELF_ENERGY_SCALE"] = cl.doubleToString(ONE_4PI_EPS0*alpha/sqrt(M_PI)); paramsDefines["EWALD_SELF_ENERGY_SCALE"] = cl.doubleToString(ONE_4PI_EPS0*alpha/sqrt(M_PI));
for (int i = 0; i < numParticles; i++)
ewaldSelfEnergy += baseParticleParamVec[i].x*baseParticleParamVec[i].x*ONE_4PI_EPS0*alpha/sqrt(M_PI);
// Create the reciprocal space kernels. // Create the reciprocal space kernels.
...@@ -1675,9 +1677,13 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb ...@@ -1675,9 +1677,13 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
if (cl.getContextIndex() == 0) { if (cl.getContextIndex() == 0) {
paramsDefines["INCLUDE_EWALD"] = "1"; paramsDefines["INCLUDE_EWALD"] = "1";
paramsDefines["EWALD_SELF_ENERGY_SCALE"] = cl.doubleToString(ONE_4PI_EPS0*alpha/sqrt(M_PI)); paramsDefines["EWALD_SELF_ENERGY_SCALE"] = cl.doubleToString(ONE_4PI_EPS0*alpha/sqrt(M_PI));
for (int i = 0; i < numParticles; i++)
ewaldSelfEnergy += baseParticleParamVec[i].x*baseParticleParamVec[i].x*ONE_4PI_EPS0*alpha/sqrt(M_PI);
if (doLJPME) { if (doLJPME) {
paramsDefines["INCLUDE_LJPME"] = "1"; paramsDefines["INCLUDE_LJPME"] = "1";
paramsDefines["LJPME_SELF_ENERGY_SCALE"] = cl.doubleToString(pow(dispersionAlpha, 6)/3.0); paramsDefines["LJPME_SELF_ENERGY_SCALE"] = cl.doubleToString(pow(dispersionAlpha, 6)/3.0);
for (int i = 0; i < numParticles; i++)
ewaldSelfEnergy += baseParticleParamVec[i].z*pow(baseParticleParamVec[i].y*dispersionAlpha, 6)/3.0;
} }
pmeDefines["PME_ORDER"] = cl.intToString(PmeOrder); pmeDefines["PME_ORDER"] = cl.intToString(PmeOrder);
pmeDefines["NUM_ATOMS"] = cl.intToString(numParticles); pmeDefines["NUM_ATOMS"] = cl.intToString(numParticles);
...@@ -1973,6 +1979,7 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb ...@@ -1973,6 +1979,7 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
exceptionOffsetIndices.upload(exceptionOffsetIndicesVec); exceptionOffsetIndices.upload(exceptionOffsetIndicesVec);
} }
globalParams.initialize(cl, max((int) paramValues.size(), 1), cl.getUseDoublePrecision() ? sizeof(double) : sizeof(float), "globalParams"); globalParams.initialize(cl, max((int) paramValues.size(), 1), cl.getUseDoublePrecision() ? sizeof(double) : sizeof(float), "globalParams");
recomputeParams = true;
// Initialize the kernel for updating parameters. // Initialize the kernel for updating parameters.
...@@ -2143,6 +2150,7 @@ double OpenCLCalcNonbondedForceKernel::execute(ContextImpl& context, bool includ ...@@ -2143,6 +2150,7 @@ double OpenCLCalcNonbondedForceKernel::execute(ContextImpl& context, bool includ
} }
} }
if (paramChanged) { if (paramChanged) {
recomputeParams = true;
if (cl.getUseDoublePrecision()) if (cl.getUseDoublePrecision())
globalParams.upload(paramValues); globalParams.upload(paramValues);
else { else {
...@@ -2152,8 +2160,12 @@ double OpenCLCalcNonbondedForceKernel::execute(ContextImpl& context, bool includ ...@@ -2152,8 +2160,12 @@ double OpenCLCalcNonbondedForceKernel::execute(ContextImpl& context, bool includ
globalParams.upload(v); globalParams.upload(v);
} }
} }
double energy = (includeReciprocal ? ewaldSelfEnergy : 0.0);
if (recomputeParams || hasOffsets) {
computeParamsKernel.setArg<cl_int>(1, includeEnergy && includeReciprocal); computeParamsKernel.setArg<cl_int>(1, includeEnergy && includeReciprocal);
cl.executeKernel(computeParamsKernel, cl.getPaddedNumAtoms()); cl.executeKernel(computeParamsKernel, cl.getPaddedNumAtoms());
energy = 0.0; // The Ewald self energy was computed in the kernel.
}
// Do reciprocal space calculations. // Do reciprocal space calculations.
...@@ -2391,7 +2403,6 @@ double OpenCLCalcNonbondedForceKernel::execute(ContextImpl& context, bool includ ...@@ -2391,7 +2403,6 @@ double OpenCLCalcNonbondedForceKernel::execute(ContextImpl& context, bool includ
cl.restoreDefaultQueue(); cl.restoreDefaultQueue();
} }
} }
double energy = (includeReciprocal ? ewaldSelfEnergy : 0.0);
if (dispersionCoefficient != 0.0 && includeDirect) { if (dispersionCoefficient != 0.0 && includeDirect) {
mm_double4 boxSize = cl.getPeriodicBoxSizeDouble(); mm_double4 boxSize = cl.getPeriodicBoxSizeDouble();
energy += dispersionCoefficient/(boxSize.x*boxSize.y*boxSize.z); energy += dispersionCoefficient/(boxSize.x*boxSize.y*boxSize.z);
...@@ -2457,6 +2468,7 @@ void OpenCLCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& contex ...@@ -2457,6 +2468,7 @@ void OpenCLCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& contex
if (force.getUseDispersionCorrection() && cl.getContextIndex() == 0 && (nonbondedMethod == CutoffPeriodic || nonbondedMethod == Ewald || nonbondedMethod == PME)) if (force.getUseDispersionCorrection() && cl.getContextIndex() == 0 && (nonbondedMethod == CutoffPeriodic || nonbondedMethod == Ewald || nonbondedMethod == PME))
dispersionCoefficient = NonbondedForceImpl::calcDispersionCorrection(context.getSystem(), force); dispersionCoefficient = NonbondedForceImpl::calcDispersionCorrection(context.getSystem(), force);
cl.invalidateMolecules(info); cl.invalidateMolecules(info);
recomputeParams = true;
} }
void OpenCLCalcNonbondedForceKernel::getPMEParameters(double& alpha, int& nx, int& ny, int& nz) const { void OpenCLCalcNonbondedForceKernel::getPMEParameters(double& alpha, int& nx, int& ny, int& nz) const {
......
...@@ -26,17 +26,17 @@ __kernel void computeParameters(__global mixed* restrict energyBuffer, int inclu ...@@ -26,17 +26,17 @@ __kernel void computeParameters(__global mixed* restrict energyBuffer, int inclu
} }
#endif #endif
#ifdef USE_POSQ_CHARGES #ifdef USE_POSQ_CHARGES
posq[i].w = params[0]; posq[i].w = params.x;
#else #else
charge[i] = params[0]; charge[i] = params.x;
#endif #endif
sigmaEpsilon[i] = (float2) (0.5f*params[1], 2*SQRT(params[2])); sigmaEpsilon[i] = (float2) (0.5f*params.y, 2*SQRT(params.z));
#ifdef INCLUDE_EWALD #ifdef INCLUDE_EWALD
energy -= EWALD_SELF_ENERGY_SCALE*params[0]*params[0]; energy -= EWALD_SELF_ENERGY_SCALE*params.x*params.x;
#endif #endif
#ifdef INCLUDE_LJPME #ifdef INCLUDE_LJPME
real sig3 = params[1]*params[1]*params[1]; real sig3 = params.y*params.y*params.y;
energy += LJPME_SELF_ENERGY_SCALE*sig3*sig3*params[2]; energy += LJPME_SELF_ENERGY_SCALE*sig3*sig3*params.z;
#endif #endif
} }
...@@ -55,7 +55,7 @@ __kernel void computeParameters(__global mixed* restrict energyBuffer, int inclu ...@@ -55,7 +55,7 @@ __kernel void computeParameters(__global mixed* restrict energyBuffer, int inclu
params.z += value*offset.z; params.z += value*offset.z;
} }
#endif #endif
exceptionParams[i] = (float4) ((float) (138.935456f*params[0]), (float) params[1], (float) (4*params[2]), 0); exceptionParams[i] = (float4) ((float) (138.935456f*params.x), (float) params.y, (float) (4*params.z), 0);
} }
#endif #endif
if (includeSelfEnergy) if (includeSelfEnergy)
......
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