Skip to content

GitLab

Unverified Commit ed9df876 authored by Peter Eastman's avatar Peter Eastman Committed by GitHub
Browse files

Fixed errors when running on multiple GPUs (#3340)

parent 768963ae
Hide whitespace changes
Inline Side-by-side
Showing with 93 additions and 47 deletions
platforms/common/src/kernels/nonbondedParameters.cc
View file @ ed9df876
...@@ -76,7 +76,7 @@ KERNEL void computeExclusionParameters(GLOBAL real4* RESTRICT posq, GLOBAL real* ...@@ -76,7 +76,7 @@ KERNEL void computeExclusionParameters(GLOBAL real4* RESTRICT posq, GLOBAL real*
#else #else
real chargeProd = charge[atoms.x]*charge[atoms.y]; real chargeProd = charge[atoms.x]*charge[atoms.y];
#endif #endif
#ifdef INCLUDE_LJPME #ifdef INCLUDE_LJPME_EXCEPTIONS
float2 sigEps1 = sigmaEpsilon[atoms.x]; float2 sigEps1 = sigmaEpsilon[atoms.x];
float2 sigEps2 = sigmaEpsilon[atoms.y]; float2 sigEps2 = sigmaEpsilon[atoms.y];
float sigma = sigEps1.x*sigEps2.x; float sigma = sigEps1.x*sigEps2.x;
......
platforms/cuda/src/CudaKernels.cpp
View file @ ed9df876
...@@ -672,6 +672,8 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon ...@@ -672,6 +672,8 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
paramsDefines["HAS_OFFSETS"] = "1"; paramsDefines["HAS_OFFSETS"] = "1";
if (usePosqCharges) if (usePosqCharges)
paramsDefines["USE_POSQ_CHARGES"] = "1"; paramsDefines["USE_POSQ_CHARGES"] = "1";
if (doLJPME)
paramsDefines["INCLUDE_LJPME_EXCEPTIONS"] = "1";
if (nonbondedMethod == Ewald) { if (nonbondedMethod == Ewald) {
// Compute the Ewald parameters. // Compute the Ewald parameters.
...@@ -723,8 +725,16 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon ...@@ -723,8 +725,16 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
defines["TWO_OVER_SQRT_PI"] = cu.doubleToString(2.0/sqrt(M_PI)); defines["TWO_OVER_SQRT_PI"] = cu.doubleToString(2.0/sqrt(M_PI));
defines["USE_EWALD"] = "1"; defines["USE_EWALD"] = "1";
defines["DO_LJPME"] = doLJPME ? "1" : "0"; defines["DO_LJPME"] = doLJPME ? "1" : "0";
if (doLJPME) if (doLJPME) {
defines["EWALD_DISPERSION_ALPHA"] = cu.doubleToString(dispersionAlpha); defines["EWALD_DISPERSION_ALPHA"] = cu.doubleToString(dispersionAlpha);
double invRCut6 = pow(force.getCutoffDistance(), -6);
double dalphaR = dispersionAlpha * force.getCutoffDistance();
double dar2 = dalphaR*dalphaR;
double dar4 = dar2*dar2;
double multShift6 = -invRCut6*(1.0 - exp(-dar2) * (1.0 + dar2 + 0.5*dar4));
defines["INVCUT6"] = cu.doubleToString(invRCut6);
defines["MULTSHIFT6"] = cu.doubleToString(multShift6);
}
if (cu.getContextIndex() == 0) { if (cu.getContextIndex() == 0) {
paramsDefines["INCLUDE_EWALD"] = "1"; paramsDefines["INCLUDE_EWALD"] = "1";
paramsDefines["EWALD_SELF_ENERGY_SCALE"] = cu.doubleToString(ONE_4PI_EPS0*alpha/sqrt(M_PI)); paramsDefines["EWALD_SELF_ENERGY_SCALE"] = cu.doubleToString(ONE_4PI_EPS0*alpha/sqrt(M_PI));
...@@ -790,13 +800,6 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon ...@@ -790,13 +800,6 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
pmeDefines["CHARGE_FROM_SIGEPS"] = "1"; pmeDefines["CHARGE_FROM_SIGEPS"] = "1";
if (cu.getUseDoublePrecision() || cu.getPlatformData().deterministicForces) if (cu.getUseDoublePrecision() || cu.getPlatformData().deterministicForces)
pmeDefines["USE_FIXED_POINT_CHARGE_SPREADING"] = "1"; pmeDefines["USE_FIXED_POINT_CHARGE_SPREADING"] = "1";
double invRCut6 = pow(force.getCutoffDistance(), -6);
double dalphaR = dispersionAlpha * force.getCutoffDistance();
double dar2 = dalphaR*dalphaR;
double dar4 = dar2*dar2;
double multShift6 = -invRCut6*(1.0 - exp(-dar2) * (1.0 + dar2 + 0.5*dar4));
defines["INVCUT6"] = cu.doubleToString(invRCut6);
defines["MULTSHIFT6"] = cu.doubleToString(multShift6);
module = cu.createModule(CudaKernelSources::vectorOps+CommonKernelSources::pme, pmeDefines); module = cu.createModule(CudaKernelSources::vectorOps+CommonKernelSources::pme, pmeDefines);
pmeDispersionFinishSpreadChargeKernel = cu.getKernel(module, "finishSpreadCharge"); pmeDispersionFinishSpreadChargeKernel = cu.getKernel(module, "finishSpreadCharge");
pmeDispersionGridIndexKernel = cu.getKernel(module, "findAtomGridIndex"); pmeDispersionGridIndexKernel = cu.getKernel(module, "findAtomGridIndex");
...@@ -1055,7 +1058,7 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon ...@@ -1055,7 +1058,7 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
// Initialize parameter offsets. // Initialize parameter offsets.
vector<vector<float4> > particleOffsetVec(force.getNumParticles()); vector<vector<float4> > particleOffsetVec(force.getNumParticles());
vector<vector<float4> > exceptionOffsetVec(force.getNumExceptions()); vector<vector<float4> > exceptionOffsetVec(numExceptions);
for (int i = 0; i < force.getNumParticleParameterOffsets(); i++) { for (int i = 0; i < force.getNumParticleParameterOffsets(); i++) {
string param; string param;
int particle; int particle;
...@@ -1076,6 +1079,9 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon ...@@ -1076,6 +1079,9 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
int exception; int exception;
double charge, sigma, epsilon; double charge, sigma, epsilon;
force.getExceptionParameterOffset(i, param, exception, charge, sigma, epsilon); force.getExceptionParameterOffset(i, param, exception, charge, sigma, epsilon);
int index = exceptionIndex[exception];
if (index < startIndex || index >= endIndex)
continue;
auto paramPos = find(paramNames.begin(), paramNames.end(), param); auto paramPos = find(paramNames.begin(), paramNames.end(), param);
int paramIndex; int paramIndex;
if (paramPos == paramNames.end()) { if (paramPos == paramNames.end()) {
...@@ -1084,13 +1090,11 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon ...@@ -1084,13 +1090,11 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
} }
else else
paramIndex = paramPos-paramNames.begin(); paramIndex = paramPos-paramNames.begin();
exceptionOffsetVec[exceptionIndex[exception]].push_back(make_float4(charge, sigma, epsilon, paramIndex)); exceptionOffsetVec[index-startIndex].push_back(make_float4(charge, sigma, epsilon, paramIndex));
} }
paramValues.resize(paramNames.size(), 0.0); paramValues.resize(paramNames.size(), 0.0);
particleParamOffsets.initialize<float4>(cu, max(force.getNumParticleParameterOffsets(), 1), "particleParamOffsets"); particleParamOffsets.initialize<float4>(cu, max(force.getNumParticleParameterOffsets(), 1), "particleParamOffsets");
exceptionParamOffsets.initialize<float4>(cu, max(force.getNumExceptionParameterOffsets(), 1), "exceptionParamOffsets");
particleOffsetIndices.initialize<int>(cu, cu.getPaddedNumAtoms()+1, "particleOffsetIndices"); particleOffsetIndices.initialize<int>(cu, cu.getPaddedNumAtoms()+1, "particleOffsetIndices");
exceptionOffsetIndices.initialize<int>(cu, force.getNumExceptions()+1, "exceptionOffsetIndices");
vector<int> particleOffsetIndicesVec, exceptionOffsetIndicesVec; vector<int> particleOffsetIndicesVec, exceptionOffsetIndicesVec;
vector<float4> p, e; vector<float4> p, e;
for (int i = 0; i < particleOffsetVec.size(); i++) { for (int i = 0; i < particleOffsetVec.size(); i++) {
...@@ -1110,7 +1114,9 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon ...@@ -1110,7 +1114,9 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
particleParamOffsets.upload(p); particleParamOffsets.upload(p);
particleOffsetIndices.upload(particleOffsetIndicesVec); particleOffsetIndices.upload(particleOffsetIndicesVec);
} }
if (force.getNumExceptionParameterOffsets() > 0) { exceptionParamOffsets.initialize<float4>(cu, max((int) e.size(), 1), "exceptionParamOffsets");
exceptionOffsetIndices.initialize<int>(cu, exceptionOffsetIndicesVec.size(), "exceptionOffsetIndices");
if (e.size() > 0) {
exceptionParamOffsets.upload(e); exceptionParamOffsets.upload(e);
exceptionOffsetIndices.upload(exceptionOffsetIndicesVec); exceptionOffsetIndices.upload(exceptionOffsetIndicesVec);
} }
...@@ -1380,20 +1386,28 @@ void CudaCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& context, ...@@ -1380,20 +1386,28 @@ void CudaCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& context,
throw OpenMMException("updateParametersInContext: The nonbonded force kernel does not include Lennard-Jones interactions, because all epsilons were originally 0"); throw OpenMMException("updateParametersInContext: The nonbonded force kernel does not include Lennard-Jones interactions, because all epsilons were originally 0");
} }
} }
set<int> exceptionsWithOffsets;
for (int i = 0; i < force.getNumExceptionParameterOffsets(); i++) {
string param;
int exception;
double charge, sigma, epsilon;
force.getExceptionParameterOffset(i, param, exception, charge, sigma, epsilon);
exceptionsWithOffsets.insert(exception);
}
vector<int> exceptions; vector<int> exceptions;
for (int i = 0; i < force.getNumExceptions(); i++) { for (int i = 0; i < force.getNumExceptions(); i++) {
int particle1, particle2; int particle1, particle2;
double chargeProd, sigma, epsilon; double chargeProd, sigma, epsilon;
force.getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon); force.getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
if (exceptionAtoms.size() > exceptions.size() && make_pair(particle1, particle2) == exceptionAtoms[exceptions.size()]) if (chargeProd != 0.0 || epsilon != 0.0 || exceptionsWithOffsets.find(i) != exceptionsWithOffsets.end())
exceptions.push_back(i); exceptions.push_back(i);
else if (chargeProd != 0.0 || epsilon != 0.0)
throw OpenMMException("updateParametersInContext: The set of non-excluded exceptions has changed");
} }
int numContexts = cu.getPlatformData().contexts.size(); int numContexts = cu.getPlatformData().contexts.size();
int startIndex = cu.getContextIndex()*exceptions.size()/numContexts; int startIndex = cu.getContextIndex()*exceptions.size()/numContexts;
int endIndex = (cu.getContextIndex()+1)*exceptions.size()/numContexts; int endIndex = (cu.getContextIndex()+1)*exceptions.size()/numContexts;
int numExceptions = endIndex-startIndex; int numExceptions = endIndex-startIndex;
if (numExceptions != exceptionAtoms.size())
throw OpenMMException("updateParametersInContext: The set of non-excluded exceptions has changed");
// Record the per-particle parameters. // Record the per-particle parameters.
...@@ -1409,11 +1423,13 @@ void CudaCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& context, ...@@ -1409,11 +1423,13 @@ void CudaCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& context,
// Record the exceptions. // Record the exceptions.
if (numExceptions > 0) { if (numExceptions > 0) {
vector<vector<int> > atoms(numExceptions, vector<int>(2));
vector<float4> baseExceptionParamsVec(numExceptions); vector<float4> baseExceptionParamsVec(numExceptions);
for (int i = 0; i < numExceptions; i++) { for (int i = 0; i < numExceptions; i++) {
int particle1, particle2;
double chargeProd, sigma, epsilon; double chargeProd, sigma, epsilon;
force.getExceptionParameters(exceptions[startIndex+i], atoms[i][0], atoms[i][1], chargeProd, sigma, epsilon); force.getExceptionParameters(exceptions[startIndex+i], particle1, particle2, chargeProd, sigma, epsilon);
if (make_pair(particle1, particle2) != exceptionAtoms[i])
throw OpenMMException("updateParametersInContext: The set of non-excluded exceptions has changed");
baseExceptionParamsVec[i] = make_float4(chargeProd, sigma, epsilon, 0); baseExceptionParamsVec[i] = make_float4(chargeProd, sigma, epsilon, 0);
} }
baseExceptionParams.upload(baseExceptionParamsVec); baseExceptionParams.upload(baseExceptionParamsVec);
......
platforms/cuda/tests/CudaTests.h
View file @ ed9df876
...@@ -6,7 +6,7 @@ ...@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for * * Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. * * Medical Research, grant U54 GM072970. See https://simtk.org. *
* * * *
* Portions copyright (c) 2015-2016 Stanford University and the Authors. * * Portions copyright (c) 2015-2021 Stanford University and the Authors. *
* Authors: Peter Eastman * * Authors: Peter Eastman *
* Contributors: * * Contributors: *
* * * *
...@@ -40,4 +40,6 @@ OpenMM::CudaPlatform platform; ...@@ -40,4 +40,6 @@ OpenMM::CudaPlatform platform;
void initializeTests(int argc, char* argv[]) { void initializeTests(int argc, char* argv[]) {
if (argc > 1) if (argc > 1)
platform.setPropertyDefaultValue("Precision", std::string(argv[1])); platform.setPropertyDefaultValue("Precision", std::string(argv[1]));
if (argc > 2)
platform.setPropertyDefaultValue("DeviceIndex", std::string(argv[2]));
} }
platforms/cuda/tests/TestCudaNonbondedForce.cpp
View file @ ed9df876
...@@ -6,7 +6,7 @@ ...@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for * * Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. * * Medical Research, grant U54 GM072970. See https://simtk.org. *
* * * *
* Portions copyright (c) 2008-2020 Stanford University and the Authors. * * Portions copyright (c) 2008-2021 Stanford University and the Authors. *
* Authors: Peter Eastman * * Authors: Peter Eastman *
* Contributors: * * Contributors: *
* * * *
...@@ -50,11 +50,17 @@ void testParallelComputation(NonbondedForce::NonbondedMethod method) { ...@@ -50,11 +50,17 @@ void testParallelComputation(NonbondedForce::NonbondedMethod method) {
vector<Vec3> positions(numParticles); vector<Vec3> positions(numParticles);
for (int i = 0; i < numParticles; i++) for (int i = 0; i < numParticles; i++)
positions[i] = Vec3(5*genrand_real2(sfmt), 5*genrand_real2(sfmt), 5*genrand_real2(sfmt)); positions[i] = Vec3(5*genrand_real2(sfmt), 5*genrand_real2(sfmt), 5*genrand_real2(sfmt));
force->addGlobalParameter("scale", 0.5);
for (int i = 0; i < numParticles; ++i) for (int i = 0; i < numParticles; ++i)
for (int j = 0; j < i; ++j) { for (int j = 0; j < i; ++j) {
Vec3 delta = positions[i]-positions[j]; Vec3 delta = positions[i]-positions[j];
if (delta.dot(delta) < 0.1) if (delta.dot(delta) < 0.1) {
force->addException(i, j, 0, 1, 0); force->addException(i, j, 0, 1, 0);
}
else if (delta.dot(delta) < 0.2) {
int index = force->addException(i, j, 0.5, 1, 1.0);
force->addExceptionParameterOffset("scale", index, 0.5, 0.4, 0.3);
}
} }
// Create two contexts, one with a single device and one with two devices. // Create two contexts, one with a single device and one with two devices.
...@@ -179,6 +185,7 @@ void runPlatformTests() { ...@@ -179,6 +185,7 @@ void runPlatformTests() {
testParallelComputation(NonbondedForce::NoCutoff); testParallelComputation(NonbondedForce::NoCutoff);
testParallelComputation(NonbondedForce::Ewald); testParallelComputation(NonbondedForce::Ewald);
testParallelComputation(NonbondedForce::PME); testParallelComputation(NonbondedForce::PME);
testParallelComputation(NonbondedForce::LJPME);
testReordering(); testReordering();
testDeterministicForces(); testDeterministicForces();
if (canRunHugeTest()) if (canRunHugeTest())
......
platforms/opencl/src/OpenCLKernels.cpp
View file @ ed9df876
...@@ -677,6 +677,8 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb ...@@ -677,6 +677,8 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
paramsDefines["HAS_OFFSETS"] = "1"; paramsDefines["HAS_OFFSETS"] = "1";
if (usePosqCharges) if (usePosqCharges)
paramsDefines["USE_POSQ_CHARGES"] = "1"; paramsDefines["USE_POSQ_CHARGES"] = "1";
if (doLJPME)
paramsDefines["INCLUDE_LJPME_EXCEPTIONS"] = "1";
if (nonbondedMethod == Ewald) { if (nonbondedMethod == Ewald) {
// Compute the Ewald parameters. // Compute the Ewald parameters.
...@@ -727,8 +729,16 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb ...@@ -727,8 +729,16 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
defines["TWO_OVER_SQRT_PI"] = cl.doubleToString(2.0/sqrt(M_PI)); defines["TWO_OVER_SQRT_PI"] = cl.doubleToString(2.0/sqrt(M_PI));
defines["USE_EWALD"] = "1"; defines["USE_EWALD"] = "1";
defines["DO_LJPME"] = doLJPME ? "1" : "0"; defines["DO_LJPME"] = doLJPME ? "1" : "0";
if (doLJPME) if (doLJPME) {
defines["EWALD_DISPERSION_ALPHA"] = cl.doubleToString(dispersionAlpha); defines["EWALD_DISPERSION_ALPHA"] = cl.doubleToString(dispersionAlpha);
double invRCut6 = pow(force.getCutoffDistance(), -6);
double dalphaR = dispersionAlpha * force.getCutoffDistance();
double dar2 = dalphaR*dalphaR;
double dar4 = dar2*dar2;
double multShift6 = -invRCut6*(1.0 - exp(-dar2) * (1.0 + dar2 + 0.5*dar4));
defines["INVCUT6"] = cl.doubleToString(invRCut6);
defines["MULTSHIFT6"] = cl.doubleToString(multShift6);
}
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));
...@@ -772,15 +782,6 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb ...@@ -772,15 +782,6 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
if (pmeio == NULL) { if (pmeio == NULL) {
// Create required data structures. // Create required data structures.
if (doLJPME) {
double invRCut6 = pow(force.getCutoffDistance(), -6);
double dalphaR = dispersionAlpha * force.getCutoffDistance();
double dar2 = dalphaR*dalphaR;
double dar4 = dar2*dar2;
double multShift6 = -invRCut6*(1.0 - exp(-dar2) * (1.0 + dar2 + 0.5*dar4));
defines["INVCUT6"] = cl.doubleToString(invRCut6);
defines["MULTSHIFT6"] = cl.doubleToString(multShift6);
}
int elementSize = (cl.getUseDoublePrecision() ? sizeof(double) : sizeof(float)); int elementSize = (cl.getUseDoublePrecision() ? sizeof(double) : sizeof(float));
int roundedZSize = PmeOrder*(int) ceil(gridSizeZ/(double) PmeOrder); int roundedZSize = PmeOrder*(int) ceil(gridSizeZ/(double) PmeOrder);
int gridElements = gridSizeX*gridSizeY*roundedZSize; int gridElements = gridSizeX*gridSizeY*roundedZSize;
...@@ -1001,7 +1002,7 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb ...@@ -1001,7 +1002,7 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
// Initialize parameter offsets. // Initialize parameter offsets.
vector<vector<mm_float4> > particleOffsetVec(force.getNumParticles()); vector<vector<mm_float4> > particleOffsetVec(force.getNumParticles());
vector<vector<mm_float4> > exceptionOffsetVec(force.getNumExceptions()); vector<vector<mm_float4> > exceptionOffsetVec(numExceptions);
for (int i = 0; i < force.getNumParticleParameterOffsets(); i++) { for (int i = 0; i < force.getNumParticleParameterOffsets(); i++) {
string param; string param;
int particle; int particle;
...@@ -1022,6 +1023,9 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb ...@@ -1022,6 +1023,9 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
int exception; int exception;
double charge, sigma, epsilon; double charge, sigma, epsilon;
force.getExceptionParameterOffset(i, param, exception, charge, sigma, epsilon); force.getExceptionParameterOffset(i, param, exception, charge, sigma, epsilon);
int index = exceptionIndex[exception];
if (index < startIndex || index >= endIndex)
continue;
auto paramPos = find(paramNames.begin(), paramNames.end(), param); auto paramPos = find(paramNames.begin(), paramNames.end(), param);
int paramIndex; int paramIndex;
if (paramPos == paramNames.end()) { if (paramPos == paramNames.end()) {
...@@ -1030,13 +1034,11 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb ...@@ -1030,13 +1034,11 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
} }
else else
paramIndex = paramPos-paramNames.begin(); paramIndex = paramPos-paramNames.begin();
exceptionOffsetVec[exceptionIndex[exception]].push_back(mm_float4(charge, sigma, epsilon, paramIndex)); exceptionOffsetVec[index-startIndex].push_back(mm_float4(charge, sigma, epsilon, paramIndex));
} }
paramValues.resize(paramNames.size(), 0.0); paramValues.resize(paramNames.size(), 0.0);
particleParamOffsets.initialize<mm_float4>(cl, max(force.getNumParticleParameterOffsets(), 1), "particleParamOffsets"); particleParamOffsets.initialize<mm_float4>(cl, max(force.getNumParticleParameterOffsets(), 1), "particleParamOffsets");
exceptionParamOffsets.initialize<mm_float4>(cl, max(force.getNumExceptionParameterOffsets(), 1), "exceptionParamOffsets");
particleOffsetIndices.initialize<cl_int>(cl, cl.getPaddedNumAtoms()+1, "particleOffsetIndices"); particleOffsetIndices.initialize<cl_int>(cl, cl.getPaddedNumAtoms()+1, "particleOffsetIndices");
exceptionOffsetIndices.initialize<cl_int>(cl, force.getNumExceptions()+1, "exceptionOffsetIndices");
vector<cl_int> particleOffsetIndicesVec, exceptionOffsetIndicesVec; vector<cl_int> particleOffsetIndicesVec, exceptionOffsetIndicesVec;
vector<mm_float4> p, e; vector<mm_float4> p, e;
for (int i = 0; i < particleOffsetVec.size(); i++) { for (int i = 0; i < particleOffsetVec.size(); i++) {
...@@ -1056,7 +1058,9 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb ...@@ -1056,7 +1058,9 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
particleParamOffsets.upload(p); particleParamOffsets.upload(p);
particleOffsetIndices.upload(particleOffsetIndicesVec); particleOffsetIndices.upload(particleOffsetIndicesVec);
} }
if (force.getNumExceptionParameterOffsets() > 0) { exceptionParamOffsets.initialize<mm_float4>(cl, max((int) e.size(), 1), "exceptionParamOffsets");
exceptionOffsetIndices.initialize<cl_int>(cl, exceptionOffsetIndicesVec.size(), "exceptionOffsetIndices");
if (e.size() > 0) {
exceptionParamOffsets.upload(e); exceptionParamOffsets.upload(e);
exceptionOffsetIndices.upload(exceptionOffsetIndicesVec); exceptionOffsetIndices.upload(exceptionOffsetIndicesVec);
} }
...@@ -1521,7 +1525,7 @@ double OpenCLCalcNonbondedForceKernel::execute(ContextImpl& context, bool includ ...@@ -1521,7 +1525,7 @@ double OpenCLCalcNonbondedForceKernel::execute(ContextImpl& context, bool includ
void OpenCLCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& context, const NonbondedForce& force) { void OpenCLCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& context, const NonbondedForce& force) {
// Make sure the new parameters are acceptable. // Make sure the new parameters are acceptable.
if (force.getNumParticles() != cl.getNumAtoms()) if (force.getNumParticles() != cl.getNumAtoms())
throw OpenMMException("updateParametersInContext: The number of particles has changed"); throw OpenMMException("updateParametersInContext: The number of particles has changed");
if (!hasCoulomb || !hasLJ) { if (!hasCoulomb || !hasLJ) {
...@@ -1534,23 +1538,31 @@ void OpenCLCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& contex ...@@ -1534,23 +1538,31 @@ void OpenCLCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& contex
throw OpenMMException("updateParametersInContext: The nonbonded force kernel does not include Lennard-Jones interactions, because all epsilons were originally 0"); throw OpenMMException("updateParametersInContext: The nonbonded force kernel does not include Lennard-Jones interactions, because all epsilons were originally 0");
} }
} }
set<int> exceptionsWithOffsets;
for (int i = 0; i < force.getNumExceptionParameterOffsets(); i++) {
string param;
int exception;
double charge, sigma, epsilon;
force.getExceptionParameterOffset(i, param, exception, charge, sigma, epsilon);
exceptionsWithOffsets.insert(exception);
}
vector<int> exceptions; vector<int> exceptions;
for (int i = 0; i < force.getNumExceptions(); i++) { for (int i = 0; i < force.getNumExceptions(); i++) {
int particle1, particle2; int particle1, particle2;
double chargeProd, sigma, epsilon; double chargeProd, sigma, epsilon;
force.getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon); force.getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
if (exceptionAtoms.size() > exceptions.size() && make_pair(particle1, particle2) == exceptionAtoms[exceptions.size()]) if (chargeProd != 0.0 || epsilon != 0.0 || exceptionsWithOffsets.find(i) != exceptionsWithOffsets.end())
exceptions.push_back(i); exceptions.push_back(i);
else if (chargeProd != 0.0 || epsilon != 0.0)
throw OpenMMException("updateParametersInContext: The set of non-excluded exceptions has changed");
} }
int numContexts = cl.getPlatformData().contexts.size(); int numContexts = cl.getPlatformData().contexts.size();
int startIndex = cl.getContextIndex()*exceptions.size()/numContexts; int startIndex = cl.getContextIndex()*exceptions.size()/numContexts;
int endIndex = (cl.getContextIndex()+1)*exceptions.size()/numContexts; int endIndex = (cl.getContextIndex()+1)*exceptions.size()/numContexts;
int numExceptions = endIndex-startIndex; int numExceptions = endIndex-startIndex;
if (numExceptions != exceptionAtoms.size())
throw OpenMMException("updateParametersInContext: The set of non-excluded exceptions has changed");
// Record the per-particle parameters. // Record the per-particle parameters.
vector<mm_float4> baseParticleParamVec(cl.getPaddedNumAtoms(), mm_float4(0, 0, 0, 0)); vector<mm_float4> baseParticleParamVec(cl.getPaddedNumAtoms(), mm_float4(0, 0, 0, 0));
for (int i = 0; i < force.getNumParticles(); i++) { for (int i = 0; i < force.getNumParticles(); i++) {
double charge, sigma, epsilon; double charge, sigma, epsilon;
...@@ -1562,11 +1574,13 @@ void OpenCLCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& contex ...@@ -1562,11 +1574,13 @@ void OpenCLCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& contex
// Record the exceptions. // Record the exceptions.
if (numExceptions > 0) { if (numExceptions > 0) {
vector<vector<int> > atoms(numExceptions, vector<int>(2));
vector<mm_float4> baseExceptionParamsVec(numExceptions); vector<mm_float4> baseExceptionParamsVec(numExceptions);
for (int i = 0; i < numExceptions; i++) { for (int i = 0; i < numExceptions; i++) {
int particle1, particle2;
double chargeProd, sigma, epsilon; double chargeProd, sigma, epsilon;
force.getExceptionParameters(exceptions[startIndex+i], atoms[i][0], atoms[i][1], chargeProd, sigma, epsilon); force.getExceptionParameters(exceptions[startIndex+i], particle1, particle2, chargeProd, sigma, epsilon);
if (make_pair(particle1, particle2) != exceptionAtoms[i])
throw OpenMMException("updateParametersInContext: The set of non-excluded exceptions has changed");
baseExceptionParamsVec[i] = mm_float4(chargeProd, sigma, epsilon, 0); baseExceptionParamsVec[i] = mm_float4(chargeProd, sigma, epsilon, 0);
} }
baseExceptionParams.upload(baseExceptionParamsVec); baseExceptionParams.upload(baseExceptionParamsVec);
......
platforms/opencl/tests/TestOpenCLNonbondedForce.cpp
View file @ ed9df876
...@@ -6,7 +6,7 @@ ...@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for * * Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. * * Medical Research, grant U54 GM072970. See https://simtk.org. *
* * * *
* Portions copyright (c) 2008-2020 Stanford University and the Authors. * * Portions copyright (c) 2008-2021 Stanford University and the Authors. *
* Authors: Peter Eastman * * Authors: Peter Eastman *
* Contributors: * * Contributors: *
* * * *
...@@ -53,11 +53,17 @@ void testParallelComputation(NonbondedForce::NonbondedMethod method) { ...@@ -53,11 +53,17 @@ void testParallelComputation(NonbondedForce::NonbondedMethod method) {
vector<Vec3> positions(numParticles); vector<Vec3> positions(numParticles);
for (int i = 0; i < numParticles; i++) for (int i = 0; i < numParticles; i++)
positions[i] = Vec3(5*genrand_real2(sfmt), 5*genrand_real2(sfmt), 5*genrand_real2(sfmt)); positions[i] = Vec3(5*genrand_real2(sfmt), 5*genrand_real2(sfmt), 5*genrand_real2(sfmt));
force->addGlobalParameter("scale", 0.5);
for (int i = 0; i < numParticles; ++i) for (int i = 0; i < numParticles; ++i)
for (int j = 0; j < i; ++j) { for (int j = 0; j < i; ++j) {
Vec3 delta = positions[i]-positions[j]; Vec3 delta = positions[i]-positions[j];
if (delta.dot(delta) < 0.1) if (delta.dot(delta) < 0.1) {
force->addException(i, j, 0, 1, 0); force->addException(i, j, 0, 1, 0);
}
else if (delta.dot(delta) < 0.2) {
int index = force->addException(i, j, 0.5, 1, 1.0);
force->addExceptionParameterOffset("scale", index, 0.5, 0.4, 0.3);
}
} }
// Create two contexts, one with a single device and one with two devices. // Create two contexts, one with a single device and one with two devices.
...@@ -154,6 +160,7 @@ void runPlatformTests() { ...@@ -154,6 +160,7 @@ void runPlatformTests() {
testParallelComputation(NonbondedForce::NoCutoff); testParallelComputation(NonbondedForce::NoCutoff);
testParallelComputation(NonbondedForce::Ewald); testParallelComputation(NonbondedForce::Ewald);
testParallelComputation(NonbondedForce::PME); testParallelComputation(NonbondedForce::PME);
testParallelComputation(NonbondedForce::LJPME);
testReordering(); testReordering();
if (canRunHugeTest()) { if (canRunHugeTest()) {
double tol = (platform.getPropertyDefaultValue("Precision") == "single" ? 1e-4 : 1e-5); double tol = (platform.getPropertyDefaultValue("Precision") == "single" ? 1e-4 : 1e-5);
......
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