Skip to content

GitLab

Commit f1603534 authored by peastman's avatar peastman
Browse files

Created API for CustomCentroidBondForce

parent 5866a98f
Show whitespace changes
Inline Side-by-side
Showing with 948 additions and 3 deletions
olla/include/openmm/kernels.h
View file @ f1603534
......@@ -38,6 +38,7 @@
#include "openmm/CMMotionRemover.h"
#include "openmm/CustomAngleForce.h"
#include "openmm/CustomBondForce.h"
#include "openmm/CustomCentroidBondForce.h"
#include "openmm/CustomCompoundBondForce.h"
#include "openmm/CustomExternalForce.h"
#include "openmm/CustomGBForce.h"
......@@ -801,6 +802,41 @@ public:
virtual void copyParametersToContext(ContextImpl& context, const CustomHbondForce& force) = 0;
};
/**
* This kernel is invoked by CustomCentroidBondForce to calculate the forces acting on the system and the energy of the system.
*/
class CalcCustomCentroidBondForceKernel : public KernelImpl {
public:
static std::string Name() {
return "CalcCustomCentroidBondForce";
}
CalcCustomCentroidBondForceKernel(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 CustomCentroidBondForce this kernel will be used for
*/
virtual void initialize(const System& system, const CustomCentroidBondForce& 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;
/**
* Copy changed parameters over to a context.
*
* @param context the context to copy parameters to
* @param force the CustomCentroidBondForce to copy the parameters from
*/
virtual void copyParametersToContext(ContextImpl& context, const CustomCentroidBondForce& force) = 0;
};
/**
* This kernel is invoked by CustomCompoundBondForce to calculate the forces acting on the system and the energy of the system.
*/
......
openmmapi/include/openmm/CustomCentroidBondForce.h 0 → 100644
View file @ f1603534
#ifndef OPENMM_CUSTOMCENTROIDBONDFORCE_H_
#define OPENMM_CUSTOMCENTROIDBONDFORCE_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-2015 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 "TabulatedFunction.h"
#include "Force.h"
#include "Vec3.h"
#include <vector>
#include "internal/windowsExport.h"
namespace OpenMM {
/**
* This class is similar to CustomCompoundBondForce, but instead of applying forces between individual particles,
* it applies them between the centers of groups of particles. This is useful for a variety of purposes, such as
* restraints to keep two molecules from moving too far apart.
*
* When using this class, you define groups of particles, and the center of each group is calculated as a weighted
* average of the particle positions. By default, the particle masses are used as weights, so the center position
* is the center of mass. You can optionally specify different weights to use. You then add bonds just as with
* CustomCompoundBondForce, but instead of specifying the particles that make up a bond, you specify the groups.
*
* When creating a CustomCentroidBondForce, you specify the number of groups involved in a bond, and an expression
* for the energy of each bond. It may depend on the center positions of individual groups, the distances between
* the centers of pairs of groups, the angles formed by sets of three groups, and the dihedral angles formed by
* sets of four groups.
*
* We refer to the groups in a bond as g1, g2, g3, etc. For each bond, CustomCentroidBondForce 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 centers of the groups. For example, x1
* is the x coordinate of the center of group g1, and y3 is the y coordinate of the center of group g3.</li>
* <li>distance(g1, g2): the distance between the centers of groups g1 and g2 (where "g1" and "g2" may be replaced
* by the names of whichever groups you want to calculate the distance between).</li>
* <li>angle(g1, g2, g3): the angle formed by the centers of the three specified groups.</li>
* <li>dihedral(g1, g2, g3, g4): the dihedral angle formed by the centers of the four specified groups.</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 CustomCentroidBondForce 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 addGroup() to define the particle groups. Each group is specified by the particles it contains, and
* the weights to use when computing the center position.
*
* Then call addBond() to define bonds and specify their parameter values. After a bond has been added, you can
* modify its parameters by calling setBondParameters(). This will have no effect on Contexts that already exist unless
* you call updateParametersInContext().
*
* As an example, the following code creates a CustomCentroidBondForce that implements a harmonic force between the
* centers of mass of two groups of particles.
*
* <tt><pre>
* CustomCentroidBondForce* force = new CustomCentroidBondForce(2, "0.5*k*distance(g1,g2)^2");
* force->addPerBondParameter("k");
* force->addGroup(particles1);
* force->addGroup(particles2);
* vector<int> bondGroups;
* bondGroups.push_back(0);
* bondGroups.push_back(1);
* vector<double> bondParameters;
* bondParameters.push_back(k);
* force->addBond(bondGroups, bondParameters);
* </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, floor, ceil, step, delta, select. 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.
* select(x,y,z) = z if x = 0, y otherwise.
*
* In addition, you can call addTabulatedFunction() to define a new function based on tabulated values. You specify the function by
* creating a TabulatedFunction object. That function can then appear in the expression.
*/
class OPENMM_EXPORT CustomCentroidBondForce : public Force {
public:
/**
* Create a CustomCentroidBondForce.
*
* @param numGroups the number of groups 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 CustomCentroidBondForce(int numGroups, const std::string& energy);
~CustomCentroidBondForce();
/**
* Get the number of groups used to define each bond.
*/
int getNumGroupsPerBond() const {
return groupsPerBond;
}
/**
* Get the number of particle groups that have been defined.
*/
int getNumGroups() const {
return groups.size();
}
/**
* 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 getNumTabulatedFunctions() const {
return functions.size();
}
/**
* Get the number of tabulated functions that have been defined.
*
* @deprecated This method exists only for backward compatibility. Use getNumTabulatedFunctions() instead.
*/
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 particle group.
*
* @param particles the indices of the particles to include in the group
* @param weights the weight to use for each particle when computing the center position.
* If this is omitted, then particle masses will be used as weights.
* @return the index of the group that was added
*/
int addGroup(const std::vector<int>& particles, const std::vector<double>& weights = std::vector<double>());
/**
* Get the properties of a group.
*
* @param index the index of the group to get
* @param particles the indices of the particles in the group
* @param weights the weight used for each particle when computing the center position.
* If no weights were specified, this vector will be empty indicating that particle
* masses should be used as weights.
*/
void getGroupParameters(int index, std::vector<int>& particles, std::vector<double>& weights) const;
/**
* Set the properties of a group.
*
* @param index the index of the group to set
* @param particles the indices of the particles in the group
* @param weights the weight to use for each particle when computing the center position.
* If this is omitted, then particle masses will be used as weights.
*/
void setGroupParameters(int index, const std::vector<int>& particles, const std::vector<double>& weights = std::vector<double>());
/**
* Add a bond to the force
*
* @param groups the indices of the groups 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>& groups, const std::vector<double>& parameters);
/**
* Get the properties of a bond.
*
* @param index the index of the bond to get
* @param groups the indices of the groups in the bond
* @param parameters the list of per-bond parameter values for the bond
*/
void getBondParameters(int index, std::vector<int>& groups, std::vector<double>& parameters) const;
/**
* Set the properties of a bond.
*
* @param index the index of the bond to set
* @param groups the indices of the groups in the bond
* @param parameters the list of per-bond parameter values for the bond
*/
void setBondParameters(int index, const std::vector<int>& groups, 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 function a TabulatedFunction object defining the function. The TabulatedFunction
* should have been created on the heap with the "new" operator. The
* Force takes over ownership of it, and deletes it when the Force itself is deleted.
* @return the index of the function that was added
*/
int addTabulatedFunction(const std::string& name, TabulatedFunction* function);
/**
* Get a const reference to a tabulated function that may appear in the energy expression.
*
* @param index the index of the function to get
* @return the TabulatedFunction object defining the function
*/
const TabulatedFunction& getTabulatedFunction(int index) const;
/**
* Get a reference to a tabulated function that may appear in the energy expression.
*
* @param index the index of the function to get
* @return the TabulatedFunction object defining the function
*/
TabulatedFunction& getTabulatedFunction(int index);
/**
* Get the name of a tabulated function that may appear in the energy expression.
*
* @param index the index of the function to get
* @return the name of the function as it appears in expressions
*/
const std::string& getTabulatedFunctionName(int index) const;
/**
* Update the per-bond parameters in a Context to match those stored in this Force object. This method provides
* an efficient method to update certain parameters in an existing Context without needing to reinitialize it.
* Simply call setBondParameters() to modify this object's parameters, then call updateParametersInContext()
* to copy them over to the Context.
*
* This method has several limitations. The only information it updates is the values of per-bond parameters.
* All other aspects of the Force (such as the energy function) are unaffected and can only be changed by reinitializing
* the Context. The set of groups involved in a bond cannot be changed, nor can new bonds be added.
*/
void updateParametersInContext(Context& context);
/**
* Returns whether or not this force makes use of periodic boundary
* conditions.
*
* @returns false
*/
bool usesPeriodicBoundaryConditions() const {
return false;
}
protected:
ForceImpl* createImpl() const;
private:
class GroupInfo;
class BondInfo;
class BondParameterInfo;
class GlobalParameterInfo;
class FunctionInfo;
int groupsPerBond;
std::string energyExpression;
std::vector<BondParameterInfo> bondParameters;
std::vector<GlobalParameterInfo> globalParameters;
std::vector<GroupInfo> groups;
std::vector<BondInfo> bonds;
std::vector<FunctionInfo> functions;
};
/**
* This is an internal class used to record information about a group.
* @private
*/
class CustomCentroidBondForce::GroupInfo {
public:
std::vector<int> particles;
std::vector<double> weights;
GroupInfo() {
}
GroupInfo(const std::vector<int>& particles, const std::vector<double>& weights) :
particles(particles), weights(weights) {
}
};
/**
* This is an internal class used to record information about a bond.
* @private
*/
class CustomCentroidBondForce::BondInfo {
public:
std::vector<int> groups;
std::vector<double> parameters;
BondInfo() {
}
BondInfo(const std::vector<int>& groups, const std::vector<double>& parameters) :
groups(groups), parameters(parameters) {
}
};
/**
* This is an internal class used to record information about a per-bond parameter.
* @private
*/
class CustomCentroidBondForce::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 CustomCentroidBondForce::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 CustomCentroidBondForce::FunctionInfo {
public:
std::string name;
TabulatedFunction* function;
FunctionInfo() {
}
FunctionInfo(const std::string& name, TabulatedFunction* function) : name(name), function(function) {
}
};
} // namespace OpenMM
#endif /*OPENMM_CUSTOMCENTROIDBONDFORCE_H_*/
openmmapi/include/openmm/CustomCompoundBondForce.h
View file @ f1603534
......@@ -231,7 +231,7 @@ public:
/**
* Set the properties of a bond.
*
* @param index the index of the bond group to set
* @param index the index of the bond to set
* @param particles the indices of the particles in the bond
* @param parameters the list of per-bond parameter values for the bond
*/
......@@ -323,7 +323,7 @@ private:
};
/**
* This is an internal class used to record information about a bond or acceptor.
* This is an internal class used to record information about a bond.
* @private
*/
class CustomCompoundBondForce::BondInfo {
......@@ -338,7 +338,7 @@ public:
};
/**
* This is an internal class used to record information about a per-bond or per-acceptor parameter.
* This is an internal class used to record information about a per-bond parameter.
* @private
*/
class CustomCompoundBondForce::BondParameterInfo {
......
openmmapi/include/openmm/internal/CustomCentroidBondForceImpl.h 0 → 100644
View file @ f1603534
#ifndef OPENMM_CUSTOMCENTROIDBONDFORCEIMPL_H_
#define OPENMM_CUSTOMCENTROIDBONDFORCEIMPL_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-2015 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/CustomCentroidBondForce.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 CustomCentroidBondForce.
*/
class OPENMM_EXPORT CustomCentroidBondForceImpl : public ForceImpl {
public:
CustomCentroidBondForceImpl(const CustomCentroidBondForce& owner);
~CustomCentroidBondForceImpl();
void initialize(ContextImpl& context);
const CustomCentroidBondForce& getOwner() const {
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();
void updateParametersInContext(ContextImpl& context);
/**
* 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 CustomCentroidBondForce 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 CustomCentroidBondForce& 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);
const CustomCentroidBondForce& owner;
Kernel kernel;
};
} // namespace OpenMM
#endif /*OPENMM_CUSTOMCENTROIDBONDFORCEIMPL_H_*/
openmmapi/src/CustomCentroidBondForce.cpp 0 → 100644
View file @ f1603534
/* -------------------------------------------------------------------------- *
* 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) 2015 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/CustomCentroidBondForce.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/internal/CustomCentroidBondForceImpl.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;
CustomCentroidBondForce::CustomCentroidBondForce(int numGroups, const string& energy) : groupsPerBond(numGroups), energyExpression(energy) {
}
CustomCentroidBondForce::~CustomCentroidBondForce() {
for (int i = 0; i < (int) functions.size(); i++)
delete functions[i].function;
}
const string& CustomCentroidBondForce::getEnergyFunction() const {
return energyExpression;
}
void CustomCentroidBondForce::setEnergyFunction(const std::string& energy) {
energyExpression = energy;
}
int CustomCentroidBondForce::addPerBondParameter(const string& name) {
bondParameters.push_back(BondParameterInfo(name));
return bondParameters.size()-1;
}
const string& CustomCentroidBondForce::getPerBondParameterName(int index) const {
ASSERT_VALID_INDEX(index, bondParameters);
return bondParameters[index].name;
}
void CustomCentroidBondForce::setPerBondParameterName(int index, const string& name) {
ASSERT_VALID_INDEX(index, bondParameters);
bondParameters[index].name = name;
}
int CustomCentroidBondForce::addGlobalParameter(const string& name, double defaultValue) {
globalParameters.push_back(GlobalParameterInfo(name, defaultValue));
return globalParameters.size()-1;
}
const string& CustomCentroidBondForce::getGlobalParameterName(int index) const {
ASSERT_VALID_INDEX(index, globalParameters);
return globalParameters[index].name;
}
void CustomCentroidBondForce::setGlobalParameterName(int index, const string& name) {
ASSERT_VALID_INDEX(index, globalParameters);
globalParameters[index].name = name;
}
double CustomCentroidBondForce::getGlobalParameterDefaultValue(int index) const {
ASSERT_VALID_INDEX(index, globalParameters);
return globalParameters[index].defaultValue;
}
void CustomCentroidBondForce::setGlobalParameterDefaultValue(int index, double defaultValue) {
ASSERT_VALID_INDEX(index, globalParameters);
globalParameters[index].defaultValue = defaultValue;
}
int CustomCentroidBondForce::addGroup(const vector<int>& particles, const vector<double>& weights) {
if (particles.size() != weights.size() && weights.size() > 0)
throw OpenMMException("CustomCentroidBondForce: wrong number of weights specified for a group.");
groups.push_back(GroupInfo(particles, weights));
return groups.size()-1;
}
void CustomCentroidBondForce::getGroupParameters(int index, vector<int>& particles, std::vector<double>& weights) const {
ASSERT_VALID_INDEX(index, groups);
particles = groups[index].particles;
weights = groups[index].weights;
}
void CustomCentroidBondForce::setGroupParameters(int index, const vector<int>& particles, const vector<double>& weights) {
ASSERT_VALID_INDEX(index, groups);
if (particles.size() != weights.size() && weights.size() > 0)
throw OpenMMException("CustomCentroidBondForce: wrong number of weights specified for a group.");
groups[index].particles = particles;
groups[index].weights = weights;
}
int CustomCentroidBondForce::addBond(const vector<int>& groups, const vector<double>& parameters) {
if (groups.size() != groupsPerBond)
throw OpenMMException("CustomCentroidBondForce: wrong number of groups specified for a bond.");
bonds.push_back(BondInfo(groups, parameters));
return bonds.size()-1;
}
void CustomCentroidBondForce::getBondParameters(int index, vector<int>& groups, std::vector<double>& parameters) const {
ASSERT_VALID_INDEX(index, bonds);
groups = bonds[index].groups;
parameters = bonds[index].parameters;
}
void CustomCentroidBondForce::setBondParameters(int index, const vector<int>& groups, const vector<double>& parameters) {
ASSERT_VALID_INDEX(index, bonds);
if (groups.size() != groupsPerBond)
throw OpenMMException("CustomCentroidBondForce: wrong number of groups specified for a bond.");
bonds[index].groups = groups;
bonds[index].parameters = parameters;
}
int CustomCentroidBondForce::addTabulatedFunction(const std::string& name, TabulatedFunction* function) {
functions.push_back(FunctionInfo(name, function));
return functions.size()-1;
}
const TabulatedFunction& CustomCentroidBondForce::getTabulatedFunction(int index) const {
ASSERT_VALID_INDEX(index, functions);
return *functions[index].function;
}
TabulatedFunction& CustomCentroidBondForce::getTabulatedFunction(int index) {
ASSERT_VALID_INDEX(index, functions);
return *functions[index].function;
}
const string& CustomCentroidBondForce::getTabulatedFunctionName(int index) const {
ASSERT_VALID_INDEX(index, functions);
return functions[index].name;
}
ForceImpl* CustomCentroidBondForce::createImpl() const {
return new CustomCentroidBondForceImpl(*this);
}
void CustomCentroidBondForce::updateParametersInContext(Context& context) {
dynamic_cast<CustomCentroidBondForceImpl&>(getImplInContext(context)).updateParametersInContext(getContextImpl(context));
}
openmmapi/src/CustomCentroidBondForceImpl.cpp 0 → 100644
View file @ f1603534
/* -------------------------------------------------------------------------- *
* 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-2015 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/CustomCentroidBondForceImpl.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 CustomCentroidBondForceImpl::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);
}
};
CustomCentroidBondForceImpl::CustomCentroidBondForceImpl(const CustomCentroidBondForce& owner) : owner(owner) {
}
CustomCentroidBondForceImpl::~CustomCentroidBondForceImpl() {
}
void CustomCentroidBondForceImpl::initialize(ContextImpl& context) {
kernel = context.getPlatform().createKernel(CalcCustomCentroidBondForceKernel::Name(), context);
// Check for errors in the specification of parameters and exclusions.
const System& system = context.getSystem();
vector<int> particles;
vector<double> weights;
for (int i = 0; i < owner.getNumGroups(); i++) {
owner.getGroupParameters(i, particles, weights);
for (int j = 0; j < (int) particles.size(); j++)
if (particles[j] < 0 || particles[j] >= system.getNumParticles()) {
stringstream msg;
msg << "CustomCentroidBondForce: Illegal particle index for a group: ";
msg << particles[j];
throw OpenMMException(msg.str());
}
if (weights.size() != particles.size() && weights.size() > 0) {
stringstream msg;
msg << "CustomCentroidBondForce: Wrong number of weights for group ";
msg << i;
throw OpenMMException(msg.str());
}
}
vector<int> groups;
vector<double> parameters;
int numBondParameters = owner.getNumPerBondParameters();
for (int i = 0; i < owner.getNumBonds(); i++) {
owner.getBondParameters(i, groups, parameters);
for (int j = 0; j < (int) groups.size(); j++)
if (groups[j] < 0 || groups[j] >= owner.getNumGroups()) {
stringstream msg;
msg << "CustomCentroidBondForce: Illegal group index for a bond: ";
msg << groups[j];
throw OpenMMException(msg.str());
}
if (parameters.size() != numBondParameters) {
stringstream msg;
msg << "CustomCentroidBondForce: Wrong number of parameters for bond ";
msg << i;
throw OpenMMException(msg.str());
}
}
kernel.getAs<CalcCustomCentroidBondForceKernel>().initialize(context.getSystem(), owner);
}
double CustomCentroidBondForceImpl::calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups) {
if ((groups&(1<<owner.getForceGroup())) != 0)
return kernel.getAs<CalcCustomCentroidBondForceKernel>().execute(context, includeForces, includeEnergy);
return 0.0;
}
vector<string> CustomCentroidBondForceImpl::getKernelNames() {
vector<string> names;
names.push_back(CalcCustomCentroidBondForceKernel::Name());
return names;
}
map<string, double> CustomCentroidBondForceImpl::getDefaultParameters() {
map<string, double> parameters;
for (int i = 0; i < owner.getNumGlobalParameters(); i++)
parameters[owner.getGlobalParameterName(i)] = owner.getGlobalParameterDefaultValue(i);
return parameters;
}
ParsedExpression CustomCentroidBondForceImpl::prepareExpression(const CustomCentroidBondForce& force, const map<string, CustomFunction*>& customFunctions, map<string, vector<int> >& distances,
map<string, vector<int> >& angles, map<string, vector<int> >& dihedrals) {
CustomCentroidBondForceImpl::FunctionPlaceholder custom(1);
CustomCentroidBondForceImpl::FunctionPlaceholder distance(2);
CustomCentroidBondForceImpl::FunctionPlaceholder angle(3);
CustomCentroidBondForceImpl::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> groups;
for (int i = 0; i < force.getNumGroupsPerBond(); i++) {
stringstream name;
name << 'g' << (i+1);
groups[name.str()] = i;
}
return ParsedExpression(replaceFunctions(expression.getRootNode(), groups, distances, angles, dihedrals)).optimize();
}
ExpressionTreeNode CustomCentroidBondForceImpl::replaceFunctions(const ExpressionTreeNode& node, map<string, int> groups,
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.getName() != "distance" && op.getName() != "angle" && op.getName() != "dihedral"))
{
// 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], groups, distances, angles, dihedrals));
return ExpressionTreeNode(op.clone(), children);
}
const Operation::Custom& custom = static_cast<const Operation::Custom&>(op);
// Identify the groups 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 = groups.find(node.getChildren()[i].getOperation().getName());
if (iter == groups.end())
throw OpenMMException("CustomCentroidBondForce: Unknown group '"+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));
}
void CustomCentroidBondForceImpl::updateParametersInContext(ContextImpl& context) {
kernel.getAs<CalcCustomCentroidBondForceKernel>().copyParametersToContext(context, owner);
}
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