Skip to content

GitLab

Commit 18587bd4 authored by peastman's avatar peastman
Browse files

Merge pull request #173 from peastman/master 

Performance improvements to CPU implementation of custom forces
parents 54bfc138 54e47d2d
Hide whitespace changes
Inline Side-by-side
Showing with 139 additions and 63 deletions
platforms/reference/src/SimTKReference/ReferenceCustomExternalIxn.cpp
View file @ 18587bd4
/* Portions copyright (c) 2009 Stanford University and Simbios.
/* Portions copyright (c) 2009-2013 Stanford University and Simbios.
* Contributors: Peter Eastman
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -40,18 +40,37 @@ using namespace OpenMM;
--------------------------------------------------------------------------------------- */
ReferenceCustomExternalIxn::ReferenceCustomExternalIxn(const Lepton::ExpressionProgram& energyExpression,
const Lepton::ExpressionProgram& forceExpressionX, const Lepton::ExpressionProgram& forceExpressionY,
const Lepton::ExpressionProgram& forceExpressionZ, const vector<string>& parameterNames, map<string, double> globalParameters) :
ReferenceCustomExternalIxn::ReferenceCustomExternalIxn(const Lepton::CompiledExpression& energyExpression,
const Lepton::CompiledExpression& forceExpressionX, const Lepton::CompiledExpression& forceExpressionY,
const Lepton::CompiledExpression& forceExpressionZ, const vector<string>& parameterNames, map<string, double> globalParameters) :
energyExpression(energyExpression), forceExpressionX(forceExpressionX), forceExpressionY(forceExpressionY),
forceExpressionZ(forceExpressionZ), paramNames(parameterNames), globalParameters(globalParameters) {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceCustomExternalIxn::ReferenceCustomExternalIxn";
// ---------------------------------------------------------------------------------------
forceExpressionZ(forceExpressionZ) {
energyX = ReferenceForce::getVariablePointer(this->energyExpression, "x");
energyY = ReferenceForce::getVariablePointer(this->energyExpression, "y");
energyZ = ReferenceForce::getVariablePointer(this->energyExpression, "z");
forceXX = ReferenceForce::getVariablePointer(this->forceExpressionX, "x");
forceXY = ReferenceForce::getVariablePointer(this->forceExpressionX, "y");
forceXZ = ReferenceForce::getVariablePointer(this->forceExpressionX, "z");
forceYX = ReferenceForce::getVariablePointer(this->forceExpressionY, "x");
forceYY = ReferenceForce::getVariablePointer(this->forceExpressionY, "y");
forceYZ = ReferenceForce::getVariablePointer(this->forceExpressionY, "z");
forceZX = ReferenceForce::getVariablePointer(this->forceExpressionZ, "x");
forceZY = ReferenceForce::getVariablePointer(this->forceExpressionZ, "y");
forceZZ = ReferenceForce::getVariablePointer(this->forceExpressionZ, "z");
numParameters = parameterNames.size();
for (int i = 0; i < (int) numParameters; i++) {
energyParams.push_back(ReferenceForce::getVariablePointer(this->energyExpression, parameterNames[i]));
forceXParams.push_back(ReferenceForce::getVariablePointer(this->forceExpressionX, parameterNames[i]));
forceYParams.push_back(ReferenceForce::getVariablePointer(this->forceExpressionY, parameterNames[i]));
forceZParams.push_back(ReferenceForce::getVariablePointer(this->forceExpressionZ, parameterNames[i]));
}
for (map<string, double>::const_iterator iter = globalParameters.begin(); iter != globalParameters.end(); ++iter) {
ReferenceForce::setVariable(ReferenceForce::getVariablePointer(this->energyExpression, iter->first), iter->second);
ReferenceForce::setVariable(ReferenceForce::getVariablePointer(this->forceExpressionX, iter->first), iter->second);
ReferenceForce::setVariable(ReferenceForce::getVariablePointer(this->forceExpressionY, iter->first), iter->second);
ReferenceForce::setVariable(ReferenceForce::getVariablePointer(this->forceExpressionZ, iter->first), iter->second);
}
}
/**---------------------------------------------------------------------------------------
......@@ -90,18 +109,30 @@ void ReferenceCustomExternalIxn::calculateForce( int atomIndex,
static const std::string methodName = "\nReferenceCustomExternalIxn::calculateBondIxn";
map<string, double> variables = globalParameters;
for (int i = 0; i < (int) paramNames.size(); ++i)
variables[paramNames[i]] = parameters[i];
variables["x"] = atomCoordinates[atomIndex][0];
variables["y"] = atomCoordinates[atomIndex][1];
variables["z"] = atomCoordinates[atomIndex][2];
for (int i = 0; i < numParameters; i++) {
ReferenceForce::setVariable(energyParams[i], parameters[i]);
ReferenceForce::setVariable(forceXParams[i], parameters[i]);
ReferenceForce::setVariable(forceYParams[i], parameters[i]);
ReferenceForce::setVariable(forceZParams[i], parameters[i]);
}
ReferenceForce::setVariable(energyX, atomCoordinates[atomIndex][0]);
ReferenceForce::setVariable(energyY, atomCoordinates[atomIndex][1]);
ReferenceForce::setVariable(energyZ, atomCoordinates[atomIndex][2]);
ReferenceForce::setVariable(forceXX, atomCoordinates[atomIndex][0]);
ReferenceForce::setVariable(forceXY, atomCoordinates[atomIndex][1]);
ReferenceForce::setVariable(forceXZ, atomCoordinates[atomIndex][2]);
ReferenceForce::setVariable(forceYX, atomCoordinates[atomIndex][0]);
ReferenceForce::setVariable(forceYY, atomCoordinates[atomIndex][1]);
ReferenceForce::setVariable(forceYZ, atomCoordinates[atomIndex][2]);
ReferenceForce::setVariable(forceZX, atomCoordinates[atomIndex][0]);
ReferenceForce::setVariable(forceZY, atomCoordinates[atomIndex][1]);
ReferenceForce::setVariable(forceZZ, atomCoordinates[atomIndex][2]);
// ---------------------------------------------------------------------------------------
forces[atomIndex][0] -= (RealOpenMM) forceExpressionX.evaluate(variables);
forces[atomIndex][1] -= (RealOpenMM) forceExpressionY.evaluate(variables);
forces[atomIndex][2] -= (RealOpenMM) forceExpressionZ.evaluate(variables);
forces[atomIndex][0] -= (RealOpenMM) forceExpressionX.evaluate();
forces[atomIndex][1] -= (RealOpenMM) forceExpressionY.evaluate();
forces[atomIndex][2] -= (RealOpenMM) forceExpressionZ.evaluate();
if (energy != NULL)
*energy += (RealOpenMM) energyExpression.evaluate(variables);
*energy += (RealOpenMM) energyExpression.evaluate();
}
platforms/reference/src/SimTKReference/ReferenceCustomNonbondedIxn.cpp
View file @ 18587bd4
/* Portions copyright (c) 2009 Stanford University and Simbios.
/* Portions copyright (c) 2009-2013 Stanford University and Simbios.
* Contributors: Peter Eastman
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -45,8 +45,8 @@ using OpenMM::RealVec;
--------------------------------------------------------------------------------------- */
ReferenceCustomNonbondedIxn::ReferenceCustomNonbondedIxn(const Lepton::ExpressionProgram& energyExpression,
const Lepton::ExpressionProgram& forceExpression, const vector<string>& parameterNames) :
ReferenceCustomNonbondedIxn::ReferenceCustomNonbondedIxn(const Lepton::CompiledExpression& energyExpression,
const Lepton::CompiledExpression& forceExpression, const vector<string>& parameterNames) :
cutoff(false), useSwitch(false), periodic(false), energyExpression(energyExpression), forceExpression(forceExpression), paramNames(parameterNames) {
// ---------------------------------------------------------------------------------------
......@@ -55,11 +55,14 @@ ReferenceCustomNonbondedIxn::ReferenceCustomNonbondedIxn(const Lepton::Expressio
// ---------------------------------------------------------------------------------------
energyR = ReferenceForce::getVariablePointer(this->energyExpression, "r");
forceR = ReferenceForce::getVariablePointer(this->forceExpression, "r");
for (int i = 0; i < (int) paramNames.size(); i++) {
for (int j = 1; j < 3; j++) {
stringstream name;
name << paramNames[i] << j;
particleParamNames.push_back(name.str());
energyParticleParams.push_back(ReferenceForce::getVariablePointer(this->energyExpression, name.str()));
forceParticleParams.push_back(ReferenceForce::getVariablePointer(this->forceExpression, name.str()));
}
}
}
......@@ -166,9 +169,12 @@ void ReferenceCustomNonbondedIxn::setUseSwitchingFunction( RealOpenMM distance )
void ReferenceCustomNonbondedIxn::calculatePairIxn( int numberOfAtoms, vector<RealVec>& atomCoordinates,
RealOpenMM** atomParameters, vector<set<int> >& exclusions,
RealOpenMM* fixedParameters, const map<string, double>& globalParameters, vector<RealVec>& forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const {
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) {
map<string, double> variables = globalParameters;
for (map<string, double>::const_iterator iter = globalParameters.begin(); iter != globalParameters.end(); ++iter) {
ReferenceForce::setVariable(ReferenceForce::getVariablePointer(energyExpression, iter->first), iter->second);
ReferenceForce::setVariable(ReferenceForce::getVariablePointer(forceExpression, iter->first), iter->second);
}
if (interactionGroups.size() > 0) {
// The user has specified interaction groups, so compute only the requested interactions.
......@@ -182,10 +188,12 @@ void ReferenceCustomNonbondedIxn::calculatePairIxn( int numberOfAtoms, vector<Re
if (*atom1 > *atom2 && set1.find(*atom2) != set1.end() && set2.find(*atom1) != set2.end())
continue; // Both atoms are in both sets, so skip duplicate interactions.
for (int j = 0; j < (int) paramNames.size(); j++) {
variables[particleParamNames[j*2]] = atomParameters[*atom1][j];
variables[particleParamNames[j*2+1]] = atomParameters[*atom2][j];
ReferenceForce::setVariable(energyParticleParams[j*2], atomParameters[*atom1][j]);
ReferenceForce::setVariable(energyParticleParams[j*2+1], atomParameters[*atom2][j]);
ReferenceForce::setVariable(forceParticleParams[j*2], atomParameters[*atom1][j]);
ReferenceForce::setVariable(forceParticleParams[j*2+1], atomParameters[*atom2][j]);
}
calculateOneIxn(*atom1, *atom2, atomCoordinates, variables, forces, energyByAtom, totalEnergy);
calculateOneIxn(*atom1, *atom2, atomCoordinates, forces, energyByAtom, totalEnergy);
}
}
}
......@@ -196,10 +204,12 @@ void ReferenceCustomNonbondedIxn::calculatePairIxn( int numberOfAtoms, vector<Re
for (int i = 0; i < (int) neighborList->size(); i++) {
OpenMM::AtomPair pair = (*neighborList)[i];
for (int j = 0; j < (int) paramNames.size(); j++) {
variables[particleParamNames[j*2]] = atomParameters[pair.first][j];
variables[particleParamNames[j*2+1]] = atomParameters[pair.second][j];
ReferenceForce::setVariable(energyParticleParams[j*2], atomParameters[pair.first][j]);
ReferenceForce::setVariable(energyParticleParams[j*2+1], atomParameters[pair.second][j]);
ReferenceForce::setVariable(forceParticleParams[j*2], atomParameters[pair.first][j]);
ReferenceForce::setVariable(forceParticleParams[j*2+1], atomParameters[pair.second][j]);
}
calculateOneIxn(pair.first, pair.second, atomCoordinates, variables, forces, energyByAtom, totalEnergy);
calculateOneIxn(pair.first, pair.second, atomCoordinates, forces, energyByAtom, totalEnergy);
}
}
else {
......@@ -209,10 +219,12 @@ void ReferenceCustomNonbondedIxn::calculatePairIxn( int numberOfAtoms, vector<Re
for (int jj = ii+1; jj < numberOfAtoms; jj++) {
if (exclusions[jj].find(ii) == exclusions[jj].end()) {
for (int j = 0; j < (int) paramNames.size(); j++) {
variables[particleParamNames[j*2]] = atomParameters[ii][j];
variables[particleParamNames[j*2+1]] = atomParameters[jj][j];
ReferenceForce::setVariable(energyParticleParams[j*2], atomParameters[ii][j]);
ReferenceForce::setVariable(energyParticleParams[j*2+1], atomParameters[jj][j]);
ReferenceForce::setVariable(forceParticleParams[j*2], atomParameters[ii][j]);
ReferenceForce::setVariable(forceParticleParams[j*2+1], atomParameters[jj][j]);
}
calculateOneIxn(ii, jj, atomCoordinates, variables, forces, energyByAtom, totalEnergy);
calculateOneIxn(ii, jj, atomCoordinates, forces, energyByAtom, totalEnergy);
}
}
}
......@@ -233,9 +245,8 @@ void ReferenceCustomNonbondedIxn::calculatePairIxn( int numberOfAtoms, vector<Re
--------------------------------------------------------------------------------------- */
void ReferenceCustomNonbondedIxn::calculateOneIxn( int ii, int jj, vector<RealVec>& atomCoordinates,
map<string, double>& variables, vector<RealVec>& forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const {
void ReferenceCustomNonbondedIxn::calculateOneIxn( int ii, int jj, vector<RealVec>& atomCoordinates, vector<RealVec>& forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) {
// ---------------------------------------------------------------------------------------
......@@ -266,9 +277,10 @@ void ReferenceCustomNonbondedIxn::calculateOneIxn( int ii, int jj, vector<RealVe
// accumulate forces
variables["r"] = r;
RealOpenMM dEdR = (RealOpenMM) (forceExpression.evaluate(variables)/(deltaR[ReferenceForce::RIndex]));
RealOpenMM energy = (RealOpenMM) energyExpression.evaluate(variables);
ReferenceForce::setVariable(energyR, r);
ReferenceForce::setVariable(forceR, r);
RealOpenMM dEdR = (RealOpenMM) (forceExpression.evaluate()/(deltaR[ReferenceForce::RIndex]));
RealOpenMM energy = (RealOpenMM) energyExpression.evaluate();
if (useSwitch) {
if (r > switchingDistance) {
RealOpenMM t = (r-switchingDistance)/(cutoffDistance-switchingDistance);
......
platforms/reference/src/SimTKReference/ReferenceCustomTorsionIxn.cpp
View file @ 18587bd4
/* Portions copyright (c) 2010 Stanford University and Simbios.
/* Portions copyright (c) 2010-2013 Stanford University and Simbios.
* Contributors: Peter Eastman
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -39,16 +39,21 @@ using namespace OpenMM;
--------------------------------------------------------------------------------------- */
ReferenceCustomTorsionIxn::ReferenceCustomTorsionIxn(const Lepton::ExpressionProgram& energyExpression,
const Lepton::ExpressionProgram& forceExpression, const vector<string>& parameterNames, map<string, double> globalParameters) :
energyExpression(energyExpression), forceExpression(forceExpression), paramNames(parameterNames), globalParameters(globalParameters) {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceCustomTorsionIxn::ReferenceCustomTorsionIxn";
// ---------------------------------------------------------------------------------------
ReferenceCustomTorsionIxn::ReferenceCustomTorsionIxn(const Lepton::CompiledExpression& energyExpression,
const Lepton::CompiledExpression& forceExpression, const vector<string>& parameterNames, map<string, double> globalParameters) :
energyExpression(energyExpression), forceExpression(forceExpression) {
energyTheta = ReferenceForce::getVariablePointer(this->energyExpression, "theta");
forceTheta = ReferenceForce::getVariablePointer(this->forceExpression, "theta");
numParameters = parameterNames.size();
for (int i = 0; i < (int) numParameters; i++) {
energyParams.push_back(ReferenceForce::getVariablePointer(this->energyExpression, parameterNames[i]));
forceParams.push_back(ReferenceForce::getVariablePointer(this->forceExpression, parameterNames[i]));
}
for (map<string, double>::const_iterator iter = globalParameters.begin(); iter != globalParameters.end(); ++iter) {
ReferenceForce::setVariable(ReferenceForce::getVariablePointer(this->energyExpression, iter->first), iter->second);
ReferenceForce::setVariable(ReferenceForce::getVariablePointer(this->forceExpression, iter->first), iter->second);
}
}
/**---------------------------------------------------------------------------------------
......@@ -91,9 +96,10 @@ void ReferenceCustomTorsionIxn::calculateBondIxn( int* atomIndices,
static const RealOpenMM one = 1.0;
RealOpenMM deltaR[3][ReferenceForce::LastDeltaRIndex];
map<string, double> variables = globalParameters;
for (int i = 0; i < (int) paramNames.size(); ++i)
variables[paramNames[i]] = parameters[i];
for (int i = 0; i < numParameters; i++) {
ReferenceForce::setVariable(energyParams[i], parameters[i]);
ReferenceForce::setVariable(forceParams[i], parameters[i]);
}
// ---------------------------------------------------------------------------------------
......@@ -118,13 +124,13 @@ void ReferenceCustomTorsionIxn::calculateBondIxn( int* atomIndices,
RealOpenMM dotDihedral;
RealOpenMM signOfAngle;
variables["theta"] = getDihedralAngleBetweenThreeVectors(deltaR[0], deltaR[1], deltaR[2],
crossProduct, &dotDihedral, deltaR[0],
&signOfAngle, 1);
RealOpenMM angle = getDihedralAngleBetweenThreeVectors(deltaR[0], deltaR[1], deltaR[2], crossProduct, &dotDihedral, deltaR[0], &signOfAngle, 1);
ReferenceForce::setVariable(energyTheta, angle);
ReferenceForce::setVariable(forceTheta, angle);
// evaluate delta angle, dE/d(angle)
RealOpenMM dEdAngle = (RealOpenMM) forceExpression.evaluate(variables);
RealOpenMM dEdAngle = (RealOpenMM) forceExpression.evaluate();
// compute force
......@@ -166,6 +172,6 @@ void ReferenceCustomTorsionIxn::calculateBondIxn( int* atomIndices,
// accumulate energies
if (totalEnergy != NULL)
*totalEnergy += (RealOpenMM) energyExpression.evaluate(variables);
*totalEnergy += (RealOpenMM) energyExpression.evaluate();
}
platforms/reference/src/SimTKReference/ReferenceForce.cpp
View file @ 18587bd4
/* Portions copyright (c) 2006 Stanford University and Simbios.
/* Portions copyright (c) 2006-2013 Stanford University and Simbios.
* Contributors: Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
......@@ -161,3 +161,14 @@ void ReferenceForce::getDeltaROnly( const RealOpenMM* atomCoordinatesI,
deltaR[YIndex] = atomCoordinatesJ[1] - atomCoordinatesI[1];
deltaR[ZIndex] = atomCoordinatesJ[2] - atomCoordinatesI[2];
}
double* ReferenceForce::getVariablePointer(Lepton::CompiledExpression& expression, const std::string& name) {
if (expression.getVariables().find(name) == expression.getVariables().end())
return NULL;
return &expression.getVariableReference(name);
}
void ReferenceForce::setVariable(double* pointer, double value) {
if (pointer != NULL)
*pointer = value;
}
tests/TestParser.cpp
View file @ 18587bd4
......@@ -56,6 +56,12 @@ void verifyEvaluation(const string& expression, double expectedValue) {
ExpressionProgram program = parsed.createProgram();
value = program.evaluate();
ASSERT_EQUAL_TOL(expectedValue, value, 1e-10);
// Create a CompiledExpression and see if that also gives the same result.
CompiledExpression compiled = parsed.createCompiledExpression();
value = compiled.evaluate();
ASSERT_EQUAL_TOL(expectedValue, value, 1e-10);
}
/**
......@@ -86,6 +92,16 @@ void verifyEvaluation(const string& expression, double x, double y, double expec
value = program.evaluate(variables);
ASSERT_EQUAL_TOL(expectedValue, value, 1e-10);
// Create a CompiledExpression and see if that also gives the same result.
CompiledExpression compiled = parsed.createCompiledExpression();
if (compiled.getVariables().find("x") != compiled.getVariables().end())
compiled.getVariableReference("x") = x;
if (compiled.getVariables().find("y") != compiled.getVariables().end())
compiled.getVariableReference("y") = y;
value = compiled.evaluate();
ASSERT_EQUAL_TOL(expectedValue, value, 1e-10);
// Make sure that variable renaming works.
variables.clear();
......
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