Skip to content

GitLab

Commit be137f98 authored by Peter Eastman's avatar Peter Eastman
Browse files

Redesigned CustomHbondForce to be more general

parent 45c23132
Show whitespace changes
Inline Side-by-side
Showing with 447 additions and 255 deletions
openmmapi/include/openmm/CustomHbondForce.h
View file @ be137f98
......@@ -44,30 +44,25 @@ namespace OpenMM {
/**
* This class supports a wide variety of energy functions used to represent hydrogen bonding. It computes
* interactions between "donor" particle pairs and "acceptor" particle pairs (so each interaction involves four
* particles in all). Typically a donor pair consists of a hydrogen atom and the atom it is bonded to,
* and an acceptor pair consists of a negatively charged atom and the atom it is bonded to. Within each pair,
* one particle is designated as the "primary particle" (usually the hydrogen atom in a donor, and the negatively
* charged atom in an acceptor), and the other particle is called the "secondary particle".
*
* We refer to the primary and secondary donor particles as D1 and D2, and the primary and secondary acceptor
* particles as A1 and A2. For each donor and each acceptor, CustomHbondForce computes the following values:
*
* r: the distance between D1 and A1
* theta: the angle formed by D2-D1-A1, measured in radians
* psi: the angle formed by A2-A1-D1, measured in radians
* chi: the dihedral angle formed by D2-D1-A1-A2, measured in radians. The angle is defined to be zero when
* D2 and A2 are on the same side of the bond formed by D1 and A1.
*
* The user then specifies an algebraic expression which gives the energy of each interaction as a function
* of r, theta, psi, and chi. The expression also may involve tabulated functions, and may depend on arbitrary
* interactions between "donor" particle groups and "acceptor" particle groups, where each group may include
* up to three particles. Typically a donor group consists of a hydrogen atom and the atoms it is bonded to,
* and an acceptor group consists of a negatively charged atom and the atoms it is bonded to.
*
* We refer to the particles in a donor group as d1, d2 and d3, and the particles in an acceptor group as
* a1, a2, and a3. For each donor and each acceptor, CustomHbondForce evaluates a user supplied algebraic
* expression to determine the interaction energy. The expression may depend on arbitrary distances, angles,
* and dihedral angles defined by any of the six particles involved. The function distance(p1, p2) is the distance
* between the particles p1 and p2 (where "p1" and "p2" should be replaced by the names of the actual particles
* to calculate the distance between), angle(p1, p2, p3) is the angle formed by the three specified particles,
* and dihedral(p1, p2, p3, p4) is the dihedral angle formed by the four specified particles.
*
* The expression also may involve tabulated functions, and may depend on arbitrary
* global, per-donor, and per-acceptor parameters. It also optionally supports periodic boundary conditions
* and cutoffs for long range interactions.
*
* To use this class, create a CustomHbondForce object, passing an algebraic expression to the constructor
* that defines the interaction energy between each donor and acceptor. The expression may depend on r, theta,
* psi, and chi, as well as on any parameters you choose. Then call addPerDonorParameter() to define per-donor
* parameters, addPerAcceptorParameter to define per-acceptor parameters, and addGlobalParameter() to define
* that defines the interaction energy between each donor and acceptor. Then call addPerDonorParameter() to define per-donor
* parameters, addPerAcceptorParameter() to define per-acceptor parameters, and addGlobalParameter() to define
* global parameters. The values of per-donor and per-acceptor parameters are specified as part of the system
* definition, while values of global parameters may be modified during a simulation by calling Context::setParameter().
*
......@@ -80,9 +75,9 @@ namespace OpenMM {
* that are bonded to each other.
*
* As an example, the following code creates a CustomHbondForce that implements a simple harmonic potential
* to keep r and theta near ideal values:
* to keep the distance between a1 and d1, and the angle formed by a1-d1-d2, near ideal values:
*
* <tt>CustomHbondForce* force = new CustomHbondForce("k*(r-r0)^2*(theta-theta0)^2");</tt>
* <tt>CustomHbondForce* force = new CustomHbondForce("k*(distance(a1,d1)-r0)^2*(angle(a1,d1,d2)-theta0)^2");</tt>
*
* This force depends on three parameters: k, r0, and theta0. The following code defines these as per-donor parameters:
*
......@@ -125,7 +120,7 @@ public:
* Create a CustomHbondForce.
*
* @param energy an algebraic expression giving the interaction energy between a donor as a function
* of the distance r and the angles theta, psi, and chi, as well as any global, per-donor, and
* of inter-particle distances, angles, and dihedrals, as well as any global, per-donor, and
* per-acceptor parameters
*/
CustomHbondForce(const std::string& energy);
......@@ -142,7 +137,7 @@ public:
return acceptors.size();
}
/**
* Get the number of particle pairs whose interactions should be excluded.
* Get the number of donor-acceptor pairs whose interactions should be excluded.
*/
int getNumExclusions() const {
return exclusions.size();
......@@ -188,13 +183,15 @@ public:
*/
void setNonbondedMethod(NonbondedMethod method);
/**
* Get the cutoff distance (in nm) being used. If the NonbondedMethod in use is NoCutoff, this value will have no effect.
* Get the cutoff distance (in nm) being used. All interactions for which the distance between d1 and a1
* is greater than the cutoff will be ignored. If the NonbondedMethod in use is NoCutoff, this value will have no effect.
*
* @return the cutoff distance, measured in nm
*/
double getCutoffDistance() const;
/**
* Set the cutoff distance (in nm) being used. If the NonbondedMethod in use is NoCutoff, this value will have no effect.
* Set the cutoff distance (in nm) being used. All interactions for which the distance between d1 and a1
* is greater than the cutoff will be ignored. If the NonbondedMethod in use is NoCutoff, this value will have no effect.
*
* @param distance the cutoff distance, measured in nm
*/
......@@ -278,59 +275,77 @@ public:
*/
void setGlobalParameterDefaultValue(int index, double defaultValue);
/**
* Add a donor pair to the force
* Add a donor group to the force
*
* @param primaryParticle the index of the primary particle for this donor pair
* @param secondaryParticle the index of the secondary particle for this donor pair
* @param d1 the index of the first particle for this donor group
* @param d2 the index of the second particle for this donor group. If the group only
* includes one particle, this must be -1.
* @param d3 the index of the third particle for this donor group. If the group includes
* less than three particles, this must be -1.
* @param parameters the list of per-donor parameter values for the new donor
* @return the index of the donor that was added
*/
int addDonor(int primaryParticle, int secondaryParticle, const std::vector<double>& parameters);
int addDonor(int d1, int d2, int d3, const std::vector<double>& parameters);
/**
* Get the properties of a donor pair.
* Get the properties of a donor group.
*
* @param index the index of the donor pair to get
* @param primaryParticle the index of the primary particle for this donor pair
* @param secondaryParticle the index of the secondary particle for this donor pair
* @param index the index of the donor group to get
* @param d1 the index of the first particle for this donor group
* @param d2 the index of the second particle for this donor group. If the group only
* includes one particle, this will be -1.
* @param d3 the index of the third particle for this donor group. If the group includes
* less than three particles, this will be -1.
* @param parameters the list of per-donor parameter values for the new donor
*/
void getDonorParameters(int index, int& primaryParticle, int& secondaryParticle, std::vector<double>& parameters) const;
void getDonorParameters(int index, int& d1, int& d2, int& d3, std::vector<double>& parameters) const;
/**
* Set the properties of a donor pair.
* Set the properties of a donor group.
*
* @param index the index of the donor pair to get
* @param primaryParticle the index of the primary particle for this donor pair
* @param secondaryParticle the index of the secondary particle for this donor pair
* @param index the index of the donor group to set
* @param d1 the index of the first particle for this donor group
* @param d2 the index of the second particle for this donor group. If the group only
* includes one particle, this must be -1.
* @param d3 the index of the third particle for this donor group. If the group includes
* less than three particles, this must be -1.
* @param parameters the list of per-donor parameter values for the new donor
*/
void setDonorParameters(int index, int primaryParticle, int secondaryParticle, const std::vector<double>& parameters);
void setDonorParameters(int index, int d1, int d2, int d3, const std::vector<double>& parameters);
/**
* Add an acceptor pair to the force
* Add an acceptor group to the force
*
* @param primaryParticle the index of the primary particle for this acceptor pair
* @param secondaryParticle the index of the secondary particle for this acceptor pair
* @param a1 the index of the first particle for this acceptor group
* @param a2 the index of the second particle for this acceptor group. If the group only
* includes one particle, this must be -1.
* @param a3 the index of the third particle for this acceptor group. If the group includes
* less than three particles, this must be -1.
* @param parameters the list of per-acceptor parameter values for the new acceptor
* @return the index of the acceptor that was added
*/
int addAcceptor(int primaryParticle, int secondaryParticle, const std::vector<double>& parameters);
int addAcceptor(int a1, int a2, int a3, const std::vector<double>& parameters);
/**
* Get the properties of an acceptor pair.
* Get the properties of an acceptor group.
*
* @param index the index of the acceptor pair to get
* @param primaryParticle the index of the primary particle for this acceptor pair
* @param secondaryParticle the index of the secondary particle for this acceptor pair
* @param index the index of the acceptor group to get
* @param a1 the index of the first particle for this acceptor group
* @param a2 the index of the second particle for this acceptor group. If the group only
* includes one particle, this will be -1.
* @param a3 the index of the third particle for this acceptor group. If the group includes
* less than three particles, this will be -1.
* @param parameters the list of per-acceptor parameter values for the new acceptor
*/
void getAcceptorParameters(int index, int& primaryParticle, int& secondaryParticle, std::vector<double>& parameters) const;
void getAcceptorParameters(int index, int& a1, int& a2, int& a3, std::vector<double>& parameters) const;
/**
* Set the properties of an acceptor pair.
* Set the properties of an acceptor group.
*
* @param index the index of the acceptor pair to get
* @param primaryParticle the index of the primary particle for this acceptor pair
* @param secondaryParticle the index of the secondary particle for this acceptor pair
* @param index the index of the acceptor group to set
* @param a1 the index of the first particle for this acceptor group
* @param a2 the index of the second particle for this acceptor group. If the group only
* includes one particle, this must be -1.
* @param a3 the index of the third particle for this acceptor group. If the group includes
* less than three particles, this must be -1.
* @param parameters the list of per-acceptor parameter values for the new acceptor
*/
void setAcceptorParameters(int index, int primaryParticle, int secondaryParticle, const std::vector<double>& parameters);
void setAcceptorParameters(int index, int a1, int a2, int a3, const std::vector<double>& parameters);
/**
* Add a donor-acceptor pair to the list of interactions that should be excluded.
*
......@@ -397,7 +412,7 @@ public:
protected:
ForceImpl* createImpl();
private:
class PairInfo;
class GroupInfo;
class PerPairParameterInfo;
class GlobalParameterInfo;
class ExclusionInfo;
......@@ -408,8 +423,8 @@ private:
std::vector<PerPairParameterInfo> donorParameters;
std::vector<PerPairParameterInfo> acceptorParameters;
std::vector<GlobalParameterInfo> globalParameters;
std::vector<PairInfo> donors;
std::vector<PairInfo> acceptors;
std::vector<GroupInfo> donors;
std::vector<GroupInfo> acceptors;
std::vector<ExclusionInfo> exclusions;
std::vector<FunctionInfo> functions;
};
......@@ -418,14 +433,14 @@ private:
* This is an internal class used to record information about a donor or acceptor.
* @private
*/
class CustomHbondForce::PairInfo {
class CustomHbondForce::GroupInfo {
public:
std::vector<double> parameters;
int primaryParticle, secondaryParticle;
PairInfo() : primaryParticle(-1), secondaryParticle(-1) {
int p1, p2, p3;
GroupInfo() : p1(-1), p2(-1), p3(-1) {
}
PairInfo(int primaryParticle, int secondaryParticle, const std::vector<double>& parameters) :
primaryParticle(primaryParticle), secondaryParticle(secondaryParticle), parameters(parameters) {
GroupInfo(int p1, int p2, int p3, const std::vector<double>& parameters) :
p1(p1), p2(p2), p3(p3), parameters(parameters) {
}
};
......
openmmapi/include/openmm/internal/CustomHbondForceImpl.h
View file @ be137f98
......@@ -35,6 +35,8 @@
#include "ForceImpl.h"
#include "openmm/CustomHbondForce.h"
#include "openmm/Kernel.h"
#include "lepton/ExpressionTreeNode.h"
#include "lepton/ParsedExpression.h"
#include <utility>
#include <map>
#include <string>
......@@ -60,7 +62,27 @@ public:
double calcEnergy(ContextImpl& context);
std::map<std::string, double> getDefaultParameters();
std::vector<std::string> getKernelNames();
/**
* This is a utility routine that parses the energy expression, identifies the angles and dihedrals
* in it, and replaces them with variables. The particle indices returned in the maps are defined
* as follows: 0=a1, 1=a2, 2=a3, 3=d1, 4=d2, 5=d3.
*
* @param force the CustomHbondForce to process
* @param distances on exist, this will contain an entry for each distance used in the expression. The key is the name
* of the corresponding variable, and the value is the list of particle indices.
* @param angles on exist, this will contain an entry for each angle used in the expression. The key is the name
* of the corresponding variable, and the value is the list of particle indices.
* @param dihedrals on exist, this will contain an entry for each dihedral used in the expression. The key is the name
* of the corresponding variable, and the value is the list of particle indices.
* @return a Parsed expression for the energy
*/
static Lepton::ParsedExpression prepareExpression(const CustomHbondForce& force, std::map<std::string, std::vector<int> >& distances,
std::map<std::string, std::vector<int> >& angles, std::map<std::string, std::vector<int> >& dihedrals);
private:
class FunctionPlaceholder;
static Lepton::ExpressionTreeNode replaceFunctions(const Lepton::ExpressionTreeNode& node, std::map<std::string, int> atoms,
std::map<std::string, std::vector<int> >& distances, std::map<std::string, std::vector<int> >& angles,
std::map<std::string, std::vector<int> >& dihedrals);
CustomHbondForce& owner;
Kernel kernel;
};
......
openmmapi/src/CustomHbondForce.cpp
View file @ be137f98
......@@ -120,37 +120,41 @@ void CustomHbondForce::setGlobalParameterDefaultValue(int index, double defaultV
globalParameters[index].defaultValue = defaultValue;
}
int CustomHbondForce::addDonor(int primaryParticle, int secondaryParticle, const vector<double>& parameters) {
donors.push_back(PairInfo(primaryParticle, secondaryParticle, parameters));
int CustomHbondForce::addDonor(int d1, int d2, int d3, const vector<double>& parameters) {
donors.push_back(GroupInfo(d1, d2, d3, parameters));
return donors.size()-1;
}
void CustomHbondForce::getDonorParameters(int index, int& primaryParticle, int& secondaryParticle, std::vector<double>& parameters) const {
primaryParticle = donors[index].primaryParticle;
secondaryParticle = donors[index].secondaryParticle;
void CustomHbondForce::getDonorParameters(int index, int& d1, int& d2, int& d3, std::vector<double>& parameters) const {
d1 = donors[index].p1;
d2 = donors[index].p2;
d3 = donors[index].p3;
parameters = donors[index].parameters;
}
void CustomHbondForce::setDonorParameters(int index, int primaryParticle, int secondaryParticle, const vector<double>& parameters) {
donors[index].primaryParticle = primaryParticle;
donors[index].secondaryParticle = secondaryParticle;
void CustomHbondForce::setDonorParameters(int index, int d1, int d2, int d3, const vector<double>& parameters) {
donors[index].p1 = d1;
donors[index].p2 = d2;
donors[index].p3 = d3;
donors[index].parameters = parameters;
}
int CustomHbondForce::addAcceptor(int primaryParticle, int secondaryParticle, const vector<double>& parameters) {
acceptors.push_back(PairInfo(primaryParticle, secondaryParticle, parameters));
int CustomHbondForce::addAcceptor(int a1, int a2, int a3, const vector<double>& parameters) {
acceptors.push_back(GroupInfo(a1, a2, a3, parameters));
return acceptors.size()-1;
}
void CustomHbondForce::getAcceptorParameters(int index, int& primaryParticle, int& secondaryParticle, std::vector<double>& parameters) const {
primaryParticle = acceptors[index].primaryParticle;
secondaryParticle = acceptors[index].secondaryParticle;
void CustomHbondForce::getAcceptorParameters(int index, int& a1, int& a2, int& a3, std::vector<double>& parameters) const {
a1 = acceptors[index].p1;
a2 = acceptors[index].p2;
a3 = acceptors[index].p3;
parameters = acceptors[index].parameters;
}
void CustomHbondForce::setAcceptorParameters(int index, int primaryParticle, int secondaryParticle, const vector<double>& parameters) {
acceptors[index].primaryParticle = primaryParticle;
acceptors[index].secondaryParticle = secondaryParticle;
void CustomHbondForce::setAcceptorParameters(int index, int a1, int a2, int a3, const vector<double>& parameters) {
acceptors[index].p1 = a1;
acceptors[index].p2 = a2;
acceptors[index].p3 = a3;
acceptors[index].parameters = parameters;
}
......
openmmapi/src/CustomHbondForceImpl.cpp
View file @ be137f98
......@@ -33,9 +33,15 @@
#include "openmm/internal/ContextImpl.h"
#include "openmm/internal/CustomHbondForceImpl.h"
#include "openmm/kernels.h"
#include "lepton/Operation.h"
#include "lepton/Parser.h"
#include <sstream>
using namespace OpenMM;
using Lepton::CustomFunction;
using Lepton::ExpressionTreeNode;
using Lepton::Operation;
using Lepton::ParsedExpression;
using std::map;
using std::pair;
using std::vector;
......@@ -43,6 +49,28 @@ using std::set;
using std::string;
using std::stringstream;
/**
* This class serves as a placeholder for angles and dihedrals in expressions.
*/
class CustomHbondForceImpl::FunctionPlaceholder : public CustomFunction {
public:
int numArguments;
FunctionPlaceholder(int numArguments) : numArguments(numArguments) {
}
int getNumArguments() const {
return numArguments;
}
double evaluate(const double* arguments) const {
return 0.0;
}
double evaluateDerivative(const double* arguments, const int* derivOrder) const {
return 0.0;
}
CustomFunction* clone() const {
return new FunctionPlaceholder(numArguments);
}
};
CustomHbondForceImpl::CustomHbondForceImpl(CustomHbondForce& owner) : owner(owner) {
}
......@@ -59,18 +87,24 @@ void CustomHbondForceImpl::initialize(ContextImpl& context) {
vector<double> parameters;
int numDonorParameters = owner.getNumPerDonorParameters();
for (int i = 0; i < owner.getNumDonors(); i++) {
int primary, secondary;
owner.getDonorParameters(i, primary, secondary, parameters);
if (primary < 0 || primary >= system.getNumParticles()) {
int d1, d2, d3;
owner.getDonorParameters(i, d1, d2, d3, parameters);
if (d1 < 0 || d1 >= system.getNumParticles()) {
stringstream msg;
msg << "CustomHbondForce: Illegal particle index for a donor: ";
msg << d1;
throw OpenMMException(msg.str());
}
if (d2 < -1 || d2 >= system.getNumParticles()) {
stringstream msg;
msg << "CustomHbondForce: Illegal particle index for a donor: ";
msg << primary;
msg << d2;
throw OpenMMException(msg.str());
}
if (secondary < 0 || secondary >= system.getNumParticles()) {
if (d3 < -1 || d3 >= system.getNumParticles()) {
stringstream msg;
msg << "CustomHbondForce: Illegal particle index for a donor: ";
msg << secondary;
msg << d3;
throw OpenMMException(msg.str());
}
if (parameters.size() != numDonorParameters) {
......@@ -82,18 +116,24 @@ void CustomHbondForceImpl::initialize(ContextImpl& context) {
}
int numAcceptorParameters = owner.getNumPerAcceptorParameters();
for (int i = 0; i < owner.getNumAcceptors(); i++) {
int primary, secondary;
owner.getAcceptorParameters(i, primary, secondary, parameters);
if (primary < 0 || primary >= system.getNumParticles()) {
int a1, a2, a3;
owner.getAcceptorParameters(i, a1, a2, a3, parameters);
if (a1 < 0 || a1 >= system.getNumParticles()) {
stringstream msg;
msg << "CustomHbondForce: Illegal particle index for an acceptor: ";
msg << a1;
throw OpenMMException(msg.str());
}
if (a2 < -1 || a2 >= system.getNumParticles()) {
stringstream msg;
msg << "CustomHbondForce: Illegal particle index for an acceptor: ";
msg << primary;
msg << a2;
throw OpenMMException(msg.str());
}
if (secondary < 0 || secondary >= system.getNumParticles()) {
if (a3 < -1 || a3 >= system.getNumParticles()) {
stringstream msg;
msg << "CustomHbondForce: Illegal particle index for an acceptor: ";
msg << secondary;
msg << a3;
throw OpenMMException(msg.str());
}
if (parameters.size() != numAcceptorParameters) {
......@@ -158,3 +198,81 @@ map<string, double> CustomHbondForceImpl::getDefaultParameters() {
parameters[owner.getGlobalParameterName(i)] = owner.getGlobalParameterDefaultValue(i);
return parameters;
}
ParsedExpression CustomHbondForceImpl::prepareExpression(const CustomHbondForce& force, map<string, vector<int> >& distances,
map<string, vector<int> >& angles, map<string, vector<int> >& dihedrals) {
CustomHbondForceImpl::FunctionPlaceholder custom(1);
CustomHbondForceImpl::FunctionPlaceholder distance(2);
CustomHbondForceImpl::FunctionPlaceholder angle(3);
CustomHbondForceImpl::FunctionPlaceholder dihedral(4);
map<string, CustomFunction*> functions;
functions["distance"] = &distance;
functions["angle"] = &angle;
functions["dihedral"] = &dihedral;
for (int i = 0; i < force.getNumFunctions(); i++) {
string name;
vector<double> values;
double min, max;
bool interpolating;
force.getFunctionParameters(i, name, values, min, max, interpolating);
functions[name] = &custom;
}
ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction(), functions);
map<string, int> atoms;
atoms["a1"] = 0;
atoms["a2"] = 1;
atoms["a3"] = 2;
atoms["d1"] = 3;
atoms["d2"] = 4;
atoms["d3"] = 5;
return ParsedExpression(replaceFunctions(expression.getRootNode(), atoms, distances, angles, dihedrals)).optimize();
}
ExpressionTreeNode CustomHbondForceImpl::replaceFunctions(const ExpressionTreeNode& node, map<string, int> atoms,
map<string, vector<int> >& distances, map<string, vector<int> >& angles, map<string, vector<int> >& dihedrals) {
const Operation& op = node.getOperation();
if (op.getId() != Operation::CUSTOM || op.getNumArguments() < 2)
{
// This is not an angle or dihedral, so process its children.
vector<ExpressionTreeNode> children;
for (int i = 0; i < (int) node.getChildren().size(); i++)
children.push_back(replaceFunctions(node.getChildren()[i], atoms, distances, angles, dihedrals));
return ExpressionTreeNode(op.clone(), children);
}
const Operation::Custom& custom = static_cast<const Operation::Custom&>(op);
// Identify the atoms this term is based on.
int numArgs = custom.getNumArguments();
vector<int> indices(numArgs);
for (int i = 0; i < numArgs; i++) {
map<string, int>::const_iterator iter = atoms.find(node.getChildren()[i].getOperation().getName());
if (iter == atoms.end())
throw OpenMMException("CustomHBondForce: Unknown particle '"+node.getChildren()[i].getOperation().getName()+"'");
indices[i] = iter->second;
}
// Select a name for the variable and add it to the appropriate map.
stringstream variable;
if (numArgs == 2)
variable << "distance";
else if (numArgs == 3)
variable << "angle";
else
variable << "dihedral";
for (int i = 0; i < numArgs; i++)
variable << indices[i];
string name = variable.str();
if (numArgs == 2)
distances[name] = indices;
else if (numArgs == 3)
angles[name] = indices;
else
dihedrals[name] = indices;
// Return a new node that represents it as a simple variable.
return ExpressionTreeNode(new Operation::Variable(name));
}
platforms/reference/src/ReferenceKernels.cpp
View file @ be137f98
......@@ -56,6 +56,7 @@
#include "openmm/CMMotionRemover.h"
#include "openmm/System.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/internal/CustomHbondForceImpl.h"
#include "openmm/internal/NonbondedForceImpl.h"
#include "openmm/Integrator.h"
#include "openmm/OpenMMException.h"
......@@ -1238,6 +1239,8 @@ ReferenceCalcCustomHbondForceKernel::~ReferenceCalcCustomHbondForceKernel() {
disposeRealArray(donorParamArray, numDonors);
disposeRealArray(acceptorParamArray, numAcceptors);
disposeIntArray(exclusionArray, numDonors);
if (ixn != NULL)
delete ixn;
}
void ReferenceCalcCustomHbondForceKernel::initialize(const System& system, const CustomHbondForce& force) {
......@@ -1256,21 +1259,29 @@ void ReferenceCalcCustomHbondForceKernel::initialize(const System& system, const
// Build the arrays.
donorParticles.resize(numDonors);
vector<vector<int> > donorParticles(numDonors);
int numDonorParameters = force.getNumPerDonorParameters();
donorParamArray = allocateRealArray(numDonors, numDonorParameters);
for (int i = 0; i < numDonors; ++i) {
vector<double> parameters;
force.getDonorParameters(i, donorParticles[i].first, donorParticles[i].second, parameters);
int d1, d2, d3;
force.getDonorParameters(i, d1, d2, d3, parameters);
donorParticles[i].push_back(d1);
donorParticles[i].push_back(d2);
donorParticles[i].push_back(d3);
for (int j = 0; j < numDonorParameters; j++)
donorParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
}
acceptorParticles.resize(numAcceptors);
vector<vector<int> > acceptorParticles(numAcceptors);
int numAcceptorParameters = force.getNumPerAcceptorParameters();
acceptorParamArray = allocateRealArray(numAcceptors, numAcceptorParameters);
for (int i = 0; i < numAcceptors; ++i) {
vector<double> parameters;
force.getAcceptorParameters(i, acceptorParticles[i].first, acceptorParticles[i].second, parameters);
int a1, a2, a3;
force.getAcceptorParameters(i, a1, a2, a3, parameters);
acceptorParticles[i].push_back(a1);
acceptorParticles[i].push_back(a2);
acceptorParticles[i].push_back(a3);
for (int j = 0; j < numAcceptorParameters; j++)
acceptorParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
}
......@@ -1282,7 +1293,7 @@ void ReferenceCalcCustomHbondForceKernel::initialize(const System& system, const
for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter)
exclusionArray[i][++index] = *iter;
}
nonbondedMethod = CalcCustomHbondForceKernel::NonbondedMethod(force.getNonbondedMethod());
NonbondedMethod nonbondedMethod = CalcCustomHbondForceKernel::NonbondedMethod(force.getNonbondedMethod());
nonbondedCutoff = (RealOpenMM) force.getCutoffDistance();
Vec3 boxVectors[3];
system.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
......@@ -1302,20 +1313,26 @@ void ReferenceCalcCustomHbondForceKernel::initialize(const System& system, const
functions[name] = new ReferenceTabulatedFunction(min, max, values, interpolating);
}
// Parse the various expressions used to calculate the force.
// Parse the expression and create the object used to calculate the interaction.
Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction(), functions).optimize();
energyExpression = expression.createProgram();
rForceExpression = expression.differentiate("r").optimize().createProgram();
thetaForceExpression = expression.differentiate("theta").optimize().createProgram();
psiForceExpression = expression.differentiate("psi").optimize().createProgram();
chiForceExpression = expression.differentiate("chi").optimize().createProgram();
map<string, vector<int> > distances;
map<string, vector<int> > angles;
map<string, vector<int> > dihedrals;
Lepton::ParsedExpression energyExpression = CustomHbondForceImpl::prepareExpression(force, distances, angles, dihedrals);
vector<string> donorParameterNames;
vector<string> acceptorParameterNames;
for (int i = 0; i < numDonorParameters; i++)
donorParameterNames.push_back(force.getPerDonorParameterName(i));
for (int i = 0; i < numAcceptorParameters; i++)
acceptorParameterNames.push_back(force.getPerAcceptorParameterName(i));
for (int i = 0; i < force.getNumGlobalParameters(); i++)
globalParameterNames.push_back(force.getGlobalParameterName(i));
ixn = new ReferenceCustomHbondIxn(donorParticles, acceptorParticles, energyExpression, donorParameterNames, acceptorParameterNames, distances, angles, dihedrals);
bool periodic = (nonbondedMethod == CutoffPeriodic);
if (nonbondedMethod != NoCutoff)
ixn->setUseCutoff(nonbondedCutoff);
if (periodic)
ixn->setPeriodic(periodicBoxSize);
// Delete the custom functions.
......@@ -1326,34 +1343,20 @@ void ReferenceCalcCustomHbondForceKernel::initialize(const System& system, const
void ReferenceCalcCustomHbondForceKernel::executeForces(ContextImpl& context) {
RealOpenMM** posData = extractPositions(context);
RealOpenMM** forceData = extractForces(context);
ReferenceCustomHbondIxn ixn(donorParticles, acceptorParticles, energyExpression, rForceExpression, thetaForceExpression, psiForceExpression,
chiForceExpression, donorParameterNames, acceptorParameterNames);
bool periodic = (nonbondedMethod == CutoffPeriodic);
if (nonbondedMethod != NoCutoff)
ixn.setUseCutoff(nonbondedCutoff);
if (periodic)
ixn.setPeriodic(periodicBoxSize);
map<string, double> globalParameters;
for (int i = 0; i < (int) globalParameterNames.size(); i++)
globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
ixn.calculatePairIxn(posData, donorParamArray, acceptorParamArray, exclusionArray, globalParameters, forceData, 0);
ixn->calculatePairIxn(posData, donorParamArray, acceptorParamArray, exclusionArray, globalParameters, forceData, 0);
}
double ReferenceCalcCustomHbondForceKernel::executeEnergy(ContextImpl& context) {
RealOpenMM** posData = extractPositions(context);
RealOpenMM** forceData = allocateRealArray(numParticles, 3);
RealOpenMM energy = 0;
ReferenceCustomHbondIxn ixn(donorParticles, acceptorParticles, energyExpression, rForceExpression, thetaForceExpression, psiForceExpression,
chiForceExpression, donorParameterNames, acceptorParameterNames);
bool periodic = (nonbondedMethod == CutoffPeriodic);
if (nonbondedMethod != NoCutoff)
ixn.setUseCutoff(nonbondedCutoff);
if (periodic)
ixn.setPeriodic(periodicBoxSize);
map<string, double> globalParameters;
for (int i = 0; i < (int) globalParameterNames.size(); i++)
globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
ixn.calculatePairIxn(posData, donorParamArray, acceptorParamArray, exclusionArray, globalParameters, forceData, &energy);
ixn->calculatePairIxn(posData, donorParamArray, acceptorParamArray, exclusionArray, globalParameters, forceData, &energy);
disposeRealArray(forceData, numParticles);
return energy;
}
......
platforms/reference/src/ReferenceKernels.h
View file @ be137f98
......@@ -41,6 +41,7 @@
class CpuObc;
class CpuGBVI;
class ReferenceAndersenThermostat;
class ReferenceCustomHbondIxn;
class ReferenceBrownianDynamics;
class ReferenceStochasticDynamics;
class ReferenceConstraintAlgorithm;
......@@ -633,7 +634,7 @@ private:
*/
class ReferenceCalcCustomHbondForceKernel : public CalcCustomHbondForceKernel {
public:
ReferenceCalcCustomHbondForceKernel(std::string name, const Platform& platform) : CalcCustomHbondForceKernel(name, platform) {
ReferenceCalcCustomHbondForceKernel(std::string name, const Platform& platform) : CalcCustomHbondForceKernel(name, platform), ixn(NULL) {
}
~ReferenceCalcCustomHbondForceKernel();
/**
......@@ -661,11 +662,9 @@ private:
int **exclusionArray;
RealOpenMM **donorParamArray, **acceptorParamArray;
RealOpenMM nonbondedCutoff, periodicBoxSize[3];
ReferenceCustomHbondIxn* ixn;
std::vector<std::set<int> > exclusions;
std::vector<std::pair<int, int> > donorParticles, acceptorParticles;
Lepton::ExpressionProgram energyExpression, rForceExpression, thetaForceExpression, psiForceExpression, chiForceExpression;
std::vector<std::string> donorParameterNames, acceptorParameterNames, globalParameterNames;
NonbondedMethod nonbondedMethod;
std::vector<std::string> globalParameterNames;
};
/**
......
platforms/reference/src/SimTKReference/ReferenceCustomHbondIxn.cpp
View file @ be137f98
......@@ -44,14 +44,17 @@ using std::vector;
--------------------------------------------------------------------------------------- */
ReferenceCustomHbondIxn::ReferenceCustomHbondIxn(const vector<pair<int, int> >& donorAtoms, const vector<pair<int, int> >& acceptorAtoms,
const Lepton::ExpressionProgram& energyExpression,
const Lepton::ExpressionProgram& rForceExpression, const Lepton::ExpressionProgram& thetaForceExpression,
const Lepton::ExpressionProgram& psiForceExpression, const Lepton::ExpressionProgram& chiForceExpression,
const vector<string>& donorParameterNames, const vector<string>& acceptorParameterNames) :
cutoff(false), periodic(false), donorAtoms(donorAtoms), acceptorAtoms(acceptorAtoms), energyExpression(energyExpression),
rForceExpression(rForceExpression), thetaForceExpression(thetaForceExpression), psiForceExpression(psiForceExpression),
chiForceExpression(chiForceExpression), donorParamNames(donorParameterNames), acceptorParamNames(acceptorParameterNames) {
ReferenceCustomHbondIxn::ReferenceCustomHbondIxn(const vector<vector<int> >& donorAtoms, const vector<vector<int> >& acceptorAtoms,
const Lepton::ParsedExpression& energyExpression, const vector<string>& donorParameterNames, const vector<string>& acceptorParameterNames,
const map<string, vector<int> >& distances, const map<string, vector<int> >& angles, const map<string, vector<int> >& dihedrals) :
cutoff(false), periodic(false), donorAtoms(donorAtoms), acceptorAtoms(acceptorAtoms), energyExpression(energyExpression.createProgram()),
donorParamNames(donorParameterNames), acceptorParamNames(acceptorParameterNames) {
for (map<string, vector<int> >::const_iterator iter = distances.begin(); iter != distances.end(); ++iter)
distanceTerms.push_back(ReferenceCustomHbondIxn::DistanceTermInfo(iter->first, iter->second, energyExpression.differentiate(iter->first).optimize().createProgram()));
for (map<string, vector<int> >::const_iterator iter = angles.begin(); iter != angles.end(); ++iter)
angleTerms.push_back(ReferenceCustomHbondIxn::AngleTermInfo(iter->first, iter->second, energyExpression.differentiate(iter->first).optimize().createProgram()));
for (map<string, vector<int> >::const_iterator iter = dihedrals.begin(); iter != dihedrals.end(); ++iter)
dihedralTerms.push_back(ReferenceCustomHbondIxn::DihedralTermInfo(iter->first, iter->second, energyExpression.differentiate(iter->first).optimize().createProgram()));
}
/**---------------------------------------------------------------------------------------
......@@ -180,136 +183,113 @@ void ReferenceCustomHbondIxn::calculateOneIxn(int donor, int acceptor, RealOpenM
// ---------------------------------------------------------------------------------------
// Visual Studio complains if crossProduct declared as 'crossProduct[2][3]'
RealOpenMM crossProductMemory[6];
RealOpenMM* crossProduct[2];
crossProduct[0] = crossProductMemory;
crossProduct[1] = crossProductMemory + 3;
int atoms[6];
atoms[0] = acceptorAtoms[acceptor][0];
atoms[1] = acceptorAtoms[acceptor][1];
atoms[2] = acceptorAtoms[acceptor][2];
atoms[3] = donorAtoms[donor][0];
atoms[4] = donorAtoms[donor][1];
atoms[5] = donorAtoms[donor][2];
// Compute the distance between the primary donor and acceptor atoms, and compare to the cutoff.
int d1 = donorAtoms[donor].first;
int a1 = acceptorAtoms[acceptor].first;
RealOpenMM deltaD1A1[ReferenceForce::LastDeltaRIndex];
if (periodic)
ReferenceForce::getDeltaRPeriodic(atomCoordinates[a1], atomCoordinates[d1], periodicBoxSize, deltaD1A1);
else
ReferenceForce::getDeltaR(atomCoordinates[a1], atomCoordinates[d1], deltaD1A1);
if (cutoff && deltaD1A1[ReferenceForce::RIndex] >= cutoffDistance)
if (cutoff) {
RealOpenMM delta[ReferenceForce::LastDeltaRIndex];
computeDelta(atoms[0], atoms[3], delta, atomCoordinates);
if (delta[ReferenceForce::RIndex] >= cutoffDistance)
return;
}
// Compute all of the variables the energy can depend on.
int d2 = donorAtoms[donor].second;
int a2 = acceptorAtoms[acceptor].second;
RealOpenMM deltaD1D2[ReferenceForce::LastDeltaRIndex];
RealOpenMM deltaA2A1[ReferenceForce::LastDeltaRIndex];
if (periodic) {
ReferenceForce::getDeltaRPeriodic(atomCoordinates[d2], atomCoordinates[d1], periodicBoxSize, deltaD1D2);
ReferenceForce::getDeltaRPeriodic(atomCoordinates[a1], atomCoordinates[a2], periodicBoxSize, deltaA2A1);
for (int i = 0; i < (int) distanceTerms.size(); i++) {
const DistanceTermInfo& term = distanceTerms[i];
computeDelta(atoms[term.p1], atoms[term.p2], term.delta, atomCoordinates);
variables[term.name] = term.delta[ReferenceForce::RIndex];
}
else {
ReferenceForce::getDeltaR(atomCoordinates[d2], atomCoordinates[d1], deltaD1D2);
ReferenceForce::getDeltaR(atomCoordinates[a1], atomCoordinates[a2], deltaA2A1);
for (int i = 0; i < (int) angleTerms.size(); i++) {
const AngleTermInfo& term = angleTerms[i];
computeDelta(atoms[term.p1], atoms[term.p2], term.delta1, atomCoordinates);
computeDelta(atoms[term.p3], atoms[term.p2], term.delta2, atomCoordinates);
variables[term.name] = computeAngle(term.delta1, term.delta2);
}
variables["r"] = deltaD1A1[ReferenceForce::RIndex];
variables["theta"] = computeAngle(deltaD1A1, deltaD1D2, 1);
variables["psi"] = computeAngle(deltaD1A1, deltaA2A1, 1);
for (int i = 0; i < (int) dihedralTerms.size(); i++) {
const DihedralTermInfo& term = dihedralTerms[i];
computeDelta(atoms[term.p2], atoms[term.p1], term.delta1, atomCoordinates);
computeDelta(atoms[term.p2], atoms[term.p3], term.delta2, atomCoordinates);
computeDelta(atoms[term.p4], atoms[term.p3], term.delta3, atomCoordinates);
RealOpenMM dotDihedral, signOfDihedral;
variables["chi"] = getDihedralAngleBetweenThreeVectors(deltaA2A1, deltaD1A1, deltaD1D2,
crossProduct, &dotDihedral, deltaA2A1, &signOfDihedral, 1);
RealOpenMM* crossProduct[] = {term.cross1, term.cross2};
variables[term.name] = getDihedralAngleBetweenThreeVectors(term.delta1, term.delta2, term.delta3, crossProduct, &dotDihedral, term.delta1, &signOfDihedral, 1);
}
// Apply forces based on r.
// Apply forces based on distances.
RealOpenMM dEdR = (RealOpenMM) (rForceExpression.evaluate(variables)/(deltaD1A1[ReferenceForce::RIndex]));
if (dEdR != 0) {
for (int i = 0; i < (int) distanceTerms.size(); i++) {
const DistanceTermInfo& term = distanceTerms[i];
RealOpenMM dEdR = (RealOpenMM) (term.forceExpression.evaluate(variables)/(term.delta[ReferenceForce::RIndex]));
for (int i = 0; i < 3; i++) {
RealOpenMM force = -dEdR*deltaD1A1[i];
forces[d1][i] += force;
forces[a1][i] -= force;
RealOpenMM force = -dEdR*term.delta[i];
forces[atoms[term.p1]][i] -= force;
forces[atoms[term.p2]][i] += force;
}
}
// Apply forces based on theta.
// Apply forces based on angles.
RealOpenMM dEdTheta = (RealOpenMM) thetaForceExpression.evaluate(variables);
if (dEdTheta != 0) {
for (int i = 0; i < (int) angleTerms.size(); i++) {
const AngleTermInfo& term = angleTerms[i];
RealOpenMM dEdTheta = (RealOpenMM) term.forceExpression.evaluate(variables);
RealOpenMM thetaCross[ReferenceForce::LastDeltaRIndex];
SimTKOpenMMUtilities::crossProductVector3(deltaD1D2, deltaD1A1, thetaCross);
SimTKOpenMMUtilities::crossProductVector3(term.delta1, term.delta2, thetaCross);
RealOpenMM lengthThetaCross = SQRT(DOT3(thetaCross, thetaCross));
if (lengthThetaCross < 1.0e-06)
lengthThetaCross = (RealOpenMM) 1.0e-06;
RealOpenMM termA = dEdTheta/(deltaD1D2[ReferenceForce::R2Index]*lengthThetaCross);
RealOpenMM termC = -dEdTheta/(deltaD1A1[ReferenceForce::R2Index]*lengthThetaCross);
RealOpenMM termA = dEdTheta/(term.delta1[ReferenceForce::R2Index]*lengthThetaCross);
RealOpenMM termC = -dEdTheta/(term.delta2[ReferenceForce::R2Index]*lengthThetaCross);
RealOpenMM deltaCrossP[3][3];
SimTKOpenMMUtilities::crossProductVector3(deltaD1D2, thetaCross, deltaCrossP[0]);
SimTKOpenMMUtilities::crossProductVector3(deltaD1A1, thetaCross, deltaCrossP[2]);
SimTKOpenMMUtilities::crossProductVector3(term.delta1, thetaCross, deltaCrossP[0]);
SimTKOpenMMUtilities::crossProductVector3(term.delta2, thetaCross, deltaCrossP[2]);
for (int i = 0; i < 3; i++) {
deltaCrossP[0][i] *= termA;
deltaCrossP[2][i] *= termC;
deltaCrossP[1][i] = -(deltaCrossP[0][i]+deltaCrossP[2][i]);
}
for (int i = 0; i < 3; i++) {
forces[d2][i] += deltaCrossP[0][i];
forces[d1][i] += deltaCrossP[1][i];
forces[a1][i] += deltaCrossP[2][i];
forces[atoms[term.p1]][i] += deltaCrossP[0][i];
forces[atoms[term.p2]][i] += deltaCrossP[1][i];
forces[atoms[term.p3]][i] += deltaCrossP[2][i];
}
}
// Apply forces based on psi.
RealOpenMM dEdPsi = (RealOpenMM) psiForceExpression.evaluate(variables);
if (dEdPsi != 0) {
RealOpenMM psiCross[ReferenceForce::LastDeltaRIndex];
SimTKOpenMMUtilities::crossProductVector3(deltaA2A1, deltaD1A1, psiCross);
RealOpenMM lengthPsiCross = SQRT(DOT3(psiCross, psiCross));
if (lengthPsiCross < 1.0e-06)
lengthPsiCross = (RealOpenMM) 1.0e-06;
RealOpenMM termA = dEdPsi/(deltaD1A1[ReferenceForce::R2Index]*lengthPsiCross);
RealOpenMM termC = -dEdPsi/(deltaA2A1[ReferenceForce::R2Index]*lengthPsiCross);
RealOpenMM deltaCrossP[3][3];
SimTKOpenMMUtilities::crossProductVector3(deltaD1A1, psiCross, deltaCrossP[0]);
SimTKOpenMMUtilities::crossProductVector3(deltaA2A1, psiCross, deltaCrossP[2]);
for (int i = 0; i < 3; i++) {
deltaCrossP[0][i] *= termA;
deltaCrossP[2][i] *= termC;
deltaCrossP[1][i] = -(deltaCrossP[0][i]+deltaCrossP[2][i]);
}
for (int i = 0; i < 3; i++) {
forces[d1][i] += deltaCrossP[0][i];
forces[a1][i] += deltaCrossP[1][i];
forces[a2][i] += deltaCrossP[2][i];
}
}
// Apply forces based on chi.
// Apply forces based on dihedrals.
RealOpenMM dEdChi = (RealOpenMM) chiForceExpression.evaluate(variables);
if (dEdChi != 0) {
for (int i = 0; i < (int) dihedralTerms.size(); i++) {
const DihedralTermInfo& term = dihedralTerms[i];
RealOpenMM dEdTheta = (RealOpenMM) term.forceExpression.evaluate(variables);
RealOpenMM internalF[4][3];
RealOpenMM forceFactors[4];
RealOpenMM normCross1 = DOT3(crossProduct[0], crossProduct[0]);
RealOpenMM normBC = deltaD1A1[ReferenceForce::RIndex];
forceFactors[0] = (-dEdChi*normBC)/normCross1;
RealOpenMM normCross2 = DOT3(crossProduct[1], crossProduct[1]);
forceFactors[3] = (dEdChi*normBC)/normCross2;
forceFactors[1] = DOT3(deltaA2A1, deltaD1A1);
forceFactors[1] /= deltaD1A1[ReferenceForce::R2Index];
forceFactors[2] = DOT3(deltaD1D2, deltaD1A1);
forceFactors[2] /= deltaD1A1[ReferenceForce::R2Index];
RealOpenMM normCross1 = DOT3(term.cross1, term.cross1);
RealOpenMM normBC = term.delta2[ReferenceForce::RIndex];
forceFactors[0] = (-dEdTheta*normBC)/normCross1;
RealOpenMM normCross2 = DOT3(term.cross2, term.cross2);
forceFactors[3] = (dEdTheta*normBC)/normCross2;
forceFactors[1] = DOT3(term.delta1, term.delta2);
forceFactors[1] /= term.delta2[ReferenceForce::R2Index];
forceFactors[2] = DOT3(term.delta3, term.delta2);
forceFactors[2] /= term.delta2[ReferenceForce::R2Index];
for (int i = 0; i < 3; i++) {
internalF[0][i] = forceFactors[0]*crossProduct[0][i];
internalF[3][i] = forceFactors[3]*crossProduct[1][i];
internalF[0][i] = forceFactors[0]*term.cross1[i];
internalF[3][i] = forceFactors[3]*term.cross2[i];
RealOpenMM s = forceFactors[1]*internalF[0][i] - forceFactors[2]*internalF[3][i];
internalF[1][i] = internalF[0][i] - s;
internalF[2][i] = internalF[3][i] + s;
}
for (int i = 0; i < 3; i++) {
forces[a2][i] += internalF[0][i];
forces[a1][i] -= internalF[1][i];
forces[d1][i] -= internalF[2][i];
forces[d2][i] += internalF[3][i];
forces[atoms[term.p1]][i] += internalF[0][i];
forces[atoms[term.p2]][i] -= internalF[1][i];
forces[atoms[term.p3]][i] -= internalF[2][i];
forces[atoms[term.p4]][i] += internalF[3][i];
}
}
......@@ -319,8 +299,15 @@ void ReferenceCustomHbondIxn::calculateOneIxn(int donor, int acceptor, RealOpenM
*totalEnergy += (RealOpenMM) energyExpression.evaluate(variables);
}
RealOpenMM ReferenceCustomHbondIxn::computeAngle(RealOpenMM* vec1, RealOpenMM* vec2, RealOpenMM sign) {
RealOpenMM dot = sign*DOT3(vec1, vec2);
void ReferenceCustomHbondIxn::computeDelta(int atom1, int atom2, RealOpenMM* delta, RealOpenMM** atomCoordinates) const {
if (periodic)
ReferenceForce::getDeltaRPeriodic(atomCoordinates[atom1], atomCoordinates[atom2], periodicBoxSize, delta);
else
ReferenceForce::getDeltaR(atomCoordinates[atom1], atomCoordinates[atom2], delta);
}
RealOpenMM ReferenceCustomHbondIxn::computeAngle(RealOpenMM* vec1, RealOpenMM* vec2) {
RealOpenMM dot = DOT3(vec1, vec2);
RealOpenMM cosine = dot/SQRT((vec1[ReferenceForce::R2Index]*vec2[ReferenceForce::R2Index]));
RealOpenMM angle;
if (cosine >= 1)
......
platforms/reference/src/SimTKReference/ReferenceCustomHbondIxn.h
View file @ be137f98
......@@ -27,6 +27,7 @@
#include "ReferenceBondIxn.h"
#include "lepton/ExpressionProgram.h"
#include "lepton/ParsedExpression.h"
#include <map>
#include <vector>
......@@ -36,17 +37,19 @@ class ReferenceCustomHbondIxn : public ReferenceBondIxn {
private:
class DistanceTermInfo;
class AngleTermInfo;
class DihedralTermInfo;
bool cutoff;
bool periodic;
RealOpenMM periodicBoxSize[3];
RealOpenMM cutoffDistance;
std::vector<std::vector<int> > donorAtoms, acceptorAtoms;
Lepton::ExpressionProgram energyExpression;
Lepton::ExpressionProgram rForceExpression;
Lepton::ExpressionProgram thetaForceExpression;
Lepton::ExpressionProgram psiForceExpression;
Lepton::ExpressionProgram chiForceExpression;
std::vector<std::string> donorParamNames, acceptorParamNames;
std::vector<std::pair<int, int> > donorAtoms, acceptorAtoms;
std::vector<DistanceTermInfo> distanceTerms;
std::vector<AngleTermInfo> angleTerms;
std::vector<DihedralTermInfo> dihedralTerms;
/**---------------------------------------------------------------------------------------
......@@ -65,7 +68,9 @@ class ReferenceCustomHbondIxn : public ReferenceBondIxn {
std::map<std::string, double>& variables, RealOpenMM** forces,
RealOpenMM* totalEnergy) const;
static RealOpenMM computeAngle(RealOpenMM* vec1, RealOpenMM* vec2, RealOpenMM sign);
void computeDelta(int atom1, int atom2, RealOpenMM* delta, RealOpenMM** atomCoordinates) const;
static RealOpenMM computeAngle(RealOpenMM* vec1, RealOpenMM* vec2);
public:
......@@ -76,11 +81,10 @@ class ReferenceCustomHbondIxn : public ReferenceBondIxn {
--------------------------------------------------------------------------------------- */
ReferenceCustomHbondIxn(const std::vector<std::pair<int, int> >& donorAtoms, const std::vector<std::pair<int, int> >& acceptorAtoms,
const Lepton::ExpressionProgram& energyExpression, const Lepton::ExpressionProgram& rForceExpression,
const Lepton::ExpressionProgram& thetaForceExpression, const Lepton::ExpressionProgram& psiForceExpression,
const Lepton::ExpressionProgram& chiForceExpression, const std::vector<std::string>& donorParameterNames,
const std::vector<std::string>& acceptorParameterNames);
ReferenceCustomHbondIxn(const std::vector<std::vector<int> >& donorAtoms, const std::vector<std::vector<int> >& acceptorAtoms,
const Lepton::ParsedExpression& energyExpression, const std::vector<std::string>& donorParameterNames,
const std::vector<std::string>& acceptorParameterNames, const std::map<std::string, std::vector<int> >& distances,
const std::map<std::string, std::vector<int> >& angles, const std::map<std::string, std::vector<int> >& dihedrals);
/**---------------------------------------------------------------------------------------
......@@ -137,4 +141,42 @@ class ReferenceCustomHbondIxn : public ReferenceBondIxn {
};
class ReferenceCustomHbondIxn::DistanceTermInfo {
public:
std::string name;
int p1, p2;
Lepton::ExpressionProgram forceExpression;
mutable RealOpenMM delta[ReferenceForce::LastDeltaRIndex];
DistanceTermInfo(const std::string& name, const std::vector<int>& atoms, const Lepton::ExpressionProgram& forceExpression) :
name(name), p1(atoms[0]), p2(atoms[1]), forceExpression(forceExpression) {
}
};
class ReferenceCustomHbondIxn::AngleTermInfo {
public:
std::string name;
int p1, p2, p3;
Lepton::ExpressionProgram forceExpression;
mutable RealOpenMM delta1[ReferenceForce::LastDeltaRIndex];
mutable RealOpenMM delta2[ReferenceForce::LastDeltaRIndex];
AngleTermInfo(const std::string& name, const std::vector<int>& atoms, const Lepton::ExpressionProgram& forceExpression) :
name(name), p1(atoms[0]), p2(atoms[1]), p3(atoms[2]), forceExpression(forceExpression) {
}
};
class ReferenceCustomHbondIxn::DihedralTermInfo {
public:
std::string name;
int p1, p2, p3, p4;
Lepton::ExpressionProgram forceExpression;
mutable RealOpenMM delta1[ReferenceForce::LastDeltaRIndex];
mutable RealOpenMM delta2[ReferenceForce::LastDeltaRIndex];
mutable RealOpenMM delta3[ReferenceForce::LastDeltaRIndex];
mutable RealOpenMM cross1[3];
mutable RealOpenMM cross2[3];
DihedralTermInfo(const std::string& name, const std::vector<int>& atoms, const Lepton::ExpressionProgram& forceExpression) :
name(name), p1(atoms[0]), p2(atoms[1]), p3(atoms[2]), p4(atoms[3]), forceExpression(forceExpression) {
}
};
#endif // __ReferenceCustomHbondIxn_H__
platforms/reference/tests/TestReferenceCustomHbondForce.cpp
View file @ be137f98
......@@ -30,7 +30,7 @@
* -------------------------------------------------------------------------- */
/**
* This tests the reference implementation of CustomTorsionForce.
* This tests the reference implementation of CustomHbondForce.
*/
#include "../../../tests/AssertionUtilities.h"
......@@ -61,7 +61,8 @@ void testHbond() {
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
CustomHbondForce* custom = new CustomHbondForce("0.5*kr*(r-r0)^2 + 0.5*ktheta*(theta-theta0)^2 + 0.5*kpsi*(psi-psi0)^2 + kchi*(1+cos(n*chi-chi0))");
customSystem.addParticle(1.0);
CustomHbondForce* custom = new CustomHbondForce("0.5*kr*(distance(d1,a1)-r0)^2 + 0.5*ktheta*(angle(a1,d1,d2)-theta0)^2 + 0.5*kpsi*(angle(d1,a1,a2)-psi0)^2 + kchi*(1+cos(n*dihedral(a3,a2,a1,d1)-chi0))");
custom->addPerDonorParameter("r0");
custom->addPerDonorParameter("theta0");
custom->addPerDonorParameter("psi0");
......@@ -75,11 +76,11 @@ void testHbond() {
parameters[0] = 1.5;
parameters[1] = 1.7;
parameters[2] = 1.9;
custom->addDonor(1, 0, parameters);
custom->addDonor(1, 0, -1, parameters);
parameters.resize(2);
parameters[0] = 2.1;
parameters[1] = 2;
custom->addAcceptor(2, 3, parameters);
custom->addAcceptor(2, 3, 4, parameters);
custom->setCutoffDistance(10.0);
customSystem.addForce(custom);
......@@ -90,6 +91,7 @@ void testHbond() {
standardSystem.addParticle(1.0);
standardSystem.addParticle(1.0);
standardSystem.addParticle(1.0);
standardSystem.addParticle(1.0);
HarmonicBondForce* bond = new HarmonicBondForce();
bond->addBond(1, 2, 1.5, 0.4);
standardSystem.addForce(bond);
......@@ -98,13 +100,13 @@ void testHbond() {
angle->addAngle(1, 2, 3, 1.9, 0.6);
standardSystem.addForce(angle);
PeriodicTorsionForce* torsion = new PeriodicTorsionForce();
torsion->addTorsion(0, 1, 2, 3, 2, 2.1, 0.7);;
torsion->addTorsion(1, 2, 3, 4, 2, 2.1, 0.7);;
standardSystem.addForce(torsion);
// Set the atoms in various positions, and verify that both systems give identical forces and energy.
init_gen_rand(0);
vector<Vec3> positions(4);
vector<Vec3> positions(5);
VerletIntegrator integrator1(0.01);
VerletIntegrator integrator2(0.01);
for (int i = 0; i < 10; i++) {
......@@ -129,10 +131,10 @@ void testExclusions() {
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomHbondForce* custom = new CustomHbondForce("(r-1)^2");
custom->addDonor(0, 1, vector<double>());
custom->addDonor(1, 0, vector<double>());
custom->addAcceptor(2, 0, vector<double>());
CustomHbondForce* custom = new CustomHbondForce("(distance(d1,a1)-1)^2");
custom->addDonor(0, 1, -1, vector<double>());
custom->addDonor(1, 0, -1, vector<double>());
custom->addAcceptor(2, 0, -1, vector<double>());
custom->addExclusion(1, 0);
system.addForce(custom);
Context context(system, integrator, platform);
......@@ -156,10 +158,10 @@ void testCutoff() {
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomHbondForce* custom = new CustomHbondForce("(r-1)^2");
custom->addDonor(0, 1, vector<double>());
custom->addDonor(1, 0, vector<double>());
custom->addAcceptor(2, 0, vector<double>());
CustomHbondForce* custom = new CustomHbondForce("(distance(d1,a1)-1)^2");
custom->addDonor(0, 1, -1, vector<double>());
custom->addDonor(1, 0, -1, vector<double>());
custom->addAcceptor(2, 0, -1, vector<double>());
custom->setNonbondedMethod(CustomHbondForce::CutoffNonPeriodic);
custom->setCutoffDistance(2.5);
system.addForce(custom);
......
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