Skip to content

GitLab

Commit 4747b0f2 authored by Peter Eastman's avatar Peter Eastman
Browse files

Continuing implementation of CustomGBForce

parent bdc71cd4
Show whitespace changes
Inline Side-by-side
Showing with 523 additions and 160 deletions
openmmapi/include/openmm/CustomGBForce.h
View file @ 4747b0f2
......@@ -43,6 +43,97 @@
namespace OpenMM {
/**
* This class implements complex, multiple stage nonbonded interactions between particles. It is designed primarily
* for implementing Generalized Born implicit solvation models, although it is not strictly limited to that purpose.
* The interaction is specified as a series of computations, each defined by an arbitrary algebraic expression.
* It also allows tabulated functions to be defined and used with the computations. It optionally supports periodic boundary
* conditions and cutoffs for long range interactions.
*
* The computation consists of calculating some number of per-particle <i>computed values</i>, followed by one or more
* <i>energy terms</i>. A computed value is a scalar value that is computed for each particle in the system. It may
* depend on an arbitrary set of global and per-particle parameters, and well as on other computed values that have
* been calculating before it. Once all computed values have been calculated, the energy terms and their derivatives
* are evaluated to determine the system energy and particle forces. The energy terms may depend on global parameters,
* per-particle parameters, and per-particle computed values.
*
* When specifying a computed value or energy term, you provide an algebraic expression to evaluate and a <i>computation type</i>
* describing how is the expression is to be evaluated. There are two main types of computations:
*
* <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
* 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
* value, the value for a particular particle is calculated by pairing it with every other particle in the system, evaluating
* the expression for each pair, and summing them. For an energy term, each particle pair makes an independent contribution to
* the System energy. (Note that energy terms are assumed to be symmetric with respect to the two interacting particles, and
* therefore are evaluated only once per pair. In contrast, computed values need not be symmetric and therefore are calculated
* twice for each pair: once when calculating the value for the first particle, and again when calculating the value for the
* second particle.)</li>
* </ul>
*
* Be aware that, although this class is extremely general in the computations it can define, particular Platforms may only support
* more restricted types of computations. In particular, all currently existing Platforms require that the first computed value
* <i>must</i> be a particle pair computation, and all computed values that follow it <i>must</i> be single particle computations.
* This is sufficient for most Generalized Born models, but might not permit some other types of calculations to be implemented.
*
* This is a complicated class to use, and an example may help to clarify it. The following code implements the OBC variant
* of the GB/SA solvation model, using the ACE approximation to estimate surface area:
*
* <tt><pre>
* CustomGBForce* custom = new CustomGBForce();
* custom->addPerParticleParameter("q");
* custom->addPerParticleParameter("radius");
* custom->addPerParticleParameter("scale");
* custom->addGlobalParameter("solventDielectric", obc->getSolventDielectric());
* custom->addGlobalParameter("soluteDielectric", obc->getSoluteDielectric());
* custom->addComputedValue("I", "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r*(1/U^2-1/L^2))+0.5*log(L/U)/r+0.25*(sr2*sr2/r)*(1/L^2-1/U^2)+C);"
* "U=r+sr2;"
* "C=2*(1/or1-1/L)*step(sr2-r-or1);"
* "L=step(or1-D)*or1+step(D-or1)*D;"
* "D=step(r-sr2)*(r-sr2)+step(sr2-r)*(sr2-r);"
* "sr1 = scale1*or1; sr2 = scale2*or2;"
* "or1 = radius1-0.009; or2 = radius2-0.009", CustomGBForce::ParticlePair);
* custom->addComputedValue("B", "1/(1/or-tanh(1*psi-0.8*psi^2+4.85*psi^3)/radius);"
* "psi=I*or; or=radius-0.009", CustomGBForce::SingleParticle);
* custom->addEnergyTerm("28.3919551*(radius+0.14)^2*(radius/B)^6-0.5*138.935485*(1/soluteDielectric-1/solventDielectric)*q^2/B",
* CustomGBForce::SingleParticle);
* custom->addEnergyTerm("-138.935485*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;"
* "f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce::ParticlePair);
* </pre></tt>
*
* It begins by defining three per-particle parameters (charge, atomic radius, and scale factor) and two global parameters
* (the dielectric constants for the solute and solvent). It then defines a computed value "I" of type ParticlePair. The
* expression for evaluating it is a complicated function of the distance between each pair of particles (r), their atomic
* radii (radius1 and radius2), and their scale factors (scale1 and scale2). Very roughly speaking, it is a measure of the
* overlap between each particle and other nearby particles.
*
* Next a computation is defined for the Born Radius (B). It is computed independently for each particle, and is a function of
* that particle's atomic radius and the intermediate value I defined above.
*
* Finally, two energy terms are defined. The first one is computed for each particle and represents the surface area term,
* as well as the self interaction part of the polarization energy. The second term is calculated for each pair of particles,
* and represents the screening of electrostatic interactions by the solvent.
*
* After defining the force as shown above, you should then call addParticle() once for each particle in the System to set the
* values of its per-particle parameters (q, radius, and scale). The number of particles for which you set parameters must be
* exactly equal to the number of particles in the System, or else an exception will be thrown when you try to create a Context.
* After a particle has been added, you can modify its parameters by calling setParticleParameters().
*
* CustomNonbondedForce also lets you specify "exclusions", particular pairs of particles whose interactions should be
* omitted from calculations. This is most often used for particles that are bonded to each other. Even if you specify exclusions,
* however, you can use the computation type ParticlePairNoExclusions to indicate that exclusions should not be applied to a
* particular piece of the computation.
*
* Expressions may involve the operators + (add), - (subtract), * (multiply), / (divide), and ^ (power), and the following
* functions: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, step. All trigonometric functions
* are defined in radians, and log is the natural logarithm. step(x) = 0 if x is less than 0, 1 otherwise. The names of per-particle parameters
* have the suffix "1" or "2" appended to them to indicate the values for the two interacting particles. As seen in the above example,
* the expression may also involve intermediate quantities that are defined following the main expression, using ";" as a separator.
*
* In addition, you can call addFunction() to define a new function based on tabulated values. You specify a vector of
* values, and an interpolating or approximating spline is created from them. That function can then appear in the expression.
*/
class OPENMM_EXPORT CustomGBForce : public Force {
......@@ -66,9 +157,22 @@ public:
*/
CutoffPeriodic = 2,
};
/**
* This is an enumeration of the different ways in which a computed value or energy term can be calculated.
*/
enum ComputationType {
/**
* The value is computed independently for each particle, based only on the parameters and computed values for that particle.
*/
SingleParticle = 0,
/**
* The value is computed as a sum over all pairs of particles, except those which have been added as exclusions.
*/
ParticlePair = 1,
/**
* The value is computed as a sum over all pairs of particles. Unlike ParticlePair, the list of exclusions is ignored
* and all pairs are included in the sum, even those marked as exclusions.
*/
ParticlePairNoExclusions = 2
};
/**
......@@ -105,9 +209,15 @@ public:
int getNumFunctions() const {
return functions.size();
}
/**
* Get the number of per-particle computed values the interaction depends on.
*/
int getNumComputedValues() const {
return computedValues.size();
}
/**
* Get the number of terms in the energy computation.
*/
int getNumEnergyTerms() const {
return energyTerms.size();
}
......@@ -208,11 +318,108 @@ public:
* @param parameters the list of parameters for the specified particle
*/
void setParticleParameters(int index, const std::vector<double>& parameters);
/**
* Add a computed value to calculate for each particle.
*
* @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
* 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.
* Append "1" to a variable name to indicate the parameter for the particle whose
* value is being calculated, and "2" to indicate the particle it is interacting with.
* @param type the method to use for computing this value
*/
int addComputedValue(const std::string& name, const std::string& expression, ComputationType type);
/**
* Get the properties of a computed value.
*
* @param index the index of the computed value for which to get parameters
* @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
* 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.
* Append "1" to a variable name to indicate the parameter for the particle whose
* value is being calculated, and "2" to indicate the particle it is interacting with.
* @param type the method to use for computing this value
*/
void getComputedValueParameters(int index, std::string& name, std::string& expression, ComputationType& type) const;
/**
* Set the properties of a computed value.
*
* @param index the index of the computed value for which to set parameters
* @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
* 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.
* Append "1" to a variable name to indicate the parameter for the particle whose
* value is being calculated, and "2" to indicate the particle it is interacting with.
* @param type the method to use for computing this value
*/
void setComputedValueParameters(int index, const std::string& name, const std::string& expression, ComputationType type);
/**
* Add a term to the energy computation.
*
* @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
* 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
* the per-particle parameters and computed values for each of them.
* Append "1" to a variable name to indicate the parameter for the first particle
* in the pair and "2" to indicate the second particle in the pair.
* @param type the method to use for computing this value
*/
int addEnergyTerm(const std::string& expression, ComputationType type);
/**
* Get the properties of a term to the energy computation.
*
* @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
* 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
* the per-particle parameters and computed values for each of them.
* Append "1" to a variable name to indicate the parameter for the first particle
* in the pair and "2" to indicate the second particle in the pair.
* @param type the method to use for computing this value
*/
void getEnergyTermParameters(int index, std::string& expression, ComputationType& type) const;
/**
* Set the properties of a term to the energy computation.
*
* @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
* 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
* the per-particle parameters and computed values for each of them.
* Append "1" to a variable name to indicate the parameter for the first particle
* in the pair and "2" to indicate the second particle in the pair.
* @param type the method to use for computing this value
*/
void setEnergyTermParameters(int index, const std::string& expression, ComputationType type);
/**
* Add a particle pair to the list of interactions that should be excluded.
......
platforms/reference/src/ReferenceKernels.cpp
View file @ 4747b0f2
......@@ -55,6 +55,7 @@
#include "openmm/internal/ContextImpl.h"
#include "openmm/internal/NonbondedForceImpl.h"
#include "openmm/Integrator.h"
#include "openmm/OpenMMException.h"
#include "SimTKUtilities/SimTKOpenMMUtilities.h"
#include "lepton/CustomFunction.h"
#include "lepton/Parser.h"
......@@ -899,6 +900,18 @@ ReferenceCalcCustomGBForceKernel::~ReferenceCalcCustomGBForceKernel() {
}
void ReferenceCalcCustomGBForceKernel::initialize(const System& system, const CustomGBForce& force) {
if (force.getNumComputedValues() > 0) {
string name, expression;
CustomGBForce::ComputationType type;
force.getComputedValueParameters(0, name, expression, type);
if (type == CustomGBForce::SingleParticle)
throw OpenMMException("ReferencePlatform requires that the first computed value for a CustomGBForce be of type ParticlePair or ParticlePairNoExclusions.");
for (int i = 1; i < force.getNumComputedValues(); i++) {
force.getComputedValueParameters(i, name, expression, type);
if (type != CustomGBForce::SingleParticle)
throw OpenMMException("ReferencePlatform requires that a CustomGBForce only have one computed value of type ParticlePair or ParticlePairNoExclusions.");
}
}
// Record the exclusions.
......@@ -951,7 +964,6 @@ void ReferenceCalcCustomGBForceKernel::initialize(const System& system, const Cu
// Parse the expressions for computed values.
valueDerivExpressions.resize(force.getNumComputedValues());
for (int i = 0; i < force.getNumComputedValues(); i++) {
string name, expression;
CustomGBForce::ComputationType type;
......@@ -960,14 +972,10 @@ void ReferenceCalcCustomGBForceKernel::initialize(const System& system, const Cu
valueExpressions.push_back(ex.createProgram());
valueTypes.push_back(type);
valueNames.push_back(name);
if (type != CustomGBForce::SingleParticle)
valueDerivExpressions[i].push_back(ex.differentiate("r").optimize().createProgram());
for (int j = 0; j < i; j++) {
if (type == CustomGBForce::SingleParticle)
valueDerivExpressions[i].push_back(ex.differentiate(valueNames[j]).optimize().createProgram());
valueDerivExpressions.push_back(ex.differentiate(valueNames[i-1]).optimize().createProgram());
else
valueDerivExpressions[i].push_back(ex.differentiate(valueNames[j]+"1").optimize().createProgram());
}
valueDerivExpressions.push_back(ex.differentiate("r").optimize().createProgram());
}
// Parse the various expressions used to calculate the force.
......@@ -1013,7 +1021,7 @@ void ReferenceCalcCustomGBForceKernel::executeForces(ContextImpl& context) {
map<string, double> globalParameters;
for (int i = 0; i < (int) globalParameterNames.size(); i++)
globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
ixn.calculateIxn(numParticles, posData, particleParamArray, exclusions, globalParameters, forceData, 0, 0);
ixn.calculateIxn(numParticles, posData, particleParamArray, exclusions, globalParameters, forceData, 0);
}
double ReferenceCalcCustomGBForceKernel::executeEnergy(ContextImpl& context) {
......@@ -1032,7 +1040,7 @@ double ReferenceCalcCustomGBForceKernel::executeEnergy(ContextImpl& context) {
map<string, double> globalParameters;
for (int i = 0; i < (int) globalParameterNames.size(); i++)
globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
ixn.calculateIxn(numParticles, posData, particleParamArray, exclusions, globalParameters, forceData, 0, &energy);
ixn.calculateIxn(numParticles, posData, particleParamArray, exclusions, globalParameters, forceData, &energy);
disposeRealArray(forceData, numParticles);
return energy;
}
......
platforms/reference/src/ReferenceKernels.h
View file @ 4747b0f2
......@@ -512,7 +512,7 @@ private:
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<Lepton::ExpressionProgram> valueDerivExpressions;
std::vector<OpenMM::CustomGBForce::ComputationType> valueTypes;
std::vector<Lepton::ExpressionProgram> energyExpressions;
std::vector<std::vector<Lepton::ExpressionProgram> > energyDerivExpressions;
......
platforms/reference/src/SimTKReference/ReferenceCustomGBIxn.cpp
View file @ 4747b0f2
......@@ -44,7 +44,7 @@ using std::vector;
--------------------------------------------------------------------------------------- */
ReferenceCustomGBIxn::ReferenceCustomGBIxn(const vector<Lepton::ExpressionProgram>& valueExpressions,
const vector<vector<Lepton::ExpressionProgram> >& valueDerivExpressions,
const vector<Lepton::ExpressionProgram>& valueDerivExpressions,
const vector<string>& valueNames,
const vector<OpenMM::CustomGBForce::ComputationType>& valueTypes,
const vector<Lepton::ExpressionProgram>& energyExpressions,
......@@ -100,17 +100,13 @@ ReferenceCustomGBIxn::~ReferenceCustomGBIxn( ){
@param distance the cutoff distance
@param neighbors the neighbor list to use
@return ReferenceForce::DefaultReturn
--------------------------------------------------------------------------------------- */
int ReferenceCustomGBIxn::setUseCutoff( RealOpenMM distance, const OpenMM::NeighborList& neighbors ) {
void ReferenceCustomGBIxn::setUseCutoff( RealOpenMM distance, const OpenMM::NeighborList& neighbors ) {
cutoff = true;
cutoffDistance = distance;
neighborList = &neighbors;
return ReferenceForce::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
......@@ -121,11 +117,9 @@ ReferenceCustomGBIxn::~ReferenceCustomGBIxn( ){
@param boxSize the X, Y, and Z widths of the periodic box
@return ReferenceForce::DefaultReturn
--------------------------------------------------------------------------------------- */
int ReferenceCustomGBIxn::setPeriodic( RealOpenMM* boxSize ) {
void ReferenceCustomGBIxn::setPeriodic( RealOpenMM* boxSize ) {
assert(cutoff);
assert(boxSize[0] >= 2.0*cutoffDistance);
......@@ -135,37 +129,15 @@ ReferenceCustomGBIxn::~ReferenceCustomGBIxn( ){
periodicBoxSize[0] = boxSize[0];
periodicBoxSize[1] = boxSize[1];
periodicBoxSize[2] = boxSize[2];
return ReferenceForce::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Calculate the custom pair ixn
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param atomParameters atom parameters atomParameters[atomIndex][paramterIndex]
@param exclusions atom exclusion indices exclusions[atomIndex][atomToExcludeIndex]
exclusions[atomIndex][0] = number of exclusions
exclusions[atomIndex][1-no.] = atom indices of atoms to excluded from
interacting w/ atom atomIndex
@param fixedParameters non atom parameters (not currently used)
@param globalParameters the values of global parameters
@param forces force array (forces added)
@param energyByAtom atom energy
@param totalEnergy total energy
@return ReferenceForce::DefaultReturn
--------------------------------------------------------------------------------------- */
int ReferenceCustomGBIxn::calculateIxn(int numberOfAtoms, RealOpenMM** atomCoordinates, RealOpenMM** atomParameters,
void ReferenceCustomGBIxn::calculateIxn(int numberOfAtoms, RealOpenMM** atomCoordinates, RealOpenMM** atomParameters,
const vector<set<int> >& exclusions, map<string, double>& globalParameters, RealOpenMM** forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy) const {
RealOpenMM* totalEnergy) const {
// First calculate the computed values.
int numValues = valueTypes.size();
vector<vector<ComputedValue> > values(numValues);
vector<vector<RealOpenMM> > values(numValues);
for (int valueIndex = 0; valueIndex < numValues; valueIndex++) {
if (valueTypes[valueIndex] == OpenMM::CustomGBForce::SingleParticle)
calculateSingleParticleValue(valueIndex, numberOfAtoms, values, globalParameters, atomParameters);
......@@ -174,18 +146,25 @@ int ReferenceCustomGBIxn::calculateIxn(int numberOfAtoms, RealOpenMM** atomCoord
else
calculateParticlePairValue(valueIndex, numberOfAtoms, atomCoordinates, atomParameters, values, globalParameters, exclusions, false);
}
// Now calculate the energy and its derivates.
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, forces, totalEnergy);
calculateSingleParticleEnergyTerm(termIndex, numberOfAtoms, values, globalParameters, atomParameters, totalEnergy, dEdV);
else if (energyTypes[termIndex] == OpenMM::CustomGBForce::ParticlePair)
calculateParticlePairEnergyTerm(termIndex, numberOfAtoms, atomCoordinates, atomParameters, values, globalParameters, exclusions, true, forces, totalEnergy);
calculateParticlePairEnergyTerm(termIndex, numberOfAtoms, atomCoordinates, atomParameters, values, globalParameters, exclusions, true, forces, totalEnergy, dEdV);
else
calculateParticlePairEnergyTerm(termIndex, numberOfAtoms, atomCoordinates, atomParameters, values, globalParameters, exclusions, false, forces, totalEnergy);
calculateParticlePairEnergyTerm(termIndex, numberOfAtoms, atomCoordinates, atomParameters, values, globalParameters, exclusions, false, forces, totalEnergy, dEdV);
}
return ReferenceForce::DefaultReturn;
// Apply the chain rule to evaluate forces.
calculateChainRuleForces(numberOfAtoms, atomCoordinates, atomParameters, values, globalParameters, exclusions, forces, dEdV);
}
void ReferenceCustomGBIxn::calculateSingleParticleValue(int index, int numAtoms, vector<vector<ComputedValue> >& values,
void ReferenceCustomGBIxn::calculateSingleParticleValue(int index, int numAtoms, vector<vector<RealOpenMM> >& values,
const map<string, double>& globalParameters, RealOpenMM** atomParameters) const {
values[index].resize(numAtoms);
map<string, double> variables = globalParameters;
......@@ -193,23 +172,19 @@ void ReferenceCustomGBIxn::calculateSingleParticleValue(int index, int numAtoms,
for (int j = 0; j < (int) paramNames.size(); j++)
variables[paramNames[j]] = atomParameters[i][j];
for (int j = 0; j < index; j++)
variables[valueNames[j]] = values[j][i].value;
values[index][i].value = (RealOpenMM) valueExpressions[index].evaluate(variables);
for (int j = 0; j < index; j++) {
RealOpenMM scale = (RealOpenMM) valueDerivExpressions[index][j].evaluate(variables);
values[index][i].gradient[0] += scale*values[j][i].gradient[0];
values[index][i].gradient[1] += scale*values[j][i].gradient[1];
values[index][i].gradient[2] += scale*values[j][i].gradient[2];
}
variables[valueNames[j]] = values[j][i];
values[index][i] = (RealOpenMM) valueExpressions[index].evaluate(variables);
}
}
void ReferenceCustomGBIxn::calculateParticlePairValue(int index, int numAtoms, RealOpenMM** atomCoordinates, RealOpenMM** atomParameters,
vector<vector<ComputedValue> >& values, const map<string, double>& globalParameters, const vector<set<int> >& exclusions, bool useExclusions) const {
vector<vector<RealOpenMM> >& values, const map<string, double>& globalParameters, const vector<set<int> >& exclusions, bool useExclusions) const {
values[index].resize(numAtoms);
for (int i = 0; i < numAtoms; i++)
values[index][i].value = (RealOpenMM) 0.0;
values[index][i] = (RealOpenMM) 0.0;
if (cutoff) {
// Loop over all pairs in the neighbor list.
for (int i = 0; i < (int) neighborList->size(); i++) {
OpenMM::AtomPair pair = (*neighborList)[i];
if (useExclusions && exclusions[pair.first].find(pair.second) != exclusions[pair.first].end())
......@@ -219,6 +194,8 @@ void ReferenceCustomGBIxn::calculateParticlePairValue(int index, int numAtoms, R
}
}
else {
// Perform an O(N^2) loop over all atom pairs.
for (int i = 0; i < numAtoms; i++){
for (int j = i+1; j < numAtoms; j++ ){
if (useExclusions && exclusions[i].find(j) != exclusions[i].end())
......@@ -231,7 +208,7 @@ void ReferenceCustomGBIxn::calculateParticlePairValue(int index, int numAtoms, R
}
void ReferenceCustomGBIxn::calculateOnePairValue(int index, int atom1, int atom2, RealOpenMM** atomCoordinates, RealOpenMM** atomParameters,
const map<string, double>& globalParameters, vector<vector<ComputedValue> >& values) const {
const map<string, double>& globalParameters, vector<vector<RealOpenMM> >& values) const {
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (periodic)
ReferenceForce::getDeltaRPeriodic(atomCoordinates[atom2], atomCoordinates[atom1], periodicBoxSize, deltaR);
......@@ -247,65 +224,59 @@ void ReferenceCustomGBIxn::calculateOnePairValue(int index, int atom1, int atom2
}
variables["r"] = r;
for (int i = 0; i < index; i++) {
variables[particleValueNames[i*2]] = values[i][atom1].value;
variables[particleValueNames[i*2+1]] = values[i][atom2].value;
}
values[index][atom1].value += (RealOpenMM) valueExpressions[index].evaluate(variables);
RealOpenMM dVdR = (RealOpenMM) valueDerivExpressions[index][0].evaluate(variables);
RealOpenMM rinv = 1/r;
for (int i = 0; i < 3; i++)
values[index][atom1].gradient[i] += dVdR*deltaR[i]*rinv;
for (int i = 0; i < index; i++) {
RealOpenMM scale = (RealOpenMM) valueDerivExpressions[index][i+1].evaluate(variables);
values[index][atom1].gradient[0] += scale*values[i][atom1].gradient[0];
values[index][atom1].gradient[1] += scale*values[i][atom1].gradient[1];
values[index][atom1].gradient[2] += scale*values[i][atom1].gradient[2];
variables[particleValueNames[i*2]] = values[i][atom1];
variables[particleValueNames[i*2+1]] = values[i][atom2];
}
values[index][atom1] += (RealOpenMM) valueExpressions[index].evaluate(variables);
}
void ReferenceCustomGBIxn::calculateSingleParticleEnergyTerm(int index, int numAtoms, const vector<vector<ComputedValue> >& values,
const map<string, double>& globalParameters, RealOpenMM** atomParameters, RealOpenMM** forces, RealOpenMM* totalEnergy) const {
void ReferenceCustomGBIxn::calculateSingleParticleEnergyTerm(int index, int numAtoms, const vector<vector<RealOpenMM> >& values,
const map<string, double>& globalParameters, RealOpenMM** atomParameters, RealOpenMM* totalEnergy,
vector<vector<RealOpenMM> >& dEdV) const {
map<string, double> variables = globalParameters;
for (int i = 0; i < numAtoms; i++) {
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].value;
variables[valueNames[j]] = values[j][i];
if (totalEnergy != NULL)
*totalEnergy += (RealOpenMM) energyExpressions[index].evaluate(variables);
for (int j = 0; j < (int) valueNames.size(); j++) {
RealOpenMM scale = (RealOpenMM) energyDerivExpressions[index][j].evaluate(variables);
forces[i][0] -= scale*values[j][i].gradient[0];
forces[i][1] -= scale*values[j][i].gradient[1];
forces[i][2] -= scale*values[j][i].gradient[2];
}
for (int j = 0; j < (int) valueNames.size(); j++)
dEdV[j][i] += (RealOpenMM) energyDerivExpressions[index][j].evaluate(variables);
}
}
void ReferenceCustomGBIxn::calculateParticlePairEnergyTerm(int index, int numAtoms, RealOpenMM** atomCoordinates, RealOpenMM** atomParameters,
const vector<vector<ComputedValue> >& values, const map<string, double>& globalParameters, const vector<set<int> >& exclusions, bool useExclusions,
RealOpenMM** forces, RealOpenMM* totalEnergy) const {
const vector<vector<RealOpenMM> >& values, const map<string, double>& globalParameters, const vector<set<int> >& exclusions, bool useExclusions,
RealOpenMM** forces, RealOpenMM* totalEnergy, vector<vector<RealOpenMM> >& dEdV) const {
if (cutoff) {
// Loop over all pairs in the neighbor list.
for (int i = 0; i < (int) neighborList->size(); i++) {
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);
calculateOnePairEnergyTerm(index, pair.first, pair.second, atomCoordinates, atomParameters, globalParameters, values, forces, totalEnergy, dEdV);
}
}
else {
// Perform an O(N^2) loop over all atom pairs.
for (int i = 0; i < numAtoms; i++){
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);
calculateOnePairEnergyTerm(index, i, j, atomCoordinates, atomParameters, globalParameters, values, forces, totalEnergy, dEdV);
}
}
}
}
void ReferenceCustomGBIxn::calculateOnePairEnergyTerm(int index, int atom1, int atom2, RealOpenMM** atomCoordinates, RealOpenMM** atomParameters,
const map<string, double>& globalParameters, const vector<vector<ComputedValue> >& values, RealOpenMM** forces, RealOpenMM* totalEnergy) const {
const map<string, double>& globalParameters, const vector<vector<RealOpenMM> >& values, RealOpenMM** forces, RealOpenMM* totalEnergy,
vector<vector<RealOpenMM> >& dEdV) const {
// Compute the displacement.
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (periodic)
ReferenceForce::getDeltaRPeriodic(atomCoordinates[atom2], atomCoordinates[atom1], periodicBoxSize, deltaR);
......@@ -314,6 +285,9 @@ void ReferenceCustomGBIxn::calculateOnePairEnergyTerm(int index, int atom1, int
RealOpenMM r = deltaR[ReferenceForce::RIndex];
if (cutoff && r >= cutoffDistance)
return;
// Record variables for evaluating expressions.
map<string, double> variables = globalParameters;
for (int i = 0; i < paramNames.size(); i++) {
variables[particleParamNames[i*2]] = atomParameters[atom1][i];
......@@ -321,9 +295,12 @@ void ReferenceCustomGBIxn::calculateOnePairEnergyTerm(int index, int atom1, int
}
variables["r"] = r;
for (int i = 0; i < (int) valueNames.size(); i++) {
variables[particleValueNames[i*2]] = values[i][atom1].value;
variables[particleValueNames[i*2+1]] = values[i][atom2].value;
variables[particleValueNames[i*2]] = values[i][atom1];
variables[particleValueNames[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);
......@@ -333,13 +310,85 @@ void ReferenceCustomGBIxn::calculateOnePairEnergyTerm(int index, int atom1, int
forces[atom2][i] += dEdR*deltaR[i];
}
for (int i = 0; i < (int) valueNames.size(); i++) {
RealOpenMM scale1 = (RealOpenMM) energyDerivExpressions[index][2*i+1].evaluate(variables);
RealOpenMM scale2 = (RealOpenMM) energyDerivExpressions[index][2*i+2].evaluate(variables);
forces[atom1][0] -= scale1*values[i][atom1].gradient[0];
forces[atom1][1] -= scale1*values[i][atom1].gradient[1];
forces[atom1][2] -= scale1*values[i][atom1].gradient[2];
forces[atom2][0] += scale2*values[i][atom2].gradient[0];
forces[atom2][1] += scale2*values[i][atom2].gradient[1];
forces[atom2][2] += scale2*values[i][atom2].gradient[2];
dEdV[i][atom1] += (RealOpenMM) energyDerivExpressions[index][2*i+1].evaluate(variables);
dEdV[i][atom2] += (RealOpenMM) energyDerivExpressions[index][2*i+2].evaluate(variables);
}
}
void ReferenceCustomGBIxn::calculateChainRuleForces(int numAtoms, RealOpenMM** atomCoordinates, RealOpenMM** atomParameters,
const vector<vector<RealOpenMM> >& values, const map<string, double>& globalParameters,
const vector<set<int> >& exclusions, RealOpenMM** forces, vector<vector<RealOpenMM> >& dEdV) const {
bool useExclusions = (energyTypes[0] == OpenMM::CustomGBForce::ParticlePair);
if (cutoff) {
// Loop over all pairs in the neighbor list.
for (int i = 0; i < (int) neighborList->size(); i++) {
OpenMM::AtomPair pair = (*neighborList)[i];
if (useExclusions && exclusions[pair.first].find(pair.second) != exclusions[pair.first].end())
continue;
calculateOnePairChainRule(pair.first, pair.second, atomCoordinates, atomParameters, globalParameters, values, forces, dEdV);
calculateOnePairChainRule(pair.second, pair.first, atomCoordinates, atomParameters, globalParameters, values, forces, dEdV);
}
}
else {
// Perform an O(N^2) loop over all atom pairs.
for (int i = 0; i < numAtoms; i++){
for (int j = i+1; j < numAtoms; j++ ){
if (useExclusions && exclusions[i].find(j) != exclusions[i].end())
continue;
calculateOnePairChainRule(i, j, atomCoordinates, atomParameters, globalParameters, values, forces, dEdV);
calculateOnePairChainRule(j, i, atomCoordinates, atomParameters, globalParameters, values, forces, dEdV);
}
}
}
}
void ReferenceCustomGBIxn::calculateOnePairChainRule(int atom1, int atom2, RealOpenMM** atomCoordinates, RealOpenMM** atomParameters,
const map<string, double>& globalParameters, const vector<vector<RealOpenMM> >& values, RealOpenMM** forces,
vector<vector<RealOpenMM> >& dEdV) const {
// Compute the displacement.
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (periodic)
ReferenceForce::getDeltaRPeriodic(atomCoordinates[atom2], atomCoordinates[atom1], periodicBoxSize, deltaR);
else
ReferenceForce::getDeltaR(atomCoordinates[atom2], atomCoordinates[atom1], deltaR);
RealOpenMM r = deltaR[ReferenceForce::RIndex];
if (cutoff && r >= cutoffDistance)
return;
// Record variables for evaluating expressions.
map<string, double> variables = globalParameters;
for (int i = 0; i < paramNames.size(); i++) {
variables[particleParamNames[i*2]] = atomParameters[atom1][i];
variables[particleParamNames[i*2+1]] = atomParameters[atom2][i];
}
variables["r"] = r;
variables[particleValueNames[0]] = values[0][atom1];
variables[particleValueNames[1]] = values[0][atom2];
// Evaluate the derivative of each parameter with respect to position and apply forces.
RealOpenMM dVdR = (RealOpenMM) valueDerivExpressions[0].evaluate(variables);
RealOpenMM rinv = 1/r;
for (int i = 0; i < 3; i++) {
RealOpenMM f = dEdV[0][atom1]*dVdR*deltaR[i]*rinv;
forces[atom1][i] -= f;
forces[atom2][i] += f;
}
variables = globalParameters;
for (int i = 0; i < 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];
dVdR *= (RealOpenMM) valueDerivExpressions[i].evaluate(variables);
for (int j = 0; j < 3; j++) {
RealOpenMM f = dEdV[i][atom1]*dVdR*deltaR[j]*rinv;
forces[atom1][j] -= f;
forces[atom2][j] += f;
}
}
}
platforms/reference/src/SimTKReference/ReferenceCustomGBIxn.h
View file @ 4747b0f2
......@@ -44,7 +44,7 @@ class ReferenceCustomGBIxn {
RealOpenMM periodicBoxSize[3];
RealOpenMM cutoffDistance;
std::vector<Lepton::ExpressionProgram> valueExpressions;
std::vector<std::vector<Lepton::ExpressionProgram> > valueDerivExpressions;
std::vector<Lepton::ExpressionProgram> valueDerivExpressions;
std::vector<std::string> valueNames;
std::vector<OpenMM::CustomGBForce::ComputationType> valueTypes;
std::vector<Lepton::ExpressionProgram> energyExpressions;
......@@ -53,57 +53,164 @@ class ReferenceCustomGBIxn {
std::vector<OpenMM::CustomGBForce::ComputationType> energyTypes;
std::vector<std::string> particleParamNames;
std::vector<std::string> particleValueNames;
struct ComputedValue {
RealOpenMM value;
RealOpenMM gradient[3];
ComputedValue() {
value = (RealOpenMM) 0;
gradient[0] = (RealOpenMM) 0;
gradient[1] = (RealOpenMM) 0;
gradient[2] = (RealOpenMM) 0;
}
};
/**---------------------------------------------------------------------------------------
Calculate custom pair ixn between two atoms
Calculate a computed value of type SingleParticle
@param atom1 the index of the first atom
@param atom2 the index of the second atom
@param atomCoordinates atom coordinates
@param atomParameters atomParameters[atomIndex][parameterIndex]
@param forces force array (forces added)
@param energyByAtom atom energy
@param totalEnergy total energy
@param index the index of the value to compute
@param numAtoms number of atoms
@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<ComputedValue> >& values,
void calculateSingleParticleValue(int index, int numAtoms, std::vector<std::vector<RealOpenMM> >& values,
const std::map<std::string, double>& globalParameters, RealOpenMM** atomParameters) const;
/**---------------------------------------------------------------------------------------
Calculate a computed value that is based on particle pairs
@param index the index of the value to compute
@param numAtoms number of atoms
@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
--------------------------------------------------------------------------------------- */
void calculateParticlePairValue(int index, int numAtoms, RealOpenMM** atomCoordinates, RealOpenMM** atomParameters,
std::vector<std::vector<ComputedValue> >& values,
std::vector<std::vector<RealOpenMM> >& values,
const std::map<std::string, double>& globalParameters,
const std::vector<std::set<int> >& exclusions, bool useExclusions) const;
/**---------------------------------------------------------------------------------------
Evaluate a single atom pair as part of calculating a computed value
@param index the index of the value to compute
@param atom1 the index of the first atom in the pair
@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, RealOpenMM** atomCoordinates, RealOpenMM** atomParameters,
const std::map<std::string, double>& globalParameters,
std::vector<std::vector<ComputedValue> >& values) const;
std::vector<std::vector<RealOpenMM> >& values) const;
/**---------------------------------------------------------------------------------------
Calculate an energy term of type SingleParticle
@param index the index of the value to compute
@param numAtoms number of atoms
@param values the vector containing computed values
@param globalParameters the values of global parameters
@param atomParameters atomParameters[atomIndex][paramterIndex]
@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<ComputedValue> >& values,
--------------------------------------------------------------------------------------- */
void calculateSingleParticleEnergyTerm(int index, int numAtoms, const std::vector<std::vector<RealOpenMM> >& values,
const std::map<std::string, double>& globalParameters, RealOpenMM** atomParameters,
RealOpenMM** forces, RealOpenMM* totalEnergy) const;
RealOpenMM* totalEnergy, std::vector<std::vector<RealOpenMM> >& dEdV) const;
/**---------------------------------------------------------------------------------------
Calculate an energy term that is based on particle pairs
@param index the index of the term to compute
@param numAtoms number of atoms
@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
@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 calculateParticlePairEnergyTerm(int index, int numAtoms, RealOpenMM** atomCoordinates, RealOpenMM** atomParameters,
const std::vector<std::vector<ComputedValue> >& values,
const std::vector<std::vector<RealOpenMM> >& values,
const std::map<std::string, double>& globalParameters,
const std::vector<std::set<int> >& exclusions, bool useExclusions,
RealOpenMM** forces, RealOpenMM* totalEnergy) const;
RealOpenMM** forces, RealOpenMM* totalEnergy, std::vector<std::vector<RealOpenMM> >& dEdV) const;
/**---------------------------------------------------------------------------------------
Evaluate a single atom pair as part of calculating an energy term
@param index the index of the term to compute
@param atom1 the index of the first atom in the pair
@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
@param dEdV the derivative of energy with respect to computed values is stored in this
--------------------------------------------------------------------------------------- */
void calculateOnePairEnergyTerm(int index, int atom1, int atom2, RealOpenMM** atomCoordinates, RealOpenMM** atomParameters,
const std::map<std::string, double>& globalParameters,
const std::vector<std::vector<ComputedValue> >& values,
RealOpenMM** forces, RealOpenMM* totalEnergy) const;
const std::vector<std::vector<RealOpenMM> >& values,
RealOpenMM** forces, RealOpenMM* totalEnergy, std::vector<std::vector<RealOpenMM> >& dEdV) const;
/**---------------------------------------------------------------------------------------
Apply the chain rule to compute forces on atoms
@param numAtoms number of atoms
@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
--------------------------------------------------------------------------------------- */
void calculateChainRuleForces(int numAtoms, RealOpenMM** atomCoordinates, RealOpenMM** atomParameters,
const std::vector<std::vector<RealOpenMM> >& values,
const std::map<std::string, double>& globalParameters,
const std::vector<std::set<int> >& exclusions,
RealOpenMM** forces, std::vector<std::vector<RealOpenMM> >& dEdV) const;
/**---------------------------------------------------------------------------------------
Evaluate a single atom pair as part of applying the chain rule
@param atom1 the index of the first atom in the pair
@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
--------------------------------------------------------------------------------------- */
void calculateOnePairChainRule(int atom1, int atom2, RealOpenMM** atomCoordinates, RealOpenMM** atomParameters,
const std::map<std::string, double>& globalParameters,
const std::vector<std::vector<RealOpenMM> >& values,
RealOpenMM** forces, std::vector<std::vector<RealOpenMM> >& dEdV) const;
public:
......@@ -114,7 +221,7 @@ class ReferenceCustomGBIxn {
--------------------------------------------------------------------------------------- */
ReferenceCustomGBIxn(const std::vector<Lepton::ExpressionProgram>& valueExpressions,
const std::vector<std::vector<Lepton::ExpressionProgram> >& valueDerivExpressions,
const std::vector<Lepton::ExpressionProgram>& valueDerivExpressions,
const std::vector<std::string>& valueNames,
const std::vector<OpenMM::CustomGBForce::ComputationType>& valueTypes,
const std::vector<Lepton::ExpressionProgram>& energyExpressions,
......@@ -137,11 +244,9 @@ class ReferenceCustomGBIxn {
@param distance the cutoff distance
@param neighbors the neighbor list to use
@return ReferenceForce::DefaultReturn
--------------------------------------------------------------------------------------- */
int setUseCutoff( RealOpenMM distance, const OpenMM::NeighborList& neighbors );
void setUseCutoff( RealOpenMM distance, const OpenMM::NeighborList& neighbors );
/**---------------------------------------------------------------------------------------
......@@ -151,36 +256,26 @@ class ReferenceCustomGBIxn {
@param boxSize the X, Y, and Z widths of the periodic box
@return ReferenceForce::DefaultReturn
--------------------------------------------------------------------------------------- */
int setPeriodic( RealOpenMM* boxSize );
void setPeriodic( RealOpenMM* boxSize );
/**---------------------------------------------------------------------------------------
Calculate custom pair ixn
Calculate custom GB ixn
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param atomParameters atom parameters (charges, c6, c12, ...) atomParameters[atomIndex][paramterIndex]
@param exclusions atom exclusion indices exclusions[atomIndex][atomToExcludeIndex]
exclusions[atomIndex][0] = number of exclusions
exclusions[atomIndex][1-no.] = atom indices of atoms to excluded from
interacting w/ atom atomIndex
@param fixedParameters non atom parameters (not currently used)
@param atomParameters atomParameters[atomIndex][paramterIndex]
@param exclusions exclusions[i] is the set of excluded indices for atom i
@param globalParameters the values of global parameters
@param forces force array (forces added)
@param energyByAtom atom energy
@param totalEnergy total energy
@return ReferenceForce::DefaultReturn
--------------------------------------------------------------------------------------- */
int calculateIxn(int numberOfAtoms, RealOpenMM** atomCoordinates, RealOpenMM** atomParameters, const std::vector<std::set<int> >& exclusions,
std::map<std::string, double>& globalParameters, RealOpenMM** forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy) const;
void calculateIxn(int numberOfAtoms, RealOpenMM** atomCoordinates, RealOpenMM** atomParameters, const std::vector<std::set<int> >& exclusions,
std::map<std::string, double>& globalParameters, RealOpenMM** forces, RealOpenMM* totalEnergy) const;
// ---------------------------------------------------------------------------------------
......
platforms/reference/tests/TestReferenceCustomGBForce.cpp
View file @ 4747b0f2
......@@ -50,8 +50,8 @@ using namespace std;
const double TOL = 1e-5;
void testOBC() {
const int numMolecules = 100;
void testOBC(GBSAOBCForce::NonbondedMethod obcMethod, CustomGBForce::NonbondedMethod customMethod) {
const int numMolecules = 70;
const int numParticles = numMolecules*2;
const double boxSize = 10.0;
ReferencePlatform platform;
......@@ -64,8 +64,12 @@ void testOBC() {
standardSystem.addParticle(1.0);
customSystem.addParticle(1.0);
}
standardSystem.setPeriodicBoxVectors(Vec3(boxSize, 0.0, 0.0), Vec3(0.0, boxSize, 0.0), Vec3(0.0, 0.0, boxSize));
customSystem.setPeriodicBoxVectors(Vec3(boxSize, 0.0, 0.0), Vec3(0.0, boxSize, 0.0), Vec3(0.0, 0.0, boxSize));
GBSAOBCForce* obc = new GBSAOBCForce();
CustomGBForce* custom = new CustomGBForce();
obc->setCutoffDistance(2.0);
custom->setCutoffDistance(2.0);
custom->addPerParticleParameter("q");
custom->addPerParticleParameter("radius");
custom->addPerParticleParameter("scale");
......@@ -115,31 +119,31 @@ void testOBC() {
velocities[2*i] = Vec3(genrand_real2(), genrand_real2(), genrand_real2());
velocities[2*i+1] = Vec3(genrand_real2(), genrand_real2(), genrand_real2());
}
obc->setNonbondedMethod(GBSAOBCForce::NoCutoff);
custom->setNonbondedMethod(CustomGBForce::NoCutoff);
obc->setNonbondedMethod(obcMethod);
custom->setNonbondedMethod(customMethod);
standardSystem.addForce(obc);
customSystem.addForce(custom);
VerletIntegrator integrator1(0.01);
VerletIntegrator integrator2(0.01);
Context context1(standardSystem, integrator1, platform);
Context context2(customSystem, integrator2, platform);
context1.setPositions(positions);
context2.setPositions(positions);
context1.setVelocities(velocities);
context2.setVelocities(velocities);
State state1 = context1.getState(State::Forces | State::Energy);
Context context2(customSystem, integrator2, platform);
context2.setPositions(positions);
context2.setVelocities(velocities);
State state2 = context2.getState(State::Forces | State::Energy);
for (int i = 0; i < numParticles; i++)
std::cout << state1.getForces()[i]<< state2.getForces()[i]<< std::endl;
ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-4);
for (int i = 0; i < numParticles; i++) {
ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-3);
ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-4);
}
}
int main() {
try {
testOBC();
testOBC(GBSAOBCForce::NoCutoff, CustomGBForce::NoCutoff);
testOBC(GBSAOBCForce::CutoffNonPeriodic, CustomGBForce::CutoffNonPeriodic);
testOBC(GBSAOBCForce::CutoffPeriodic, CustomGBForce::CutoffPeriodic);
}
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