Skip to content

GitLab

Commit 35aef079 authored by Peter Eastman's avatar Peter Eastman
Browse files

Created CustomCompoundBondForce, including reference implementation

parent 36fced55
Hide whitespace changes
Inline Side-by-side
Showing with 1491 additions and 6 deletions
olla/include/openmm/kernels.h
View file @ 35aef079
...@@ -38,6 +38,7 @@ ...@@ -38,6 +38,7 @@
#include "openmm/CMMotionRemover.h" #include "openmm/CMMotionRemover.h"
#include "openmm/CustomAngleForce.h" #include "openmm/CustomAngleForce.h"
#include "openmm/CustomBondForce.h" #include "openmm/CustomBondForce.h"
#include "openmm/CustomCompoundBondForce.h"
#include "openmm/CustomExternalForce.h" #include "openmm/CustomExternalForce.h"
#include "openmm/CustomGBForce.h" #include "openmm/CustomGBForce.h"
#include "openmm/CustomHbondForce.h" #include "openmm/CustomHbondForce.h"
...@@ -679,6 +680,34 @@ public: ...@@ -679,6 +680,34 @@ public:
virtual double execute(ContextImpl& context, bool includeForces, bool includeEnergy) = 0; virtual double execute(ContextImpl& context, bool includeForces, bool includeEnergy) = 0;
}; };
/**
* This kernel is invoked by CustomCompoundBondForce to calculate the forces acting on the system and the energy of the system.
*/
class CalcCustomCompoundBondForceKernel : public KernelImpl {
public:
static std::string Name() {
return "CalcCustomCompoundBondForce";
}
CalcCustomCompoundBondForceKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
}
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param force the CustomCompoundBondForce this kernel will be used for
*/
virtual void initialize(const System& system, const CustomCompoundBondForce& force) = 0;
/**
* Execute the kernel to calculate the forces and/or energy.
*
* @param context the context in which to execute this kernel
* @param includeForces true if forces should be calculated
* @param includeEnergy true if the energy should be calculated
* @return the potential energy due to the force
*/
virtual double execute(ContextImpl& context, bool includeForces, bool includeEnergy) = 0;
};
/** /**
* This kernel is invoked by VerletIntegrator to take one time step. * This kernel is invoked by VerletIntegrator to take one time step.
*/ */
......
openmmapi/include/openmm/CustomCompoundBondForce.h 0 → 100644
View file @ 35aef079
#ifndef OPENMM_CUSTOMCOMPOUNDBONDFORCE_H_
#define OPENMM_CUSTOMCOMPOUNDBONDFORCE_H_
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2012 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "Force.h"
#include "Vec3.h"
#include <vector>
#include "internal/windowsExport.h"
namespace OpenMM {
/**
* This class supports a wide variety of bonded interactions. It defines a "bond" as a single energy term
* that depends on the positions of a fixed set of particles. The number of particles involved in a bond, and how
* the energy depends on their positions, is configurable. It may depend on the positions of individual particles,
* the distances between pairs of particles, the angles formed by sets of three particles, and the dihedral
* angles formed by sets of four particles.
*
* We refer to the particles in a bond as p1, p2, p3, etc. For each bond, CustomCompoundBondForce evaluates a
* user supplied algebraic expression to determine the interaction energy. The expression may depend on the
* following variables and functions:
*
* <ul>
* <li>x1, y1, z1, x2, y2, z2, etc.: The x, y, and z coordinates of the particle positions. For example, x1
* is the x coordinate of particle p1, and y3 is the y coordinate of particle p3.</li>
* <li>distance(p1, p2): the distance between particles p1 and p2 (where "p1" and "p2" may be replaced by the names
* of whichever particles you want to calculate the distance between).</li>
* <li>angle(p1, p2, p3): the angle formed by the three specified particles.</li>
* <li>dihedral(p1, p2, p3, p4): the dihedral angle formed by the four specified particles.</li>
* </ul>
*
* The expression also may involve tabulated functions, and may depend on arbitrary
* global and per-bond parameters.
*
* To use this class, create a CustomCompoundBondForce object, passing an algebraic expression to the constructor
* that defines the interaction energy of each bond. Then call addPerBondParameter() to define per-bond
* parameters and addGlobalParameter() to define global parameters. The values of per-bond parameters are specified
* as part of the system definition, while values of global parameters may be modified during a simulation by calling
* Context::setParameter().
*
* Next, call addBond() to define bonds and specify their parameter values. After a bond has been added, you can
* modify its parameters by calling setBondParameters().
*
* As an example, the following code creates a CustomCompoundBondForce that implements a Urey-Bradley potential. This
* is an interaction between three particles that depends on the angle formed by p1-p2-p3, and on the distance between
* p1 and p3.
*
* <tt>CustomCompoundBondForce* force = new CustomCompoundBondForce(3, "0.5*(kangle*(angle(p1,p2,p3)-theta0)^2+kbond*(distance(p1,p3)-r0)^2)");</tt>
*
* This force depends on four parameters: kangle, kbond, theta0, and r0. The following code defines these as per-bond parameters:
*
* <tt><pre>
* force->addPerBondParameter("kangle");
* force->addPerBondParameter("kbond");
* force->addPerBondParameter("theta0");
* force->addPerBondParameter("r0");
* </pre></tt>
*
* 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, erf, erfc, min, max, abs, step, delta. All trigonometric functions
* are defined in radians, and log is the natural logarithm. step(x) = 0 if x is less than 0, 1 otherwise. delta(x) = 1 if x is 0, 0 otherwise.
*
* In addition, you can call addFunction() to define a new function based on tabulated values. You specify a vector of
* values, and a natural spline is created from them. That function can then appear in the expression.
*/
class OPENMM_EXPORT CustomCompoundBondForce : public Force {
public:
/**
* Create a CustomCompoundBondForce.
*
* @param numParticles the number of particles used to define each bond
* @param energy an algebraic expression giving the interaction energy of each bond as a function
* of particle positions, inter-particle distances, angles, and dihedrals, and any global
* and per-bond parameters
*/
explicit CustomCompoundBondForce(int numParticles, const std::string& energy);
/**
* Get the number of particles used to define each bond.
*/
int getNumParticlesPerBond() const {
return particlesPerBond;
}
/**
* Get the number of bonds for which force field parameters have been defined.
*/
int getNumBonds() const {
return bonds.size();
}
/**
* Get the number of per-bond parameters that the interaction depends on.
*/
int getNumPerBondParameters() const {
return bondParameters.size();
}
/**
* Get the number of global parameters that the interaction depends on.
*/
int getNumGlobalParameters() const {
return globalParameters.size();
}
/**
* Get the number of tabulated functions that have been defined.
*/
int getNumFunctions() const {
return functions.size();
}
/**
* Get the algebraic expression that gives the interaction energy of each bond
*/
const std::string& getEnergyFunction() const;
/**
* Set the algebraic expression that gives the interaction energy of each bond
*/
void setEnergyFunction(const std::string& energy);
/**
* Add a new per-bond parameter that the interaction may depend on.
*
* @param name the name of the parameter
* @return the index of the parameter that was added
*/
int addPerBondParameter(const std::string& name);
/**
* Get the name of a per-bond parameter.
*
* @param index the index of the parameter for which to get the name
* @return the parameter name
*/
const std::string& getPerBondParameterName(int index) const;
/**
* Set the name of a per-bond parameter.
*
* @param index the index of the parameter for which to set the name
* @param name the name of the parameter
*/
void setPerBondParameterName(int index, const std::string& name);
/**
* Add a new global parameter that the interaction may depend on.
*
* @param name the name of the parameter
* @param defaultValue the default value of the parameter
* @return the index of the parameter that was added
*/
int addGlobalParameter(const std::string& name, double defaultValue);
/**
* Get the name of a global parameter.
*
* @param index the index of the parameter for which to get the name
* @return the parameter name
*/
const std::string& getGlobalParameterName(int index) const;
/**
* Set the name of a global parameter.
*
* @param index the index of the parameter for which to set the name
* @param name the name of the parameter
*/
void setGlobalParameterName(int index, const std::string& name);
/**
* Get the default value of a global parameter.
*
* @param index the index of the parameter for which to get the default value
* @return the parameter default value
*/
double getGlobalParameterDefaultValue(int index) const;
/**
* Set the default value of a global parameter.
*
* @param index the index of the parameter for which to set the default value
* @param name the default value of the parameter
*/
void setGlobalParameterDefaultValue(int index, double defaultValue);
/**
* Add a bond to the force
*
* @param particles the indices of the particles the bond depends on
* @param parameters the list of per-bond parameter values for the new bond
* @return the index of the bond that was added
*/
int addBond(const std::vector<int>& particles, const std::vector<double>& parameters);
/**
* Get the properties of a bond.
*
* @param index the index of the bond to get
* @param particles the indices of the particles in the bond
* @param parameters the list of per-bond parameter values for the bond
*/
void getBondParameters(int index, std::vector<int>& particles, std::vector<double>& parameters) const;
/**
* Set the properties of a bond.
*
* @param index the index of the bond group to set
* @param particles the indices of the particles in the bond
* @param parameters the list of per-bond parameter values for the bond
*/
void setBondParameters(int index, const std::vector<int>& particles, const std::vector<double>& parameters);
/**
* Add a tabulated function that may appear in the energy expression.
*
* @param name the name of the function as it appears in expressions
* @param values the tabulated values of the function f(x) at uniformly spaced values of x between min and max.
* The function is assumed to be zero for x &lt; min or x &gt; max.
* @param min the value of the independent variable corresponding to the first element of values
* @param max the value of the independent variable corresponding to the last element of values
* @return the index of the function that was added
*/
int addFunction(const std::string& name, const std::vector<double>& values, double min, double max);
/**
* Get the parameters for a tabulated function that may appear in the energy expression.
*
* @param index the index of the function for which to get parameters
* @param name the name of the function as it appears in expressions
* @param values the tabulated values of the function f(x) at uniformly spaced values of x between min and max.
* The function is assumed to be zero for x &lt; min or x &gt; max.
* @param min the value of the independent variable corresponding to the first element of values
* @param max the value of the independent variable corresponding to the last element of values
*/
void getFunctionParameters(int index, std::string& name, std::vector<double>& values, double& min, double& max) const;
/**
* Set the parameters for a tabulated function that may appear in algebraic expressions.
*
* @param index the index of the function for which to set parameters
* @param name the name of the function as it appears in expressions
* @param values the tabulated values of the function f(x) at uniformly spaced values of x between min and max.
* The function is assumed to be zero for x &lt; min or x &gt; max.
* @param min the value of the independent variable corresponding to the first element of values
* @param max the value of the independent variable corresponding to the last element of values
*/
void setFunctionParameters(int index, const std::string& name, const std::vector<double>& values, double min, double max);
protected:
ForceImpl* createImpl();
private:
class BondInfo;
class BondParameterInfo;
class GlobalParameterInfo;
class FunctionInfo;
int particlesPerBond;
std::string energyExpression;
std::vector<BondParameterInfo> bondParameters;
std::vector<GlobalParameterInfo> globalParameters;
std::vector<BondInfo> bonds;
std::vector<FunctionInfo> functions;
};
/**
* This is an internal class used to record information about a bond or acceptor.
* @private
*/
class CustomCompoundBondForce::BondInfo {
public:
std::vector<int> particles;
std::vector<double> parameters;
BondInfo() {
}
BondInfo(const std::vector<int>& particles, const std::vector<double>& parameters) :
particles(particles), parameters(parameters) {
}
};
/**
* This is an internal class used to record information about a per-bond or per-acceptor parameter.
* @private
*/
class CustomCompoundBondForce::BondParameterInfo {
public:
std::string name;
BondParameterInfo() {
}
BondParameterInfo(const std::string& name) : name(name) {
}
};
/**
* This is an internal class used to record information about a global parameter.
* @private
*/
class CustomCompoundBondForce::GlobalParameterInfo {
public:
std::string name;
double defaultValue;
GlobalParameterInfo() {
}
GlobalParameterInfo(const std::string& name, double defaultValue) : name(name), defaultValue(defaultValue) {
}
};
/**
* This is an internal class used to record information about a tabulated function.
* @private
*/
class CustomCompoundBondForce::FunctionInfo {
public:
std::string name;
std::vector<double> values;
double min, max;
FunctionInfo() {
}
FunctionInfo(const std::string& name, const std::vector<double>& values, double min, double max) :
name(name), values(values), min(min), max(max) {
}
};
} // namespace OpenMM
#endif /*OPENMM_CUSTOMCOMPOUNDBONDFORCE_H_*/
openmmapi/include/openmm/CustomHbondForce.h
View file @ 35aef079
...@@ -9,7 +9,7 @@ ...@@ -9,7 +9,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for * * Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. * * Medical Research, grant U54 GM072970. See https://simtk.org. *
* * * *
* Portions copyright (c) 2008-2010 Stanford University and the Authors. * * Portions copyright (c) 2008-2012 Stanford University and the Authors. *
* Authors: Peter Eastman * * Authors: Peter Eastman *
* Contributors: * * Contributors: *
* * * *
...@@ -119,7 +119,7 @@ public: ...@@ -119,7 +119,7 @@ public:
/** /**
* Create a CustomHbondForce. * Create a CustomHbondForce.
* *
* @param energy an algebraic expression giving the interaction energy between a donor as a function * @param energy an algebraic expression giving the interaction energy between a donor and an acceptor as a function
* of inter-particle distances, angles, and dihedrals, 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 * per-acceptor parameters
*/ */
...@@ -295,7 +295,7 @@ public: ...@@ -295,7 +295,7 @@ public:
* includes one particle, this will be -1. * includes one particle, this will be -1.
* @param d3 the index of the third particle for this donor group. If the group includes * @param d3 the index of the third particle for this donor group. If the group includes
* less than three particles, this will be -1. * less than three particles, this will be -1.
* @param parameters the list of per-donor parameter values for the new donor * @param parameters the list of per-donor parameter values for the donor
*/ */
void getDonorParameters(int index, int& d1, int& d2, int& d3, std::vector<double>& parameters) const; void getDonorParameters(int index, int& d1, int& d2, int& d3, std::vector<double>& parameters) const;
/** /**
...@@ -307,7 +307,7 @@ public: ...@@ -307,7 +307,7 @@ public:
* includes one particle, this must be -1. * includes one particle, this must be -1.
* @param d3 the index of the third particle for this donor group. If the group includes * @param d3 the index of the third particle for this donor group. If the group includes
* less than three particles, this must be -1. * less than three particles, this must be -1.
* @param parameters the list of per-donor parameter values for the new donor * @param parameters the list of per-donor parameter values for the donor
*/ */
void setDonorParameters(int index, int d1, int d2, int d3, const std::vector<double>& parameters); void setDonorParameters(int index, int d1, int d2, int d3, const std::vector<double>& parameters);
/** /**
...@@ -331,7 +331,7 @@ public: ...@@ -331,7 +331,7 @@ public:
* includes one particle, this will be -1. * includes one particle, this will be -1.
* @param a3 the index of the third particle for this acceptor group. If the group includes * @param a3 the index of the third particle for this acceptor group. If the group includes
* less than three particles, this will be -1. * less than three particles, this will be -1.
* @param parameters the list of per-acceptor parameter values for the new acceptor * @param parameters the list of per-acceptor parameter values for the acceptor
*/ */
void getAcceptorParameters(int index, int& a1, int& a2, int& a3, std::vector<double>& parameters) const; void getAcceptorParameters(int index, int& a1, int& a2, int& a3, std::vector<double>& parameters) const;
/** /**
...@@ -343,7 +343,7 @@ public: ...@@ -343,7 +343,7 @@ public:
* includes one particle, this must be -1. * includes one particle, this must be -1.
* @param a3 the index of the third particle for this acceptor group. If the group includes * @param a3 the index of the third particle for this acceptor group. If the group includes
* less than three particles, this must be -1. * less than three particles, this must be -1.
* @param parameters the list of per-acceptor parameter values for the new acceptor * @param parameters the list of per-acceptor parameter values for the acceptor
*/ */
void setAcceptorParameters(int index, int a1, int a2, int a3, const std::vector<double>& parameters); void setAcceptorParameters(int index, int a1, int a2, int a3, const std::vector<double>& parameters);
/** /**
......
openmmapi/include/openmm/internal/CustomCompoundBondForceImpl.h 0 → 100644
View file @ 35aef079
#ifndef OPENMM_CUSTOMCOMPOUNDBONDFORCEIMPL_H_
#define OPENMM_CUSTOMCOMPOUNDBONDFORCEIMPL_H_
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2012 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "ForceImpl.h"
#include "openmm/CustomCompoundBondForce.h"
#include "openmm/Kernel.h"
#include "lepton/CustomFunction.h"
#include "lepton/ExpressionTreeNode.h"
#include "lepton/ParsedExpression.h"
#include <utility>
#include <map>
#include <string>
namespace OpenMM {
/**
* This is the internal implementation of CustomNonbondedForce.
*/
class OPENMM_EXPORT CustomCompoundBondForceImpl : public ForceImpl {
public:
CustomCompoundBondForceImpl(CustomCompoundBondForce& owner);
~CustomCompoundBondForceImpl();
void initialize(ContextImpl& context);
CustomCompoundBondForce& getOwner() {
return owner;
}
void updateContextState(ContextImpl& context) {
// This force field doesn't update the state directly.
}
double calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups);
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.
*
* @param force the CustomCompoundBondForce to process
* @param functions definitions of custom function that may appear in the expression
* @param distances on exit, 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 exit, 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 exit, 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 CustomCompoundBondForce& force, const std::map<std::string, Lepton::CustomFunction*>& functions, 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);
CustomCompoundBondForce& owner;
Kernel kernel;
};
} // namespace OpenMM
#endif /*OPENMM_CUSTOMCOMPOUNDBONDFORCEIMPL_H_*/
openmmapi/src/CustomCompoundBondForce.cpp 0 → 100644
View file @ 35aef079
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2012 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/Force.h"
#include "openmm/OpenMMException.h"
#include "openmm/CustomCompoundBondForce.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/internal/CustomCompoundBondForceImpl.h"
#include <cmath>
#include <map>
#include <set>
#include <sstream>
#include <utility>
using namespace OpenMM;
using std::map;
using std::pair;
using std::set;
using std::string;
using std::stringstream;
using std::vector;
CustomCompoundBondForce::CustomCompoundBondForce(int numParticles, const string& energy) : particlesPerBond(numParticles), energyExpression(energy) {
}
const string& CustomCompoundBondForce::getEnergyFunction() const {
return energyExpression;
}
void CustomCompoundBondForce::setEnergyFunction(const std::string& energy) {
energyExpression = energy;
}
int CustomCompoundBondForce::addPerBondParameter(const string& name) {
bondParameters.push_back(BondParameterInfo(name));
return bondParameters.size()-1;
}
const string& CustomCompoundBondForce::getPerBondParameterName(int index) const {
ASSERT_VALID_INDEX(index, bondParameters);
return bondParameters[index].name;
}
void CustomCompoundBondForce::setPerBondParameterName(int index, const string& name) {
ASSERT_VALID_INDEX(index, bondParameters);
bondParameters[index].name = name;
}
int CustomCompoundBondForce::addGlobalParameter(const string& name, double defaultValue) {
globalParameters.push_back(GlobalParameterInfo(name, defaultValue));
return globalParameters.size()-1;
}
const string& CustomCompoundBondForce::getGlobalParameterName(int index) const {
ASSERT_VALID_INDEX(index, globalParameters);
return globalParameters[index].name;
}
void CustomCompoundBondForce::setGlobalParameterName(int index, const string& name) {
ASSERT_VALID_INDEX(index, globalParameters);
globalParameters[index].name = name;
}
double CustomCompoundBondForce::getGlobalParameterDefaultValue(int index) const {
ASSERT_VALID_INDEX(index, globalParameters);
return globalParameters[index].defaultValue;
}
void CustomCompoundBondForce::setGlobalParameterDefaultValue(int index, double defaultValue) {
ASSERT_VALID_INDEX(index, globalParameters);
globalParameters[index].defaultValue = defaultValue;
}
int CustomCompoundBondForce::addBond(const vector<int>& particles, const vector<double>& parameters) {
if (particles.size() != particlesPerBond)
throw OpenMMException("CustomCompoundBondForce: wrong number of particles specified for a bond.");
bonds.push_back(BondInfo(particles, parameters));
return bonds.size()-1;
}
void CustomCompoundBondForce::getBondParameters(int index, vector<int>& particles, std::vector<double>& parameters) const {
ASSERT_VALID_INDEX(index, bonds);
particles = bonds[index].particles;
parameters = bonds[index].parameters;
}
void CustomCompoundBondForce::setBondParameters(int index, const vector<int>& particles, const vector<double>& parameters) {
ASSERT_VALID_INDEX(index, bonds);
if (particles.size() != particlesPerBond)
throw OpenMMException("CustomCompoundBondForce: wrong number of particles specified for a bond.");
bonds[index].particles = particles;
bonds[index].parameters = parameters;
}
int CustomCompoundBondForce::addFunction(const std::string& name, const std::vector<double>& values, double min, double max) {
if (max <= min)
throw OpenMMException("CustomCompoundBondForce: max <= min for a tabulated function.");
if (values.size() < 2)
throw OpenMMException("CustomCompoundBondForce: a tabulated function must have at least two points");
functions.push_back(FunctionInfo(name, values, min, max));
return functions.size()-1;
}
void CustomCompoundBondForce::getFunctionParameters(int index, std::string& name, std::vector<double>& values, double& min, double& max) const {
ASSERT_VALID_INDEX(index, functions);
name = functions[index].name;
values = functions[index].values;
min = functions[index].min;
max = functions[index].max;
}
void CustomCompoundBondForce::setFunctionParameters(int index, const std::string& name, const std::vector<double>& values, double min, double max) {
if (max <= min)
throw OpenMMException("CustomCompoundBondForce: max <= min for a tabulated function.");
if (values.size() < 2)
throw OpenMMException("CustomCompoundBondForce: a tabulated function must have at least two points");
ASSERT_VALID_INDEX(index, functions);
functions[index].name = name;
functions[index].values = values;
functions[index].min = min;
functions[index].max = max;
}
ForceImpl* CustomCompoundBondForce::createImpl() {
return new CustomCompoundBondForceImpl(*this);
}
openmmapi/src/CustomCompoundBondForceImpl.cpp 0 → 100644
View file @ 35aef079
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2012 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/OpenMMException.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/internal/CustomCompoundBondForceImpl.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;
using std::set;
using std::string;
using std::stringstream;
/**
* This class serves as a placeholder for angles and dihedrals in expressions.
*/
class CustomCompoundBondForceImpl::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);
}
};
CustomCompoundBondForceImpl::CustomCompoundBondForceImpl(CustomCompoundBondForce& owner) : owner(owner) {
}
CustomCompoundBondForceImpl::~CustomCompoundBondForceImpl() {
}
void CustomCompoundBondForceImpl::initialize(ContextImpl& context) {
kernel = context.getPlatform().createKernel(CalcCustomCompoundBondForceKernel::Name(), context);
// Check for errors in the specification of parameters and exclusions.
System& system = context.getSystem();
vector<int> particles;
vector<double> parameters;
int numBondParameters = owner.getNumPerBondParameters();
for (int i = 0; i < owner.getNumBonds(); i++) {
owner.getBondParameters(i, particles, parameters);
for (int j = 0; j < (int) particles.size(); j++)
if (particles[j] < 0 || particles[j] >= system.getNumParticles()) {
stringstream msg;
msg << "CustomCompoundBondForce: Illegal particle index for a bond: ";
msg << particles[j];
throw OpenMMException(msg.str());
}
if (parameters.size() != numBondParameters) {
stringstream msg;
msg << "CustomCompoundBondForce: Wrong number of parameters for bond ";
msg << i;
throw OpenMMException(msg.str());
}
}
dynamic_cast<CalcCustomCompoundBondForceKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner);
}
double CustomCompoundBondForceImpl::calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups) {
if ((groups&(1<<owner.getForceGroup())) != 0)
return dynamic_cast<CalcCustomCompoundBondForceKernel&>(kernel.getImpl()).execute(context, includeForces, includeEnergy);
return 0.0;
}
vector<string> CustomCompoundBondForceImpl::getKernelNames() {
vector<string> names;
names.push_back(CalcCustomCompoundBondForceKernel::Name());
return names;
}
map<string, double> CustomCompoundBondForceImpl::getDefaultParameters() {
map<string, double> parameters;
for (int i = 0; i < owner.getNumGlobalParameters(); i++)
parameters[owner.getGlobalParameterName(i)] = owner.getGlobalParameterDefaultValue(i);
return parameters;
}
ParsedExpression CustomCompoundBondForceImpl::prepareExpression(const CustomCompoundBondForce& force, const map<string, CustomFunction*>& customFunctions, map<string, vector<int> >& distances,
map<string, vector<int> >& angles, map<string, vector<int> >& dihedrals) {
CustomCompoundBondForceImpl::FunctionPlaceholder custom(1);
CustomCompoundBondForceImpl::FunctionPlaceholder distance(2);
CustomCompoundBondForceImpl::FunctionPlaceholder angle(3);
CustomCompoundBondForceImpl::FunctionPlaceholder dihedral(4);
map<string, CustomFunction*> functions = customFunctions;
functions["distance"] = &distance;
functions["angle"] = &angle;
functions["dihedral"] = &dihedral;
ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction(), functions);
map<string, int> atoms;
for (int i = 0; i < force.getNumParticlesPerBond(); i++) {
stringstream name;
name << 'p' << (i+1);
atoms[name.str()] = i;
}
return ParsedExpression(replaceFunctions(expression.getRootNode(), atoms, distances, angles, dihedrals)).optimize();
}
ExpressionTreeNode CustomCompoundBondForceImpl::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("CustomCompoundBondForce: 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/ReferenceKernelFactory.cpp
View file @ 35aef079
...@@ -78,6 +78,8 @@ KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Pla ...@@ -78,6 +78,8 @@ KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Pla
return new ReferenceCalcCustomExternalForceKernel(name, platform); return new ReferenceCalcCustomExternalForceKernel(name, platform);
if (name == CalcCustomHbondForceKernel::Name()) if (name == CalcCustomHbondForceKernel::Name())
return new ReferenceCalcCustomHbondForceKernel(name, platform); return new ReferenceCalcCustomHbondForceKernel(name, platform);
if (name == CalcCustomCompoundBondForceKernel::Name())
return new ReferenceCalcCustomCompoundBondForceKernel(name, platform);
if (name == IntegrateVerletStepKernel::Name()) if (name == IntegrateVerletStepKernel::Name())
return new ReferenceIntegrateVerletStepKernel(name, platform, data); return new ReferenceIntegrateVerletStepKernel(name, platform, data);
if (name == IntegrateLangevinStepKernel::Name()) if (name == IntegrateLangevinStepKernel::Name())
......
platforms/reference/src/ReferenceKernels.cpp
View file @ 35aef079
...@@ -40,6 +40,7 @@ ...@@ -40,6 +40,7 @@
#include "SimTKReference/ReferenceCMAPTorsionIxn.h" #include "SimTKReference/ReferenceCMAPTorsionIxn.h"
#include "SimTKReference/ReferenceCustomAngleIxn.h" #include "SimTKReference/ReferenceCustomAngleIxn.h"
#include "SimTKReference/ReferenceCustomBondIxn.h" #include "SimTKReference/ReferenceCustomBondIxn.h"
#include "SimTKReference/ReferenceCustomCompoundBondIxn.h"
#include "SimTKReference/ReferenceCustomDynamics.h" #include "SimTKReference/ReferenceCustomDynamics.h"
#include "SimTKReference/ReferenceCustomExternalIxn.h" #include "SimTKReference/ReferenceCustomExternalIxn.h"
#include "SimTKReference/ReferenceCustomGBIxn.h" #include "SimTKReference/ReferenceCustomGBIxn.h"
...@@ -62,6 +63,7 @@ ...@@ -62,6 +63,7 @@
#include "openmm/System.h" #include "openmm/System.h"
#include "openmm/internal/AndersenThermostatImpl.h" #include "openmm/internal/AndersenThermostatImpl.h"
#include "openmm/internal/ContextImpl.h" #include "openmm/internal/ContextImpl.h"
#include "openmm/internal/CustomCompoundBondForceImpl.h"
#include "openmm/internal/CustomHbondForceImpl.h" #include "openmm/internal/CustomHbondForceImpl.h"
#include "openmm/internal/CMAPTorsionForceImpl.h" #include "openmm/internal/CMAPTorsionForceImpl.h"
#include "openmm/internal/NonbondedForceImpl.h" #include "openmm/internal/NonbondedForceImpl.h"
...@@ -1236,6 +1238,69 @@ double ReferenceCalcCustomHbondForceKernel::execute(ContextImpl& context, bool i ...@@ -1236,6 +1238,69 @@ double ReferenceCalcCustomHbondForceKernel::execute(ContextImpl& context, bool i
return energy; return energy;
} }
ReferenceCalcCustomCompoundBondForceKernel::~ReferenceCalcCustomCompoundBondForceKernel() {
disposeRealArray(bondParamArray, numBonds);
if (ixn != NULL)
delete ixn;
}
void ReferenceCalcCustomCompoundBondForceKernel::initialize(const System& system, const CustomCompoundBondForce& force) {
// Build the arrays.
numBonds = force.getNumBonds();
numParticles = system.getNumParticles();
vector<vector<int> > bondParticles(numBonds);
int numBondParameters = force.getNumPerBondParameters();
bondParamArray = allocateRealArray(numBonds, numBondParameters);
for (int i = 0; i < numBonds; ++i) {
vector<double> parameters;
force.getBondParameters(i, bondParticles[i], parameters);
for (int j = 0; j < numBondParameters; j++)
bondParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
}
// Create custom functions for the tabulated functions.
map<string, Lepton::CustomFunction*> functions;
for (int i = 0; i < force.getNumFunctions(); i++) {
string name;
vector<double> values;
double min, max;
force.getFunctionParameters(i, name, values, min, max);
functions[name] = new ReferenceTabulatedFunction(min, max, values);
}
// Parse the expression and create the object used to calculate the interaction.
map<string, vector<int> > distances;
map<string, vector<int> > angles;
map<string, vector<int> > dihedrals;
Lepton::ParsedExpression energyExpression = CustomCompoundBondForceImpl::prepareExpression(force, functions, distances, angles, dihedrals);
vector<string> bondParameterNames;
for (int i = 0; i < numBondParameters; i++)
bondParameterNames.push_back(force.getPerBondParameterName(i));
for (int i = 0; i < force.getNumGlobalParameters(); i++)
globalParameterNames.push_back(force.getGlobalParameterName(i));
ixn = new ReferenceCustomCompoundBondIxn(force.getNumParticlesPerBond(), bondParticles, energyExpression, bondParameterNames, distances, angles, dihedrals);
// Delete the custom functions.
for (map<string, Lepton::CustomFunction*>::iterator iter = functions.begin(); iter != functions.end(); iter++)
delete iter->second;
}
double ReferenceCalcCustomCompoundBondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
vector<RealVec>& posData = extractPositions(context);
vector<RealVec>& forceData = extractForces(context);
RealOpenMM energy = 0;
map<string, double> globalParameters;
for (int i = 0; i < (int) globalParameterNames.size(); i++)
globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
ixn->calculatePairIxn(posData, bondParamArray, globalParameters, forceData, includeEnergy ? &energy : NULL);
return energy;
}
ReferenceIntegrateVerletStepKernel::~ReferenceIntegrateVerletStepKernel() { ReferenceIntegrateVerletStepKernel::~ReferenceIntegrateVerletStepKernel() {
if (dynamics) if (dynamics)
delete dynamics; delete dynamics;
......
platforms/reference/src/ReferenceKernels.h
View file @ 35aef079
...@@ -41,6 +41,7 @@ ...@@ -41,6 +41,7 @@
class CpuObc; class CpuObc;
class CpuGBVI; class CpuGBVI;
class ReferenceAndersenThermostat; class ReferenceAndersenThermostat;
class ReferenceCustomCompoundBondIxn;
class ReferenceCustomHbondIxn; class ReferenceCustomHbondIxn;
class ReferenceBrownianDynamics; class ReferenceBrownianDynamics;
class ReferenceStochasticDynamics; class ReferenceStochasticDynamics;
...@@ -712,6 +713,37 @@ private: ...@@ -712,6 +713,37 @@ private:
std::vector<std::string> globalParameterNames; std::vector<std::string> globalParameterNames;
}; };
/**
* This kernel is invoked by CustomCompoundBondForce to calculate the forces acting on the system.
*/
class ReferenceCalcCustomCompoundBondForceKernel : public CalcCustomCompoundBondForceKernel {
public:
ReferenceCalcCustomCompoundBondForceKernel(std::string name, const Platform& platform) : CalcCustomCompoundBondForceKernel(name, platform), ixn(NULL) {
}
~ReferenceCalcCustomCompoundBondForceKernel();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param force the CustomCompoundBondForce this kernel will be used for
*/
void initialize(const System& system, const CustomCompoundBondForce& force);
/**
* Execute the kernel to calculate the forces and/or energy.
*
* @param context the context in which to execute this kernel
* @param includeForces true if forces should be calculated
* @param includeEnergy true if the energy should be calculated
* @return the potential energy due to the force
*/
double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
private:
int numBonds, numParticles;
RealOpenMM **bondParamArray;
ReferenceCustomCompoundBondIxn* ixn;
std::vector<std::string> globalParameterNames;
};
/** /**
* This kernel is invoked by VerletIntegrator to take one time step. * This kernel is invoked by VerletIntegrator to take one time step.
*/ */
......
platforms/reference/src/ReferencePlatform.cpp
View file @ 35aef079
...@@ -61,6 +61,7 @@ ReferencePlatform::ReferencePlatform() { ...@@ -61,6 +61,7 @@ ReferencePlatform::ReferencePlatform() {
registerKernelFactory(CalcCustomGBForceKernel::Name(), factory); registerKernelFactory(CalcCustomGBForceKernel::Name(), factory);
registerKernelFactory(CalcCustomExternalForceKernel::Name(), factory); registerKernelFactory(CalcCustomExternalForceKernel::Name(), factory);
registerKernelFactory(CalcCustomHbondForceKernel::Name(), factory); registerKernelFactory(CalcCustomHbondForceKernel::Name(), factory);
registerKernelFactory(CalcCustomCompoundBondForceKernel::Name(), factory);
registerKernelFactory(IntegrateVerletStepKernel::Name(), factory); registerKernelFactory(IntegrateVerletStepKernel::Name(), factory);
registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory); registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);
registerKernelFactory(IntegrateBrownianStepKernel::Name(), factory); registerKernelFactory(IntegrateBrownianStepKernel::Name(), factory);
......
platforms/reference/src/SimTKReference/ReferenceCustomCompoundBondIxn.cpp 0 → 100644
View file @ 35aef079
/* Portions copyright (c) 2009-2010 Stanford University and Simbios.
* Contributors: Peter Eastman
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject
* to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#include <string.h>
#include <sstream>
#include <utility>
#include "../SimTKUtilities/SimTKOpenMMCommon.h"
#include "../SimTKUtilities/SimTKOpenMMLog.h"
#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
#include "ReferenceForce.h"
#include "ReferenceCustomCompoundBondIxn.h"
using std::map;
using std::pair;
using std::string;
using std::stringstream;
using std::vector;
using OpenMM::RealVec;
/**---------------------------------------------------------------------------------------
ReferenceCustomCompoundBondIxn constructor
--------------------------------------------------------------------------------------- */
ReferenceCustomCompoundBondIxn::ReferenceCustomCompoundBondIxn(int numParticlesPerBond, const vector<vector<int> >& bondAtoms,
const Lepton::ParsedExpression& energyExpression, const vector<string>& bondParameterNames,
const map<string, vector<int> >& distances, const map<string, vector<int> >& angles, const map<string, vector<int> >& dihedrals) :
bondAtoms(bondAtoms), energyExpression(energyExpression.createProgram()), bondParamNames(bondParameterNames) {
for (int i = 0; i < numParticlesPerBond; i++) {
stringstream xname, yname, zname;
xname << 'x' << (i+1);
yname << 'y' << (i+1);
zname << 'z' << (i+1);
particleTerms.push_back(ReferenceCustomCompoundBondIxn::ParticleTermInfo(xname.str(), i, 0, energyExpression.differentiate(xname.str()).optimize().createProgram()));
particleTerms.push_back(ReferenceCustomCompoundBondIxn::ParticleTermInfo(yname.str(), i, 1, energyExpression.differentiate(yname.str()).optimize().createProgram()));
particleTerms.push_back(ReferenceCustomCompoundBondIxn::ParticleTermInfo(zname.str(), i, 2, energyExpression.differentiate(zname.str()).optimize().createProgram()));
}
for (map<string, vector<int> >::const_iterator iter = distances.begin(); iter != distances.end(); ++iter)
distanceTerms.push_back(ReferenceCustomCompoundBondIxn::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(ReferenceCustomCompoundBondIxn::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(ReferenceCustomCompoundBondIxn::DihedralTermInfo(iter->first, iter->second, energyExpression.differentiate(iter->first).optimize().createProgram()));
}
/**---------------------------------------------------------------------------------------
ReferenceCustomCompoundBondIxn destructor
--------------------------------------------------------------------------------------- */
ReferenceCustomCompoundBondIxn::~ReferenceCustomCompoundBondIxn( ){
}
/**---------------------------------------------------------------------------------------
Calculate custom hbond interaction
@param atomCoordinates atom coordinates
@param bondParameters bond parameters values bondParameters[bondIndex][parameterIndex]
@param globalParameters the values of global parameters
@param forces force array (forces added)
@param totalEnergy total energy
--------------------------------------------------------------------------------------- */
void ReferenceCustomCompoundBondIxn::calculatePairIxn(vector<RealVec>& atomCoordinates, RealOpenMM** bondParameters,
const map<string, double>& globalParameters, vector<RealVec>& forces,
RealOpenMM* totalEnergy) const {
map<string, double> variables = globalParameters;
int numBonds = bondAtoms.size();
for (int bond = 0; bond < numBonds; bond++){
for (int j = 0; j < (int) bondParamNames.size(); j++)
variables[bondParamNames[j]] = bondParameters[bond][j];
calculateOneIxn(bond, atomCoordinates, variables, forces, totalEnergy);
}
}
/**---------------------------------------------------------------------------------------
Calculate interaction for one bond
@param bond the index of the bond
@param atomCoordinates atom coordinates
@param variables the values of variables that may appear in expressions
@param forces force array (forces added)
@param energyByAtom atom energy
@param totalEnergy total energy
--------------------------------------------------------------------------------------- */
void ReferenceCustomCompoundBondIxn::calculateOneIxn(int bond, vector<RealVec>& atomCoordinates,
map<string, double>& variables, vector<RealVec>& forces, RealOpenMM* totalEnergy) const {
// ---------------------------------------------------------------------------------------
static const std::string methodName = "\nReferenceCustomCompoundBondIxn::calculateOneIxn";
// ---------------------------------------------------------------------------------------
// Compute all of the variables the energy can depend on.
const vector<int>& atoms = bondAtoms[0];
for (int i = 0; i < (int) particleTerms.size(); i++) {
const ParticleTermInfo& term = particleTerms[i];
variables[term.name] = atomCoordinates[term.atom][term.component];
}
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];
}
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);
}
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;
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 individual particle coordinates.
for (int i = 0; i < (int) particleTerms.size(); i++) {
const ParticleTermInfo& term = particleTerms[i];
forces[atoms[term.atom]][term.component] -= term.forceExpression.evaluate(variables);
}
// Apply forces based on distances.
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*term.delta[i];
forces[atoms[term.p1]][i] -= force;
forces[atoms[term.p2]][i] += force;
}
}
// Apply forces based on angles.
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(term.delta1, term.delta2, thetaCross);
RealOpenMM lengthThetaCross = SQRT(DOT3(thetaCross, thetaCross));
if (lengthThetaCross < 1.0e-06)
lengthThetaCross = (RealOpenMM) 1.0e-06;
RealOpenMM termA = dEdTheta/(term.delta1[ReferenceForce::R2Index]*lengthThetaCross);
RealOpenMM termC = -dEdTheta/(term.delta2[ReferenceForce::R2Index]*lengthThetaCross);
RealOpenMM deltaCrossP[3][3];
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[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 dihedrals.
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(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]*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[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];
}
}
// Add the energy
if (totalEnergy)
*totalEnergy += (RealOpenMM) energyExpression.evaluate(variables);
}
void ReferenceCustomCompoundBondIxn::computeDelta(int atom1, int atom2, RealOpenMM* delta, vector<RealVec>& atomCoordinates) const {
ReferenceForce::getDeltaR(atomCoordinates[atom1], atomCoordinates[atom2], delta);
}
RealOpenMM ReferenceCustomCompoundBondIxn::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)
angle = 0;
else if (cosine <= -1)
angle = PI_M;
else
angle = ACOS(cosine);
return angle;
}
platforms/reference/src/SimTKReference/ReferenceCustomCompoundBondIxn.h 0 → 100644
View file @ 35aef079
/* Portions copyright (c) 2009-2012 Stanford University and Simbios.
* Contributors: Peter Eastman
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject
* to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
* LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
* OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
* WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#ifndef __ReferenceCustomCompoundBondIxn_H__
#define __ReferenceCustomCompoundBondIxn_H__
#include "ReferenceBondIxn.h"
#include "lepton/ExpressionProgram.h"
#include "lepton/ParsedExpression.h"
#include <map>
#include <vector>
// ---------------------------------------------------------------------------------------
class ReferenceCustomCompoundBondIxn : public ReferenceBondIxn {
private:
class ParticleTermInfo;
class DistanceTermInfo;
class AngleTermInfo;
class DihedralTermInfo;
std::vector<std::vector<int> > bondAtoms;
Lepton::ExpressionProgram energyExpression;
std::vector<std::string> bondParamNames;
std::vector<ParticleTermInfo> particleTerms;
std::vector<DistanceTermInfo> distanceTerms;
std::vector<AngleTermInfo> angleTerms;
std::vector<DihedralTermInfo> dihedralTerms;
/**---------------------------------------------------------------------------------------
Calculate custom interaction for one bond
@param bond the index of the bond
@param atomCoordinates atom coordinates
@param variables the values of variables that may appear in expressions
@param forces force array (forces added)
@param totalEnergy total energy
--------------------------------------------------------------------------------------- */
void calculateOneIxn(int bond, std::vector<OpenMM::RealVec>& atomCoordinates,
std::map<std::string, double>& variables, std::vector<OpenMM::RealVec>& forces,
RealOpenMM* totalEnergy) const;
void computeDelta(int atom1, int atom2, RealOpenMM* delta, std::vector<OpenMM::RealVec>& atomCoordinates) const;
static RealOpenMM computeAngle(RealOpenMM* vec1, RealOpenMM* vec2);
public:
/**---------------------------------------------------------------------------------------
Constructor
--------------------------------------------------------------------------------------- */
ReferenceCustomCompoundBondIxn(int numParticlesPerBond, const std::vector<std::vector<int> >& bondAtoms, const Lepton::ParsedExpression& energyExpression,
const std::vector<std::string>& bondParameterNames, 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);
/**---------------------------------------------------------------------------------------
Destructor
--------------------------------------------------------------------------------------- */
~ReferenceCustomCompoundBondIxn();
/**---------------------------------------------------------------------------------------
Calculate custom compound bond interaction
@param atomCoordinates atom coordinates
@param bondParameters bond parameters values bondParameters[bondIndex][parameterIndex]
@param globalParameters the values of global parameters
@param forces force array (forces added)
@param totalEnergy total energy
--------------------------------------------------------------------------------------- */
void calculatePairIxn(std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** bondParameters,
const std::map<std::string, double>& globalParameters,
std::vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy) const;
// ---------------------------------------------------------------------------------------
};
class ReferenceCustomCompoundBondIxn::ParticleTermInfo {
public:
std::string name;
int atom, component;
Lepton::ExpressionProgram forceExpression;
ParticleTermInfo(const std::string& name, int atom, int component, const Lepton::ExpressionProgram& forceExpression) :
name(name), atom(atom), component(component), forceExpression(forceExpression) {
}
};
class ReferenceCustomCompoundBondIxn::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 ReferenceCustomCompoundBondIxn::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 ReferenceCustomCompoundBondIxn::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 // __ReferenceCustomCompoundBondIxn_H__
platforms/reference/tests/TestReferenceCustomCompoundBondForce.cpp 0 → 100644
View file @ 35aef079
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2012 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the reference implementation of CustomCompoundBondForce.
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "ReferencePlatform.h"
#include "openmm/CustomCompoundBondForce.h"
#include "openmm/HarmonicAngleForce.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/PeriodicTorsionForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testBond() {
ReferencePlatform platform;
// Create a system using a CustomCompoundBondForce.
System customSystem;
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
CustomCompoundBondForce* custom = new CustomCompoundBondForce(4, "0.5*kb*((distance(p1,p2)-b0)^2+(distance(p2,p3)-b0)^2)+0.5*ka*(angle(p2,p3,p4)-a0)^2+kt*(1+cos(dihedral(p1,p2,p3,p4)-t0))");
custom->addPerBondParameter("kb");
custom->addPerBondParameter("ka");
custom->addPerBondParameter("kt");
custom->addPerBondParameter("b0");
custom->addPerBondParameter("a0");
custom->addPerBondParameter("t0");
vector<int> particles(4);
particles[0] = 0;
particles[1] = 1;
particles[2] = 3;
particles[3] = 2;
vector<double> parameters(6);
parameters[0] = 1.5;
parameters[1] = 0.8;
parameters[2] = 0.6;
parameters[3] = 1.1;
parameters[4] = 2.9;
parameters[5] = 1.3;
custom->addBond(particles, parameters);
customSystem.addForce(custom);
// Create an identical system using standard forces.
System standardSystem;
standardSystem.addParticle(1.0);
standardSystem.addParticle(1.0);
standardSystem.addParticle(1.0);
standardSystem.addParticle(1.0);
HarmonicBondForce* bonds = new HarmonicBondForce();
bonds->addBond(0, 1, 1.1, 1.5);
bonds->addBond(1, 3, 1.1, 1.5);
standardSystem.addForce(bonds);
HarmonicAngleForce* angles = new HarmonicAngleForce();
angles->addAngle(1, 3, 2, 2.9, 0.8);
standardSystem.addForce(angles);
PeriodicTorsionForce* torsions = new PeriodicTorsionForce();
torsions->addTorsion(0, 1, 3, 2, 1, 1.3, 0.6);
standardSystem.addForce(torsions);
// Set the atoms in various positions, and verify that both systems give identical forces and energy.
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<Vec3> positions(4);
for (int i = 0; i < 10; i++) {
VerletIntegrator integrator1(0.01);
VerletIntegrator integrator2(0.01);
Context c1(customSystem, integrator1, platform);
Context c2(standardSystem, integrator2, platform);
for (int j = 0; j < (int) positions.size(); j++)
positions[j] = Vec3(5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt));
c1.setPositions(positions);
c2.setPositions(positions);
State s1 = c1.getState(State::Forces | State::Energy);
State s2 = c2.getState(State::Forces | State::Energy);
for (int i = 0; i < customSystem.getNumParticles(); i++)
ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], TOL);
ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), TOL);
}
}
void testPositionDependence() {
ReferencePlatform platform;
System customSystem;
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
CustomCompoundBondForce* custom = new CustomCompoundBondForce(2, "scale1*distance(p1,p2)+scale2*x1+2*y2");
custom->addGlobalParameter("scale1", 0.3);
custom->addGlobalParameter("scale2", 0.2);
vector<int> particles(2);
particles[0] = 0;
particles[1] = 1;
vector<double> parameters;
custom->addBond(particles, parameters);
customSystem.addForce(custom);
vector<Vec3> positions(2);
positions[0] = Vec3(0.5, 1, 0);
positions[1] = Vec3(1.5, 1, 0);
VerletIntegrator integrator(0.01);
Context context(customSystem, integrator, platform);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
ASSERT_EQUAL_TOL(0.3*1.0+0.2*0.5+2*1, state.getPotentialEnergy(), 1e-5);
ASSERT_EQUAL_VEC(Vec3(0.3-0.2, 0, 0), state.getForces()[0], 1e-5);
ASSERT_EQUAL_VEC(Vec3(-0.3, -2, 0), state.getForces()[1], 1e-5);
}
int main() {
try {
testBond();
testPositionDependence();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
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