Skip to content

GitLab

Commit 0b1f03db authored by peastman's avatar peastman
Browse files

Use CompiledExpressionSet for reference CustomGBForce

parent f54a0079
Show whitespace changes
Inline Side-by-side
Showing with 183 additions and 192 deletions
platforms/reference/include/ReferenceCustomGBIxn.h
View file @ 0b1f03db
......@@ -26,7 +26,7 @@
#define __ReferenceCustomGBIxn_H__
#include "ReferenceNeighborList.h"
#include "lepton/ExpressionProgram.h"
#include "openmm/internal/CompiledExpressionSet.h"
#include "openmm/CustomGBForce.h"
#include <map>
#include <set>
......@@ -43,18 +43,20 @@ class ReferenceCustomGBIxn {
const OpenMM::NeighborList* neighborList;
OpenMM::RealVec periodicBoxVectors[3];
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;
CompiledExpressionSet expressionSet;
std::vector<Lepton::CompiledExpression> valueExpressions;
std::vector<std::vector<Lepton::CompiledExpression> > valueDerivExpressions;
std::vector<std::vector<Lepton::CompiledExpression> > 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<std::string> paramNames;
std::vector<Lepton::CompiledExpression> energyExpressions;
std::vector<std::vector<Lepton::CompiledExpression> > energyDerivExpressions;
std::vector<std::vector<Lepton::CompiledExpression> > energyGradientExpressions;
std::vector<OpenMM::CustomGBForce::ComputationType> energyTypes;
std::vector<std::string> particleParamNames;
std::vector<std::string> particleValueNames;
std::vector<int> paramIndex;
std::vector<int> valueIndex;
std::vector<int> particleParamIndex;
std::vector<int> particleValueIndex;
int rIndex, xIndex, yIndex, zIndex;
/**---------------------------------------------------------------------------------------
......@@ -64,13 +66,12 @@ class ReferenceCustomGBIxn {
@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<OpenMM::RealVec>& atomCoordinates, std::vector<std::vector<RealOpenMM> >& values,
const std::map<std::string, double>& globalParameters, RealOpenMM** atomParameters) const;
RealOpenMM** atomParameters);
/**---------------------------------------------------------------------------------------
......@@ -81,7 +82,6 @@ class ReferenceCustomGBIxn {
@param atomCoordinates atom coordinates
@param atomParameters atomParameters[atomIndex][paramterIndex]
@param values the vector to store computed values into
@param globalParameters the values of global parameters
@param exclusions exclusions[i] is the set of excluded indices for atom i
@param useExclusions specifies whether to use exclusions
......@@ -89,8 +89,7 @@ class ReferenceCustomGBIxn {
void calculateParticlePairValue(int index, int numAtoms, std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** atomParameters,
std::vector<std::vector<RealOpenMM> >& values,
const std::map<std::string, double>& globalParameters,
const std::vector<std::set<int> >& exclusions, bool useExclusions) const;
const std::vector<std::set<int> >& exclusions, bool useExclusions);
/**---------------------------------------------------------------------------------------
......@@ -101,14 +100,12 @@ class ReferenceCustomGBIxn {
@param atom2 the index of the second atom in the pair
@param atomCoordinates atom coordinates
@param atomParameters atomParameters[atomIndex][paramterIndex]
@param globalParameters the values of global parameters
@param values the vector to store computed values into
--------------------------------------------------------------------------------------- */
void calculateOnePairValue(int index, int atom1, int atom2, std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** atomParameters,
const std::map<std::string, double>& globalParameters,
std::vector<std::vector<RealOpenMM> >& values) const;
std::vector<std::vector<RealOpenMM> >& values);
/**---------------------------------------------------------------------------------------
......@@ -118,7 +115,6 @@ class ReferenceCustomGBIxn {
@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
......@@ -127,8 +123,8 @@ class ReferenceCustomGBIxn {
--------------------------------------------------------------------------------------- */
void calculateSingleParticleEnergyTerm(int index, int numAtoms, std::vector<OpenMM::RealVec>& atomCoordinates, const std::vector<std::vector<RealOpenMM> >& values,
const std::map<std::string, double>& globalParameters, RealOpenMM** atomParameters, std::vector<OpenMM::RealVec>& forces,
RealOpenMM* totalEnergy, std::vector<std::vector<RealOpenMM> >& dEdV) const;
RealOpenMM** atomParameters, std::vector<OpenMM::RealVec>& forces,
RealOpenMM* totalEnergy, std::vector<std::vector<RealOpenMM> >& dEdV);
/**---------------------------------------------------------------------------------------
......@@ -139,7 +135,6 @@ class ReferenceCustomGBIxn {
@param atomCoordinates atom coordinates
@param atomParameters atomParameters[atomIndex][paramterIndex]
@param values the vector containing computed values
@param globalParameters the values of global parameters
@param exclusions exclusions[i] is the set of excluded indices for atom i
@param useExclusions specifies whether to use exclusions
@param forces forces on atoms are added to this
......@@ -150,9 +145,8 @@ class ReferenceCustomGBIxn {
void calculateParticlePairEnergyTerm(int index, int numAtoms, std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** atomParameters,
const std::vector<std::vector<RealOpenMM> >& values,
const std::map<std::string, double>& globalParameters,
const std::vector<std::set<int> >& exclusions, bool useExclusions,
std::vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy, std::vector<std::vector<RealOpenMM> >& dEdV) const;
std::vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy, std::vector<std::vector<RealOpenMM> >& dEdV);
/**---------------------------------------------------------------------------------------
......@@ -163,7 +157,6 @@ class ReferenceCustomGBIxn {
@param atom2 the index of the second atom in the pair
@param atomCoordinates atom coordinates
@param atomParameters atomParameters[atomIndex][paramterIndex]
@param globalParameters the values of global parameters
@param values the vector containing computed values
@param forces forces on atoms are added to this
@param totalEnergy the energy contribution is added to this
......@@ -172,9 +165,8 @@ class ReferenceCustomGBIxn {
--------------------------------------------------------------------------------------- */
void calculateOnePairEnergyTerm(int index, int atom1, int atom2, std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** atomParameters,
const std::map<std::string, double>& globalParameters,
const std::vector<std::vector<RealOpenMM> >& values,
std::vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy, std::vector<std::vector<RealOpenMM> >& dEdV) const;
std::vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy, std::vector<std::vector<RealOpenMM> >& dEdV);
/**---------------------------------------------------------------------------------------
......@@ -184,7 +176,6 @@ class ReferenceCustomGBIxn {
@param atomCoordinates atom coordinates
@param atomParameters atomParameters[atomIndex][paramterIndex]
@param values the vector containing computed values
@param globalParameters the values of global parameters
@param exclusions exclusions[i] is the set of excluded indices for atom i
@param forces forces on atoms are added to this
@param dEdV the derivative of energy with respect to computed values is stored in this
......@@ -193,9 +184,8 @@ class ReferenceCustomGBIxn {
void calculateChainRuleForces(int numAtoms, std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** atomParameters,
const std::vector<std::vector<RealOpenMM> >& values,
const std::map<std::string, double>& globalParameters,
const std::vector<std::set<int> >& exclusions,
std::vector<OpenMM::RealVec>& forces, std::vector<std::vector<RealOpenMM> >& dEdV) const;
std::vector<OpenMM::RealVec>& forces, std::vector<std::vector<RealOpenMM> >& dEdV);
/**---------------------------------------------------------------------------------------
......@@ -205,7 +195,6 @@ class ReferenceCustomGBIxn {
@param atom2 the index of the second atom in the pair
@param atomCoordinates atom coordinates
@param atomParameters atomParameters[atomIndex][paramterIndex]
@param globalParameters the values of global parameters
@param values the vector containing computed values
@param forces forces on atoms are added to this
@param dEdV the derivative of energy with respect to computed values is stored in this
......@@ -214,10 +203,9 @@ class ReferenceCustomGBIxn {
--------------------------------------------------------------------------------------- */
void calculateOnePairChainRule(int atom1, int atom2, std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** atomParameters,
const std::map<std::string, double>& globalParameters,
const std::vector<std::vector<RealOpenMM> >& values,
std::vector<OpenMM::RealVec>& forces, std::vector<std::vector<RealOpenMM> >& dEdV,
bool isExcluded) const;
bool isExcluded);
public:
......@@ -227,14 +215,14 @@ 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,
ReferenceCustomGBIxn(const std::vector<Lepton::CompiledExpression>& valueExpressions,
const std::vector<std::vector<Lepton::CompiledExpression> > valueDerivExpressions,
const std::vector<std::vector<Lepton::CompiledExpression> > 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<Lepton::CompiledExpression>& energyExpressions,
const std::vector<std::vector<Lepton::CompiledExpression> > energyDerivExpressions,
const std::vector<std::vector<Lepton::CompiledExpression> > energyGradientExpressions,
const std::vector<OpenMM::CustomGBForce::ComputationType>& energyTypes,
const std::vector<std::string>& parameterNames);
......@@ -284,7 +272,7 @@ class ReferenceCustomGBIxn {
--------------------------------------------------------------------------------------- */
void calculateIxn(int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** atomParameters, const std::vector<std::set<int> >& exclusions,
std::map<std::string, double>& globalParameters, std::vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy) const;
std::map<std::string, double>& globalParameters, std::vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy);
// ---------------------------------------------------------------------------------------
......
platforms/reference/include/ReferenceKernels.h
View file @ 0b1f03db
......@@ -38,7 +38,6 @@
#include "ReferenceNeighborList.h"
#include "lepton/CompiledExpression.h"
#include "lepton/CustomFunction.h"
#include "lepton/ExpressionProgram.h"
namespace OpenMM {
......@@ -739,13 +738,13 @@ private:
RealOpenMM nonbondedCutoff;
std::vector<std::set<int> > exclusions;
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<Lepton::CompiledExpression> valueExpressions;
std::vector<std::vector<Lepton::CompiledExpression> > valueDerivExpressions;
std::vector<std::vector<Lepton::CompiledExpression> > 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<Lepton::CompiledExpression> energyExpressions;
std::vector<std::vector<Lepton::CompiledExpression> > energyDerivExpressions;
std::vector<std::vector<Lepton::CompiledExpression> > energyGradientExpressions;
std::vector<OpenMM::CustomGBForce::ComputationType> energyTypes;
NonbondedMethod nonbondedMethod;
NeighborList* neighborList;
......
platforms/reference/src/ReferenceKernels.cpp
View file @ 0b1f03db
......@@ -1376,19 +1376,19 @@ void ReferenceCalcCustomGBForceKernel::initialize(const System& system, const Cu
CustomGBForce::ComputationType type;
force.getComputedValueParameters(i, name, expression, type);
Lepton::ParsedExpression ex = Lepton::Parser::parse(expression, functions).optimize();
valueExpressions.push_back(ex.createProgram());
valueExpressions.push_back(ex.createCompiledExpression());
valueTypes.push_back(type);
valueNames.push_back(name);
if (i == 0) {
valueDerivExpressions[i].push_back(ex.differentiate("r").optimize().createProgram());
valueDerivExpressions[i].push_back(ex.differentiate("r").optimize().createCompiledExpression());
validateVariables(ex.getRootNode(), pairVariables);
}
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());
valueGradientExpressions[i].push_back(ex.differentiate("x").optimize().createCompiledExpression());
valueGradientExpressions[i].push_back(ex.differentiate("y").optimize().createCompiledExpression());
valueGradientExpressions[i].push_back(ex.differentiate("z").optimize().createCompiledExpression());
for (int j = 0; j < i; j++)
valueDerivExpressions[i].push_back(ex.differentiate(valueNames[j]).optimize().createProgram());
valueDerivExpressions[i].push_back(ex.differentiate(valueNames[j]).optimize().createCompiledExpression());
validateVariables(ex.getRootNode(), particleVariables);
}
particleVariables.insert(name);
......@@ -1405,21 +1405,21 @@ void ReferenceCalcCustomGBForceKernel::initialize(const System& system, const Cu
CustomGBForce::ComputationType type;
force.getEnergyTermParameters(i, expression, type);
Lepton::ParsedExpression ex = Lepton::Parser::parse(expression, functions).optimize();
energyExpressions.push_back(ex.createProgram());
energyExpressions.push_back(ex.createCompiledExpression());
energyTypes.push_back(type);
if (type != CustomGBForce::SingleParticle)
energyDerivExpressions[i].push_back(ex.differentiate("r").optimize().createProgram());
energyDerivExpressions[i].push_back(ex.differentiate("r").optimize().createCompiledExpression());
for (int j = 0; j < force.getNumComputedValues(); j++) {
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());
energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]).optimize().createCompiledExpression());
energyGradientExpressions[i].push_back(ex.differentiate("x").optimize().createCompiledExpression());
energyGradientExpressions[i].push_back(ex.differentiate("y").optimize().createCompiledExpression());
energyGradientExpressions[i].push_back(ex.differentiate("z").optimize().createCompiledExpression());
validateVariables(ex.getRootNode(), particleVariables);
}
else {
energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]+"1").optimize().createProgram());
energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]+"2").optimize().createProgram());
energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]+"1").optimize().createCompiledExpression());
energyDerivExpressions[i].push_back(ex.differentiate(valueNames[j]+"2").optimize().createCompiledExpression());
validateVariables(ex.getRootNode(), pairVariables);
}
}
......
platforms/reference/src/SimTKReference/ReferenceCustomGBIxn.cpp
View file @ 0b1f03db
......@@ -42,38 +42,54 @@ using namespace OpenMM;
--------------------------------------------------------------------------------------- */
ReferenceCustomGBIxn::ReferenceCustomGBIxn(const vector<Lepton::ExpressionProgram>& valueExpressions,
const vector<vector<Lepton::ExpressionProgram> > valueDerivExpressions,
const vector<vector<Lepton::ExpressionProgram> > valueGradientExpressions,
ReferenceCustomGBIxn::ReferenceCustomGBIxn(const vector<Lepton::CompiledExpression>& valueExpressions,
const vector<vector<Lepton::CompiledExpression> > valueDerivExpressions,
const vector<vector<Lepton::CompiledExpression> > 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<Lepton::CompiledExpression>& energyExpressions,
const vector<vector<Lepton::CompiledExpression> > energyDerivExpressions,
const vector<vector<Lepton::CompiledExpression> > energyGradientExpressions,
const vector<OpenMM::CustomGBForce::ComputationType>& energyTypes,
const vector<string>& parameterNames) :
cutoff(false), periodic(false), valueExpressions(valueExpressions), valueDerivExpressions(valueDerivExpressions), valueGradientExpressions(valueGradientExpressions),
valueNames(valueNames), valueTypes(valueTypes), energyExpressions(energyExpressions), energyDerivExpressions(energyDerivExpressions), energyGradientExpressions(energyGradientExpressions),
energyTypes(energyTypes), paramNames(parameterNames) {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceCustomGBIxn::ReferenceCustomGBIxn";
// ---------------------------------------------------------------------------------------
for (int i = 0; i < (int) paramNames.size(); i++) {
valueTypes(valueTypes), energyExpressions(energyExpressions), energyDerivExpressions(energyDerivExpressions), energyGradientExpressions(energyGradientExpressions),
energyTypes(energyTypes) {
for (int i = 0; i < this->valueExpressions.size(); i++)
expressionSet.registerExpression(this->valueExpressions[i]);
for (int i = 0; i < this->valueDerivExpressions.size(); i++)
for (int j = 0; j < this->valueDerivExpressions[i].size(); j++)
expressionSet.registerExpression(this->valueDerivExpressions[i][j]);
for (int i = 0; i < this->valueGradientExpressions.size(); i++)
for (int j = 0; j < this->valueGradientExpressions[i].size(); j++)
expressionSet.registerExpression(this->valueGradientExpressions[i][j]);
for (int i = 0; i < this->energyExpressions.size(); i++)
expressionSet.registerExpression(this->energyExpressions[i]);
for (int i = 0; i < this->energyDerivExpressions.size(); i++)
for (int j = 0; j < this->energyDerivExpressions[i].size(); j++)
expressionSet.registerExpression(this->energyDerivExpressions[i][j]);
for (int i = 0; i < this->energyGradientExpressions.size(); i++)
for (int j = 0; j < this->energyGradientExpressions[i].size(); j++)
expressionSet.registerExpression(this->energyGradientExpressions[i][j]);
rIndex = expressionSet.getVariableIndex("r");
xIndex = expressionSet.getVariableIndex("x");
yIndex = expressionSet.getVariableIndex("y");
zIndex = expressionSet.getVariableIndex("z");
for (int i = 0; i < (int) parameterNames.size(); i++) {
paramIndex.push_back(expressionSet.getVariableIndex(parameterNames[i]));
for (int j = 1; j < 3; j++) {
stringstream name;
name << paramNames[i] << j;
particleParamNames.push_back(name.str());
name << parameterNames[i] << j;
particleParamIndex.push_back(expressionSet.getVariableIndex(name.str()));
}
}
for (int i = 0; i < (int) valueNames.size(); i++) {
valueIndex.push_back(expressionSet.getVariableIndex(valueNames[i]));
for (int j = 1; j < 3; j++) {
stringstream name;
name << valueNames[i] << j;
particleValueNames.push_back(name.str());
particleValueIndex.push_back(expressionSet.getVariableIndex(name.str()));
}
}
}
......@@ -85,13 +101,6 @@ ReferenceCustomGBIxn::ReferenceCustomGBIxn(const vector<Lepton::ExpressionProgra
--------------------------------------------------------------------------------------- */
ReferenceCustomGBIxn::~ReferenceCustomGBIxn() {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceCustomGBIxn::~ReferenceCustomGBIxn";
// ---------------------------------------------------------------------------------------
}
/**---------------------------------------------------------------------------------------
......@@ -135,18 +144,21 @@ ReferenceCustomGBIxn::~ReferenceCustomGBIxn() {
void ReferenceCustomGBIxn::calculateIxn(int numberOfAtoms, vector<RealVec>& atomCoordinates, RealOpenMM** atomParameters,
const vector<set<int> >& exclusions, map<string, double>& globalParameters, vector<RealVec>& forces,
RealOpenMM* totalEnergy) const {
RealOpenMM* totalEnergy) {
for (map<string, double>::const_iterator iter = globalParameters.begin(); iter != globalParameters.end(); ++iter)
expressionSet.setVariable(expressionSet.getVariableIndex(iter->first), iter->second);
// First calculate the computed values.
int numValues = valueTypes.size();
vector<vector<RealOpenMM> > values(numValues);
for (int valueIndex = 0; valueIndex < numValues; valueIndex++) {
if (valueTypes[valueIndex] == OpenMM::CustomGBForce::SingleParticle)
calculateSingleParticleValue(valueIndex, numberOfAtoms, atomCoordinates, values, globalParameters, atomParameters);
calculateSingleParticleValue(valueIndex, numberOfAtoms, atomCoordinates, values, atomParameters);
else if (valueTypes[valueIndex] == OpenMM::CustomGBForce::ParticlePair)
calculateParticlePairValue(valueIndex, numberOfAtoms, atomCoordinates, atomParameters, values, globalParameters, exclusions, true);
calculateParticlePairValue(valueIndex, numberOfAtoms, atomCoordinates, atomParameters, values, exclusions, true);
else
calculateParticlePairValue(valueIndex, numberOfAtoms, atomCoordinates, atomParameters, values, globalParameters, exclusions, false);
calculateParticlePairValue(valueIndex, numberOfAtoms, atomCoordinates, atomParameters, values, exclusions, false);
}
// Now calculate the energy and its derivates.
......@@ -154,36 +166,35 @@ void ReferenceCustomGBIxn::calculateIxn(int numberOfAtoms, vector<RealVec>& atom
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, atomCoordinates, values, globalParameters, atomParameters, forces, totalEnergy, dEdV);
calculateSingleParticleEnergyTerm(termIndex, numberOfAtoms, atomCoordinates, values, atomParameters, forces, totalEnergy, dEdV);
else if (energyTypes[termIndex] == OpenMM::CustomGBForce::ParticlePair)
calculateParticlePairEnergyTerm(termIndex, numberOfAtoms, atomCoordinates, atomParameters, values, globalParameters, exclusions, true, forces, totalEnergy, dEdV);
calculateParticlePairEnergyTerm(termIndex, numberOfAtoms, atomCoordinates, atomParameters, values, exclusions, true, forces, totalEnergy, dEdV);
else
calculateParticlePairEnergyTerm(termIndex, numberOfAtoms, atomCoordinates, atomParameters, values, globalParameters, exclusions, false, forces, totalEnergy, dEdV);
calculateParticlePairEnergyTerm(termIndex, numberOfAtoms, atomCoordinates, atomParameters, values, exclusions, false, forces, totalEnergy, dEdV);
}
// Apply the chain rule to evaluate forces.
calculateChainRuleForces(numberOfAtoms, atomCoordinates, atomParameters, values, globalParameters, exclusions, forces, dEdV);
calculateChainRuleForces(numberOfAtoms, atomCoordinates, atomParameters, values, exclusions, forces, dEdV);
}
void ReferenceCustomGBIxn::calculateSingleParticleValue(int index, int numAtoms, vector<RealVec>& atomCoordinates, vector<vector<RealOpenMM> >& values,
const map<string, double>& globalParameters, RealOpenMM** atomParameters) const {
RealOpenMM** atomParameters) {
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];
expressionSet.setVariable(xIndex, atomCoordinates[i][0]);
expressionSet.setVariable(yIndex, atomCoordinates[i][1]);
expressionSet.setVariable(zIndex, atomCoordinates[i][2]);
for (int j = 0; j < (int) paramIndex.size(); j++)
expressionSet.setVariable(paramIndex[j], atomParameters[i][j]);
for (int j = 0; j < index; j++)
variables[valueNames[j]] = values[j][i];
values[index][i] = (RealOpenMM) valueExpressions[index].evaluate(variables);
expressionSet.setVariable(valueIndex[j], values[j][i]);
values[index][i] = (RealOpenMM) valueExpressions[index].evaluate();
}
}
void ReferenceCustomGBIxn::calculateParticlePairValue(int index, int numAtoms, vector<RealVec>& atomCoordinates, RealOpenMM** atomParameters,
vector<vector<RealOpenMM> >& values, const map<string, double>& globalParameters, const vector<set<int> >& exclusions, bool useExclusions) const {
vector<vector<RealOpenMM> >& values, const vector<set<int> >& exclusions, bool useExclusions) {
values[index].resize(numAtoms);
for (int i = 0; i < numAtoms; i++)
values[index][i] = (RealOpenMM) 0.0;
......@@ -194,8 +205,8 @@ void ReferenceCustomGBIxn::calculateParticlePairValue(int index, int numAtoms, v
OpenMM::AtomPair pair = (*neighborList)[i];
if (useExclusions && exclusions[pair.first].find(pair.second) != exclusions[pair.first].end())
continue;
calculateOnePairValue(index, pair.first, pair.second, atomCoordinates, atomParameters, globalParameters, values);
calculateOnePairValue(index, pair.second, pair.first, atomCoordinates, atomParameters, globalParameters, values);
calculateOnePairValue(index, pair.first, pair.second, atomCoordinates, atomParameters, values);
calculateOnePairValue(index, pair.second, pair.first, atomCoordinates, atomParameters, values);
}
}
else {
......@@ -205,15 +216,15 @@ void ReferenceCustomGBIxn::calculateParticlePairValue(int index, int numAtoms, v
for (int j = i+1; j < numAtoms; j++) {
if (useExclusions && exclusions[i].find(j) != exclusions[i].end())
continue;
calculateOnePairValue(index, i, j, atomCoordinates, atomParameters, globalParameters, values);
calculateOnePairValue(index, j, i, atomCoordinates, atomParameters, globalParameters, values);
calculateOnePairValue(index, i, j, atomCoordinates, atomParameters, values);
calculateOnePairValue(index, j, i, atomCoordinates, atomParameters, values);
}
}
}
}
void ReferenceCustomGBIxn::calculateOnePairValue(int index, int atom1, int atom2, vector<RealVec>& atomCoordinates, RealOpenMM** atomParameters,
const map<string, double>& globalParameters, vector<vector<RealOpenMM> >& values) const {
vector<vector<RealOpenMM> >& values) {
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (periodic)
ReferenceForce::getDeltaRPeriodic(atomCoordinates[atom2], atomCoordinates[atom1], periodicBoxVectors, deltaR);
......@@ -222,44 +233,42 @@ void ReferenceCustomGBIxn::calculateOnePairValue(int index, int atom1, int atom2
RealOpenMM r = deltaR[ReferenceForce::RIndex];
if (cutoff && r >= cutoffDistance)
return;
map<string, double> variables = globalParameters;
for (int i = 0; i < (int) paramNames.size(); i++) {
variables[particleParamNames[i*2]] = atomParameters[atom1][i];
variables[particleParamNames[i*2+1]] = atomParameters[atom2][i];
for (int i = 0; i < (int) paramIndex.size(); i++) {
expressionSet.setVariable(particleParamIndex[i*2], atomParameters[atom1][i]);
expressionSet.setVariable(particleParamIndex[i*2+1], atomParameters[atom2][i]);
}
variables["r"] = r;
expressionSet.setVariable(rIndex, r);
for (int i = 0; i < index; i++) {
variables[particleValueNames[i*2]] = values[i][atom1];
variables[particleValueNames[i*2+1]] = values[i][atom2];
expressionSet.setVariable(particleValueIndex[i*2], values[i][atom1]);
expressionSet.setVariable(particleValueIndex[i*2+1], values[i][atom2]);
}
values[index][atom1] += (RealOpenMM) valueExpressions[index].evaluate(variables);
values[index][atom1] += (RealOpenMM) valueExpressions[index].evaluate();
}
void ReferenceCustomGBIxn::calculateSingleParticleEnergyTerm(int index, int numAtoms, vector<RealVec>& atomCoordinates, const vector<vector<RealOpenMM> >& values,
const map<string, double>& globalParameters, RealOpenMM** atomParameters, vector<RealVec>& forces, RealOpenMM* totalEnergy,
vector<vector<RealOpenMM> >& dEdV) const {
map<string, double> variables = globalParameters;
RealOpenMM** atomParameters, vector<RealVec>& forces, RealOpenMM* totalEnergy,
vector<vector<RealOpenMM> >& dEdV) {
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++)
variables[valueNames[j]] = values[j][i];
expressionSet.setVariable(xIndex, atomCoordinates[i][0]);
expressionSet.setVariable(yIndex, atomCoordinates[i][1]);
expressionSet.setVariable(zIndex, atomCoordinates[i][2]);
for (int j = 0; j < (int) paramIndex.size(); j++)
expressionSet.setVariable(paramIndex[j], atomParameters[i][j]);
for (int j = 0; j < valueIndex.size(); j++)
expressionSet.setVariable(valueIndex[j], values[j][i]);
if (totalEnergy != NULL)
*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);
*totalEnergy += (RealOpenMM) energyExpressions[index].evaluate();
for (int j = 0; j < (int) valueIndex.size(); j++)
dEdV[j][i] += (RealOpenMM) energyDerivExpressions[index][j].evaluate();
forces[i][0] -= (RealOpenMM) energyGradientExpressions[index][0].evaluate();
forces[i][1] -= (RealOpenMM) energyGradientExpressions[index][1].evaluate();
forces[i][2] -= (RealOpenMM) energyGradientExpressions[index][2].evaluate();
}
}
void ReferenceCustomGBIxn::calculateParticlePairEnergyTerm(int index, int numAtoms, vector<RealVec>& atomCoordinates, RealOpenMM** atomParameters,
const vector<vector<RealOpenMM> >& values, const map<string, double>& globalParameters, const vector<set<int> >& exclusions, bool useExclusions,
vector<RealVec>& forces, RealOpenMM* totalEnergy, vector<vector<RealOpenMM> >& dEdV) const {
const vector<vector<RealOpenMM> >& values, const vector<set<int> >& exclusions, bool useExclusions,
vector<RealVec>& forces, RealOpenMM* totalEnergy, vector<vector<RealOpenMM> >& dEdV) {
if (cutoff) {
// Loop over all pairs in the neighbor list.
......@@ -267,7 +276,7 @@ void ReferenceCustomGBIxn::calculateParticlePairEnergyTerm(int index, int numAto
OpenMM::AtomPair pair = (*neighborList)[i];
if (useExclusions && exclusions[pair.first].find(pair.second) != exclusions[pair.first].end())
continue;
calculateOnePairEnergyTerm(index, pair.first, pair.second, atomCoordinates, atomParameters, globalParameters, values, forces, totalEnergy, dEdV);
calculateOnePairEnergyTerm(index, pair.first, pair.second, atomCoordinates, atomParameters, values, forces, totalEnergy, dEdV);
}
}
else {
......@@ -277,15 +286,15 @@ void ReferenceCustomGBIxn::calculateParticlePairEnergyTerm(int index, int numAto
for (int j = i+1; j < numAtoms; j++) {
if (useExclusions && exclusions[i].find(j) != exclusions[i].end())
continue;
calculateOnePairEnergyTerm(index, i, j, atomCoordinates, atomParameters, globalParameters, values, forces, totalEnergy, dEdV);
calculateOnePairEnergyTerm(index, i, j, atomCoordinates, atomParameters, values, forces, totalEnergy, dEdV);
}
}
}
}
void ReferenceCustomGBIxn::calculateOnePairEnergyTerm(int index, int atom1, int atom2, vector<RealVec>& atomCoordinates, RealOpenMM** atomParameters,
const map<string, double>& globalParameters, const vector<vector<RealOpenMM> >& values, vector<RealVec>& forces, RealOpenMM* totalEnergy,
vector<vector<RealOpenMM> >& dEdV) const {
const vector<vector<RealOpenMM> >& values, vector<RealVec>& forces, RealOpenMM* totalEnergy,
vector<vector<RealOpenMM> >& dEdV) {
// Compute the displacement.
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
......@@ -299,44 +308,42 @@ void ReferenceCustomGBIxn::calculateOnePairEnergyTerm(int index, int atom1, int
// Record variables for evaluating expressions.
map<string, double> variables = globalParameters;
for (int i = 0; i < (int) paramNames.size(); i++) {
variables[particleParamNames[i*2]] = atomParameters[atom1][i];
variables[particleParamNames[i*2+1]] = atomParameters[atom2][i];
for (int i = 0; i < (int) paramIndex.size(); i++) {
expressionSet.setVariable(particleParamIndex[i*2], atomParameters[atom1][i]);
expressionSet.setVariable(particleParamIndex[i*2+1], atomParameters[atom2][i]);
}
variables["r"] = r;
for (int i = 0; i < (int) valueNames.size(); i++) {
variables[particleValueNames[i*2]] = values[i][atom1];
variables[particleValueNames[i*2+1]] = values[i][atom2];
expressionSet.setVariable(rIndex, r);
for (int i = 0; i < (int) valueIndex.size(); i++) {
expressionSet.setVariable(particleValueIndex[i*2], values[i][atom1]);
expressionSet.setVariable(particleValueIndex[i*2+1], values[i][atom2]);
}
// Evaluate the energy and its derivatives.
if (totalEnergy != NULL)
*totalEnergy += (RealOpenMM) energyExpressions[index].evaluate(variables);
RealOpenMM dEdR = (RealOpenMM) energyDerivExpressions[index][0].evaluate(variables);
*totalEnergy += (RealOpenMM) energyExpressions[index].evaluate();
RealOpenMM dEdR = (RealOpenMM) energyDerivExpressions[index][0].evaluate();
dEdR *= 1/r;
for (int i = 0; i < 3; i++) {
forces[atom1][i] -= dEdR*deltaR[i];
forces[atom2][i] += dEdR*deltaR[i];
}
for (int i = 0; i < (int) valueNames.size(); i++) {
dEdV[i][atom1] += (RealOpenMM) energyDerivExpressions[index][2*i+1].evaluate(variables);
dEdV[i][atom2] += (RealOpenMM) energyDerivExpressions[index][2*i+2].evaluate(variables);
for (int i = 0; i < (int) valueIndex.size(); i++) {
dEdV[i][atom1] += (RealOpenMM) energyDerivExpressions[index][2*i+1].evaluate();
dEdV[i][atom2] += (RealOpenMM) energyDerivExpressions[index][2*i+2].evaluate();
}
}
void ReferenceCustomGBIxn::calculateChainRuleForces(int numAtoms, vector<RealVec>& atomCoordinates, RealOpenMM** atomParameters,
const vector<vector<RealOpenMM> >& values, const map<string, double>& globalParameters,
const vector<set<int> >& exclusions, vector<RealVec>& forces, vector<vector<RealOpenMM> >& dEdV) const {
const vector<vector<RealOpenMM> >& values, const vector<set<int> >& exclusions, vector<RealVec>& forces, vector<vector<RealOpenMM> >& dEdV) {
if (cutoff) {
// Loop over all pairs in the neighbor list.
for (int i = 0; i < (int) neighborList->size(); i++) {
OpenMM::AtomPair pair = (*neighborList)[i];
bool isExcluded = (exclusions[pair.first].find(pair.second) != exclusions[pair.first].end());
calculateOnePairChainRule(pair.first, pair.second, atomCoordinates, atomParameters, globalParameters, values, forces, dEdV, isExcluded);
calculateOnePairChainRule(pair.second, pair.first, atomCoordinates, atomParameters, globalParameters, values, forces, dEdV, isExcluded);
calculateOnePairChainRule(pair.first, pair.second, atomCoordinates, atomParameters, values, forces, dEdV, isExcluded);
calculateOnePairChainRule(pair.second, pair.first, atomCoordinates, atomParameters, values, forces, dEdV, isExcluded);
}
}
else {
......@@ -345,35 +352,34 @@ void ReferenceCustomGBIxn::calculateChainRuleForces(int numAtoms, vector<RealVec
for (int i = 0; i < numAtoms; i++) {
for (int j = i+1; j < numAtoms; j++) {
bool isExcluded = (exclusions[i].find(j) != exclusions[i].end());
calculateOnePairChainRule(i, j, atomCoordinates, atomParameters, globalParameters, values, forces, dEdV, isExcluded);
calculateOnePairChainRule(j, i, atomCoordinates, atomParameters, globalParameters, values, forces, dEdV, isExcluded);
calculateOnePairChainRule(i, j, atomCoordinates, atomParameters, values, forces, dEdV, isExcluded);
calculateOnePairChainRule(j, i, atomCoordinates, atomParameters, values, forces, dEdV, isExcluded);
}
}
}
// Compute chain rule terms for computed values that depend explicitly on particle coordinates.
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];
expressionSet.setVariable(xIndex, atomCoordinates[i][0]);
expressionSet.setVariable(yIndex, atomCoordinates[i][1]);
expressionSet.setVariable(zIndex, atomCoordinates[i][2]);
vector<RealOpenMM> dVdX(valueDerivExpressions.size(), 0.0);
vector<RealOpenMM> dVdY(valueDerivExpressions.size(), 0.0);
vector<RealOpenMM> dVdZ(valueDerivExpressions.size(), 0.0);
for (int j = 0; j < (int) paramNames.size(); j++)
variables[paramNames[j]] = atomParameters[i][j];
for (int j = 1; j < (int) valueNames.size(); j++) {
variables[valueNames[j-1]] = values[j-1][i];
for (int j = 0; j < (int) paramIndex.size(); j++)
expressionSet.setVariable(paramIndex[j], atomParameters[i][j]);
for (int j = 1; j < (int) valueIndex.size(); j++) {
expressionSet.setVariable(valueIndex[j-1], values[j-1][i]);
for (int k = 1; k < j; k++) {
RealOpenMM dVdV = (RealOpenMM) valueDerivExpressions[j][k].evaluate(variables);
RealOpenMM dVdV = (RealOpenMM) valueDerivExpressions[j][k].evaluate();
dVdX[j] += dVdV*dVdX[k];
dVdY[j] += dVdV*dVdY[k];
dVdZ[j] += dVdV*dVdZ[k];
}
dVdX[j] += (RealOpenMM) valueGradientExpressions[j][0].evaluate(variables);
dVdY[j] += (RealOpenMM) valueGradientExpressions[j][1].evaluate(variables);
dVdZ[j] += (RealOpenMM) valueGradientExpressions[j][2].evaluate(variables);
dVdX[j] += (RealOpenMM) valueGradientExpressions[j][0].evaluate();
dVdY[j] += (RealOpenMM) valueGradientExpressions[j][1].evaluate();
dVdZ[j] += (RealOpenMM) valueGradientExpressions[j][2].evaluate();
forces[i][0] -= dEdV[j][i]*dVdX[j];
forces[i][1] -= dEdV[j][i]*dVdY[j];
forces[i][2] -= dEdV[j][i]*dVdZ[j];
......@@ -382,8 +388,8 @@ void ReferenceCustomGBIxn::calculateChainRuleForces(int numAtoms, vector<RealVec
}
void ReferenceCustomGBIxn::calculateOnePairChainRule(int atom1, int atom2, vector<RealVec>& atomCoordinates, RealOpenMM** atomParameters,
const map<string, double>& globalParameters, const vector<vector<RealOpenMM> >& values, vector<RealVec>& forces,
vector<vector<RealOpenMM> >& dEdV, bool isExcluded) const {
const vector<vector<RealOpenMM> >& values, vector<RealVec>& forces,
vector<vector<RealOpenMM> >& dEdV, bool isExcluded) {
// Compute the displacement.
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
......@@ -397,14 +403,13 @@ void ReferenceCustomGBIxn::calculateOnePairChainRule(int atom1, int atom2, vecto
// Record variables for evaluating expressions.
map<string, double> variables = globalParameters;
for (int i = 0; i < (int) paramNames.size(); i++) {
variables[particleParamNames[i*2]] = atomParameters[atom1][i];
variables[particleParamNames[i*2+1]] = atomParameters[atom2][i];
for (int i = 0; i < (int) paramIndex.size(); i++) {
expressionSet.setVariable(particleParamIndex[i*2], atomParameters[atom1][i]);
expressionSet.setVariable(particleParamIndex[i*2+1], atomParameters[atom2][i]);
}
variables["r"] = r;
variables[particleValueNames[0]] = values[0][atom1];
variables[particleValueNames[1]] = values[0][atom2];
expressionSet.setVariable(rIndex, r);
expressionSet.setVariable(particleValueIndex[0], values[0][atom1]);
expressionSet.setVariable(particleValueIndex[1], values[0][atom2]);
// Evaluate the derivative of each parameter with respect to position and apply forces.
......@@ -415,24 +420,23 @@ void ReferenceCustomGBIxn::calculateOnePairChainRule(int atom1, int atom2, vecto
vector<RealOpenMM> dVdR1(valueDerivExpressions.size(), 0.0);
vector<RealOpenMM> dVdR2(valueDerivExpressions.size(), 0.0);
if (!isExcluded || valueTypes[0] != OpenMM::CustomGBForce::ParticlePair) {
dVdR1[0] = (RealOpenMM) valueDerivExpressions[0][0].evaluate(variables);;
dVdR1[0] = (RealOpenMM) valueDerivExpressions[0][0].evaluate();
dVdR2[0] = -dVdR1[0];
for (int i = 0; i < 3; i++) {
forces[atom1][i] -= dEdV[0][atom1]*dVdR1[0]*deltaR[i];
forces[atom2][i] -= dEdV[0][atom1]*dVdR2[0]*deltaR[i];
}
}
variables = globalParameters;
for (int i = 0; i < (int) paramNames.size(); i++)
variables[paramNames[i]] = atomParameters[atom1][i];
variables[valueNames[0]] = values[0][atom1];
for (int i = 1; i < (int) valueNames.size(); i++) {
variables[valueNames[i]] = values[i][atom1];
variables["x"] = atomCoordinates[atom1][0];
variables["y"] = atomCoordinates[atom1][1];
variables["z"] = atomCoordinates[atom1][2];
for (int i = 0; i < (int) paramIndex.size(); i++)
expressionSet.setVariable(paramIndex[i], atomParameters[atom1][i]);
expressionSet.setVariable(valueIndex[0], values[0][atom1]);
for (int i = 1; i < (int) valueIndex.size(); i++) {
expressionSet.setVariable(valueIndex[i], values[i][atom1]);
expressionSet.setVariable(xIndex, atomCoordinates[atom1][0]);
expressionSet.setVariable(yIndex, atomCoordinates[atom1][1]);
expressionSet.setVariable(zIndex, atomCoordinates[atom1][2]);
for (int j = 0; j < i; j++) {
RealOpenMM dVdV = (RealOpenMM) valueDerivExpressions[i][j].evaluate(variables);
RealOpenMM dVdV = (RealOpenMM) valueDerivExpressions[i][j].evaluate();
dVdR1[i] += dVdV*dVdR1[j];
dVdR2[i] += dVdV*dVdR2[j];
}
......
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