Skip to content

GitLab

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

Created CustomTorsionForce, including reference implementation

parent 4af689c5
Hide whitespace changes
Inline Side-by-side
Showing with 1089 additions and 13 deletions
olla/include/openmm/kernels.h
View file @ bce05131
...@@ -40,6 +40,7 @@ ...@@ -40,6 +40,7 @@
#include "openmm/CustomExternalForce.h" #include "openmm/CustomExternalForce.h"
#include "openmm/CustomGBForce.h" #include "openmm/CustomGBForce.h"
#include "openmm/CustomNonbondedForce.h" #include "openmm/CustomNonbondedForce.h"
#include "openmm/CustomTorsionForce.h"
#include "openmm/GBSAOBCForce.h" #include "openmm/GBSAOBCForce.h"
#include "openmm/GBVIForce.h" #include "openmm/GBVIForce.h"
#include "openmm/HarmonicAngleForce.h" #include "openmm/HarmonicAngleForce.h"
...@@ -364,6 +365,38 @@ public: ...@@ -364,6 +365,38 @@ public:
virtual double executeEnergy(ContextImpl& context) = 0; virtual double executeEnergy(ContextImpl& context) = 0;
}; };
/**
* This kernel is invoked by CustomTorsionForce to calculate the forces acting on the system and the energy of the system.
*/
class CalcCustomTorsionForceKernel : public KernelImpl {
public:
static std::string Name() {
return "CalcCustomTorsionForce";
}
CalcCustomTorsionForceKernel(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 CustomTorsionForce this kernel will be used for
*/
virtual void initialize(const System& system, const CustomTorsionForce& force) = 0;
/**
* Execute the kernel to calculate the forces.
*
* @param context the context in which to execute this kernel
*/
virtual void executeForces(ContextImpl& context) = 0;
/**
* Execute the kernel to calculate the energy.
*
* @param context the context in which to execute this kernel
* @return the potential energy due to the CustomTorsionForce
*/
virtual double executeEnergy(ContextImpl& context) = 0;
};
/** /**
* This kernel is invoked by NonbondedForce to calculate the forces acting on the system and the energy of the system. * This kernel is invoked by NonbondedForce to calculate the forces acting on the system and the energy of the system.
*/ */
......
openmmapi/include/openmm/CustomTorsionForce.h 0 → 100644
View file @ bce05131
#ifndef OPENMM_CUSTOMTORSIONFORCE_H_
#define OPENMM_CUSTOMTORSIONFORCE_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) 2010 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 implements interactions between sets of four particles that depend on the torsion angle between them.
* Unlike PeriodicTorsionForce, the functional form of the interaction is completely customizable, and may
* involve arbitrary algebraic expressions. In addition to the angle formed by the particles, it may depend
* on arbitrary global and per-torsion parameters.
*
* To use this class, create a CustomTorsionForce object, passing an algebraic expression to the constructor
* that defines the interaction energy between each set of particles. The expression may depend on theta, the torsion angle
* formed by the particles, as well as on any parameters you choose. Then call addPerTorsionParameter() to define per-torsion
* parameters, and addGlobalParameter() to define global parameters. The values of per-torsion parameters are specified as
* part of the system definition, while values of global parameters may be modified during a simulation by calling Context::setParameter().
* Finally, call addTorsion() once for each torsion. After an torsion has been added, you can modify its parameters by calling setTorsionParameters().
*
* As an example, the following code creates a CustomTorsionForce that implements a harmonic potential:
*
* <tt>CustomTorsionForce* force = new CustomTorsionForce("0.5*k*(theta-theta0)^2");</tt>
*
* This force depends on two parameters: the spring constant k and equilibrium angle theta0. The following code defines these parameters:
*
* <tt><pre>
* force->addPerTorsionParameter("k");
* force->addPerTorsionParameter("theta0");
* </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, step. All trigonometric functions
* are defined in radians, and log is the natural logarithm. step(x) = 0 if x is less than 0, 1 otherwise.
*/
class OPENMM_EXPORT CustomTorsionForce : public Force {
public:
/**
* Create a CustomTorsionForce.
*
* @param energy an algebraic expression giving the interaction energy between three particles as a function
* of theta, the torsion angle between them
*/
CustomTorsionForce(const std::string& energy);
/**
* Get the number of torsions for which force field parameters have been defined.
*/
int getNumTorsions() const {
return torsions.size();
}
/**
* Get the number of per-torsion parameters that the interaction depends on.
*/
int getNumPerTorsionParameters() const {
return parameters.size();
}
/**
* Get the number of global parameters that the interaction depends on.
*/
int getNumGlobalParameters() const {
return globalParameters.size();
}
/**
* Get the algebraic expression that gives the interaction energy for each torsion
*/
const std::string& getEnergyFunction() const;
/**
* Set the algebraic expression that gives the interaction energy for each torsion
*/
void setEnergyFunction(const std::string& energy);
/**
* Add a new per-torsion parameter that the interaction may depend on.
*
* @param name the name of the parameter
* @return the index of the parameter that was added
*/
int addPerTorsionParameter(const std::string& name);
/**
* Get the name of a per-torsion parameter.
*
* @param index the index of the parameter for which to get the name
* @return the parameter name
*/
const std::string& getPerTorsionParameterName(int index) const;
/**
* Set the name of a per-torsion parameter.
*
* @param index the index of the parameter for which to set the name
* @param name the name of the parameter
*/
void setPerTorsionParameterName(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 torsion term to the force field.
*
* @param particle1 the index of the first particle connected by the torsion
* @param particle2 the index of the second particle connected by the torsion
* @param particle3 the index of the third particle connected by the torsion
* @param particle4 the index of the fourth particle connected by the torsion
* @param parameters the list of parameters for the new torsion
* @return the index of the torsion that was added
*/
int addTorsion(int particle1, int particle2, int particle3, int particle4, const std::vector<double>& parameters);
/**
* Get the force field parameters for a torsion term.
*
* @param index the index of the torsion for which to get parameters
* @param particle1 the index of the first particle connected by the torsion
* @param particle2 the index of the second particle connected by the torsion
* @param particle3 the index of the third particle connected by the torsion
* @param particle4 the index of the fourth particle connected by the torsion
* @param parameters the list of parameters for the torsion
*/
void getTorsionParameters(int index, int& particle1, int& particle2, int& particle3, int& particle4, std::vector<double>& parameters) const;
/**
* Set the force field parameters for a torsion term.
*
* @param index the index of the torsion for which to set parameters
* @param particle1 the index of the first particle connected by the torsion
* @param particle2 the index of the second particle connected by the torsion
* @param particle3 the index of the third particle connected by the torsion
* @param particle4 the index of the fourth particle connected by the torsion
* @param parameters the list of parameters for the torsion
*/
void setTorsionParameters(int index, int particle1, int particle2, int particle3, int particle4, const std::vector<double>& parameters);
protected:
ForceImpl* createImpl();
private:
class TorsionInfo;
class TorsionParameterInfo;
class GlobalParameterInfo;
std::string energyExpression;
std::vector<TorsionParameterInfo> parameters;
std::vector<GlobalParameterInfo> globalParameters;
std::vector<TorsionInfo> torsions;
};
/**
* This is an internal class used to record information about a torsion.
* @private
*/
class CustomTorsionForce::TorsionInfo {
public:
int particle1, particle2, particle3, particle4;
std::vector<double> parameters;
TorsionInfo() : particle1(-1), particle2(-1), particle3(-1), particle4(-1) {
}
TorsionInfo(int particle1, int particle2, int particle3, int particle4, const std::vector<double>& parameters) :
particle1(particle1), particle2(particle2), particle3(particle3), particle4(particle4), parameters(parameters) {
}
};
/**
* This is an internal class used to record information about a per-torsion parameter.
* @private
*/
class CustomTorsionForce::TorsionParameterInfo {
public:
std::string name;
TorsionParameterInfo() {
}
TorsionParameterInfo(const std::string& name) : name(name) {
}
};
/**
* This is an internal class used to record information about a global parameter.
* @private
*/
class CustomTorsionForce::GlobalParameterInfo {
public:
std::string name;
double defaultValue;
GlobalParameterInfo() {
}
GlobalParameterInfo(const std::string& name, double defaultValue) : name(name), defaultValue(defaultValue) {
}
};
} // namespace OpenMM
#endif /*OPENMM_CUSTOMTORSIONFORCE_H_*/
openmmapi/include/openmm/internal/CustomTorsionForceImpl.h 0 → 100644
View file @ bce05131
#ifndef OPENMM_CUSTOMTORSIONFORCEIMPL_H_
#define OPENMM_CUSTOMTORSIONFORCEIMPL_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) 2010 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/CustomTorsionForce.h"
#include "openmm/Kernel.h"
#include <utility>
#include <map>
#include <string>
namespace OpenMM {
/**
* This is the internal implementation of CustomTorsionForce.
*/
class CustomTorsionForceImpl : public ForceImpl {
public:
CustomTorsionForceImpl(CustomTorsionForce& owner);
~CustomTorsionForceImpl();
void initialize(ContextImpl& context);
CustomTorsionForce& getOwner() {
return owner;
}
void updateContextState(ContextImpl& context) {
// This force field doesn't update the state directly.
}
void calcForces(ContextImpl& context);
double calcEnergy(ContextImpl& context);
std::map<std::string, double> getDefaultParameters();
std::vector<std::string> getKernelNames();
private:
CustomTorsionForce& owner;
Kernel kernel;
};
} // namespace OpenMM
#endif /*OPENMM_CUSTOMTORSIONFORCEIMPL_H_*/
openmmapi/src/CustomTorsionForce.cpp 0 → 100644
View file @ bce05131
/* -------------------------------------------------------------------------- *
* 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) 2010 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/CustomTorsionForce.h"
#include "openmm/internal/CustomTorsionForceImpl.h"
#include <cmath>
#include <map>
#include <sstream>
#include <utility>
using namespace OpenMM;
using std::string;
using std::stringstream;
using std::vector;
CustomTorsionForce::CustomTorsionForce(const string& energy) : energyExpression(energy) {
}
const string& CustomTorsionForce::getEnergyFunction() const {
return energyExpression;
}
void CustomTorsionForce::setEnergyFunction(const std::string& energy) {
energyExpression = energy;
}
int CustomTorsionForce::addPerTorsionParameter(const string& name) {
parameters.push_back(TorsionParameterInfo(name));
return parameters.size()-1;
}
const string& CustomTorsionForce::getPerTorsionParameterName(int index) const {
return parameters[index].name;
}
void CustomTorsionForce::setPerTorsionParameterName(int index, const string& name) {
parameters[index].name = name;
}
int CustomTorsionForce::addGlobalParameter(const string& name, double defaultValue) {
globalParameters.push_back(GlobalParameterInfo(name, defaultValue));
return globalParameters.size()-1;
}
const string& CustomTorsionForce::getGlobalParameterName(int index) const {
return globalParameters[index].name;
}
void CustomTorsionForce::setGlobalParameterName(int index, const string& name) {
globalParameters[index].name = name;
}
double CustomTorsionForce::getGlobalParameterDefaultValue(int index) const {
return globalParameters[index].defaultValue;
}
void CustomTorsionForce::setGlobalParameterDefaultValue(int index, double defaultValue) {
globalParameters[index].defaultValue = defaultValue;
}
int CustomTorsionForce::addTorsion(int particle1, int particle2, int particle3, int particle4, const vector<double>& parameters) {
torsions.push_back(TorsionInfo(particle1, particle2, particle3, particle4, parameters));
return torsions.size()-1;
}
void CustomTorsionForce::getTorsionParameters(int index, int& particle1, int& particle2, int& particle3, int& particle4, std::vector<double>& parameters) const {
particle1 = torsions[index].particle1;
particle2 = torsions[index].particle2;
particle3 = torsions[index].particle3;
particle4 = torsions[index].particle4;
parameters = torsions[index].parameters;
}
void CustomTorsionForce::setTorsionParameters(int index, int particle1, int particle2, int particle3, int particle4, const vector<double>& parameters) {
torsions[index].parameters = parameters;
torsions[index].particle1 = particle1;
torsions[index].particle2 = particle2;
torsions[index].particle3 = particle3;
torsions[index].particle4 = particle4;
}
ForceImpl* CustomTorsionForce::createImpl() {
return new CustomTorsionForceImpl(*this);
}
openmmapi/src/CustomTorsionForceImpl.cpp 0 → 100644
View file @ bce05131
/* -------------------------------------------------------------------------- *
* 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) 2010 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/CustomTorsionForceImpl.h"
#include "openmm/kernels.h"
#include <sstream>
using namespace OpenMM;
using std::map;
using std::pair;
using std::vector;
using std::set;
using std::string;
using std::stringstream;
CustomTorsionForceImpl::CustomTorsionForceImpl(CustomTorsionForce& owner) : owner(owner) {
}
CustomTorsionForceImpl::~CustomTorsionForceImpl() {
}
void CustomTorsionForceImpl::initialize(ContextImpl& context) {
kernel = context.getPlatform().createKernel(CalcCustomTorsionForceKernel::Name(), context);
// Check for errors in the specification of bonds.
System& system = context.getSystem();
vector<double> parameters;
int numParameters = owner.getNumPerTorsionParameters();
for (int i = 0; i < owner.getNumTorsions(); i++) {
int particle1, particle2, particle3, particle4;
owner.getTorsionParameters(i, particle1, particle2, particle3, particle4, parameters);
if (particle1 < 0 || particle1 >= system.getNumParticles()) {
stringstream msg;
msg << "CustomTorsionForce: Illegal particle index for an torsion: ";
msg << particle1;
throw OpenMMException(msg.str());
}
if (particle2 < 0 || particle2 >= system.getNumParticles()) {
stringstream msg;
msg << "CustomTorsionForce: Illegal particle index for an torsion: ";
msg << particle2;
throw OpenMMException(msg.str());
}
if (particle3 < 0 || particle3 >= system.getNumParticles()) {
stringstream msg;
msg << "CustomTorsionForce: Illegal particle index for an torsion: ";
msg << particle3;
throw OpenMMException(msg.str());
}
if (particle4 < 0 || particle4 >= system.getNumParticles()) {
stringstream msg;
msg << "CustomTorsionForce: Illegal particle index for an torsion: ";
msg << particle4;
throw OpenMMException(msg.str());
}
if (parameters.size() != numParameters) {
stringstream msg;
msg << "CustomTorsionForce: Wrong number of parameters for torsion ";
msg << i;
throw OpenMMException(msg.str());
}
}
dynamic_cast<CalcCustomTorsionForceKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner);
}
void CustomTorsionForceImpl::calcForces(ContextImpl& context) {
dynamic_cast<CalcCustomTorsionForceKernel&>(kernel.getImpl()).executeForces(context);
}
double CustomTorsionForceImpl::calcEnergy(ContextImpl& context) {
return dynamic_cast<CalcCustomTorsionForceKernel&>(kernel.getImpl()).executeEnergy(context);
}
vector<string> CustomTorsionForceImpl::getKernelNames() {
vector<string> names;
names.push_back(CalcCustomTorsionForceKernel::Name());
return names;
}
map<string, double> CustomTorsionForceImpl::getDefaultParameters() {
map<string, double> parameters;
for (int i = 0; i < owner.getNumGlobalParameters(); i++)
parameters[owner.getGlobalParameterName(i)] = owner.getGlobalParameterDefaultValue(i);
return parameters;
}
platforms/cuda/tests/TestCudaCustomAngleForce.cpp
View file @ bce05131
...@@ -110,9 +110,8 @@ void testAngles() { ...@@ -110,9 +110,8 @@ void testAngles() {
energy2 = s.getPotentialEnergy(); energy2 = s.getPotentialEnergy();
forces2 = s.getForces(); forces2 = s.getForces();
} }
ASSERT_EQUAL_VEC(forces2[0], forces1[0], TOL); for (int i = 0; i < customSystem.getNumParticles(); i++)
ASSERT_EQUAL_VEC(forces2[1], forces1[1], TOL); ASSERT_EQUAL_VEC(forces2[i], forces1[i], TOL);
ASSERT_EQUAL_VEC(forces2[2], forces1[2], TOL);
ASSERT_EQUAL_TOL(energy2, energy1, TOL); ASSERT_EQUAL_TOL(energy2, energy1, TOL);
} }
} }
......
platforms/opencl/tests/TestOpenCLCustomAngleForce.cpp
View file @ bce05131
...@@ -101,9 +101,8 @@ void testAngles() { ...@@ -101,9 +101,8 @@ void testAngles() {
State s1 = c1.getState(State::Forces | State::Energy); State s1 = c1.getState(State::Forces | State::Energy);
State s2 = c2.getState(State::Forces | State::Energy); State s2 = c2.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = s1.getForces(); const vector<Vec3>& forces = s1.getForces();
ASSERT_EQUAL_VEC(s1.getForces()[0], s2.getForces()[0], TOL); for (int i = 0; i < customSystem.getNumParticles(); i++)
ASSERT_EQUAL_VEC(s1.getForces()[1], s2.getForces()[1], TOL); ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], TOL);
ASSERT_EQUAL_VEC(s1.getForces()[2], s2.getForces()[2], TOL);
ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), TOL); ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), TOL);
} }
} }
......
platforms/reference/src/ReferenceKernelFactory.cpp
View file @ bce05131
...@@ -60,6 +60,8 @@ KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Pla ...@@ -60,6 +60,8 @@ KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Pla
return new ReferenceCalcPeriodicTorsionForceKernel(name, platform); return new ReferenceCalcPeriodicTorsionForceKernel(name, platform);
if (name == CalcRBTorsionForceKernel::Name()) if (name == CalcRBTorsionForceKernel::Name())
return new ReferenceCalcRBTorsionForceKernel(name, platform); return new ReferenceCalcRBTorsionForceKernel(name, platform);
if (name == CalcCustomTorsionForceKernel::Name())
return new ReferenceCalcCustomTorsionForceKernel(name, platform);
if (name == CalcGBSAOBCForceKernel::Name()) if (name == CalcGBSAOBCForceKernel::Name())
return new ReferenceCalcGBSAOBCForceKernel(name, platform); return new ReferenceCalcGBSAOBCForceKernel(name, platform);
if (name == CalcGBVIForceKernel::Name()) if (name == CalcGBVIForceKernel::Name())
......
platforms/reference/src/ReferenceKernels.cpp
View file @ bce05131
...@@ -42,6 +42,7 @@ ...@@ -42,6 +42,7 @@
#include "SimTKReference/ReferenceCustomExternalIxn.h" #include "SimTKReference/ReferenceCustomExternalIxn.h"
#include "SimTKReference/ReferenceCustomGBIxn.h" #include "SimTKReference/ReferenceCustomGBIxn.h"
#include "SimTKReference/ReferenceCustomNonbondedIxn.h" #include "SimTKReference/ReferenceCustomNonbondedIxn.h"
#include "SimTKReference/ReferenceCustomTorsionIxn.h"
#include "SimTKReference/ReferenceHarmonicBondIxn.h" #include "SimTKReference/ReferenceHarmonicBondIxn.h"
#include "SimTKReference/ReferenceLJCoulomb14.h" #include "SimTKReference/ReferenceLJCoulomb14.h"
#include "SimTKReference/ReferenceLJCoulombIxn.h" #include "SimTKReference/ReferenceLJCoulombIxn.h"
...@@ -416,8 +417,8 @@ void ReferenceCalcCustomAngleForceKernel::executeForces(ContextImpl& context) { ...@@ -416,8 +417,8 @@ void ReferenceCalcCustomAngleForceKernel::executeForces(ContextImpl& context) {
for (int i = 0; i < (int) globalParameterNames.size(); i++) for (int i = 0; i < (int) globalParameterNames.size(); i++)
globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]); globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
ReferenceBondForce refBondForce; ReferenceBondForce refBondForce;
ReferenceCustomAngleIxn harmonicAngle(energyExpression, forceExpression, parameterNames, globalParameters); ReferenceCustomAngleIxn customAngle(energyExpression, forceExpression, parameterNames, globalParameters);
refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, 0, 0, 0, harmonicAngle); refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, 0, 0, 0, customAngle);
} }
double ReferenceCalcCustomAngleForceKernel::executeEnergy(ContextImpl& context) { double ReferenceCalcCustomAngleForceKernel::executeEnergy(ContextImpl& context) {
...@@ -429,10 +430,10 @@ double ReferenceCalcCustomAngleForceKernel::executeEnergy(ContextImpl& context) ...@@ -429,10 +430,10 @@ double ReferenceCalcCustomAngleForceKernel::executeEnergy(ContextImpl& context)
for (int i = 0; i < (int) globalParameterNames.size(); i++) for (int i = 0; i < (int) globalParameterNames.size(); i++)
globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]); globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
ReferenceBondForce refBondForce; ReferenceBondForce refBondForce;
ReferenceCustomAngleIxn harmonicAngle(energyExpression, forceExpression, parameterNames, globalParameters); ReferenceCustomAngleIxn customAngle(energyExpression, forceExpression, parameterNames, globalParameters);
for (int i = 0; i < numAngles; ++i) for (int i = 0; i < numAngles; ++i)
energyArray[i] = 0; energyArray[i] = 0;
refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, energyArray, 0, &energy, harmonicAngle); refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, energyArray, 0, &energy, customAngle);
disposeRealArray(forceData, context.getSystem().getNumParticles()); disposeRealArray(forceData, context.getSystem().getNumParticles());
delete[] energyArray; delete[] energyArray;
return energy; return energy;
...@@ -533,6 +534,71 @@ double ReferenceCalcRBTorsionForceKernel::executeEnergy(ContextImpl& context) { ...@@ -533,6 +534,71 @@ double ReferenceCalcRBTorsionForceKernel::executeEnergy(ContextImpl& context) {
return energy; return energy;
} }
ReferenceCalcCustomTorsionForceKernel::~ReferenceCalcCustomTorsionForceKernel() {
disposeIntArray(torsionIndexArray, numTorsions);
disposeRealArray(torsionParamArray, numTorsions);
}
void ReferenceCalcCustomTorsionForceKernel::initialize(const System& system, const CustomTorsionForce& force) {
numTorsions = force.getNumTorsions();
int numParameters = force.getNumPerTorsionParameters();
// Build the arrays.
torsionIndexArray = allocateIntArray(numTorsions, 4);
torsionParamArray = allocateRealArray(numTorsions, numParameters);
vector<double> params;
for (int i = 0; i < force.getNumTorsions(); ++i) {
int particle1, particle2, particle3, particle4;
force.getTorsionParameters(i, particle1, particle2, particle3, particle4, params);
torsionIndexArray[i][0] = particle1;
torsionIndexArray[i][1] = particle2;
torsionIndexArray[i][2] = particle3;
torsionIndexArray[i][3] = particle4;
for (int j = 0; j < numParameters; j++)
torsionParamArray[i][j] = (RealOpenMM) params[j];
}
// Parse the expression used to calculate the force.
Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction()).optimize();
energyExpression = expression.createProgram();
forceExpression = expression.differentiate("theta").optimize().createProgram();
for (int i = 0; i < numParameters; i++)
parameterNames.push_back(force.getPerTorsionParameterName(i));
for (int i = 0; i < force.getNumGlobalParameters(); i++)
globalParameterNames.push_back(force.getGlobalParameterName(i));
}
void ReferenceCalcCustomTorsionForceKernel::executeForces(ContextImpl& context) {
RealOpenMM** posData = extractPositions(context);
RealOpenMM** forceData = extractForces(context);
map<string, double> globalParameters;
for (int i = 0; i < (int) globalParameterNames.size(); i++)
globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
ReferenceBondForce refBondForce;
ReferenceCustomTorsionIxn customTorsion(energyExpression, forceExpression, parameterNames, globalParameters);
refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, 0, 0, 0, customTorsion);
}
double ReferenceCalcCustomTorsionForceKernel::executeEnergy(ContextImpl& context) {
RealOpenMM** posData = extractPositions(context);
RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
RealOpenMM* energyArray = new RealOpenMM[numTorsions];
RealOpenMM energy = 0;
map<string, double> globalParameters;
for (int i = 0; i < (int) globalParameterNames.size(); i++)
globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
ReferenceBondForce refBondForce;
ReferenceCustomTorsionIxn customTorsion(energyExpression, forceExpression, parameterNames, globalParameters);
for (int i = 0; i < numTorsions; ++i)
energyArray[i] = 0;
refBondForce.calculateForce(numTorsions, torsionIndexArray, posData, torsionParamArray, forceData, energyArray, 0, &energy, customTorsion);
disposeRealArray(forceData, context.getSystem().getNumParticles());
delete[] energyArray;
return energy;
}
ReferenceCalcNonbondedForceKernel::~ReferenceCalcNonbondedForceKernel() { ReferenceCalcNonbondedForceKernel::~ReferenceCalcNonbondedForceKernel() {
disposeRealArray(particleParamArray, numParticles); disposeRealArray(particleParamArray, numParticles);
disposeIntArray(exclusionArray, numParticles); disposeIntArray(exclusionArray, numParticles);
......
platforms/reference/src/ReferenceKernels.h
View file @ bce05131
...@@ -366,6 +366,42 @@ private: ...@@ -366,6 +366,42 @@ private:
RealOpenMM **torsionParamArray; RealOpenMM **torsionParamArray;
}; };
/**
* This kernel is invoked by CustomTorsionForce to calculate the forces acting on the system and the energy of the system.
*/
class ReferenceCalcCustomTorsionForceKernel : public CalcCustomTorsionForceKernel {
public:
ReferenceCalcCustomTorsionForceKernel(std::string name, const Platform& platform) : CalcCustomTorsionForceKernel(name, platform) {
}
~ReferenceCalcCustomTorsionForceKernel();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param force the CustomTorsionForce this kernel will be used for
*/
void initialize(const System& system, const CustomTorsionForce& force);
/**
* Execute the kernel to calculate the forces.
*
* @param context the context in which to execute this kernel
*/
void executeForces(ContextImpl& context);
/**
* Execute the kernel to calculate the energy.
*
* @param context the context in which to execute this kernel
* @return the potential energy due to the CustomTorsionForce
*/
double executeEnergy(ContextImpl& context);
private:
int numTorsions;
int **torsionIndexArray;
RealOpenMM **torsionParamArray;
Lepton::ExpressionProgram energyExpression, forceExpression;
std::vector<std::string> parameterNames, globalParameterNames;
};
/** /**
* This kernel is invoked by NonbondedForce to calculate the forces acting on the system. * This kernel is invoked by NonbondedForce to calculate the forces acting on the system.
*/ */
......
platforms/reference/src/ReferencePlatform.cpp
View file @ bce05131
...@@ -48,6 +48,7 @@ ReferencePlatform::ReferencePlatform() { ...@@ -48,6 +48,7 @@ ReferencePlatform::ReferencePlatform() {
registerKernelFactory(CalcCustomAngleForceKernel::Name(), factory); registerKernelFactory(CalcCustomAngleForceKernel::Name(), factory);
registerKernelFactory(CalcPeriodicTorsionForceKernel::Name(), factory); registerKernelFactory(CalcPeriodicTorsionForceKernel::Name(), factory);
registerKernelFactory(CalcRBTorsionForceKernel::Name(), factory); registerKernelFactory(CalcRBTorsionForceKernel::Name(), factory);
registerKernelFactory(CalcCustomTorsionForceKernel::Name(), factory);
registerKernelFactory(CalcNonbondedForceKernel::Name(), factory); registerKernelFactory(CalcNonbondedForceKernel::Name(), factory);
registerKernelFactory(CalcCustomNonbondedForceKernel::Name(), factory); registerKernelFactory(CalcCustomNonbondedForceKernel::Name(), factory);
registerKernelFactory(CalcGBSAOBCForceKernel::Name(), factory); registerKernelFactory(CalcGBSAOBCForceKernel::Name(), factory);
......
platforms/reference/src/SimTKReference/ReferenceCustomTorsionIxn.cpp 0 → 100644
View file @ bce05131
/* Portions copyright (c) 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 "../SimTKUtilities/SimTKOpenMMCommon.h"
#include "../SimTKUtilities/SimTKOpenMMLog.h"
#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
#include "ReferenceCustomTorsionIxn.h"
#include "ReferenceForce.h"
using namespace std;
/**---------------------------------------------------------------------------------------
ReferenceCustomTorsionIxn constructor
--------------------------------------------------------------------------------------- */
ReferenceCustomTorsionIxn::ReferenceCustomTorsionIxn(const Lepton::ExpressionProgram& energyExpression,
const Lepton::ExpressionProgram& forceExpression, const vector<string>& parameterNames, map<string, double> globalParameters) :
energyExpression(energyExpression), forceExpression(forceExpression), paramNames(parameterNames), globalParameters(globalParameters) {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceCustomTorsionIxn::ReferenceCustomTorsionIxn";
// ---------------------------------------------------------------------------------------
}
/**---------------------------------------------------------------------------------------
ReferenceCustomTorsionIxn destructor
--------------------------------------------------------------------------------------- */
ReferenceCustomTorsionIxn::~ReferenceCustomTorsionIxn( ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceCustomTorsionIxn::~ReferenceCustomTorsionIxn";
// ---------------------------------------------------------------------------------------
}
/**---------------------------------------------------------------------------------------
Calculate Custom Torsion Ixn
@param atomIndices atom indices of atom participating in bond
@param atomCoordinates atom coordinates
@param parameters parameters values
@param forces force array (forces added to input values)
@param energiesByBond energies by bond: energiesByBond[bondIndex]
@param energiesByAtom energies by atom: energiesByAtom[atomIndex]
@return ReferenceForce::DefaultReturn;
--------------------------------------------------------------------------------------- */
int ReferenceCustomTorsionIxn::calculateBondIxn( int* atomIndices,
RealOpenMM** atomCoordinates,
RealOpenMM* parameters,
RealOpenMM** forces,
RealOpenMM* energiesByBond,
RealOpenMM* energiesByAtom ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceCustomTorsionIxn::calculateTorsionIxn";
// ---------------------------------------------------------------------------------------
static const std::string methodName = "\nReferenceCustomTorsionIxn::calculateTorsionIxn";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
RealOpenMM deltaR[3][ReferenceForce::LastDeltaRIndex];
map<string, double> variables = globalParameters;
for (int i = 0; i < (int) paramNames.size(); ++i)
variables[paramNames[i]] = parameters[i];
// ---------------------------------------------------------------------------------------
// get deltaR, R2, and R between three pairs of atoms: [j,i], [j,k], [l,k]
int atomAIndex = atomIndices[0];
int atomBIndex = atomIndices[1];
int atomCIndex = atomIndices[2];
int atomDIndex = atomIndices[3];
ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], deltaR[0]);
ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomCIndex], deltaR[1]);
ReferenceForce::getDeltaR(atomCoordinates[atomDIndex], atomCoordinates[atomCIndex], deltaR[2]);
// Visual Studio complains if crossProduct declared as 'crossProduct[2][3]'
RealOpenMM crossProductMemory[6];
RealOpenMM* crossProduct[2];
crossProduct[0] = crossProductMemory;
crossProduct[1] = crossProductMemory + 3;
// get dihedral angle
RealOpenMM dotDihedral;
RealOpenMM signOfAngle;
variables["theta"] = getDihedralAngleBetweenThreeVectors(deltaR[0], deltaR[1], deltaR[2],
crossProduct, &dotDihedral, deltaR[0],
&signOfAngle, 1);
// evaluate delta angle, dE/d(angle)
RealOpenMM energy = (RealOpenMM) energyExpression.evaluate(variables);
RealOpenMM dEdAngle = (RealOpenMM) forceExpression.evaluate(variables);
// compute force
RealOpenMM internalF[4][3];
RealOpenMM forceFactors[4];
RealOpenMM normCross1 = DOT3( crossProduct[0], crossProduct[0] );
RealOpenMM normBC = deltaR[1][ReferenceForce::RIndex];
forceFactors[0] = (-dEdAngle*normBC)/normCross1;
RealOpenMM normCross2 = DOT3( crossProduct[1], crossProduct[1] );
forceFactors[3] = (dEdAngle*normBC)/normCross2;
forceFactors[1] = DOT3( deltaR[0], deltaR[1] );
forceFactors[1] /= deltaR[1][ReferenceForce::R2Index];
forceFactors[2] = DOT3( deltaR[2], deltaR[1] );
forceFactors[2] /= deltaR[1][ReferenceForce::R2Index];
for( int ii = 0; ii < 3; ii++ ){
internalF[0][ii] = forceFactors[0]*crossProduct[0][ii];
internalF[3][ii] = forceFactors[3]*crossProduct[1][ii];
RealOpenMM s = forceFactors[1]*internalF[0][ii] - forceFactors[2]*internalF[3][ii];
internalF[1][ii] = internalF[0][ii] - s;
internalF[2][ii] = internalF[3][ii] + s;
}
// accumulate forces
for( int ii = 0; ii < 3; ii++ ){
forces[atomAIndex][ii] += internalF[0][ii];
forces[atomBIndex][ii] -= internalF[1][ii];
forces[atomCIndex][ii] -= internalF[2][ii];
forces[atomDIndex][ii] += internalF[3][ii];
}
// accumulate energies
updateEnergy( energy, energiesByBond, 4, atomIndices, energiesByAtom );
return ReferenceForce::DefaultReturn;
}
platforms/reference/src/SimTKReference/ReferenceCustomTorsionIxn.h 0 → 100644
View file @ bce05131
/* Portions copyright (c) 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.
*/
#ifndef __ReferenceCustomTorsionIxn_H__
#define __ReferenceCustomTorsionIxn_H__
#include "ReferenceBondIxn.h"
#include "lepton/ExpressionProgram.h"
// ---------------------------------------------------------------------------------------
class ReferenceCustomTorsionIxn : public ReferenceBondIxn {
private:
Lepton::ExpressionProgram energyExpression;
Lepton::ExpressionProgram forceExpression;
std::vector<std::string> paramNames;
std::map<std::string, double> globalParameters;
public:
/**---------------------------------------------------------------------------------------
Constructor
--------------------------------------------------------------------------------------- */
ReferenceCustomTorsionIxn(const Lepton::ExpressionProgram& energyExpression, const Lepton::ExpressionProgram& forceExpression,
const std::vector<std::string>& parameterNames, std::map<std::string, double> globalParameters);
/**---------------------------------------------------------------------------------------
Destructor
--------------------------------------------------------------------------------------- */
~ReferenceCustomTorsionIxn( );
/**---------------------------------------------------------------------------------------
Calculate Custom Torsion Ixn
@param atomIndices two bond indices
@param atomCoordinates atom coordinates
@param parameters parameter values
@param forces force array (forces added)
@param energiesByBond energies by bond: energiesByBond[bondIndex]
@param energiesByAtom energies by atom: energiesByAtom[atomIndex]
--------------------------------------------------------------------------------------- */
int calculateBondIxn( int* atomIndices, RealOpenMM** atomCoordinates,
RealOpenMM* parameters, RealOpenMM** forces,
RealOpenMM* energiesByBond, RealOpenMM* energiesByAtom ) const;
};
// ---------------------------------------------------------------------------------------
#endif // _ReferenceCustomTorsionIxn___
platforms/reference/tests/TestReferenceCustomAngleForce.cpp
View file @ bce05131
...@@ -101,9 +101,8 @@ void testAngles() { ...@@ -101,9 +101,8 @@ void testAngles() {
State s1 = c1.getState(State::Forces | State::Energy); State s1 = c1.getState(State::Forces | State::Energy);
State s2 = c2.getState(State::Forces | State::Energy); State s2 = c2.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = s1.getForces(); const vector<Vec3>& forces = s1.getForces();
ASSERT_EQUAL_VEC(s1.getForces()[0], s2.getForces()[0], TOL); for (int i = 0; i < customSystem.getNumParticles(); i++)
ASSERT_EQUAL_VEC(s1.getForces()[1], s2.getForces()[1], TOL); ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], TOL);
ASSERT_EQUAL_VEC(s1.getForces()[2], s2.getForces()[2], TOL);
ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), TOL); ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), TOL);
} }
} }
......
platforms/reference/tests/TestReferenceCustomTorsionForce.cpp 0 → 100644
View file @ bce05131
/* -------------------------------------------------------------------------- *
* 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-2010 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 CustomTorsionForce.
*/
#include "../../../tests/AssertionUtilities.h"
#include "openmm/Context.h"
#include "ReferencePlatform.h"
#include "openmm/CustomTorsionForce.h"
#include "openmm/PeriodicTorsionForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "../src/sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testTorsions() {
ReferencePlatform platform;
// Create a system using a CustomTorsionForce.
System customSystem;
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
customSystem.addParticle(1.0);
CustomTorsionForce* custom = new CustomTorsionForce("k*(1+cos(n*theta-theta0))");
custom->addPerTorsionParameter("theta0");
custom->addPerTorsionParameter("n");
custom->addGlobalParameter("k", 0.5);
vector<double> parameters(2);
parameters[0] = 1.5;
parameters[1] = 1;
custom->addTorsion(0, 1, 2, 3, parameters);
parameters[0] = 2.0;
parameters[1] = 2;
custom->addTorsion(1, 2, 3, 4, parameters);
customSystem.addForce(custom);
// Create an identical system using a PeriodicTorsionForce.
System harmonicSystem;
harmonicSystem.addParticle(1.0);
harmonicSystem.addParticle(1.0);
harmonicSystem.addParticle(1.0);
harmonicSystem.addParticle(1.0);
harmonicSystem.addParticle(1.0);
VerletIntegrator integrator(0.01);
PeriodicTorsionForce* periodic = new PeriodicTorsionForce();
periodic->addTorsion(0, 1, 2, 3, 1, 1.5, 0.5);
periodic->addTorsion(1, 2, 3, 4, 2, 2.0, 0.5);
harmonicSystem.addForce(periodic);
// Set the atoms in various positions, and verify that both systems give identical forces and energy.
init_gen_rand(0);
vector<Vec3> positions(5);
VerletIntegrator integrator1(0.01);
VerletIntegrator integrator2(0.01);
for (int i = 0; i < 10; i++) {
Context c1(customSystem, integrator1, platform);
Context c2(harmonicSystem, integrator2, platform);
for (int j = 0; j < (int) positions.size(); j++)
positions[j] = Vec3(5.0*genrand_real2(), 5.0*genrand_real2(), 5.0*genrand_real2());
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);
}
}
int main() {
try {
testTorsions();
}
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