Skip to content

GitLab

Commit 43e6e68c authored by Peter Eastman's avatar Peter Eastman
Browse files

CustomGBForce allows single-particle values to depend on position

parent 74ef687d
Hide whitespace changes
Inline Side-by-side
Showing with 401 additions and 72 deletions
openmmapi/include/openmm/CustomGBForce.h
View file @ 43e6e68c
......@@ -61,7 +61,7 @@ namespace OpenMM {
*
* <ul>
* <li><b>Single Particle</b>: The expression is evaluated once for each particle in the System. In the case of a computed
* value, this means the value for a particle depends only on other properties of that particle (its parameters and other
* value, this means the value for a particle depends only on other properties of that particle (its position, parameters, and other
* computed values). In the case of an energy term, it means each particle makes an independent contribution to the System
* energy.</li>
* <li><b>Particle Pairs</b>: The expression is evaluated for every pair of particles in the system. In the case of a computed
......@@ -329,9 +329,9 @@ public:
* @param name the name of the value
* @param expression an algebraic expression to evaluate when calculating the computed value. If the
* ComputationType is SingleParticle, the expression is evaluated independently
* for each particle, and may depend on the per-particle parameters and previous
* computed values for that particle. If the ComputationType is ParticlePair or
* ParticlePairNoExclusions, the expression is evaluated once for every other
* for each particle, and may depend on its x, y, and z coordinates, as well as the per-particle
* parameters and previous computed values for that particle. If the ComputationType is ParticlePair
* or ParticlePairNoExclusions, the expression is evaluated once for every other
* particle in the system and summed to get the final value. In the latter case,
* the expression may depend on the distance r between the two particles, and on
* the per-particle parameters and previous computed values for each of them.
......@@ -347,9 +347,9 @@ public:
* @param name the name of the value
* @param expression an algebraic expression to evaluate when calculating the computed value. If the
* ComputationType is SingleParticle, the expression is evaluated independently
* for each particle, and may depend on the per-particle parameters and previous
* computed values for that particle. If the ComputationType is ParticlePair or
* ParticlePairNoExclusions, the expression is evaluated once for every other
* for each particle, and may depend on its x, y, and z coordinates, as well as the per-particle
* parameters and previous computed values for that particle. If the ComputationType is ParticlePair
* or ParticlePairNoExclusions, the expression is evaluated once for every other
* particle in the system and summed to get the final value. In the latter case,
* the expression may depend on the distance r between the two particles, and on
* the per-particle parameters and previous computed values for each of them.
......@@ -365,9 +365,9 @@ public:
* @param name the name of the value
* @param expression an algebraic expression to evaluate when calculating the computed value. If the
* ComputationType is SingleParticle, the expression is evaluated independently
* for each particle, and may depend on the per-particle parameters and previous
* computed values for that particle. If the ComputationType is ParticlePair or
* ParticlePairNoExclusions, the expression is evaluated once for every other
* for each particle, and may depend on its x, y, and z coordinates, as well as the per-particle
* parameters and previous computed values for that particle. If the ComputationType is ParticlePair
* or ParticlePairNoExclusions, the expression is evaluated once for every other
* particle in the system and summed to get the final value. In the latter case,
* the expression may depend on the distance r between the two particles, and on
* the per-particle parameters and previous computed values for each of them.
......@@ -381,8 +381,8 @@ public:
*
* @param expression an algebraic expression to evaluate when calculating the energy. If the
* ComputationType is SingleParticle, the expression is evaluated once
* for each particle, and may depend on the per-particle parameters and
* computed values for that particle. If the ComputationType is ParticlePair or
* for each particle, and may depend on its x, y, and z coordinates, as well as the per-particle
* parameters and computed values for that particle. If the ComputationType is ParticlePair or
* ParticlePairNoExclusions, the expression is evaluated once for every pair of
* particles in the system. In the latter case,
* the expression may depend on the distance r between the two particles, and on
......@@ -398,8 +398,8 @@ public:
* @param index the index of the term for which to get parameters
* @param expression an algebraic expression to evaluate when calculating the energy. If the
* ComputationType is SingleParticle, the expression is evaluated once
* for each particle, and may depend on the per-particle parameters and
* computed values for that particle. If the ComputationType is ParticlePair or
* for each particle, and may depend on its x, y, and z coordinates, as well as the per-particle
* parameters and computed values for that particle. If the ComputationType is ParticlePair or
* ParticlePairNoExclusions, the expression is evaluated once for every pair of
* particles in the system. In the latter case,
* the expression may depend on the distance r between the two particles, and on
......@@ -415,8 +415,8 @@ public:
* @param index the index of the term for which to set parameters
* @param expression an algebraic expression to evaluate when calculating the energy. If the
* ComputationType is SingleParticle, the expression is evaluated once
* for each particle, and may depend on the per-particle parameters and
* computed values for that particle. If the ComputationType is ParticlePair or
* for each particle, and may depend on its x, y, and z coordinates, as well as the per-particle
* parameters and computed values for that particle. If the ComputationType is ParticlePair or
* ParticlePairNoExclusions, the expression is evaluated once for every pair of
* particles in the system. In the latter case,
* the expression may depend on the distance r between the two particles, and on
......
platforms/opencl/src/OpenCLKernels.cpp
View file @ 43e6e68c
......@@ -59,6 +59,13 @@ static string intToString(int value) {
return s.str();
}
static bool isZeroExpression(const Lepton::ParsedExpression& expression) {
const Lepton::Operation& op = expression.getRootNode().getOperation();
if (op.getId() != Lepton::Operation::CONSTANT)
return false;
return (dynamic_cast<const Lepton::Operation::Constant&>(op).getValue() == 0.0);
}
void OpenCLCalcForcesAndEnergyKernel::initialize(const System& system) {
}
......@@ -1793,6 +1800,16 @@ void OpenCLCalcCustomGBForceKernel::initialize(const System& system, const Custo
// Record derivatives of expressions needed for the chain rule terms.
vector<vector<Lepton::ParsedExpression> > valueGradientExpressions(force.getNumComputedValues());
bool needParameterGradient = false;
for (int i = 1; i < force.getNumComputedValues(); i++) {
Lepton::ParsedExpression ex = Lepton::Parser::parse(computedValueExpressions[i], functions).optimize();
valueGradientExpressions[i].push_back(ex.differentiate("x").optimize());
valueGradientExpressions[i].push_back(ex.differentiate("y").optimize());
valueGradientExpressions[i].push_back(ex.differentiate("z").optimize());
if (!isZeroExpression(valueGradientExpressions[i][0]) || !isZeroExpression(valueGradientExpressions[i][1]) || !isZeroExpression(valueGradientExpressions[i][2]))
needParameterGradient = true;
}
vector<vector<Lepton::ParsedExpression> > energyDerivExpressions(force.getNumEnergyTerms());
for (int i = 0; i < force.getNumEnergyTerms(); i++) {
string expression;
......@@ -1900,6 +1917,9 @@ void OpenCLCalcCustomGBForceKernel::initialize(const System& system, const Custo
}
reductionSource << "local_values" << computedValues->getParameterSuffix(0) << " = sum;\n";
map<string, string> variables;
variables["x"] = "pos.x";
variables["y"] = "pos.y";
variables["z"] = "pos.z";
for (int i = 0; i < force.getNumPerParticleParameters(); i++)
variables[force.getPerParticleParameterName(i)] = "params"+params->getParameterSuffix(i, "[index]");
for (int i = 0; i < force.getNumGlobalParameters(); i++)
......@@ -2036,7 +2056,6 @@ void OpenCLCalcCustomGBForceKernel::initialize(const System& system, const Custo
// Create the kernel to reduce the derivatives and calculate per-particle energy terms.
stringstream compute, extraArgs, reduce;
map<string, Lepton::ParsedExpression> energyExpressions;
if (force.getNumGlobalParameters() > 0)
extraArgs << ", __constant float* globals";
for (int i = 0; i < (int) params->getBuffers().size(); i++) {
......@@ -2057,27 +2076,42 @@ void OpenCLCalcCustomGBForceKernel::initialize(const System& system, const Custo
compute << buffer.getType() << " deriv" << index << " = derivBuffers" << index << "[index];\n";
}
map<string, string> variables;
variables["x"] = "pos.x";
variables["y"] = "pos.y";
variables["z"] = "pos.z";
for (int i = 0; i < force.getNumPerParticleParameters(); i++)
variables[force.getPerParticleParameterName(i)] = "params"+params->getParameterSuffix(i, "[index]");
for (int i = 0; i < force.getNumGlobalParameters(); i++)
variables[force.getGlobalParameterName(i)] = "globals["+intToString(i)+"]";
for (int i = 0; i < force.getNumComputedValues(); i++)
variables[computedValueNames[i]] = "values"+computedValues->getParameterSuffix(i, "[index]");
map<string, Lepton::ParsedExpression> energyExpressions;
for (int i = 0; i < force.getNumEnergyTerms(); i++) {
string expression;
CustomGBForce::ComputationType type;
force.getEnergyTermParameters(i, expression, type);
if (type != CustomGBForce::SingleParticle)
continue;
energyExpressions["/*"+intToString(i+1)+"*/ energy += "] = Lepton::Parser::parse(expression, functions).optimize();
Lepton::ParsedExpression parsed = Lepton::Parser::parse(expression, functions).optimize();
energyExpressions["/*"+intToString(i+1)+"*/ energy += "] = parsed;
for (int j = 0; j < force.getNumComputedValues(); j++)
energyExpressions["/*"+intToString(i+1)+"*/ deriv"+energyDerivs->getParameterSuffix(j)+" += "] = energyDerivExpressions[i][j];
Lepton::ParsedExpression gradx = parsed.differentiate("x").optimize();
Lepton::ParsedExpression grady = parsed.differentiate("y").optimize();
Lepton::ParsedExpression gradz = parsed.differentiate("z").optimize();
if (!isZeroExpression(gradx))
energyExpressions["/*"+intToString(i+1)+"*/ force.x -= "] = gradx;
if (!isZeroExpression(grady))
energyExpressions["/*"+intToString(i+1)+"*/ force.y -= "] = grady;
if (!isZeroExpression(gradz))
energyExpressions["/*"+intToString(i+1)+"*/ force.z -= "] = gradz;
}
compute << OpenCLExpressionUtilities::createExpressions(energyExpressions, variables, functionDefinitions, "temp", prefix+"functionParams");
for (int i = 0; i < (int) energyDerivs->getBuffers().size(); i++) {
string index = intToString(i+1);
compute << "derivBuffers" << index << "[index] = deriv" << index << ";\n";
}
compute << "forceBuffers[index] = forceBuffers[index]+force;\n";
map<string, string> replacements;
replacements["PARAMETER_ARGUMENTS"] = extraArgs.str()+tableArgs.str();
replacements["REDUCE_DERIVATIVES"] = reduce.str();
......@@ -2088,7 +2122,7 @@ void OpenCLCalcCustomGBForceKernel::initialize(const System& system, const Custo
perParticleEnergyKernel = cl::Kernel(program, "computePerParticleEnergy");
}
{
// Create the code to calculate chain rules terms (as part of the default nonbonded kernel).
// Create the code to calculate chain rules terms (possibly as part of the default nonbonded kernel).
map<string, string> globalVariables;
for (int i = 0; i < force.getNumGlobalParameters(); i++) {
......@@ -2112,8 +2146,20 @@ void OpenCLCalcCustomGBForceKernel::initialize(const System& system, const Custo
derivExpressions["float dV0dR1 = "] = dVdR;
derivExpressions["float dV0dR2 = "] = dVdR.renameVariables(rename);
chainSource << OpenCLExpressionUtilities::createExpressions(derivExpressions, variables, functionDefinitions, prefix+"temp0_", prefix+"functionParams");
chainSource << "tempForce -= dV0dR1*" << prefix << "dEdV" << energyDerivs->getParameterSuffix(0, "1") << ";\n";
chainSource << "tempForce -= dV0dR2*" << prefix << "dEdV" << energyDerivs->getParameterSuffix(0, "2") << ";\n";
if (needParameterGradient) {
chainSource << "float4 grad1_0_1 = dV0dR1*delta*invR;\n";
chainSource << "float4 grad1_0_2 = dV0dR2*delta*invR;\n";
chainSource << "float4 grad2_0_1 = -grad1_0_1;\n";
chainSource << "float4 grad2_0_2 = -grad1_0_2;\n";
chainSource << "tempForce1 -= grad1_0_1*" << prefix << "dEdV" << energyDerivs->getParameterSuffix(0, "1") << ";\n";
chainSource << "tempForce1 -= grad1_0_2*" << prefix << "dEdV" << energyDerivs->getParameterSuffix(0, "2") << ";\n";
chainSource << "tempForce2 -= grad2_0_1*" << prefix << "dEdV" << energyDerivs->getParameterSuffix(0, "1") << ";\n";
chainSource << "tempForce2 -= grad2_0_2*" << prefix << "dEdV" << energyDerivs->getParameterSuffix(0, "2") << ";\n";
}
else {
chainSource << "tempForce -= dV0dR1*" << prefix << "dEdV" << energyDerivs->getParameterSuffix(0, "1") << ";\n";
chainSource << "tempForce -= dV0dR2*" << prefix << "dEdV" << energyDerivs->getParameterSuffix(0, "2") << ";\n";
}
variables = globalVariables;
map<string, string> rename1;
map<string, string> rename2;
......@@ -2132,40 +2178,96 @@ void OpenCLCalcCustomGBForceKernel::initialize(const System& system, const Custo
rename2[name] = name+"2";
if (i == 0)
continue;
chainSource << "float dV"+intToString(i)+"dR1 = 0;\n";
chainSource << "float dV"+intToString(i)+"dR2 = 0;\n";
for (int j = 0; j < i; j++) {
Lepton::ParsedExpression dVdV = Lepton::Parser::parse(computedValueExpressions[i], functions).differentiate(computedValueNames[j]).optimize();
string is = intToString(i);
if (needParameterGradient) {
chainSource << "float4 grad1_"+is+"_1 = 0;\n";
chainSource << "float4 grad1_"+is+"_2 = 0;\n";
chainSource << "float4 grad2_"+is+"_1 = 0;\n";
chainSource << "float4 grad2_"+is+"_2 = 0;\n";
for (int j = 0; j < i; j++) {
string js = intToString(j);
Lepton::ParsedExpression dVdV = Lepton::Parser::parse(computedValueExpressions[i], functions).differentiate(computedValueNames[j]).optimize();
derivExpressions.clear();
derivExpressions["float dV"+is+"dV"+js+"_1 = "] = dVdV.renameVariables(rename1);
derivExpressions["float dV"+is+"dV"+js+"_2 = "] = dVdV.renameVariables(rename2);
chainSource << OpenCLExpressionUtilities::createExpressions(derivExpressions, variables, functionDefinitions, prefix+"temp"+is+"_"+js+"_", prefix+"functionParams");
chainSource << "grad1_"+is+"_1 += dV"+is+"dV"+js+"_1*grad1_"+js+"_1;\n";
chainSource << "grad2_"+is+"_1 += dV"+is+"dV"+js+"_1*grad2_"+js+"_1;\n";
chainSource << "grad1_"+is+"_2 += dV"+is+"dV"+js+"_2*grad1_"+js+"_2;\n";
chainSource << "grad2_"+is+"_2 += dV"+is+"dV"+js+"_2*grad2_"+js+"_2;\n";
}
derivExpressions.clear();
derivExpressions["dV"+intToString(i)+"dR1 += dV"+intToString(j)+"dR1*"] = dVdV.renameVariables(rename1);
derivExpressions["dV"+intToString(i)+"dR2 += dV"+intToString(j)+"dR2*"] = dVdV.renameVariables(rename2);
chainSource << OpenCLExpressionUtilities::createExpressions(derivExpressions, variables, functionDefinitions, prefix+"temp"+intToString(i)+"_"+intToString(j)+"_", prefix+"functionParams");
if (!isZeroExpression(valueGradientExpressions[i][0])) {
derivExpressions["grad1_"+is+"_1.x -= "] = valueGradientExpressions[i][0].renameVariables(rename1);
derivExpressions["grad2_"+is+"_2.x -= "] = valueGradientExpressions[i][0].renameVariables(rename2);
}
if (!isZeroExpression(valueGradientExpressions[i][1])) {
derivExpressions["grad1_"+is+"_1.y -= "] = valueGradientExpressions[i][1].renameVariables(rename1);
derivExpressions["grad2_"+is+"_2.y -= "] = valueGradientExpressions[i][1].renameVariables(rename2);
}
if (!isZeroExpression(valueGradientExpressions[i][2])) {
derivExpressions["grad1_"+is+"_1.z -= "] = valueGradientExpressions[i][2].renameVariables(rename1);
derivExpressions["grad2_"+is+"_2.z -= "] = valueGradientExpressions[i][2].renameVariables(rename2);
}
chainSource << OpenCLExpressionUtilities::createExpressions(derivExpressions, variables, functionDefinitions, prefix+"temp"+is+"_", prefix+"functionParams");
chainSource << "tempForce1 -= grad1_"<<is<<"_1*"<<prefix<<"dEdV"<<energyDerivs->getParameterSuffix(i, "1")<<";\n";
chainSource << "tempForce2 -= grad2_"<<is<<"_1*"<<prefix<<"dEdV"<<energyDerivs->getParameterSuffix(i, "1")<<";\n";
chainSource << "tempForce1 -= grad1_"<<is<<"_2*"<<prefix<<"dEdV"<<energyDerivs->getParameterSuffix(i, "2")<<";\n";
chainSource << "tempForce2 -= grad2_"<<is<<"_2*"<<prefix<<"dEdV"<<energyDerivs->getParameterSuffix(i, "2")<<";\n";
}
else {
chainSource << "float dV"+is+"dR1 = 0;\n";
chainSource << "float dV"+is+"dR2 = 0;\n";
for (int j = 0; j < i; j++) {
string js = intToString(j);
Lepton::ParsedExpression dVdV = Lepton::Parser::parse(computedValueExpressions[i], functions).differentiate(computedValueNames[j]).optimize();
derivExpressions.clear();
derivExpressions["dV"+is+"dR1 += dV"+js+"dR1*"] = dVdV.renameVariables(rename1);
derivExpressions["dV"+is+"dR2 += dV"+js+"dR2*"] = dVdV.renameVariables(rename2);
chainSource << OpenCLExpressionUtilities::createExpressions(derivExpressions, variables, functionDefinitions, prefix+"temp"+is+"_"+js+"_", prefix+"functionParams");
}
chainSource << "tempForce -= dV"<< is << "dR1*" << prefix << "dEdV" << energyDerivs->getParameterSuffix(i, "1") << ";\n";
chainSource << "tempForce -= dV"<< is << "dR2*" << prefix << "dEdV" << energyDerivs->getParameterSuffix(i, "2") << ";\n";
}
chainSource << "tempForce -= dV"<< intToString(i) << "dR1*" << prefix << "dEdV" << energyDerivs->getParameterSuffix(i, "1") << ";\n";
chainSource << "tempForce -= dV"<< intToString(i) << "dR2*" << prefix << "dEdV" << energyDerivs->getParameterSuffix(i, "2") << ";\n";
}
map<string, string> replacements;
replacements["COMPUTE_FORCE"] = chainSource.str();
string source = cl.replaceStrings(OpenCLKernelSources::customGBChainRule, replacements);
cl.getNonbondedUtilities().addInteraction(useCutoff, usePeriodic, true, force.getCutoffDistance(), exclusionList, source);
vector<OpenCLNonbondedUtilities::ParameterInfo> parameters;
vector<OpenCLNonbondedUtilities::ParameterInfo> arguments;
for (int i = 0; i < (int) params->getBuffers().size(); i++) {
const OpenCLNonbondedUtilities::ParameterInfo& buffer = params->getBuffers()[i];
string paramName = prefix+"params"+intToString(i+1);
cl.getNonbondedUtilities().addParameter(OpenCLNonbondedUtilities::ParameterInfo(paramName, buffer.getComponentType(), buffer.getNumComponents(), buffer.getSize(), buffer.getMemory()));
parameters.push_back(OpenCLNonbondedUtilities::ParameterInfo(paramName, buffer.getComponentType(), buffer.getNumComponents(), buffer.getSize(), buffer.getMemory()));
}
for (int i = 0; i < (int) computedValues->getBuffers().size(); i++) {
const OpenCLNonbondedUtilities::ParameterInfo& buffer = computedValues->getBuffers()[i];
string paramName = prefix+"values"+intToString(i+1);
cl.getNonbondedUtilities().addParameter(OpenCLNonbondedUtilities::ParameterInfo(paramName, buffer.getComponentType(), buffer.getNumComponents(), buffer.getSize(), buffer.getMemory()));
parameters.push_back(OpenCLNonbondedUtilities::ParameterInfo(paramName, buffer.getComponentType(), buffer.getNumComponents(), buffer.getSize(), buffer.getMemory()));
}
for (int i = 0; i < (int) energyDerivs->getBuffers().size(); i++) {
const OpenCLNonbondedUtilities::ParameterInfo& buffer = energyDerivs->getBuffers()[i];
string paramName = prefix+"dEdV"+intToString(i+1);
cl.getNonbondedUtilities().addParameter(OpenCLNonbondedUtilities::ParameterInfo(paramName, buffer.getComponentType(), buffer.getNumComponents(), buffer.getSize(), buffer.getMemory()));
parameters.push_back(OpenCLNonbondedUtilities::ParameterInfo(paramName, buffer.getComponentType(), buffer.getNumComponents(), buffer.getSize(), buffer.getMemory()));
}
if (globals != NULL) {
globals->upload(globalParamValues);
cl.getNonbondedUtilities().addArgument(OpenCLNonbondedUtilities::ParameterInfo(prefix+"globals", "float", 1, sizeof(cl_float), globals->getDeviceBuffer()));
arguments.push_back(OpenCLNonbondedUtilities::ParameterInfo(prefix+"globals", "float", 1, sizeof(cl_float), globals->getDeviceBuffer()));
}
if (needParameterGradient) {
chainRuleParameters = parameters;
chainRuleArguments = arguments;
chainRuleSource = source;
separateChainRuleKernel = true;
cl.getNonbondedUtilities().addInteraction(useCutoff, usePeriodic, true, force.getCutoffDistance(), exclusionList, "");
}
else {
cl.getNonbondedUtilities().addInteraction(useCutoff, usePeriodic, true, force.getCutoffDistance(), exclusionList, source);
for (int i = 0; i < (int) parameters.size(); i++)
cl.getNonbondedUtilities().addParameter(parameters[i]);
for (int i = 0; i < (int) arguments.size(); i++)
cl.getNonbondedUtilities().addArgument(arguments[i]);
separateChainRuleKernel = false;
}
}
cl.addForce(new OpenCLCustomGBForceInfo(cl.getNonbondedUtilities().getNumForceBuffers(), force));
......@@ -2209,6 +2311,7 @@ void OpenCLCalcCustomGBForceKernel::executeForces(ContextImpl& context) {
perParticleValueKernel.setArg<cl_int>(index++, cl.getPaddedNumAtoms());
perParticleValueKernel.setArg<cl_int>(index++, nb.getNumForceBuffers());
perParticleValueKernel.setArg<cl::Buffer>(index++, valueBuffers->getDeviceBuffer());
perParticleValueKernel.setArg<cl::Buffer>(index++, cl.getPosq().getDeviceBuffer());
if (globals != NULL)
perParticleValueKernel.setArg<cl::Buffer>(index++, globals->getDeviceBuffer());
for (int i = 0; i < (int) params->getBuffers().size(); i++)
......@@ -2263,7 +2366,9 @@ void OpenCLCalcCustomGBForceKernel::executeForces(ContextImpl& context) {
index = 0;
perParticleEnergyKernel.setArg<cl_int>(index++, cl.getPaddedNumAtoms());
perParticleEnergyKernel.setArg<cl_int>(index++, nb.getNumForceBuffers());
perParticleEnergyKernel.setArg<cl::Buffer>(index++, cl.getForceBuffers().getDeviceBuffer());
perParticleEnergyKernel.setArg<cl::Buffer>(index++, cl.getEnergyBuffer().getDeviceBuffer());
perParticleEnergyKernel.setArg<cl::Buffer>(index++, cl.getPosq().getDeviceBuffer());
if (globals != NULL)
perParticleEnergyKernel.setArg<cl::Buffer>(index++, globals->getDeviceBuffer());
for (int i = 0; i < (int) params->getBuffers().size(); i++)
......@@ -2277,6 +2382,7 @@ void OpenCLCalcCustomGBForceKernel::executeForces(ContextImpl& context) {
perParticleEnergyKernel.setArg<cl::Buffer>(index++, tabulatedFunctions[i]->getDeviceBuffer());
perParticleEnergyKernel.setArg<cl::Buffer>(index++, tabulatedFunctionParams->getDeviceBuffer());
}
chainRuleKernel = nb.createInteractionKernel(chainRuleSource, chainRuleParameters, chainRuleArguments, true, false);
}
if (globals != NULL) {
bool changed = false;
......@@ -2298,6 +2404,8 @@ void OpenCLCalcCustomGBForceKernel::executeForces(ContextImpl& context) {
cl.executeKernel(perParticleValueKernel, cl.getPaddedNumAtoms());
cl.executeKernel(pairEnergyKernel, nb.getTiles().getSize()*OpenCLContext::TileSize);
cl.executeKernel(perParticleEnergyKernel, cl.getPaddedNumAtoms());
if (separateChainRuleKernel)
cl.executeKernel(chainRuleKernel, nb.getTiles().getSize()*OpenCLContext::TileSize);
}
double OpenCLCalcCustomGBForceKernel::executeEnergy(ContextImpl& context) {
......
platforms/opencl/src/OpenCLKernels.h
View file @ 43e6e68c
......@@ -608,7 +608,7 @@ public:
*/
double executeEnergy(ContextImpl& context);
private:
bool hasInitializedKernels;
bool hasInitializedKernels, separateChainRuleKernel;
OpenCLContext& cl;
OpenCLParameterSet* params;
OpenCLParameterSet* computedValues;
......@@ -619,8 +619,11 @@ private:
std::vector<std::string> globalParamNames;
std::vector<cl_float> globalParamValues;
std::vector<OpenCLArray<mm_float4>*> tabulatedFunctions;
std::vector<OpenCLNonbondedUtilities::ParameterInfo> chainRuleParameters;
std::vector<OpenCLNonbondedUtilities::ParameterInfo> chainRuleArguments;
std::string chainRuleSource;
System& system;
cl::Kernel pairValueKernel, perParticleValueKernel, pairEnergyKernel, perParticleEnergyKernel;
cl::Kernel pairValueKernel, perParticleValueKernel, pairEnergyKernel, perParticleEnergyKernel, chainRuleKernel;
};
/**
......
platforms/opencl/src/OpenCLNonbondedUtilities.cpp
View file @ 43e6e68c
......@@ -226,7 +226,7 @@ void OpenCLNonbondedUtilities::initialize(const System& system) {
// Create kernels.
forceKernel = createInteractionKernel(kernelSource, parameters, arguments, true);
forceKernel = createInteractionKernel(kernelSource, parameters, arguments, true, true);
if (useCutoff) {
map<string, string> defines;
if (forceBufferPerAtomBlock)
......@@ -279,7 +279,7 @@ void OpenCLNonbondedUtilities::computeInteractions() {
context.executeKernel(forceKernel, tiles->getSize()*OpenCLContext::TileSize);
}
cl::Kernel OpenCLNonbondedUtilities::createInteractionKernel(const string& source, const vector<ParameterInfo>& params, const vector<ParameterInfo>& arguments, bool useExclusions) const {
cl::Kernel OpenCLNonbondedUtilities::createInteractionKernel(const string& source, const vector<ParameterInfo>& params, const vector<ParameterInfo>& arguments, bool useExclusions, bool isSymmetric) const {
map<string, string> replacements;
replacements["COMPUTE_INTERACTION"] = source;
int localDataSize = 7*sizeof(cl_float);
......@@ -380,6 +380,8 @@ cl::Kernel OpenCLNonbondedUtilities::createInteractionKernel(const string& sourc
defines["USE_PERIODIC"] = "1";
if (useExclusions)
defines["USE_EXCLUSIONS"] = "1";
if (isSymmetric)
defines["USE_SYMMETRIC"] = "1";
defines["PERIODIC_BOX_SIZE_X"] = OpenCLExpressionUtilities::doubleToString(periodicBoxSize.x);
defines["PERIODIC_BOX_SIZE_Y"] = OpenCLExpressionUtilities::doubleToString(periodicBoxSize.y);
defines["PERIODIC_BOX_SIZE_Z"] = OpenCLExpressionUtilities::doubleToString(periodicBoxSize.z);
......
platforms/opencl/src/OpenCLNonbondedUtilities.h
View file @ 43e6e68c
......@@ -192,8 +192,9 @@ public:
* @param params the per-atom parameters this kernel may depend on
* @param arguments arrays (other than per-atom parameters) that should be passed as arguments to the kernel
* @param useExclusions specifies whether exclusions are applied to this interaction
* @param isSymmetric specifies whether the interaction is symmetric
*/
cl::Kernel createInteractionKernel(const std::string& source, const std::vector<ParameterInfo>& params, const std::vector<ParameterInfo>& arguments, bool useExclusions) const;
cl::Kernel createInteractionKernel(const std::string& source, const std::vector<ParameterInfo>& params, const std::vector<ParameterInfo>& arguments, bool useExclusions, bool isSymmetric) const;
private:
OpenCLContext& context;
cl::Kernel forceKernel;
......
platforms/opencl/src/kernels/customGBChainRule.cl
View file @ 43e6e68c
......@@ -3,7 +3,17 @@ if (!isExcluded && atom1 < NUM_ATOMS && atom2 < NUM_ATOMS && atom1 != atom2 && r
#else
if (!isExcluded && atom1 < NUM_ATOMS && atom2 < NUM_ATOMS && atom1 != atom2) {
#endif
#ifdef USE_SYMMETRIC
float tempForce = 0.0f;
#else
float4 tempForce1 = (float4) 0.0f;
float4 tempForce2 = (float4) 0.0f;
#endif
COMPUTE_FORCE
#ifdef USE_SYMMETRIC
dEdR += tempForce*invR;
#else
dEdR1 += tempForce1;
dEdR2 += tempForce2;
#endif
}
platforms/opencl/src/kernels/customGBEnergyPerParticle.cl
View file @ 43e6e68c
......@@ -8,7 +8,7 @@
* Reduce the derivatives computed in the N^2 energy kernel, and compute all per-particle energy terms.
*/
__kernel void computePerParticleEnergy(int bufferSize, int numBuffers, __global float* energyBuffer
__kernel void computePerParticleEnergy(int bufferSize, int numBuffers, __global float4* forceBuffers, __global float* energyBuffer, __global float4* posq
PARAMETER_ARGUMENTS) {
float energy = 0.0f;
unsigned int index = get_global_id(0);
......@@ -20,6 +20,8 @@ __kernel void computePerParticleEnergy(int bufferSize, int numBuffers, __global
// Now calculate the per-particle energy terms.
float4 pos = posq[index];
float4 force = (float4) 0.0f;
COMPUTE_ENERGY
index += get_global_size(0);
}
......
platforms/opencl/src/kernels/customGBValuePerParticle.cl
View file @ 43e6e68c
......@@ -2,7 +2,7 @@
* Reduce a pairwise computed value, and compute per-particle values.
*/
__kernel void computePerParticleValues(int bufferSize, int numBuffers, __global float* valueBuffers
__kernel void computePerParticleValues(int bufferSize, int numBuffers, __global float* valueBuffers, __global float4* posq
PARAMETER_ARGUMENTS) {
unsigned int index = get_global_id(0);
while (index < NUM_ATOMS) {
......@@ -15,6 +15,7 @@ __kernel void computePerParticleValues(int bufferSize, int numBuffers, __global
// Now calculate other values
float4 pos = posq[index];
COMPUTE_VALUES
index += get_global_size(0);
}
......
platforms/opencl/src/kernels/nonbonded_default.cl
View file @ 43e6e68c
......@@ -72,12 +72,20 @@ void computeNonbonded(__global float4* forceBuffers, __global float* energyBuffe
float r = RECIP(invR);
LOAD_ATOM2_PARAMETERS
atom2 = y+baseLocalAtom+j;
#ifdef USE_SYMMETRIC
float dEdR = 0.0f;
#else
float4 dEdR1 = (float4) 0.0f;
float4 dEdR2 = (float4) 0.0f;
#endif
float tempEnergy = 0.0f;
COMPUTE_INTERACTION
energy += 0.5f*tempEnergy;
delta.xyz *= dEdR;
force.xyz -= delta.xyz;
#ifdef USE_SYMMETRIC
force.xyz -= delta.xyz*dEdR;
#else
force.xyz -= dEdR1.xyz;
#endif
excl >>= 1;
}
......@@ -140,15 +148,27 @@ void computeNonbonded(__global float4* forceBuffers, __global float* energyBuffe
float r = RECIP(invR);
LOAD_ATOM2_PARAMETERS
atom2 = y+baseLocalAtom+tj;
#ifdef USE_SYMMETRIC
float dEdR = 0.0f;
#else
float4 dEdR1 = (float4) 0.0f;
float4 dEdR2 = (float4) 0.0f;
#endif
float tempEnergy = 0.0f;
COMPUTE_INTERACTION
energy += tempEnergy;
#ifdef USE_SYMMETRIC
delta.xyz *= dEdR;
force.xyz -= delta.xyz;
localData[baseLocalAtom+tj+forceBufferOffset].fx += delta.x;
localData[baseLocalAtom+tj+forceBufferOffset].fy += delta.y;
localData[baseLocalAtom+tj+forceBufferOffset].fz += delta.z;
#else
force.xyz -= dEdR1.xyz;
localData[baseLocalAtom+tj+forceBufferOffset].fx += dEdR2.x;
localData[baseLocalAtom+tj+forceBufferOffset].fy += dEdR2.y;
localData[baseLocalAtom+tj+forceBufferOffset].fz += dEdR2.z;
#endif
barrier(CLK_LOCAL_MEM_FENCE);
excl >>= 1;
tj = (tj+1)%(TILE_SIZE/2);
......
platforms/opencl/src/kernels/nonbonded_nvidia.cl
View file @ 43e6e68c
......@@ -72,12 +72,20 @@ void computeNonbonded(__global float4* forceBuffers, __global float* energyBuffe
float invR = 1.0f/r;
LOAD_ATOM2_PARAMETERS
atom2 = y+j;
#ifdef USE_SYMMETRIC
float dEdR = 0.0f;
#else
float4 dEdR1 = (float4) 0.0f;
float4 dEdR2 = (float4) 0.0f;
#endif
float tempEnergy = 0.0f;
COMPUTE_INTERACTION
energy += 0.5f*tempEnergy;
delta.xyz *= dEdR;
force.xyz -= delta.xyz;
#ifdef USE_SYMMETRIC
force.xyz -= delta.xyz*dEdR;
#else
force.xyz -= dEdR1.xyz;
#endif
excl >>= 1;
}
......@@ -129,13 +137,23 @@ void computeNonbonded(__global float4* forceBuffers, __global float* energyBuffe
float r = RECIP(invR);
LOAD_ATOM2_PARAMETERS
atom2 = y+j;
#ifdef USE_SYMMETRIC
float dEdR = 0.0f;
#else
float4 dEdR1 = (float4) 0.0f;
float4 dEdR2 = (float4) 0.0f;
#endif
float tempEnergy = 0.0f;
COMPUTE_INTERACTION
energy += tempEnergy;
#ifdef USE_SYMMETRIC
delta.xyz *= dEdR;
force.xyz -= delta.xyz;
tempBuffer[get_local_id(0)] = delta;
#else
force.xyz -= dEdR1.xyz;
tempBuffer[get_local_id(0)] = dEdR2;
#endif
// Sum the forces on atom2.
......@@ -186,15 +204,27 @@ void computeNonbonded(__global float4* forceBuffers, __global float* energyBuffe
float r = RECIP(invR);
LOAD_ATOM2_PARAMETERS
atom2 = y+tj;
#ifdef USE_SYMMETRIC
float dEdR = 0.0f;
#else
float4 dEdR1 = (float4) 0.0f;
float4 dEdR2 = (float4) 0.0f;
#endif
float tempEnergy = 0.0f;
COMPUTE_INTERACTION
energy += tempEnergy;
#ifdef USE_SYMMETRIC
delta.xyz *= dEdR;
force.xyz -= delta.xyz;
localData[tbx+tj].fx += delta.x;
localData[tbx+tj].fy += delta.y;
localData[tbx+tj].fz += delta.z;
#else
force.xyz -= dEdR1.xyz;
localData[tbx+tj].fx += dEdR2.x;
localData[tbx+tj].fy += dEdR2.y;
localData[tbx+tj].fz += dEdR2.z;
#endif
excl >>= 1;
tj = (tj + 1) & (TILE_SIZE - 1);
}
......
platforms/opencl/tests/TestOpenCLCustomGBForce.cpp
View file @ 43e6e68c
......@@ -200,6 +200,63 @@ void testMultipleChainRules() {
}
}
void testPositionDependence() {
OpenCLPlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomGBForce* force = new CustomGBForce();
force->addComputedValue("a", "r", CustomGBForce::ParticlePair);
force->addComputedValue("b", "a+y", CustomGBForce::SingleParticle);
force->addEnergyTerm("b*z", CustomGBForce::SingleParticle);
force->addParticle(vector<double>());
force->addParticle(vector<double>());
system.addForce(force);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
vector<Vec3> forces(2);
init_gen_rand(0);
for (int i = 0; i < 5; i++) {
positions[0] = Vec3(genrand_real2(), genrand_real2(), genrand_real2());
positions[1] = Vec3(genrand_real2(), genrand_real2(), genrand_real2());
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
Vec3 delta = positions[0]-positions[1];
double r = sqrt(delta.dot(delta));
double energy = 0;
for (int j = 0; j < 2; j++)
energy += positions[j][2]*(r+positions[j][1]);
Vec3 force1(-positions[0][2]*delta[0]/r-positions[1][2]*delta[0]/r,
-positions[0][2]*(delta[1]/r+1)-positions[1][2]*delta[1]/r,
-positions[0][2]*delta[2]/r-(r+positions[0][1])-positions[1][2]*delta[2]/r);
Vec3 force2(positions[0][2]*delta[0]/r+positions[1][2]*delta[0]/r,
positions[0][2]*delta[1]/r+positions[1][2]*(delta[1]/r-1),
positions[0][2]*delta[2]/r+positions[1][2]*delta[2]/r-(r+positions[1][1]));
ASSERT_EQUAL_VEC(force1, forces[0], 1e-4);
ASSERT_EQUAL_VEC(force2, forces[1], 1e-4);
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), 0.02);
// Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.
double norm = 0.0;
for (int i = 0; i < (int) forces.size(); ++i)
norm += forces[i].dot(forces[i]);
norm = std::sqrt(norm);
const double stepSize = 1e-3;
double step = stepSize/norm;
for (int i = 0; i < (int) positions.size(); ++i) {
Vec3 p = positions[i];
Vec3 f = forces[i];
positions[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
}
context.setPositions(positions);
State state2 = context.getState(State::Energy);
ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state.getPotentialEnergy())/stepSize, 1e-3*abs(state.getPotentialEnergy()));
}
}
int main() {
try {
testOBC(GBSAOBCForce::NoCutoff, CustomGBForce::NoCutoff);
......@@ -208,6 +265,7 @@ int main() {
testTabulatedFunction(true);
testTabulatedFunction(false);
testMultipleChainRules();
testPositionDependence();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
......
platforms/reference/src/ReferenceKernels.cpp
View file @ 43e6e68c
......@@ -1094,6 +1094,7 @@ void ReferenceCalcCustomGBForceKernel::initialize(const System& system, const Cu
// Parse the expressions for computed values.
valueDerivExpressions.resize(force.getNumComputedValues());
valueGradientExpressions.resize(force.getNumComputedValues());
for (int i = 0; i < force.getNumComputedValues(); i++) {
string name, expression;
CustomGBForce::ComputationType type;
......@@ -1105,6 +1106,9 @@ void ReferenceCalcCustomGBForceKernel::initialize(const System& system, const Cu
if (i == 0)
valueDerivExpressions[i].push_back(ex.differentiate("r").optimize().createProgram());
else {
valueGradientExpressions[i].push_back(ex.differentiate("x").optimize().createProgram());
valueGradientExpressions[i].push_back(ex.differentiate("y").optimize().createProgram());
valueGradientExpressions[i].push_back(ex.differentiate("z").optimize().createProgram());
for (int j = 0; j < i; j++)
valueDerivExpressions[i].push_back(ex.differentiate(valueNames[j]).optimize().createProgram());
}
......@@ -1113,6 +1117,7 @@ void ReferenceCalcCustomGBForceKernel::initialize(const System& system, const Cu
// Parse the expressions for energy terms.
energyDerivExpressions.resize(force.getNumEnergyTerms());
energyGradientExpressions.resize(force.getNumEnergyTerms());
for (int i = 0; i < force.getNumEnergyTerms(); i++) {
string expression;
CustomGBForce::ComputationType type;
......@@ -1123,8 +1128,12 @@ void ReferenceCalcCustomGBForceKernel::initialize(const System& system, const Cu
if (type != CustomGBForce::SingleParticle)
energyDerivExpressions[i].push_back(ex.differentiate("r").optimize().createProgram());
for (int j = 0; j < force.getNumComputedValues(); j++) {
if (type == CustomGBForce::SingleParticle)
if (type == CustomGBForce::SingleParticle) {
energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]).optimize().createProgram());
energyGradientExpressions[i].push_back(ex.differentiate("x").optimize().createProgram());
energyGradientExpressions[i].push_back(ex.differentiate("y").optimize().createProgram());
energyGradientExpressions[i].push_back(ex.differentiate("z").optimize().createProgram());
}
else {
energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]+"1").optimize().createProgram());
energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]+"2").optimize().createProgram());
......@@ -1141,8 +1150,8 @@ void ReferenceCalcCustomGBForceKernel::initialize(const System& system, const Cu
void ReferenceCalcCustomGBForceKernel::executeForces(ContextImpl& context) {
RealOpenMM** posData = extractPositions(context);
RealOpenMM** forceData = extractForces(context);
ReferenceCustomGBIxn ixn(valueExpressions, valueDerivExpressions, valueNames, valueTypes, energyExpressions,
energyDerivExpressions, energyTypes, particleParameterNames);
ReferenceCustomGBIxn ixn(valueExpressions, valueDerivExpressions, valueGradientExpressions, valueNames, valueTypes, energyExpressions,
energyDerivExpressions, energyGradientExpressions, energyTypes, particleParameterNames);
bool periodic = (nonbondedMethod == CutoffPeriodic);
if (nonbondedMethod != NoCutoff) {
computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, periodic ? periodicBoxSize : NULL, nonbondedCutoff, 0.0);
......@@ -1160,8 +1169,8 @@ double ReferenceCalcCustomGBForceKernel::executeEnergy(ContextImpl& context) {
RealOpenMM** posData = extractPositions(context);
RealOpenMM** forceData = allocateRealArray(numParticles, 3);
RealOpenMM energy = 0;
ReferenceCustomGBIxn ixn(valueExpressions, valueDerivExpressions, valueNames, valueTypes, energyExpressions,
energyDerivExpressions, energyTypes, particleParameterNames);
ReferenceCustomGBIxn ixn(valueExpressions, valueDerivExpressions, valueGradientExpressions, valueNames, valueTypes, energyExpressions,
energyDerivExpressions, energyGradientExpressions, energyTypes, particleParameterNames);
bool periodic = (nonbondedMethod == CutoffPeriodic);
if (nonbondedMethod != NoCutoff) {
computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, periodic ? periodicBoxSize : NULL, nonbondedCutoff, 0.0);
......
platforms/reference/src/ReferenceKernels.h
View file @ 43e6e68c
......@@ -585,9 +585,11 @@ private:
std::vector<std::string> particleParameterNames, globalParameterNames, valueNames;
std::vector<Lepton::ExpressionProgram> valueExpressions;
std::vector<std::vector<Lepton::ExpressionProgram> > valueDerivExpressions;
std::vector<std::vector<Lepton::ExpressionProgram> > valueGradientExpressions;
std::vector<OpenMM::CustomGBForce::ComputationType> valueTypes;
std::vector<Lepton::ExpressionProgram> energyExpressions;
std::vector<std::vector<Lepton::ExpressionProgram> > energyDerivExpressions;
std::vector<std::vector<Lepton::ExpressionProgram> > energyGradientExpressions;
std::vector<OpenMM::CustomGBForce::ComputationType> energyTypes;
NonbondedMethod nonbondedMethod;
NeighborList* neighborList;
......
platforms/reference/src/SimTKReference/ReferenceCustomGBIxn.cpp
View file @ 43e6e68c
......@@ -45,14 +45,16 @@ using std::vector;
ReferenceCustomGBIxn::ReferenceCustomGBIxn(const vector<Lepton::ExpressionProgram>& valueExpressions,
const vector<vector<Lepton::ExpressionProgram> > valueDerivExpressions,
const vector<vector<Lepton::ExpressionProgram> > valueGradientExpressions,
const vector<string>& valueNames,
const vector<OpenMM::CustomGBForce::ComputationType>& valueTypes,
const vector<Lepton::ExpressionProgram>& energyExpressions,
const vector<vector<Lepton::ExpressionProgram> > energyDerivExpressions,
const vector<vector<Lepton::ExpressionProgram> > energyGradientExpressions,
const vector<OpenMM::CustomGBForce::ComputationType>& energyTypes,
const vector<string>& parameterNames) :
cutoff(false), periodic(false), valueExpressions(valueExpressions), valueDerivExpressions(valueDerivExpressions), valueNames(valueNames),
valueTypes(valueTypes), energyExpressions(energyExpressions), energyDerivExpressions(energyDerivExpressions),
cutoff(false), periodic(false), valueExpressions(valueExpressions), valueDerivExpressions(valueDerivExpressions), valueGradientExpressions(valueGradientExpressions),
valueNames(valueNames), valueTypes(valueTypes), energyExpressions(energyExpressions), energyDerivExpressions(energyDerivExpressions), energyGradientExpressions(energyGradientExpressions),
energyTypes(energyTypes), paramNames(parameterNames) {
// ---------------------------------------------------------------------------------------
......@@ -140,7 +142,7 @@ void ReferenceCustomGBIxn::calculateIxn(int numberOfAtoms, RealOpenMM** atomCoor
vector<vector<RealOpenMM> > values(numValues);
for (int valueIndex = 0; valueIndex < numValues; valueIndex++) {
if (valueTypes[valueIndex] == OpenMM::CustomGBForce::SingleParticle)
calculateSingleParticleValue(valueIndex, numberOfAtoms, values, globalParameters, atomParameters);
calculateSingleParticleValue(valueIndex, numberOfAtoms, atomCoordinates, values, globalParameters, atomParameters);
else if (valueTypes[valueIndex] == OpenMM::CustomGBForce::ParticlePair)
calculateParticlePairValue(valueIndex, numberOfAtoms, atomCoordinates, atomParameters, values, globalParameters, exclusions, true);
else
......@@ -152,7 +154,7 @@ void ReferenceCustomGBIxn::calculateIxn(int numberOfAtoms, RealOpenMM** atomCoor
vector<vector<RealOpenMM> > dEdV(numValues, vector<RealOpenMM>(numberOfAtoms, (RealOpenMM) 0));
for (int termIndex = 0; termIndex < (int) energyExpressions.size(); termIndex++) {
if (energyTypes[termIndex] == OpenMM::CustomGBForce::SingleParticle)
calculateSingleParticleEnergyTerm(termIndex, numberOfAtoms, values, globalParameters, atomParameters, totalEnergy, dEdV);
calculateSingleParticleEnergyTerm(termIndex, numberOfAtoms, atomCoordinates, values, globalParameters, atomParameters, forces, totalEnergy, dEdV);
else if (energyTypes[termIndex] == OpenMM::CustomGBForce::ParticlePair)
calculateParticlePairEnergyTerm(termIndex, numberOfAtoms, atomCoordinates, atomParameters, values, globalParameters, exclusions, true, forces, totalEnergy, dEdV);
else
......@@ -164,11 +166,14 @@ void ReferenceCustomGBIxn::calculateIxn(int numberOfAtoms, RealOpenMM** atomCoor
calculateChainRuleForces(numberOfAtoms, atomCoordinates, atomParameters, values, globalParameters, exclusions, forces, dEdV);
}
void ReferenceCustomGBIxn::calculateSingleParticleValue(int index, int numAtoms, vector<vector<RealOpenMM> >& values,
void ReferenceCustomGBIxn::calculateSingleParticleValue(int index, int numAtoms, RealOpenMM** atomCoordinates, vector<vector<RealOpenMM> >& values,
const map<string, double>& globalParameters, RealOpenMM** atomParameters) const {
values[index].resize(numAtoms);
map<string, double> variables = globalParameters;
for (int i = 0; i < numAtoms; i++) {
variables["x"] = atomCoordinates[i][0];
variables["y"] = atomCoordinates[i][1];
variables["z"] = atomCoordinates[i][2];
for (int j = 0; j < (int) paramNames.size(); j++)
variables[paramNames[j]] = atomParameters[i][j];
for (int j = 0; j < index; j++)
......@@ -230,11 +235,14 @@ void ReferenceCustomGBIxn::calculateOnePairValue(int index, int atom1, int atom2
values[index][atom1] += (RealOpenMM) valueExpressions[index].evaluate(variables);
}
void ReferenceCustomGBIxn::calculateSingleParticleEnergyTerm(int index, int numAtoms, const vector<vector<RealOpenMM> >& values,
const map<string, double>& globalParameters, RealOpenMM** atomParameters, RealOpenMM* totalEnergy,
void ReferenceCustomGBIxn::calculateSingleParticleEnergyTerm(int index, int numAtoms, RealOpenMM** atomCoordinates, const vector<vector<RealOpenMM> >& values,
const map<string, double>& globalParameters, RealOpenMM** atomParameters, RealOpenMM** forces, RealOpenMM* totalEnergy,
vector<vector<RealOpenMM> >& dEdV) const {
map<string, double> variables = globalParameters;
for (int i = 0; i < numAtoms; i++) {
variables["x"] = atomCoordinates[i][0];
variables["y"] = atomCoordinates[i][1];
variables["z"] = atomCoordinates[i][2];
for (int j = 0; j < (int) paramNames.size(); j++)
variables[paramNames[j]] = atomParameters[i][j];
for (int j = 0; j < (int) valueNames.size(); j++)
......@@ -243,6 +251,9 @@ void ReferenceCustomGBIxn::calculateSingleParticleEnergyTerm(int index, int numA
*totalEnergy += (RealOpenMM) energyExpressions[index].evaluate(variables);
for (int j = 0; j < (int) valueNames.size(); j++)
dEdV[j][i] += (RealOpenMM) energyDerivExpressions[index][j].evaluate(variables);
forces[i][0] -= (RealOpenMM) energyGradientExpressions[index][0].evaluate(variables);
forces[i][1] -= (RealOpenMM) energyGradientExpressions[index][1].evaluate(variables);
forces[i][2] -= (RealOpenMM) energyGradientExpressions[index][2].evaluate(variables);
}
}
......@@ -371,13 +382,15 @@ void ReferenceCustomGBIxn::calculateOnePairChainRule(int atom1, int atom2, RealO
// Evaluate the derivative of each parameter with respect to position and apply forces.
vector<RealOpenMM> dVdR(valueDerivExpressions.size(), 0.0);
dVdR[0] = (RealOpenMM) valueDerivExpressions[0][0].evaluate(variables);
vector<vector<RealOpenMM> > gradient1(valueDerivExpressions.size(), vector<RealOpenMM>(3, 0.0));
vector<vector<RealOpenMM> > gradient2(valueDerivExpressions.size(), vector<RealOpenMM>(3, 0.0));
RealOpenMM dVdR = (RealOpenMM) valueDerivExpressions[0][0].evaluate(variables);
RealOpenMM rinv = 1/r;
for (int i = 0; i < 3; i++) {
RealOpenMM f = dEdV[0][atom1]*dVdR[0]*deltaR[i]*rinv;
forces[atom1][i] -= f;
forces[atom2][i] += f;
gradient1[0][i] = dVdR*deltaR[i]*rinv;
gradient2[0][i] = -gradient1[0][i];
forces[atom1][i] -= dEdV[0][atom1]*gradient1[0][i];
forces[atom2][i] -= dEdV[0][atom1]*gradient2[0][i];
}
variables = globalParameters;
for (int i = 0; i < (int) paramNames.size(); i++)
......@@ -385,12 +398,17 @@ void ReferenceCustomGBIxn::calculateOnePairChainRule(int atom1, int atom2, RealO
variables[valueNames[0]] = values[0][atom1];
for (int i = 1; i < (int) valueNames.size(); i++) {
variables[valueNames[i]] = values[i][atom1];
for (int j = 0; j < i; j++)
dVdR[i] += (RealOpenMM) (valueDerivExpressions[i][j].evaluate(variables)*dVdR[j]);
for (int j = 0; j < 3; j++) {
RealOpenMM f = dEdV[i][atom1]*dVdR[i]*deltaR[j]*rinv;
forces[atom1][j] -= f;
forces[atom2][j] += f;
for (int j = 0; j < i; j++) {
RealOpenMM dVdV = (RealOpenMM) valueDerivExpressions[i][j].evaluate(variables);
for (int k = 0; k < 3; k++) {
gradient1[i][k] += dVdV*gradient1[j][k];
gradient2[i][k] += dVdV*gradient2[j][k];
}
}
for (int k = 0; k < 3; k++) {
gradient1[i][k] += (RealOpenMM) valueGradientExpressions[i][k].evaluate(variables);
forces[atom1][k] -= dEdV[i][atom1]*gradient1[i][k];
forces[atom2][k] -= dEdV[i][atom1]*gradient2[i][k];
}
}
}
platforms/reference/src/SimTKReference/ReferenceCustomGBIxn.h
View file @ 43e6e68c
......@@ -45,10 +45,12 @@ class ReferenceCustomGBIxn {
RealOpenMM cutoffDistance;
std::vector<Lepton::ExpressionProgram> valueExpressions;
std::vector<std::vector<Lepton::ExpressionProgram> > valueDerivExpressions;
std::vector<std::vector<Lepton::ExpressionProgram> > valueGradientExpressions;
std::vector<std::string> valueNames;
std::vector<OpenMM::CustomGBForce::ComputationType> valueTypes;
std::vector<Lepton::ExpressionProgram> energyExpressions;
std::vector<std::vector<Lepton::ExpressionProgram> > energyDerivExpressions;
std::vector<std::vector<Lepton::ExpressionProgram> > energyGradientExpressions;
std::vector<std::string> paramNames;
std::vector<OpenMM::CustomGBForce::ComputationType> energyTypes;
std::vector<std::string> particleParamNames;
......@@ -60,13 +62,14 @@ class ReferenceCustomGBIxn {
@param index the index of the value to compute
@param numAtoms number of atoms
@param atomCoordinates atom coordinates
@param values the vector to store computed values into
@param globalParameters the values of global parameters
@param atomParameters atomParameters[atomIndex][paramterIndex]
--------------------------------------------------------------------------------------- */
void calculateSingleParticleValue(int index, int numAtoms, std::vector<std::vector<RealOpenMM> >& values,
void calculateSingleParticleValue(int index, int numAtoms, RealOpenMM** atomCoordinates, std::vector<std::vector<RealOpenMM> >& values,
const std::map<std::string, double>& globalParameters, RealOpenMM** atomParameters) const;
/**---------------------------------------------------------------------------------------
......@@ -113,16 +116,18 @@ class ReferenceCustomGBIxn {
@param index the index of the value to compute
@param numAtoms number of atoms
@param atomCoordinates atom coordinates
@param values the vector containing computed values
@param globalParameters the values of global parameters
@param atomParameters atomParameters[atomIndex][paramterIndex]
@param forces forces on atoms are added to this
@param totalEnergy the energy contribution is added to this
@param dEdV the derivative of energy with respect to computed values is stored in this
--------------------------------------------------------------------------------------- */
void calculateSingleParticleEnergyTerm(int index, int numAtoms, const std::vector<std::vector<RealOpenMM> >& values,
const std::map<std::string, double>& globalParameters, RealOpenMM** atomParameters,
void calculateSingleParticleEnergyTerm(int index, int numAtoms, RealOpenMM** atomCoordinates, const std::vector<std::vector<RealOpenMM> >& values,
const std::map<std::string, double>& globalParameters, RealOpenMM** atomParameters, RealOpenMM** forces,
RealOpenMM* totalEnergy, std::vector<std::vector<RealOpenMM> >& dEdV) const;
/**---------------------------------------------------------------------------------------
......@@ -222,10 +227,12 @@ class ReferenceCustomGBIxn {
ReferenceCustomGBIxn(const std::vector<Lepton::ExpressionProgram>& valueExpressions,
const std::vector<std::vector<Lepton::ExpressionProgram> > valueDerivExpressions,
const std::vector<std::vector<Lepton::ExpressionProgram> > valueGradientExpressions,
const std::vector<std::string>& valueNames,
const std::vector<OpenMM::CustomGBForce::ComputationType>& valueTypes,
const std::vector<Lepton::ExpressionProgram>& energyExpressions,
const std::vector<std::vector<Lepton::ExpressionProgram> > energyDerivExpressions,
const std::vector<std::vector<Lepton::ExpressionProgram> > energyGradientExpressions,
const std::vector<OpenMM::CustomGBForce::ComputationType>& energyTypes,
const std::vector<std::string>& parameterNames);
......
platforms/reference/tests/TestReferenceCustomGBForce.cpp
View file @ 43e6e68c
......@@ -200,6 +200,63 @@ void testMultipleChainRules() {
}
}
void testPositionDependence() {
ReferencePlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomGBForce* force = new CustomGBForce();
force->addComputedValue("a", "r", CustomGBForce::ParticlePair);
force->addComputedValue("b", "a+y", CustomGBForce::SingleParticle);
force->addEnergyTerm("b*z", CustomGBForce::SingleParticle);
force->addParticle(vector<double>());
force->addParticle(vector<double>());
system.addForce(force);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
vector<Vec3> forces(2);
init_gen_rand(0);
for (int i = 0; i < 5; i++) {
positions[0] = Vec3(genrand_real2(), genrand_real2(), genrand_real2());
positions[1] = Vec3(genrand_real2(), genrand_real2(), genrand_real2());
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
Vec3 delta = positions[0]-positions[1];
double r = sqrt(delta.dot(delta));
double energy = 0;
for (int j = 0; j < 2; j++)
energy += positions[j][2]*(r+positions[j][1]);
Vec3 force1(-positions[0][2]*delta[0]/r-positions[1][2]*delta[0]/r,
-positions[0][2]*(delta[1]/r+1)-positions[1][2]*delta[1]/r,
-positions[0][2]*delta[2]/r-(r+positions[0][1])-positions[1][2]*delta[2]/r);
Vec3 force2(positions[0][2]*delta[0]/r+positions[1][2]*delta[0]/r,
positions[0][2]*delta[1]/r+positions[1][2]*(delta[1]/r-1),
positions[0][2]*delta[2]/r+positions[1][2]*delta[2]/r-(r+positions[1][1]));
ASSERT_EQUAL_VEC(force1, forces[0], 1e-4);
ASSERT_EQUAL_VEC(force2, forces[1], 1e-4);
ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), 0.02);
// Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.
double norm = 0.0;
for (int i = 0; i < (int) forces.size(); ++i)
norm += forces[i].dot(forces[i]);
norm = std::sqrt(norm);
const double stepSize = 1e-3;
double step = stepSize/norm;
for (int i = 0; i < (int) positions.size(); ++i) {
Vec3 p = positions[i];
Vec3 f = forces[i];
positions[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
}
context.setPositions(positions);
State state2 = context.getState(State::Energy);
ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state.getPotentialEnergy())/stepSize, 1e-3*abs(state.getPotentialEnergy()));
}
}
int main() {
try {
testOBC(GBSAOBCForce::NoCutoff, CustomGBForce::NoCutoff);
......@@ -208,6 +265,7 @@ int main() {
testTabulatedFunction(true);
testTabulatedFunction(false);
testMultipleChainRules();
testPositionDependence();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
......
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