Skip to content

GitLab

Commit 70015b75 authored by Peter Eastman's avatar Peter Eastman
Browse files

Redesigned CustomNonbondedForce to eliminate combining rules and allow the... 

Redesigned CustomNonbondedForce to eliminate combining rules and allow the expression to directly use parameters for individual particles.
parent 37b2436e
Show whitespace changes
Inline Side-by-side
Showing with 215 additions and 698 deletions
openmmapi/include/OpenMM.h
View file @ 70015b75
......@@ -35,6 +35,8 @@
#include "openmm/AndersenThermostat.h"
#include "openmm/BrownianIntegrator.h"
#include "openmm/CMMotionRemover.h"
#include "openmm/CustomBondForce.h"
#include "openmm/CustomExternalForce.h"
#include "openmm/CustomNonbondedForce.h"
#include "openmm/Force.h"
#include "openmm/GBSAOBCForce.h"
......
openmmapi/include/openmm/CustomNonbondedForce.h
View file @ 70015b75
......@@ -50,40 +50,38 @@ namespace OpenMM {
*
* To use this class, create a CustomNonbondedForce object, passing an algebraic expression to the constructor
* that defines the interaction energy between each pair of particles. The expression may depend on r, the distance
* between the particles, as well as on any parameters you choose. Then call addParameter() to define per-particle
* parameters, and addGlobalParameter() to define global parameters. When defining a per-particle parameter, you
* specify an arbitrary algebraic expression which serves as the combining rule for calculating the parameter value
* based on the values for the two particles involved. The values of global parameters may be modified during a
* simulation by calling Context::setParameter().
* between the particles, as well as on any parameters you choose. Then call addPerParticleParameter() to define per-particle
* parameters, and addGlobalParameter() to define global parameters. The values of per-particle 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 addParticle() once for each particle in the System to set the values of its per-particle parameters.
* The number of particles for which you set parameters must be exactly equal to the number of particles in the
* System, or else an exception will be thrown when you try to create a Context. After a particle has been added,
* you can modify its parameters by calling setParticleParameters().
*
* CustomNonbondedForce also lets you specify "exceptions", particular pairs of particles whose interactions should be
* computed based on different parameters than those defined for the individual particles. This can be used to
* completely exclude certain interactions from the force calculation, or to alter how they interact with each other.
* CustomNonbondedForce also lets you specify "exclusions", particular pairs of particles whose interactions should be
* omitted from force and energy calculations. This is most often used for particles that are bonded to each other.
*
* As an example, the following code creates a CustomNonbondedForce that implements a 12-6 Lennard-Jones potential:
*
* <tt>CustomNonbondedForce* force = new CustomNonbondedForce("4*epsilon*((sigma/r)^12-(sigma/r)^6)");</tt>
* <tt>CustomNonbondedForce* force = new CustomNonbondedForce("4*epsilon*((sigma/r)^12-(sigma/r)^6); sigma=0.5*(sigma1*sigma2); epsilon=sqrt(epsilon1*epsilon2)");</tt>
*
* This force depends on two parameters: sigma and epsilon. The following code defines these parameters, and
* specifies combining rules for them which correspond to the standard Lorentz-Bertelot combining rules:
*
* <tt><pre>
* force->addParameter("sigma", "0.5*(sigma1*sigma2)");
* force->addParameter("epsilon", "sqrt(epsilon1*epsilon2)");
* force->addPerParticleParameter("sigma");
* force->addPerParticleParameter("epsilon");
* </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. All trigonometric functions
* are defined in radians, and log is the natural logarithm.
* are defined in radians, and log is the natural logarithm. The names of per-particle parameters have the suffix "1" or "2"
* appended to them to indicate the values for the two interacting particles. As seen in the above example, the expression may
* also involve intermediate quantities that are defined following the main expression, using ";" as a separator.
*
* In addition, you can call addFunction() to define a new function based on tabulated values. You specify a vector of
* values, and an interpolating or approximating spline is created from them. That function can then appear in expressions
* that define energy or combining rules.
* values, and an interpolating or approximating spline is created from them. That function can then appear in the expression.
*/
class OPENMM_EXPORT CustomNonbondedForce : public Force {
......@@ -111,7 +109,7 @@ public:
* Create a CustomNonbondedForce.
*
* @param energy an algebraic expression giving the interaction energy between two particles as a function
* of r, the distance between them
* of r, the distance between them, as well as any global and per-particle parameters
*/
CustomNonbondedForce(const std::string& energy);
/**
......@@ -121,15 +119,15 @@ public:
return particles.size();
}
/**
* Get the number of special interactions that should be calculated differently from other interactions.
* Get the number of particle pairs whose interactions should be excluded.
*/
int getNumExceptions() const {
return exceptions.size();
int getNumExclusions() const {
return exclusions.size();
}
/**
* Get the number of per-particle parameters that the interaction depends on.
*/
int getNumParameters() const {
int getNumPerParticleParameters() const {
return parameters.size();
}
/**
......@@ -174,38 +172,23 @@ public:
* Add a new per-particle parameter that the interaction may depend on.
*
* @param name the name of the parameter
* @param combiningRule an algebraic expression giving the combining rule for this parameter
* @return the index of the parameter that was added
*/
int addParameter(const std::string& name, const std::string& combiningRule);
int addPerParticleParameter(const std::string& name);
/**
* Get the name of a per-particle parameter.
*
* @param index the index of the parameter for which to get the name
* @return the parameter name
*/
const std::string& getParameterName(int index) const;
const std::string& getPerParticleParameterName(int index) const;
/**
* Set the name of a per-particle parameter.
*
* @param index the index of the parameter for which to set the name
* @param name the name of the parameter
*/
void setParameterName(int index, const std::string& name);
/**
* Get the combining rule for a per-particle parameter.
*
* @param index the index of the parameter for which to get the combining rule
* @return an algebraic expression giving the combining rule for the parameter
*/
const std::string& getParameterCombiningRule(int index) const;
/**
* Set the combining rule for a per-particle parameter.
*
* @param index the index of the parameter for which to set the combining rule
* @param combiningRule an algebraic expression giving the combining rule for the parameter
*/
void setParameterCombiningRule(int index, const std::string& combiningRule);
void setPerParticleParameterName(int index, const std::string& name);
/**
* Add a new global parameter that the interaction may depend on.
*
......@@ -265,39 +248,31 @@ public:
*/
void setParticleParameters(int index, const std::vector<double>& parameters);
/**
* Add an interaction to the list of exceptions that should be calculated differently from other interactions.
* Add a particle pair to the list of interactions that should be excluded.
*
* @param particle1 the index of the first particle involved in the interaction
* @param particle2 the index of the second particle involved in the interaction
* @param parameters the list of parameters for the new interaction. If this is an empty (zero length) vector, it
* will cause the interaction to be completely omitted from force and energy calculations.
* @param replace determines the behavior if there is already an exception for the same two particles. If true, the existing one is replaced. If false,
* an exception is thrown.
* @return the index of the exception that was added
* @param particle1 the index of the first particle in the pair
* @param particle2 the index of the second particle in the pair
* @return the index of the exclusion that was added
*/
int addException(int particle1, int particle2, const std::vector<double>& parameters, bool replace = false);
int addExclusion(int particle1, int particle2);
/**
* Get the force field parameters for an interaction that should be calculated differently from others.
* Get the particles in a pair whose interaction should be excluded.
*
* @param index the index of the interaction for which to get parameters
* @param particle1 the index of the first particle involved in the interaction
* @param particle2 the index of the second particle involved in the interaction
* @param parameters the list of parameters for the interaction. If this is an empty (zero length) vector, it means
* the interaction will be completely omitted from force and energy calculations.
* @param index the index of the exclusion for which to get particle indices
* @param particle1 the index of the first particle in the pair
* @param particle2 the index of the second particle in the pair
*/
void getExceptionParameters(int index, int& particle1, int& particle2, std::vector<double>& parameters) const;
void getExclusionParticles(int index, int& particle1, int& particle2) const;
/**
* Set the force field parameters for an interaction that should be calculated differently from others.
* Set the particles in a pair whose interaction should be excluded.
*
* @param index the index of the interaction for which to get parameters
* @param particle1 the index of the first particle involved in the interaction
* @param particle2 the index of the second particle involved in the interaction
* @param parameters the list of parameters for the interaction. If this is an empty (zero length) vector, it
* will cause the interaction to be completely omitted from force and energy calculations.
* @param index the index of the exclusion for which to set particle indices
* @param particle1 the index of the first particle in the pair
* @param particle2 the index of the second particle in the pair
*/
void setExceptionParameters(int index, int particle1, int particle2, const std::vector<double>& parameters);
void setExclusionParticles(int index, int particle1, int particle2);
/**
* Add a tabulated function that may appear in algebraic expressions.
* 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.
......@@ -310,7 +285,7 @@ public:
*/
int addFunction(const std::string& name, const std::vector<double>& values, double min, double max, bool interpolating);
/**
* Get the parameters for a tabulated function that may appear in algebraic expressions.
* 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
......@@ -339,19 +314,19 @@ protected:
ForceImpl* createImpl();
private:
class ParticleInfo;
class ParameterInfo;
class PerParticleParameterInfo;
class GlobalParameterInfo;
class ExceptionInfo;
class ExclusionInfo;
class FunctionInfo;
NonbondedMethod nonbondedMethod;
double cutoffDistance;
std::string energyExpression;
std::vector<ParameterInfo> parameters;
std::vector<PerParticleParameterInfo> parameters;
std::vector<GlobalParameterInfo> globalParameters;
std::vector<ParticleInfo> particles;
std::vector<ExceptionInfo> exceptions;
std::vector<ExclusionInfo> exclusions;
std::vector<FunctionInfo> functions;
std::map<std::pair<int, int>, int> exceptionMap;
std::map<std::pair<int, int>, int> exclusionMap;
};
class CustomNonbondedForce::ParticleInfo {
......@@ -363,12 +338,12 @@ public:
}
};
class CustomNonbondedForce::ParameterInfo {
class CustomNonbondedForce::PerParticleParameterInfo {
public:
std::string name, combiningRule;
ParameterInfo() {
std::string name;
PerParticleParameterInfo() {
}
ParameterInfo(const std::string& name, const std::string& combiningRule) : name(name), combiningRule(combiningRule) {
PerParticleParameterInfo(const std::string& name) : name(name) {
}
};
......@@ -382,15 +357,14 @@ public:
}
};
class CustomNonbondedForce::ExceptionInfo {
class CustomNonbondedForce::ExclusionInfo {
public:
int particle1, particle2;
std::vector<double> parameters;
ExceptionInfo() {
ExclusionInfo() {
particle1 = particle2 = -1;
}
ExceptionInfo(int particle1, int particle2, const std::vector<double>& parameters) :
particle1(particle1), particle2(particle2), parameters(parameters) {
ExclusionInfo(int particle1, int particle2) :
particle1(particle1), particle2(particle2) {
}
};
......
openmmapi/src/CustomNonbondedForce.cpp
View file @ 70015b75
......@@ -73,27 +73,19 @@ void CustomNonbondedForce::setCutoffDistance(double distance) {
cutoffDistance = distance;
}
int CustomNonbondedForce::addParameter(const string& name, const string& combiningRule) {
parameters.push_back(ParameterInfo(name, combiningRule));
int CustomNonbondedForce::addPerParticleParameter(const string& name) {
parameters.push_back(PerParticleParameterInfo(name));
return parameters.size()-1;
}
const string& CustomNonbondedForce::getParameterName(int index) const {
const string& CustomNonbondedForce::getPerParticleParameterName(int index) const {
return parameters[index].name;
}
void CustomNonbondedForce::setParameterName(int index, const string& name) {
void CustomNonbondedForce::setPerParticleParameterName(int index, const string& name) {
parameters[index].name = name;
}
const string& CustomNonbondedForce::getParameterCombiningRule(int index) const {
return parameters[index].combiningRule;
}
void CustomNonbondedForce::setParameterCombiningRule(int index, const string& combiningRule) {
parameters[index].combiningRule = combiningRule;
}
int CustomNonbondedForce::addGlobalParameter(const string& name, double defaultValue) {
globalParameters.push_back(GlobalParameterInfo(name, defaultValue));
return globalParameters.size()-1;
......@@ -128,41 +120,18 @@ void CustomNonbondedForce::setParticleParameters(int index, const vector<double>
particles[index].parameters = parameters;
}
int CustomNonbondedForce::addException(int particle1, int particle2, const vector<double>& parameters, bool replace) {
map<pair<int, int>, int>::iterator iter = exceptionMap.find(pair<int, int>(particle1, particle2));
int newIndex;
if (iter == exceptionMap.end())
iter = exceptionMap.find(pair<int, int>(particle2, particle1));
if (iter != exceptionMap.end()) {
if (!replace) {
stringstream msg;
msg << "CustomNonbondedForce: There is already an exception for particles ";
msg << particle1;
msg << " and ";
msg << particle2;
throw OpenMMException(msg.str());
}
exceptions[iter->second] = ExceptionInfo(particle1, particle2, parameters);
newIndex = iter->second;
exceptionMap.erase(iter->first);
}
else {
exceptions.push_back(ExceptionInfo(particle1, particle2, parameters));
newIndex = exceptions.size()-1;
}
exceptionMap[pair<int, int>(particle1, particle2)] = newIndex;
return newIndex;
}
void CustomNonbondedForce::getExceptionParameters(int index, int& particle1, int& particle2, vector<double>& parameters) const {
particle1 = exceptions[index].particle1;
particle2 = exceptions[index].particle2;
parameters = exceptions[index].parameters;
}
void CustomNonbondedForce::setExceptionParameters(int index, int particle1, int particle2, const vector<double>& parameters) {
exceptions[index].particle1 = particle1;
exceptions[index].particle2 = particle2;
exceptions[index].parameters = parameters;
int CustomNonbondedForce::addExclusion(int particle1, int particle2) {
exclusions.push_back(ExclusionInfo(particle1, particle2));
return exclusions.size()-1;
}
void CustomNonbondedForce::getExclusionParticles(int index, int& particle1, int& particle2) const {
particle1 = exclusions[index].particle1;
particle2 = exclusions[index].particle2;
}
void CustomNonbondedForce::setExclusionParticles(int index, int particle1, int particle2) {
exclusions[index].particle1 = particle1;
exclusions[index].particle2 = particle2;
}
int CustomNonbondedForce::addFunction(const std::string& name, const std::vector<double>& values, double min, double max, bool interpolating) {
......
openmmapi/src/CustomNonbondedForceImpl.cpp
View file @ 70015b75
......@@ -52,14 +52,14 @@ CustomNonbondedForceImpl::~CustomNonbondedForceImpl() {
void CustomNonbondedForceImpl::initialize(ContextImpl& context) {
kernel = context.getPlatform().createKernel(CalcCustomNonbondedForceKernel::Name(), context);
// Check for errors in the specification of parameters and exceptions.
// Check for errors in the specification of parameters and exclusions.
System& system = context.getSystem();
if (owner.getNumParticles() != system.getNumParticles())
throw OpenMMException("CustomNonbondedForce must have exactly as many particles as the System it belongs to.");
vector<set<int> > exceptions(owner.getNumParticles());
vector<set<int> > exclusions(owner.getNumParticles());
vector<double> parameters;
int numParameters = owner.getNumParameters();
int numParameters = owner.getNumPerParticleParameters();
for (int i = 0; i < owner.getNumParticles(); i++) {
owner.getParticleParameters(i, parameters);
if (parameters.size() != numParameters) {
......@@ -69,37 +69,31 @@ void CustomNonbondedForceImpl::initialize(ContextImpl& context) {
throw OpenMMException(msg.str());
}
}
for (int i = 0; i < owner.getNumExceptions(); i++) {
for (int i = 0; i < owner.getNumExclusions(); i++) {
int particle1, particle2;
owner.getExceptionParameters(i, particle1, particle2, parameters);
owner.getExclusionParticles(i, particle1, particle2);
if (particle1 < 0 || particle1 >= owner.getNumParticles()) {
stringstream msg;
msg << "CustomNonbondedForce: Illegal particle index for an exception: ";
msg << "CustomNonbondedForce: Illegal particle index for an exclusion: ";
msg << particle1;
throw OpenMMException(msg.str());
}
if (particle2 < 0 || particle2 >= owner.getNumParticles()) {
stringstream msg;
msg << "CustomNonbondedForce: Illegal particle index for an exception: ";
msg << "CustomNonbondedForce: Illegal particle index for an exclusion: ";
msg << particle2;
throw OpenMMException(msg.str());
}
if (exceptions[particle1].count(particle2) > 0 || exceptions[particle2].count(particle1) > 0) {
if (exclusions[particle1].count(particle2) > 0 || exclusions[particle2].count(particle1) > 0) {
stringstream msg;
msg << "CustomNonbondedForce: Multiple exceptions are specified for particles ";
msg << "CustomNonbondedForce: Multiple exclusions are specified for particles ";
msg << particle1;
msg << " and ";
msg << particle2;
throw OpenMMException(msg.str());
}
if (parameters.size() != 0 && parameters.size() != numParameters) {
stringstream msg;
msg << "CustomNonbondedForce: Wrong number of parameters for exception ";
msg << i;
throw OpenMMException(msg.str());
}
exceptions[particle1].insert(particle2);
exceptions[particle2].insert(particle1);
exclusions[particle1].insert(particle2);
exclusions[particle2].insert(particle1);
}
dynamic_cast<CalcCustomNonbondedForceKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner);
}
......
platforms/cuda/src/CudaKernels.cpp
View file @ 70015b75
......@@ -463,21 +463,6 @@ void CudaCalcCustomNonbondedForceKernel::initialize(const System& system, const
numParticles = force.getNumParticles();
_gpuContext* gpu = data.gpu;
// Identify which exceptions are actual interactions.
vector<pair<int, int> > exclusions;
vector<int> exceptions;
{
vector<double> parameters;
for (int i = 0; i < force.getNumExceptions(); i++) {
int particle1, particle2;
force.getExceptionParameters(i, particle1, particle2, parameters);
exclusions.push_back(pair<int, int>(particle1, particle2));
if (parameters.size() > 0)
exceptions.push_back(i);
}
}
// Initialize nonbonded interactions.
vector<int> particle(numParticles);
......@@ -488,9 +473,11 @@ void CudaCalcCustomNonbondedForceKernel::initialize(const System& system, const
particle[i] = i;
exclusionList[i].push_back(i);
}
for (int i = 0; i < (int)exclusions.size(); i++) {
exclusionList[exclusions[i].first].push_back(exclusions[i].second);
exclusionList[exclusions[i].second].push_back(exclusions[i].first);
for (int i = 0; i < force.getNumExclusions(); i++) {
int particle1, particle2;
force.getExclusionParticles(i, particle1, particle2);
exclusionList[particle1].push_back(particle2);
exclusionList[particle2].push_back(particle1);
}
Vec3 boxVectors[3];
system.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
......@@ -503,15 +490,6 @@ void CudaCalcCustomNonbondedForceKernel::initialize(const System& system, const
}
data.customNonbondedMethod = method;
// Initialize exceptions.
int numExceptions = exceptions.size();
vector<int> exceptionParticle1(numExceptions);
vector<int> exceptionParticle2(numExceptions);
vector<vector<double> > exceptionParams(numExceptions);
for (int i = 0; i < numExceptions; i++)
force.getExceptionParameters(exceptions[i], exceptionParticle1[i], exceptionParticle2[i], exceptionParams[i]);
// Record the tabulated functions.
for (int i = 0; i < force.getNumFunctions(); i++) {
......@@ -526,19 +504,15 @@ void CudaCalcCustomNonbondedForceKernel::initialize(const System& system, const
// Record information for the expressions.
vector<string> paramNames;
vector<string> combiningRules;
for (int i = 0; i < force.getNumParameters(); i++) {
paramNames.push_back(force.getParameterName(i));
combiningRules.push_back(force.getParameterCombiningRule(i));
}
for (int i = 0; i < force.getNumPerParticleParameters(); i++)
paramNames.push_back(force.getPerParticleParameterName(i));
globalParamNames.resize(force.getNumGlobalParameters());
globalParamValues.resize(force.getNumGlobalParameters());
for (int i = 0; i < force.getNumGlobalParameters(); i++) {
globalParamNames[i] = force.getGlobalParameterName(i);
globalParamValues[i] = (float) force.getGlobalParameterDefaultValue(i);
}
gpuSetCustomNonbondedParameters(gpu, parameters, exclusionList, exceptionParticle1, exceptionParticle2, exceptionParams, method,
(float)force.getCutoffDistance(), force.getEnergyFunction(), combiningRules, paramNames, globalParamNames);
gpuSetCustomNonbondedParameters(gpu, parameters, exclusionList, method, (float) force.getCutoffDistance(), force.getEnergyFunction(), paramNames, globalParamNames);
if (globalParamValues.size() > 0)
SetCustomNonbondedGlobalParams(&globalParamValues[0]);
}
......
platforms/cuda/src/kernels/cudaKernels.h
View file @ 70015b75
......@@ -113,5 +113,4 @@ extern void SetCustomBondEnergyExpression(const Expression<128>& expression);
extern void SetCustomBondGlobalParams(float* paramValues);
extern void SetCustomNonbondedForceExpression(const Expression<128>& expression);
extern void SetCustomNonbondedEnergyExpression(const Expression<128>& expression);
extern void SetCustomNonbondedCombiningRules(const Expression<64>* expressions);
extern void SetCustomNonbondedGlobalParams(float* paramValues);
platforms/cuda/src/kernels/cudatypes.h
View file @ 70015b75
......@@ -249,7 +249,7 @@ enum CudaNonbondedMethod
};
enum ExpressionOp {
VARIABLE0 = 0, VARIABLE1, VARIABLE2, VARIABLE3, VARIABLE4, VARIABLE5, VARIABLE6, VARIABLE7, MULTIPLY, DIVIDE, ADD, SUBTRACT, POWER, MULTIPLY_CONSTANT, POWER_CONSTANT, ADD_CONSTANT,
VARIABLE0 = 0, VARIABLE1, VARIABLE2, VARIABLE3, VARIABLE4, VARIABLE5, VARIABLE6, VARIABLE7, VARIABLE8, MULTIPLY, DIVIDE, ADD, SUBTRACT, POWER, MULTIPLY_CONSTANT, POWER_CONSTANT, ADD_CONSTANT,
GLOBAL, CONSTANT, CUSTOM, CUSTOM_DERIV, NEGATE, RECIPROCAL, SQRT, EXP, LOG, SQUARE, CUBE, SIN, COS, SEC, CSC, TAN, COT, ASIN, ACOS, ATAN, SINH, COSH, TANH
};
......
platforms/cuda/src/kernels/gpu.cpp
View file @ 70015b75
......@@ -174,6 +174,8 @@ static Expression<SIZE> createExpression(gpuContext gpu, const string& expressio
exp.op[i] = VARIABLE6;
else if (variables.size() > 7 && op.getName() == variables[7])
exp.op[i] = VARIABLE7;
else if (variables.size() > 8 && op.getName() == variables[8])
exp.op[i] = VARIABLE8;
else {
int j;
for (j = 0; j < globalParamNames.size() && op.getName() != globalParamNames[j]; j++);
......@@ -693,8 +695,7 @@ void gpuSetCustomBondParameters(gpuContext gpu, const vector<int>& bondAtom1, co
extern "C"
void gpuSetCustomNonbondedParameters(gpuContext gpu, const vector<vector<double> >& parameters, const vector<vector<int> >& exclusions,
const vector<int>& exceptionAtom1, const vector<int>& exceptionAtom2, const vector<vector<double> >& exceptionParams,
CudaNonbondedMethod method, float cutoffDistance, const string& energyExp, const vector<string>& combiningRules,
CudaNonbondedMethod method, float cutoffDistance, const string& energyExp,
const vector<string>& paramNames, const vector<string>& globalParamNames)
{
if (gpu->sim.nonbondedCutoff != 0.0f && gpu->sim.nonbondedCutoff != cutoffDistance)
......@@ -706,20 +707,10 @@ void gpuSetCustomNonbondedParameters(gpuContext gpu, const vector<vector<double>
gpu->sim.nonbondedCutoff = cutoffDistance;
gpu->sim.nonbondedCutoffSqr = cutoffDistance*cutoffDistance;
gpu->sim.customNonbondedMethod = method;
gpu->sim.customExceptions = exceptionAtom1.size();
gpu->sim.customParameters = paramNames.size();
gpu->sim.custom_exception_threads_per_block = (gpu->sim.customExceptions+gpu->sim.blocks-1)/gpu->sim.blocks;
if (gpu->sim.custom_exception_threads_per_block < 1)
gpu->sim.custom_exception_threads_per_block = 1;
if (gpu->sim.custom_exception_threads_per_block > gpu->sim.max_localForces_threads_per_block)
gpu->sim.custom_exception_threads_per_block = gpu->sim.max_localForces_threads_per_block;
setExclusions(gpu, exclusions);
gpu->psCustomParams = new CUDAStream<float4>(gpu->sim.paddedNumberOfAtoms, 1, "CustomParams");
gpu->sim.pCustomParams = gpu->psCustomParams->_pDevData;
gpu->psCustomExceptionID = new CUDAStream<int4>(gpu->sim.customExceptions, 1, "CustomExceptionId");
gpu->sim.pCustomExceptionID = gpu->psCustomExceptionID->_pDevData;
gpu->psCustomExceptionParams = new CUDAStream<float4>(gpu->sim.customExceptions, 1, "CustomExceptionParams");
gpu->sim.pCustomExceptionParams = gpu->psCustomExceptionParams->_pDevData;
for (int i = 0; i < (int) parameters.size(); i++) {
if (parameters[i].size() > 0)
(*gpu->psCustomParams)[i].x = (float) parameters[i][0];
......@@ -730,23 +721,7 @@ void gpuSetCustomNonbondedParameters(gpuContext gpu, const vector<vector<double>
if (parameters[i].size() > 3)
(*gpu->psCustomParams)[i].w = (float) parameters[i][3];
}
for (int i = 0; i < (int) exceptionAtom1.size(); i++) {
(*gpu->psCustomExceptionID)[i].x = exceptionAtom1[i];
(*gpu->psCustomExceptionID)[i].y = exceptionAtom2[i];
(*gpu->psCustomExceptionID)[i].z = gpu->pOutputBufferCounter[exceptionAtom1[i]]++;
(*gpu->psCustomExceptionID)[i].w = gpu->pOutputBufferCounter[exceptionAtom2[i]]++;
if (exceptionParams[i].size() > 0)
(*gpu->psCustomExceptionParams)[i].x = exceptionParams[i][0];
if (exceptionParams[i].size() > 1)
(*gpu->psCustomExceptionParams)[i].y = exceptionParams[i][1];
if (exceptionParams[i].size() > 2)
(*gpu->psCustomExceptionParams)[i].z = exceptionParams[i][2];
if (exceptionParams[i].size() > 3)
(*gpu->psCustomExceptionParams)[i].w = exceptionParams[i][3];
}
gpu->psCustomParams->Upload();
gpu->psCustomExceptionID->Upload();
gpu->psCustomExceptionParams->Upload();
// This class serves as a placeholder for custom functions in expressions.
......@@ -784,25 +759,18 @@ void gpuSetCustomNonbondedParameters(gpuContext gpu, const vector<vector<double>
// Create the Expressions.
vector<string> variables;
variables.push_back("r");
for (int i = 0; i < (int) paramNames.size(); i++)
variables.push_back(paramNames[i]);
SetCustomNonbondedEnergyExpression(createExpression<128>(gpu, energyExp, Lepton::Parser::parse(energyExp, functions).optimize().createProgram(), variables, globalParamNames, gpu->sim.customExpressionStackSize));
SetCustomNonbondedForceExpression(createExpression<128>(gpu, energyExp, Lepton::Parser::parse(energyExp, functions).differentiate("r").optimize().createProgram(), variables, globalParamNames, gpu->sim.customExpressionStackSize));
Expression<64> paramExpressions[4];
vector<string> combiningRuleParams;
for (int j = 1; j < 3; j++) {
for (int i = 0; i < (int) paramNames.size(); i++) {
stringstream name;
name << paramNames[i] << j;
combiningRuleParams.push_back(name.str());
variables.push_back(name.str());
}
for (int i = paramNames.size(); i < 4; i++)
combiningRuleParams.push_back("");
variables.push_back("");
}
for (int i = 0; i < (int) paramNames.size(); i++)
paramExpressions[i] = createExpression<64>(gpu, combiningRules[i], Lepton::Parser::parse(combiningRules[i], functions).optimize().createProgram(), combiningRuleParams, globalParamNames, gpu->sim.customExpressionStackSize);
SetCustomNonbondedCombiningRules(paramExpressions);
variables.push_back("r");
SetCustomNonbondedEnergyExpression(createExpression<128>(gpu, energyExp, Lepton::Parser::parse(energyExp, functions).optimize().createProgram(), variables, globalParamNames, gpu->sim.customExpressionStackSize));
SetCustomNonbondedForceExpression(createExpression<128>(gpu, energyExp, Lepton::Parser::parse(energyExp, functions).differentiate("r").optimize().createProgram(), variables, globalParamNames, gpu->sim.customExpressionStackSize));
delete fp;
}
......
platforms/cuda/src/kernels/gputypes.h
View file @ 70015b75
......@@ -214,8 +214,7 @@ void gpuSetCustomBondParameters(gpuContext gpu, const std::vector<int>& bondAtom
extern "C"
void gpuSetCustomNonbondedParameters(gpuContext gpu, const std::vector<std::vector<double> >& parameters, const std::vector<std::vector<int> >& exclusions,
const std::vector<int>& exceptionAtom1, const std::vector<int>& exceptionAtom2, const std::vector<std::vector<double> >& exceptionParams,
CudaNonbondedMethod method, float cutoffDistance, const std::string& energyExp, const std::vector<std::string>& combiningRules,
CudaNonbondedMethod method, float cutoffDistance, const std::string& energyExp,
const std::vector<std::string>& paramNames, const std::vector<std::string>& globalParamNames);
extern "C"
......
platforms/cuda/src/kernels/kCalculateCDLJForces.h
View file @ 70015b75
......@@ -36,7 +36,7 @@
__global__ void METHOD_NAME(kCalculateCDLJ, Forces_kernel)(unsigned int* workUnit)
{
extern __shared__ Atom sA[];
unsigned int totalWarps = cSim.nonbond_blocks*cSim.nonbond_threads_per_block/GRID;
unsigned int totalWarps = gridDim.x*blockDim.x/GRID;
unsigned int warp = (blockIdx.x*blockDim.x+threadIdx.x)/GRID;
unsigned int numWorkUnits = cSim.pInteractionCount[0];
unsigned int pos = warp*numWorkUnits/totalWarps;
......
platforms/cuda/src/kernels/kCalculateCustomNonbondedForces.cu
View file @ 70015b75
......@@ -52,7 +52,6 @@ struct Atom {
static __constant__ cudaGmxSimulation cSim;
static __constant__ Expression<128> forceExp;
static __constant__ Expression<128> energyExp;
static __constant__ Expression<64> combiningRules[4];
#include "kEvaluateExpression.h"
......@@ -84,13 +83,6 @@ void SetCustomNonbondedEnergyExpression(const Expression<128>& expression)
RTERROR(status, "SetCustomNonbondedEnergyExpression: cudaMemcpyToSymbol failed");
}
void SetCustomNonbondedCombiningRules(const Expression<64>* expressions)
{
cudaError_t status;
status = cudaMemcpyToSymbol(combiningRules, expressions, sizeof(combiningRules));
RTERROR(status, "SetCustomNonbondedCombiningRules: cudaMemcpyToSymbol failed");
}
void SetCustomNonbondedGlobalParams(float* paramValues)
{
cudaError_t status;
......@@ -154,7 +146,7 @@ void kCalculateCustomNonbondedForces(gpuContext gpu, bool neighborListValid)
cudaBindTexture(NULL, &texRef3, gpu->tabulatedFunctions[3].coefficients->_pDevData, &channelDesc, gpu->tabulatedFunctions[3].coefficients->_length*sizeof(float4));
gpu->tabulatedFunctionsChanged = false;
}
int sharedPerThread = sizeof(Atom)+gpu->sim.customExpressionStackSize*sizeof(float)+8*sizeof(float);
int sharedPerThread = sizeof(Atom)+gpu->sim.customExpressionStackSize*sizeof(float)+9*sizeof(float);
if (gpu->sim.customNonbondedMethod != NO_CUTOFF)
sharedPerThread += sizeof(float3);
int threads = gpu->sim.nonbond_threads_per_block;
......@@ -169,9 +161,6 @@ void kCalculateCustomNonbondedForces(gpuContext gpu, bool neighborListValid)
else
kCalculateCustomNonbondedN2Forces_kernel<<<gpu->sim.nonbond_blocks, threads, sharedPerThread*threads>>>(gpu->sim.pWorkUnit);
LAUNCHERROR("kCalculateCustomNonbondedN2Forces");
kCalculateCustomNonbondedN2Exceptions_kernel<<<gpu->sim.blocks, gpu->sim.custom_exception_threads_per_block,
gpu->sim.customExpressionStackSize*sizeof(float)*gpu->sim.custom_exception_threads_per_block>>>();
LAUNCHERROR("kCalculateCustomNonbondedN2Exceptions");
break;
case CUTOFF:
if (!neighborListValid)
......@@ -189,9 +178,6 @@ void kCalculateCustomNonbondedForces(gpuContext gpu, bool neighborListValid)
else
kCalculateCustomNonbondedCutoffForces_kernel<<<gpu->sim.nonbond_blocks, threads, sharedPerThread*threads>>>(gpu->sim.pInteractingWorkUnit);
LAUNCHERROR("kCalculateCustomNonbondedCutoffForces");
kCalculateCustomNonbondedCutoffExceptions_kernel<<<gpu->sim.blocks, gpu->sim.custom_exception_threads_per_block,
gpu->sim.customExpressionStackSize*sizeof(float)*gpu->sim.custom_exception_threads_per_block>>>();
LAUNCHERROR("kCalculateCustomNonbondedCutoffExceptions");
break;
case PERIODIC:
if (!neighborListValid)
......@@ -209,9 +195,6 @@ void kCalculateCustomNonbondedForces(gpuContext gpu, bool neighborListValid)
else
kCalculateCustomNonbondedPeriodicForces_kernel<<<gpu->sim.nonbond_blocks, threads, sharedPerThread*threads>>>(gpu->sim.pInteractingWorkUnit);
LAUNCHERROR("kCalculateCustomNonbondedPeriodicForces");
kCalculateCustomNonbondedPeriodicExceptions_kernel<<<gpu->sim.blocks, gpu->sim.custom_exception_threads_per_block,
(gpu->sim.customExpressionStackSize+8)*sizeof(float)*gpu->sim.custom_exception_threads_per_block>>>();
LAUNCHERROR("kCalculateCustomNonbondedPeriodicExceptions");
break;
}
}
platforms/cuda/src/kernels/kCalculateCustomNonbondedForces.h
View file @ 70015b75
......@@ -35,7 +35,7 @@ __global__ void METHOD_NAME(kCalculateCustomNonbonded, Forces_kernel)(unsigned i
extern __shared__ float stack[];
Atom* sA = (Atom*) &stack[cSim.customExpressionStackSize*blockDim.x];
float* variables = (float*) &sA[blockDim.x];
unsigned int totalWarps = cSim.nonbond_blocks*cSim.nonbond_threads_per_block/GRID;
unsigned int totalWarps = gridDim.x*blockDim.x/GRID;
unsigned int warp = (blockIdx.x*blockDim.x+threadIdx.x)/GRID;
unsigned int numWorkUnits = cSim.pInteractionCount[0];
unsigned int pos = warp*numWorkUnits/totalWarps;
......@@ -48,7 +48,6 @@ __global__ void METHOD_NAME(kCalculateCustomNonbonded, Forces_kernel)(unsigned i
unsigned int lasty = 0xFFFFFFFF;
while (pos < end)
{
// Extract cell coordinates from appropriate work unit
unsigned int x = workUnit[pos];
unsigned int y = ((x >> 2) & 0x7fff) << GRIDBITS;
......@@ -78,9 +77,8 @@ __global__ void METHOD_NAME(kCalculateCustomNonbonded, Forces_kernel)(unsigned i
unsigned int excl = cSim.pExclusion[cSim.pExclusionIndex[cell]+tgx];
for (unsigned int j = 0; j < GRID; j++)
{
// Apply the combining rules to the parameters.
// Record the parameters.
float combinedParams[4];
VARIABLE(0) = params.x;
VARIABLE(1) = params.y;
VARIABLE(2) = params.z;
......@@ -89,29 +87,6 @@ __global__ void METHOD_NAME(kCalculateCustomNonbonded, Forces_kernel)(unsigned i
VARIABLE(5) = psA[j].params.y;
VARIABLE(6) = psA[j].params.z;
VARIABLE(7) = psA[j].params.w;
for (int k = 0; k < cSim.customParameters; k++)
{
float value = kEvaluateExpression_kernel(&combiningRules[k], stack, variables);
switch (k)
{
case 0:
combinedParams[0] = value;
break;
case 1:
combinedParams[1] = value;
break;
case 2:
combinedParams[2] = value;
break;
case 3:
combinedParams[3] = value;
break;
}
}
VARIABLE(1) = combinedParams[0];
VARIABLE(2) = combinedParams[1];
VARIABLE(3) = combinedParams[2];
VARIABLE(4) = combinedParams[3];
// Compute the force.
......@@ -125,7 +100,7 @@ __global__ void METHOD_NAME(kCalculateCustomNonbonded, Forces_kernel)(unsigned i
#endif
float r = sqrt(dx*dx + dy*dy + dz*dz);
float invR = 1.0f/r;
VARIABLE(0) = r;
VARIABLE(8) = r;
float dEdR = -kEvaluateExpression_kernel(&forceExp, stack, variables)*invR;
float energy = kEvaluateExpression_kernel(&energyExp, stack, variables);
#ifdef USE_CUTOFF
......@@ -195,9 +170,8 @@ __global__ void METHOD_NAME(kCalculateCustomNonbonded, Forces_kernel)(unsigned i
for (unsigned int j = 0; j < GRID; j++)
{
// Apply the combining rules to the parameters.
// Record the parameters.
float combinedParams[4];
VARIABLE(0) = params.x;
VARIABLE(1) = params.y;
VARIABLE(2) = params.z;
......@@ -206,29 +180,6 @@ __global__ void METHOD_NAME(kCalculateCustomNonbonded, Forces_kernel)(unsigned i
VARIABLE(5) = psA[tj].params.y;
VARIABLE(6) = psA[tj].params.z;
VARIABLE(7) = psA[tj].params.w;
for (int k = 0; k < cSim.customParameters; k++)
{
float value = kEvaluateExpression_kernel(&combiningRules[k], stack, variables);
switch (k)
{
case 0:
combinedParams[0] = value;
break;
case 1:
combinedParams[1] = value;
break;
case 2:
combinedParams[2] = value;
break;
case 3:
combinedParams[3] = value;
break;
}
}
VARIABLE(1) = combinedParams[0];
VARIABLE(2) = combinedParams[1];
VARIABLE(3) = combinedParams[2];
VARIABLE(4) = combinedParams[3];
// Compute the force.
......@@ -242,7 +193,7 @@ __global__ void METHOD_NAME(kCalculateCustomNonbonded, Forces_kernel)(unsigned i
#endif
float r = sqrt(dx*dx + dy*dy + dz*dz);
float invR = 1.0f/r;
VARIABLE(0) = r;
VARIABLE(8) = r;
float dEdR = -kEvaluateExpression_kernel(&forceExp, stack, variables)*invR;
float energy = kEvaluateExpression_kernel(&energyExp, stack, variables);
#ifdef USE_CUTOFF
......@@ -274,9 +225,8 @@ __global__ void METHOD_NAME(kCalculateCustomNonbonded, Forces_kernel)(unsigned i
{
if ((flags&(1<<j)) != 0)
{
// Apply the combining rules to the parameters.
// Record the parameters.
float combinedParams[4];
VARIABLE(0) = params.x;
VARIABLE(1) = params.y;
VARIABLE(2) = params.z;
......@@ -285,29 +235,6 @@ __global__ void METHOD_NAME(kCalculateCustomNonbonded, Forces_kernel)(unsigned i
VARIABLE(5) = psA[j].params.y;
VARIABLE(6) = psA[j].params.z;
VARIABLE(7) = psA[j].params.w;
for (int k = 0; k < cSim.customParameters; k++)
{
float value = kEvaluateExpression_kernel(&combiningRules[k], stack, variables);
switch (k)
{
case 0:
combinedParams[0] = value;
break;
case 1:
combinedParams[1] = value;
break;
case 2:
combinedParams[2] = value;
break;
case 3:
combinedParams[3] = value;
break;
}
}
VARIABLE(1) = combinedParams[0];
VARIABLE(2) = combinedParams[1];
VARIABLE(3) = combinedParams[2];
VARIABLE(4) = combinedParams[3];
// Compute the force.
......@@ -321,7 +248,7 @@ __global__ void METHOD_NAME(kCalculateCustomNonbonded, Forces_kernel)(unsigned i
#endif
float r = sqrt(dx*dx + dy*dy + dz*dz);
float invR = 1.0f/r;
VARIABLE(0) = r;
VARIABLE(8) = r;
float dEdR = -kEvaluateExpression_kernel(&forceExp, stack, variables)*invR;
float energy = kEvaluateExpression_kernel(&energyExp, stack, variables);
#ifdef USE_CUTOFF
......@@ -389,9 +316,8 @@ __global__ void METHOD_NAME(kCalculateCustomNonbonded, Forces_kernel)(unsigned i
excl = (excl >> tgx) | (excl << (GRID - tgx));
for (unsigned int j = 0; j < GRID; j++)
{
// Apply the combining rules to the parameters.
// Record the parameters.
float combinedParams[4];
VARIABLE(0) = params.x;
VARIABLE(1) = params.y;
VARIABLE(2) = params.z;
......@@ -400,29 +326,6 @@ __global__ void METHOD_NAME(kCalculateCustomNonbonded, Forces_kernel)(unsigned i
VARIABLE(5) = psA[tj].params.y;
VARIABLE(6) = psA[tj].params.z;
VARIABLE(7) = psA[tj].params.w;
for (int k = 0; k < cSim.customParameters; k++)
{
float value = kEvaluateExpression_kernel(&combiningRules[k], stack, variables);
switch (k)
{
case 0:
combinedParams[0] = value;
break;
case 1:
combinedParams[1] = value;
break;
case 2:
combinedParams[2] = value;
break;
case 3:
combinedParams[3] = value;
break;
}
}
VARIABLE(1) = combinedParams[0];
VARIABLE(2) = combinedParams[1];
VARIABLE(3) = combinedParams[2];
VARIABLE(4) = combinedParams[3];
// Compute the force.
......@@ -436,7 +339,7 @@ __global__ void METHOD_NAME(kCalculateCustomNonbonded, Forces_kernel)(unsigned i
#endif
float r = sqrt(dx*dx + dy*dy + dz*dz);
float invR = 1.0f/r;
VARIABLE(0) = r;
VARIABLE(8) = r;
float dEdR = -kEvaluateExpression_kernel(&forceExp, stack, variables)*invR;
float energy = kEvaluateExpression_kernel(&energyExp, stack, variables);
#ifdef USE_CUTOFF
......@@ -498,62 +401,3 @@ __global__ void METHOD_NAME(kCalculateCustomNonbonded, Forces_kernel)(unsigned i
}
cSim.pEnergy[blockIdx.x*blockDim.x+threadIdx.x] += totalEnergy;
}
__global__ void METHOD_NAME(kCalculateCustomNonbonded, Exceptions_kernel)()
{
extern __shared__ float stack[];
float* variables = (float*) &stack[cSim.customExpressionStackSize*blockDim.x];
unsigned int pos = blockIdx.x * blockDim.x + threadIdx.x;
float totalEnergy = 0.0f;
while (pos < cSim.customExceptions)
{
int4 atom = cSim.pCustomExceptionID[pos];
float4 params = cSim.pCustomExceptionParams[pos];
float4 a1 = cSim.pPosq[atom.x];
float4 a2 = cSim.pPosq[atom.y];
float dx = a1.x - a2.x;
float dy = a1.y - a2.y;
float dz = a1.z - a2.z;
#ifdef USE_PERIODIC
dx -= floor(dx/cSim.periodicBoxSizeX+0.5f)*cSim.periodicBoxSizeX;
dy -= floor(dy/cSim.periodicBoxSizeY+0.5f)*cSim.periodicBoxSizeY;
dz -= floor(dz/cSim.periodicBoxSizeZ+0.5f)*cSim.periodicBoxSizeZ;
#endif
float r = sqrt(dx*dx + dy*dy + dz*dz);
float invR = 1.0f/r;
VARIABLE(0) = r;
VARIABLE(1) = params.x;
VARIABLE(2) = params.y;
VARIABLE(3) = params.z;
VARIABLE(4) = params.w;
float dEdR = -kEvaluateExpression_kernel(&forceExp, stack, variables)*invR;
float energy = kEvaluateExpression_kernel(&energyExp, stack, variables);
#ifdef USE_CUTOFF
if (r > cSim.nonbondedCutoff)
{
dEdR = 0.0f;
energy = 0.0f;
}
#endif
totalEnergy += energy;
dx *= dEdR;
dy *= dEdR;
dz *= dEdR;
unsigned int offsetA = atom.x + atom.z * cSim.stride;
unsigned int offsetB = atom.y + atom.w * cSim.stride;
float4 forceA = cSim.pForce4[offsetA];
float4 forceB = cSim.pForce4[offsetB];
forceA.x += dx;
forceA.y += dy;
forceA.z += dz;
forceB.x -= dx;
forceB.y -= dy;
forceB.z -= dz;
cSim.pForce4[offsetA] = forceA;
cSim.pForce4[offsetB] = forceB;
pos += blockDim.x * gridDim.x;
}
cSim.pEnergy[blockIdx.x * blockDim.x + threadIdx.x] += totalEnergy;
}
platforms/cuda/tests/TestCudaCustomNonbondedForce.cpp
View file @ 70015b75
......@@ -79,9 +79,9 @@ void testParameters() {
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomNonbondedForce* forceField = new CustomNonbondedForce("scale*a*(r*b)^3");
forceField->addParameter("a", "a1*a2");
forceField->addParameter("b", "c+b1+b2");
CustomNonbondedForce* forceField = new CustomNonbondedForce("scale*a*(r*b)^3; a=a1*a2; b=c+b1+b2");
forceField->addPerParticleParameter("a");
forceField->addPerParticleParameter("b");
forceField->addGlobalParameter("scale", 3.0);
forceField->addGlobalParameter("c", -1.0);
vector<double> params(2);
......@@ -115,12 +115,12 @@ void testParameters() {
ASSERT_EQUAL_TOL(1.5*3.0*(12*12*12), state.getPotentialEnergy(), TOL);
}
void testExceptions() {
void testExclusions() {
CudaPlatform platform;
System system;
VerletIntegrator integrator(0.01);
CustomNonbondedForce* nonbonded = new CustomNonbondedForce("a*r");
nonbonded->addParameter("a", "a1+a2");
CustomNonbondedForce* nonbonded = new CustomNonbondedForce("a*r; a=a1+a2");
nonbonded->addPerParticleParameter("a");
vector<double> params(1);
vector<Vec3> positions(4);
for (int i = 0; i < 4; i++) {
......@@ -129,22 +129,21 @@ void testExceptions() {
nonbonded->addParticle(params);
positions[i] = Vec3(i, 0, 0);
}
nonbonded->addException(0, 1, vector<double>());
nonbonded->addException(1, 2, vector<double>());
nonbonded->addException(2, 3, vector<double>());
params[0] = 0.5;
nonbonded->addException(0, 2, params);
nonbonded->addException(1, 3, params);
nonbonded->addExclusion(0, 1);
nonbonded->addExclusion(1, 2);
nonbonded->addExclusion(2, 3);
nonbonded->addExclusion(0, 2);
nonbonded->addExclusion(1, 3);
system.addForce(nonbonded);
Context context(system, integrator, platform);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
ASSERT_EQUAL_VEC(Vec3(0.5+1+4, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0.5, 0, 0), forces[1], TOL);
ASSERT_EQUAL_VEC(Vec3(-0.5, 0, 0), forces[2], TOL);
ASSERT_EQUAL_VEC(Vec3(-(0.5+1+4), 0, 0), forces[3], TOL);
ASSERT_EQUAL_TOL((1+4)*3+0.5*2+0.5*2, state.getPotentialEnergy(), TOL);
ASSERT_EQUAL_VEC(Vec3(1+4, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[1], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[2], TOL);
ASSERT_EQUAL_VEC(Vec3(-(1+4), 0, 0), forces[3], TOL);
ASSERT_EQUAL_TOL((1+4)*3, state.getPotentialEnergy(), TOL);
}
void testCutoff() {
......@@ -251,10 +250,10 @@ void testCoulombLennardJones() {
customSystem.addParticle(1.0);
}
NonbondedForce* standardNonbonded = new NonbondedForce();
CustomNonbondedForce* customNonbonded = new CustomNonbondedForce("4*eps*((sigma/r)^12-(sigma/r)^6) + 138.935485*q/r");
customNonbonded->addParameter("q", "q1*q2");
customNonbonded->addParameter("sigma", "0.5*(sigma1+sigma2)");
customNonbonded->addParameter("eps", "sqrt(eps1*eps2)");
CustomNonbondedForce* customNonbonded = new CustomNonbondedForce("4*eps*((sigma/r)^12-(sigma/r)^6)+138.935485*q/r; q=q1*q2; sigma=0.5*(sigma1+sigma2); eps=sqrt(eps1*eps2)");
customNonbonded->addPerParticleParameter("q");
customNonbonded->addPerParticleParameter("sigma");
customNonbonded->addPerParticleParameter("eps");
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
init_gen_rand(0);
......@@ -266,7 +265,8 @@ void testCoulombLennardJones() {
params[1] = 0.2;
params[2] = 0.1;
customNonbonded->addParticle(params);
standardNonbonded->addParticle(1.0, 0.1, 0.1);
standardNonbonded->addParticle(-1.0, 0.1, 0.1);
params[0] = -1.0;
params[1] = 0.1;
customNonbonded->addParticle(params);
}
......@@ -276,7 +276,8 @@ void testCoulombLennardJones() {
params[1] = 0.2;
params[2] = 0.2;
customNonbonded->addParticle(params);
standardNonbonded->addParticle(1.0, 0.1, 0.2);
standardNonbonded->addParticle(-1.0, 0.1, 0.2);
params[0] = -1.0;
params[1] = 0.1;
customNonbonded->addParticle(params);
}
......@@ -285,7 +286,7 @@ void testCoulombLennardJones() {
velocities[2*i] = Vec3(genrand_real2(), genrand_real2(), genrand_real2());
velocities[2*i+1] = Vec3(genrand_real2(), genrand_real2(), genrand_real2());
standardNonbonded->addException(2*i, 2*i+1, 0.0, 1.0, 0.0);
customNonbonded->addException(2*i, 2*i+1, vector<double>());
customNonbonded->addExclusion(2*i, 2*i+1);
}
standardNonbonded->setNonbondedMethod(NonbondedForce::NoCutoff);
customNonbonded->setNonbondedMethod(CustomNonbondedForce::NoCutoff);
......@@ -294,34 +295,29 @@ void testCoulombLennardJones() {
VerletIntegrator integrator1(0.01);
VerletIntegrator integrator2(0.01);
Context context1(standardSystem, integrator1, platform);
// Context context2(customSystem, integrator2, platform);
context1.setPositions(positions);
// context2.setPositions(positions);
context1.setVelocities(velocities);
// context2.setVelocities(velocities);
State state1 = context1.getState(State::Forces | State::Energy);
// State state2 = context2.getState(State::Forces | State::Energy);
for (int i = 0; i < numParticles; i++)
std::cout << i<<": "<<state1.getForces()[i]<< std::endl;
// for (int i = 0; i < numParticles; i++)
// std::cout << i<<": "<<state1.getForces()[i]<<" "<<state2.getForces()[i]<< std::endl;
// std::cout <<state1.getPotentialEnergy()<<" "<<state2.getPotentialEnergy() << std::endl;
// ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-4);
// for (int i = 0; i < numParticles; i++) {
// ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-4);
// }
Context context2(customSystem, integrator2, platform);
context2.setPositions(positions);
context2.setVelocities(velocities);
State state2 = context2.getState(State::Forces | State::Energy);
ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-4);
for (int i = 0; i < numParticles; i++) {
ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-4);
}
}
int main() {
try {
testSimpleExpression();
testParameters();
testExceptions();
testExclusions();
testCutoff();
testPeriodic();
testTabulatedFunction(true);
testTabulatedFunction(false);
// testCoulombLennardJones();
testCoulombLennardJones();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
......
platforms/opencl/src/OpenCLKernels.cpp
View file @ 70015b75
......@@ -859,25 +859,16 @@ public:
return true;
}
int getNumParticleGroups() {
return force.getNumExceptions();
return force.getNumExclusions();
}
void getParticlesInGroup(int index, std::vector<int>& particles) {
int particle1, particle2;
vector<double> params;
force.getExceptionParameters(index, particle1, particle2, params);
force.getExclusionParticles(index, particle1, particle2);
particles.resize(2);
particles[0] = particle1;
particles[1] = particle2;
}
bool areGroupsIdentical(int group1, int group2) {
int particle1, particle2;
vector<double> params1;
vector<double> params2;
force.getExceptionParameters(group1, particle1, particle2, params1);
force.getExceptionParameters(group2, particle1, particle2, params2);
for (int i = 0; i < params1.size(); i++)
if (params1[i] != params2[i])
return false;
return true;
}
private:
......@@ -889,10 +880,6 @@ OpenCLCalcCustomNonbondedForceKernel::~OpenCLCalcCustomNonbondedForceKernel() {
delete params;
if (globals != NULL)
delete globals;
if (exceptionParams != NULL)
delete exceptionParams;
if (exceptionIndices != NULL)
delete exceptionIndices;
if (tabulatedFunctionParams != NULL)
delete tabulatedFunctionParams;
for (int i = 0; i < (int) tabulatedFunctions.size(); i++)
......@@ -900,28 +887,13 @@ OpenCLCalcCustomNonbondedForceKernel::~OpenCLCalcCustomNonbondedForceKernel() {
}
void OpenCLCalcCustomNonbondedForceKernel::initialize(const System& system, const CustomNonbondedForce& force) {
if (force.getNumParameters() > 4)
if (force.getNumPerParticleParameters() > 4)
throw OpenMMException("OpenCLPlatform only supports four per-atom parameters for custom nonbonded forces");
int forceIndex;
for (forceIndex = 0; forceIndex < system.getNumForces() && &system.getForce(forceIndex) != &force; ++forceIndex)
;
string prefix = "custom"+intToString(forceIndex)+"_";
// Identify which exceptions are actual interactions.
vector<pair<int, int> > exclusions;
vector<int> exceptions;
{
vector<double> parameters;
for (int i = 0; i < force.getNumExceptions(); i++) {
int particle1, particle2;
force.getExceptionParameters(i, particle1, particle2, parameters);
exclusions.push_back(pair<int, int>(particle1, particle2));
if (parameters.size() > 0)
exceptions.push_back(i);
}
}
// Record parameters and exclusions.
int numParticles = force.getNumParticles();
......@@ -946,9 +918,11 @@ void OpenCLCalcCustomNonbondedForceKernel::initialize(const System& system, cons
paramVec[i].w = (cl_float) parameters[3];
exclusionList[i].push_back(i);
}
for (int i = 0; i < (int) exclusions.size(); i++) {
exclusionList[exclusions[i].first].push_back(exclusions[i].second);
exclusionList[exclusions[i].second].push_back(exclusions[i].first);
for (int i = 0; i < force.getNumExclusions(); i++) {
int particle1, particle2;
force.getExclusionParticles(i, particle1, particle2);
exclusionList[particle1].push_back(particle2);
exclusionList[particle2].push_back(particle1);
}
params->upload(paramVec);
......@@ -1025,11 +999,8 @@ void OpenCLCalcCustomNonbondedForceKernel::initialize(const System& system, cons
// Record information for the expressions.
vector<string> paramNames;
vector<string> combiningRules;
for (int i = 0; i < force.getNumParameters(); i++) {
paramNames.push_back(force.getParameterName(i));
combiningRules.push_back(force.getParameterCombiningRule(i));
}
for (int i = 0; i < force.getNumPerParticleParameters(); i++)
paramNames.push_back(force.getPerParticleParameterName(i));
globalParamNames.resize(force.getNumGlobalParameters());
globalParamValues.resize(force.getNumGlobalParameters());
for (int i = 0; i < force.getNumGlobalParameters(); i++) {
......@@ -1048,33 +1019,21 @@ void OpenCLCalcCustomNonbondedForceKernel::initialize(const System& system, cons
// Create the kernels.
map<string, string> paramVariables;
map<string, string> forceVariables;
map<string, string> exceptionVariables;
forceVariables["r"] = "r";
exceptionVariables["r"] = "r";
map<string, string> variables;
variables["r"] = "r";
string suffixes[] = {".x", ".y", ".z", ".w"};
for (int i = 0; i < force.getNumParameters(); i++) {
const string& name = force.getParameterName(i);
paramVariables[name+"1"] = prefix+"params1"+suffixes[i];
paramVariables[name+"2"] = prefix+"params2"+suffixes[i];
forceVariables[name] = prefix+name;
exceptionVariables[name] = "exceptionParams"+suffixes[i];
for (int i = 0; i < force.getNumPerParticleParameters(); i++) {
const string& name = force.getPerParticleParameterName(i);
variables[name+"1"] = prefix+"params1"+suffixes[i];
variables[name+"2"] = prefix+"params2"+suffixes[i];
}
for (int i = 0; i < force.getNumGlobalParameters(); i++) {
const string& name = force.getGlobalParameterName(i);
string value = "globals["+intToString(i)+"]";
paramVariables[name] = prefix+value;
forceVariables[name] = prefix+value;
exceptionVariables[name] = value;
}
map<string, Lepton::ParsedExpression> paramExpressions;
for (int i = 0; i < force.getNumParameters(); i++) {
paramExpressions["float "+prefix+force.getParameterName(i)+" = " ] = Lepton::Parser::parse(force.getParameterCombiningRule(i)).optimize();
variables[name] = prefix+value;
}
stringstream compute;
compute << OpenCLExpressionUtilities::createExpressions(paramExpressions, paramVariables, functionDefinitions, prefix+"param_temp", prefix+"functionParams");
compute << OpenCLExpressionUtilities::createExpressions(forceExpressions, forceVariables, functionDefinitions, prefix+"force_temp", prefix+"functionParams");
compute << OpenCLExpressionUtilities::createExpressions(forceExpressions, variables, functionDefinitions, prefix+"temp", prefix+"functionParams");
map<string, string> replacements;
replacements["COMPUTE_FORCE"] = compute.str();
string source = cl.loadSourceFromFile("customNonbonded.cl", replacements);
......@@ -1084,53 +1043,7 @@ void OpenCLCalcCustomNonbondedForceKernel::initialize(const System& system, cons
globals->upload(globalParamValues);
cl.getNonbondedUtilities().addArgument(OpenCLNonbondedUtilities::ParameterInfo(prefix+"globals", "float", sizeof(cl_float), globals->getDeviceBuffer()));
}
map<string, Lepton::ParsedExpression> exceptionExpressions;
stringstream computeExceptions;
exceptionExpressions["energy += "] = energyExpression;
exceptionExpressions["dEdR = "] = forceExpression;
computeExceptions << OpenCLExpressionUtilities::createExpressions(exceptionExpressions, exceptionVariables, functionDefinitions, "temp", prefix+"functionParams");
replacements["COMPUTE_FORCE"] = computeExceptions.str();
replacements["EXTRA_ARGUMENTS"] = extraArguments;
map<string, string> defines;
defines["CUTOFF_SQUARED"] = doubleToString(force.getCutoffDistance()*force.getCutoffDistance());
Vec3 boxVectors[3];
system.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
defines["PERIODIC_BOX_SIZE_X"] = doubleToString(boxVectors[0][0]);
defines["PERIODIC_BOX_SIZE_Y"] = doubleToString(boxVectors[1][1]);
defines["PERIODIC_BOX_SIZE_Z"] = doubleToString(boxVectors[2][2]);
cl::Program program = cl.createProgram(cl.loadSourceFromFile("customNonbondedExceptions.cl", replacements), defines);
exceptionsKernel = cl::Kernel(program, "computeCustomNonbondedExceptions");
// Initialize exception parameters.
int numExceptions = exceptions.size();
int maxBuffers = cl.getNonbondedUtilities().getNumForceBuffers();
if (numExceptions > 0) {
exceptionParams = new OpenCLArray<mm_float4>(cl, numExceptions, "customExceptionParams");
exceptionIndices = new OpenCLArray<mm_int4>(cl, numExceptions, "customExceptionIndices");
vector<mm_float4> exceptionParamsVector(numExceptions);
vector<mm_int4> exceptionIndicesVector(numExceptions);
vector<int> forceBufferCounter(system.getNumParticles(), 0);
for (int i = 0; i < numExceptions; i++) {
int particle1, particle2;
vector<double> parameters;
force.getExceptionParameters(exceptions[i], particle1, particle2, parameters);
if (parameters.size() > 0)
exceptionParamsVector[i].x = (cl_float) parameters[0];
if (parameters.size() > 1)
exceptionParamsVector[i].y = (cl_float) parameters[1];
if (parameters.size() > 2)
exceptionParamsVector[i].z = (cl_float) parameters[2];
if (parameters.size() > 3)
exceptionParamsVector[i].w = (cl_float) parameters[3];
exceptionIndicesVector[i] = (mm_int4) {particle1, particle2, forceBufferCounter[particle1]++, forceBufferCounter[particle2]++};
}
exceptionParams->upload(exceptionParamsVector);
exceptionIndices->upload(exceptionIndicesVector);
for (int i = 0; i < forceBufferCounter.size(); i++)
maxBuffers = max(maxBuffers, forceBufferCounter[i]);
}
cl.addForce(new OpenCLCustomNonbondedForceInfo(maxBuffers, force));
cl.addForce(new OpenCLCustomNonbondedForceInfo(cl.getNonbondedUtilities().getNumForceBuffers(), force));
delete fp;
}
......@@ -1146,21 +1059,6 @@ void OpenCLCalcCustomNonbondedForceKernel::executeForces(ContextImpl& context) {
if (changed)
globals->upload(globalParamValues);
}
if (exceptionParams != NULL) {
if (!hasInitializedKernel) {
hasInitializedKernel = true;
exceptionsKernel.setArg<cl_int>(0, cl.getPaddedNumAtoms());
exceptionsKernel.setArg<cl_int>(1, exceptionParams->getSize());
exceptionsKernel.setArg<cl::Buffer>(2, cl.getForceBuffers().getDeviceBuffer());
exceptionsKernel.setArg<cl::Buffer>(3, cl.getEnergyBuffer().getDeviceBuffer());
exceptionsKernel.setArg<cl::Buffer>(4, cl.getPosq().getDeviceBuffer());
exceptionsKernel.setArg<cl::Buffer>(5, exceptionParams->getDeviceBuffer());
exceptionsKernel.setArg<cl::Buffer>(6, exceptionIndices->getDeviceBuffer());
if (globals != NULL)
exceptionsKernel.setArg<cl::Buffer>(7, globals->getDeviceBuffer());
}
cl.executeKernel(exceptionsKernel, exceptionIndices->getSize());
}
}
double OpenCLCalcCustomNonbondedForceKernel::executeEnergy(ContextImpl& context) {
......
platforms/opencl/src/OpenCLKernels.h
View file @ 70015b75
......@@ -391,7 +391,7 @@ private:
class OpenCLCalcCustomNonbondedForceKernel : public CalcCustomNonbondedForceKernel {
public:
OpenCLCalcCustomNonbondedForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, System& system) : CalcCustomNonbondedForceKernel(name, platform),
hasInitializedKernel(false), cl(cl), params(NULL), globals(NULL), exceptionParams(NULL), exceptionIndices(NULL), tabulatedFunctionParams(NULL), system(system) {
hasInitializedKernel(false), cl(cl), params(NULL), globals(NULL), tabulatedFunctionParams(NULL), system(system) {
}
~OpenCLCalcCustomNonbondedForceKernel();
/**
......@@ -419,10 +419,7 @@ private:
OpenCLContext& cl;
OpenCLArray<mm_float4>* params;
OpenCLArray<cl_float>* globals;
OpenCLArray<mm_float4>* exceptionParams;
OpenCLArray<mm_int4>* exceptionIndices;
OpenCLArray<mm_float4>* tabulatedFunctionParams;
cl::Kernel exceptionsKernel;
std::vector<std::string> globalParamNames;
std::vector<cl_float> globalParamValues;
std::vector<OpenCLArray<mm_float4>*> tabulatedFunctions;
......
platforms/opencl/tests/TestOpenCLCustomNonbondedForce.cpp
View file @ 70015b75
......@@ -79,9 +79,9 @@ void testParameters() {
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomNonbondedForce* forceField = new CustomNonbondedForce("scale*a*(r*b)^3");
forceField->addParameter("a", "a1*a2");
forceField->addParameter("b", "c+b1+b2");
CustomNonbondedForce* forceField = new CustomNonbondedForce("scale*a*(r*b)^3; a=a1*a2; b=c+b1+b2");
forceField->addPerParticleParameter("a");
forceField->addPerParticleParameter("b");
forceField->addGlobalParameter("scale", 3.0);
forceField->addGlobalParameter("c", -1.0);
vector<double> params(2);
......@@ -115,12 +115,12 @@ void testParameters() {
ASSERT_EQUAL_TOL(1.5*3.0*(12*12*12), state.getPotentialEnergy(), TOL);
}
void testExceptions() {
void testExclusions() {
OpenCLPlatform platform;
System system;
VerletIntegrator integrator(0.01);
CustomNonbondedForce* nonbonded = new CustomNonbondedForce("a*r");
nonbonded->addParameter("a", "a1+a2");
CustomNonbondedForce* nonbonded = new CustomNonbondedForce("a*r; a=a1+a2");
nonbonded->addPerParticleParameter("a");
vector<double> params(1);
vector<Vec3> positions(4);
for (int i = 0; i < 4; i++) {
......@@ -129,22 +129,21 @@ void testExceptions() {
nonbonded->addParticle(params);
positions[i] = Vec3(i, 0, 0);
}
nonbonded->addException(0, 1, vector<double>());
nonbonded->addException(1, 2, vector<double>());
nonbonded->addException(2, 3, vector<double>());
params[0] = 0.5;
nonbonded->addException(0, 2, params);
nonbonded->addException(1, 3, params);
nonbonded->addExclusion(0, 1);
nonbonded->addExclusion(1, 2);
nonbonded->addExclusion(2, 3);
nonbonded->addExclusion(0, 2);
nonbonded->addExclusion(1, 3);
system.addForce(nonbonded);
Context context(system, integrator, platform);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
ASSERT_EQUAL_VEC(Vec3(0.5+1+4, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0.5, 0, 0), forces[1], TOL);
ASSERT_EQUAL_VEC(Vec3(-0.5, 0, 0), forces[2], TOL);
ASSERT_EQUAL_VEC(Vec3(-(0.5+1+4), 0, 0), forces[3], TOL);
ASSERT_EQUAL_TOL((1+4)*3+0.5*2+0.5*2, state.getPotentialEnergy(), TOL);
ASSERT_EQUAL_VEC(Vec3(1+4, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[1], TOL);
ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[2], TOL);
ASSERT_EQUAL_VEC(Vec3(-(1+4), 0, 0), forces[3], TOL);
ASSERT_EQUAL_TOL((1+4)*3, state.getPotentialEnergy(), TOL);
}
void testCutoff() {
......@@ -251,10 +250,10 @@ void testCoulombLennardJones() {
customSystem.addParticle(1.0);
}
NonbondedForce* standardNonbonded = new NonbondedForce();
CustomNonbondedForce* customNonbonded = new CustomNonbondedForce("4*eps*((sigma/r)^12-(sigma/r)^6) + 138.935485*q/r");
customNonbonded->addParameter("q", "q1*q2");
customNonbonded->addParameter("sigma", "0.5*(sigma1+sigma2)");
customNonbonded->addParameter("eps", "sqrt(eps1*eps2)");
CustomNonbondedForce* customNonbonded = new CustomNonbondedForce("4*eps*((sigma/r)^12-(sigma/r)^6)+138.935485*q/r; q=q1*q2; sigma=0.5*(sigma1+sigma2); eps=sqrt(eps1*eps2)");
customNonbonded->addPerParticleParameter("q");
customNonbonded->addPerParticleParameter("sigma");
customNonbonded->addPerParticleParameter("eps");
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
init_gen_rand(0);
......@@ -266,7 +265,8 @@ void testCoulombLennardJones() {
params[1] = 0.2;
params[2] = 0.1;
customNonbonded->addParticle(params);
standardNonbonded->addParticle(1.0, 0.1, 0.1);
standardNonbonded->addParticle(-1.0, 0.1, 0.1);
params[0] = -1.0;
params[1] = 0.1;
customNonbonded->addParticle(params);
}
......@@ -276,7 +276,8 @@ void testCoulombLennardJones() {
params[1] = 0.2;
params[2] = 0.2;
customNonbonded->addParticle(params);
standardNonbonded->addParticle(1.0, 0.1, 0.2);
standardNonbonded->addParticle(-1.0, 0.1, 0.2);
params[0] = -1.0;
params[1] = 0.1;
customNonbonded->addParticle(params);
}
......@@ -285,7 +286,7 @@ void testCoulombLennardJones() {
velocities[2*i] = Vec3(genrand_real2(), genrand_real2(), genrand_real2());
velocities[2*i+1] = Vec3(genrand_real2(), genrand_real2(), genrand_real2());
standardNonbonded->addException(2*i, 2*i+1, 0.0, 1.0, 0.0);
customNonbonded->addException(2*i, 2*i+1, vector<double>());
customNonbonded->addExclusion(2*i, 2*i+1);
}
standardNonbonded->setNonbondedMethod(NonbondedForce::NoCutoff);
customNonbonded->setNonbondedMethod(CustomNonbondedForce::NoCutoff);
......@@ -311,7 +312,7 @@ int main() {
try {
testSimpleExpression();
testParameters();
testExceptions();
testExclusions();
testCutoff();
testPeriodic();
testTabulatedFunction(true);
......
platforms/reference/src/ReferenceKernels.cpp
View file @ 70015b75
......@@ -668,42 +668,33 @@ public:
ReferenceCalcCustomNonbondedForceKernel::~ReferenceCalcCustomNonbondedForceKernel() {
disposeRealArray(particleParamArray, numParticles);
disposeIntArray(exclusionArray, numParticles);
disposeIntArray(exceptionIndexArray, numExceptions);
disposeRealArray(exceptionParamArray, numExceptions);
if (neighborList != NULL)
delete neighborList;
}
void ReferenceCalcCustomNonbondedForceKernel::initialize(const System& system, const CustomNonbondedForce& force) {
// Identify which exceptions are actual interactions.
// Record the exclusions.
numParticles = force.getNumParticles();
exclusions.resize(numParticles);
vector<int> exceptions;
vector<double> parameters;
for (int i = 0; i < force.getNumExceptions(); i++) {
for (int i = 0; i < force.getNumExclusions(); i++) {
int particle1, particle2;
force.getExceptionParameters(i, particle1, particle2, parameters);
force.getExclusionParticles(i, particle1, particle2);
exclusions[particle1].insert(particle2);
exclusions[particle2].insert(particle1);
if (parameters.size() > 0)
exceptions.push_back(i);
}
// Build the arrays.
numExceptions = exceptions.size();
int numParameters = force.getNumParameters();
exceptionIndexArray = allocateIntArray(numExceptions, 2);
exceptionParamArray = allocateRealArray(numExceptions, numParameters);
int numParameters = force.getNumPerParticleParameters();
particleParamArray = allocateRealArray(numParticles, numParameters);
for (int i = 0; i < numParticles; ++i) {
vector<double> parameters;
force.getParticleParameters(i, parameters);
for (int j = 0; j < numParameters; j++)
particleParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
}
this->exclusions = exclusions;
exclusionArray = new int*[numParticles];
for (int i = 0; i < numParticles; ++i) {
exclusionArray[i] = new int[exclusions[i].size()+1];
......@@ -712,14 +703,6 @@ void ReferenceCalcCustomNonbondedForceKernel::initialize(const System& system, c
for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter)
exclusionArray[i][++index] = *iter;
}
for (int i = 0; i < numExceptions; ++i) {
int particle1, particle2;
force.getExceptionParameters(exceptions[i], particle1, particle2, parameters);
exceptionIndexArray[i][0] = particle1;
exceptionIndexArray[i][1] = particle2;
for (int j = 0; j < numParameters; j++)
exceptionParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
}
nonbondedMethod = CalcCustomNonbondedForceKernel::NonbondedMethod(force.getNonbondedMethod());
nonbondedCutoff = (RealOpenMM) force.getCutoffDistance();
Vec3 boxVectors[3];
......@@ -749,10 +732,8 @@ void ReferenceCalcCustomNonbondedForceKernel::initialize(const System& system, c
Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction(), functions).optimize();
energyExpression = expression.createProgram();
forceExpression = expression.differentiate("r").optimize().createProgram();
for (int i = 0; i < numParameters; i++) {
parameterNames.push_back(force.getParameterName(i));
combiningRules.push_back(Lepton::Parser::parse(force.getParameterCombiningRule(i), functions).optimize().createProgram());
}
for (int i = 0; i < numParameters; i++)
parameterNames.push_back(force.getPerParticleParameterName(i));
for (int i = 0; i < force.getNumGlobalParameters(); i++)
globalParameterNames.push_back(force.getGlobalParameterName(i));
......@@ -765,7 +746,7 @@ void ReferenceCalcCustomNonbondedForceKernel::initialize(const System& system, c
void ReferenceCalcCustomNonbondedForceKernel::executeForces(ContextImpl& context) {
RealOpenMM** posData = extractPositions(context);
RealOpenMM** forceData = extractForces(context);
ReferenceCustomNonbondedIxn ixn(energyExpression, forceExpression, parameterNames, combiningRules);
ReferenceCustomNonbondedIxn ixn(energyExpression, forceExpression, parameterNames);
bool periodic = (nonbondedMethod == CutoffPeriodic);
if (nonbondedMethod != NoCutoff) {
computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, periodic ? periodicBoxSize : NULL, nonbondedCutoff, 0.0);
......@@ -777,14 +758,13 @@ void ReferenceCalcCustomNonbondedForceKernel::executeForces(ContextImpl& context
for (int i = 0; i < (int) globalParameterNames.size(); i++)
globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
ixn.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, globalParameters, forceData, 0, 0);
ixn.calculateExceptionIxn(numExceptions, exceptionIndexArray, posData, exceptionParamArray, globalParameters, forceData, 0, 0);
}
double ReferenceCalcCustomNonbondedForceKernel::executeEnergy(ContextImpl& context) {
RealOpenMM** posData = extractPositions(context);
RealOpenMM** forceData = allocateRealArray(numParticles, 3);
RealOpenMM energy = 0;
ReferenceCustomNonbondedIxn ixn(energyExpression, forceExpression, parameterNames, combiningRules);
ReferenceCustomNonbondedIxn ixn(energyExpression, forceExpression, parameterNames);
bool periodic = (nonbondedMethod == CutoffPeriodic);
if (nonbondedMethod != NoCutoff) {
computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, periodic ? periodicBoxSize : NULL, nonbondedCutoff, 0.0);
......@@ -796,7 +776,6 @@ double ReferenceCalcCustomNonbondedForceKernel::executeEnergy(ContextImpl& conte
for (int i = 0; i < (int) globalParameterNames.size(); i++)
globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
ixn.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, globalParameters, forceData, 0, &energy);
ixn.calculateExceptionIxn(numExceptions, exceptionIndexArray, posData, exceptionParamArray, globalParameters, forceData, 0, &energy);
disposeRealArray(forceData, numParticles);
return energy;
}
......
platforms/reference/src/ReferenceKernels.h
View file @ 70015b75
......@@ -398,14 +398,13 @@ public:
*/
double executeEnergy(ContextImpl& context);
private:
int numParticles, numExceptions;
int **exclusionArray, **exceptionIndexArray;
RealOpenMM **particleParamArray, **exceptionParamArray;
int numParticles;
int **exclusionArray;
RealOpenMM **particleParamArray;
RealOpenMM nonbondedCutoff, periodicBoxSize[3];
std::vector<std::set<int> > exclusions;
Lepton::ExpressionProgram energyExpression, forceExpression;
std::vector<std::string> parameterNames, globalParameterNames;
std::vector<Lepton::ExpressionProgram> combiningRules;
NonbondedMethod nonbondedMethod;
NeighborList* neighborList;
class TabulatedFunction;
......
platforms/reference/src/SimTKReference/ReferenceCustomNonbondedIxn.cpp
View file @ 70015b75
......@@ -43,9 +43,8 @@ using std::vector;
--------------------------------------------------------------------------------------- */
ReferenceCustomNonbondedIxn::ReferenceCustomNonbondedIxn(const Lepton::ExpressionProgram& energyExpression,
const Lepton::ExpressionProgram& forceExpression, const vector<string>& parameterNames, const vector<Lepton::ExpressionProgram>& combiningRules) :
cutoff(false), periodic(false), energyExpression(energyExpression), forceExpression(forceExpression),
paramNames(parameterNames), combiningRules(combiningRules) {
const Lepton::ExpressionProgram& forceExpression, const vector<string>& parameterNames) :
cutoff(false), periodic(false), energyExpression(energyExpression), forceExpression(forceExpression), paramNames(parameterNames) {
// ---------------------------------------------------------------------------------------
......@@ -151,21 +150,15 @@ int ReferenceCustomNonbondedIxn::calculatePairIxn( int numberOfAtoms, RealOpenMM
RealOpenMM* fixedParameters, map<string, double> globalParameters, RealOpenMM** forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const {
map<string, double> variablesForParams = globalParameters;
map<string, double> variablesForForce = globalParameters;
map<string, double> variables = globalParameters;
if (cutoff) {
for (int i = 0; i < (int) neighborList->size(); i++) {
OpenMM::AtomPair pair = (*neighborList)[i];
// Apply the combining rules to compute the force field parameters.
for (int j = 0; j < (int) combiningRules.size(); j++) {
variablesForParams[particleParamNames[j*2]] = atomParameters[pair.first][j];
variablesForParams[particleParamNames[j*2+1]] = atomParameters[pair.second][j];
for (int j = 0; j < (int) paramNames.size(); j++) {
variables[particleParamNames[j*2]] = atomParameters[pair.first][j];
variables[particleParamNames[j*2+1]] = atomParameters[pair.second][j];
}
for (int j = 0; j < (int) combiningRules.size(); j++)
variablesForForce[paramNames[j]] = combiningRules[j].evaluate(variablesForParams);
calculateOneIxn(pair.first, pair.second, atomCoordinates, variablesForForce, forces, energyByAtom, totalEnergy);
calculateOneIxn(pair.first, pair.second, atomCoordinates, variables, forces, energyByAtom, totalEnergy);
}
}
else {
......@@ -189,16 +182,11 @@ int ReferenceCustomNonbondedIxn::calculatePairIxn( int numberOfAtoms, RealOpenMM
for( int jj = ii+1; jj < numberOfAtoms; jj++ ){
if( exclusionIndices[jj] != ii ){
// Apply the combining rules to compute the force field parameters.
for (int j = 0; j < (int) combiningRules.size(); j++) {
variablesForParams[particleParamNames[j*2]] = atomParameters[ii][j];
variablesForParams[particleParamNames[j*2+1]] = atomParameters[jj][j];
for (int j = 0; j < (int) paramNames.size(); j++) {
variables[particleParamNames[j*2]] = atomParameters[ii][j];
variables[particleParamNames[j*2+1]] = atomParameters[jj][j];
}
for (int j = 0; j < (int) combiningRules.size(); j++)
variablesForForce[paramNames[j]] = combiningRules[j].evaluate(variablesForParams);
calculateOneIxn(ii, jj, atomCoordinates, variablesForForce, forces, energyByAtom, totalEnergy);
calculateOneIxn(ii, jj, atomCoordinates, variables, forces, energyByAtom, totalEnergy);
}
}
}
......@@ -208,32 +196,6 @@ int ReferenceCustomNonbondedIxn::calculatePairIxn( int numberOfAtoms, RealOpenMM
return ReferenceForce::DefaultReturn;
}
/**---------------------------------------------------------------------------------------
Calculate custom pair ixn for exceptions
@param numberOfExceptions number of exceptions
@param atomIndices indices of atoms participating in exception ixn: atomIndices[exceptionIndex][indices]
@param atomCoordinates atom coordinates: atomCoordinates[atomIndex][3]
@param parameters parameters: parameters[exceptionIndex][*]; contents of array
depend on ixn
@param globalParameters the values of global parameters
@param forces force array (forces added to current values): forces[atomIndex][3]
@param totalEnergy totalEnergy: sum over { energies[atomIndex] }
--------------------------------------------------------------------------------------- */
void ReferenceCustomNonbondedIxn::calculateExceptionIxn( int numberOfExceptions, int** atomIndices, RealOpenMM** atomCoordinates,
RealOpenMM** parameters, map<string, double> globalParameters, RealOpenMM** forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const {
map<string, double> variables = globalParameters;
for (int i = 0; i < numberOfExceptions; i++) {
for (int j = 0; j < (int) combiningRules.size(); j++)
variables[paramNames[j]] = parameters[i][j];
calculateOneIxn(atomIndices[i][0], atomIndices[i][1], atomCoordinates, variables, forces, energyByAtom, totalEnergy);
}
}
/**---------------------------------------------------------------------------------------
......
platforms/reference/src/SimTKReference/ReferenceCustomNonbondedIxn.h
View file @ 70015b75
......@@ -44,7 +44,6 @@ class ReferenceCustomNonbondedIxn {
RealOpenMM cutoffDistance;
Lepton::ExpressionProgram energyExpression;
Lepton::ExpressionProgram forceExpression;
std::vector<Lepton::ExpressionProgram> combiningRules;
std::vector<std::string> paramNames;
std::vector<std::string> particleParamNames;
......@@ -76,7 +75,7 @@ class ReferenceCustomNonbondedIxn {
--------------------------------------------------------------------------------------- */
ReferenceCustomNonbondedIxn(const Lepton::ExpressionProgram& energyExpression, const Lepton::ExpressionProgram& forceExpression,
const std::vector<std::string>& parameterNames, const std::vector<Lepton::ExpressionProgram>& combiningRules);
const std::vector<std::string>& parameterNames);
/**---------------------------------------------------------------------------------------
......@@ -139,26 +138,6 @@ class ReferenceCustomNonbondedIxn {
RealOpenMM* fixedParameters, std::map<std::string, double> globalParameters,
RealOpenMM** forces, RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const;
/**---------------------------------------------------------------------------------------
Calculate custom pair ixn for exceptions
@param numberOfExceptions number of exceptions
@param atomIndices indices of atoms participating in exception ixn: atomIndices[exceptionIndex][indices]
@param atomCoordinates atom coordinates: atomCoordinates[atomIndex][3]
@param parameters parameters: parameters[excptionIndex][*]; contents of array
depend on ixn
@param globalParameters the values of global parameters
@param forces force array (forces added to current values): forces[atomIndex][3]
@param energyByAtom atom energy
@param totalEnergy total energy
--------------------------------------------------------------------------------------- */
void calculateExceptionIxn( int numberOfExceptions, int** atomIndices, RealOpenMM** atomCoordinates,
RealOpenMM** parameters, std::map<std::string, double> globalParameters,
RealOpenMM** forces, RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const;
// ---------------------------------------------------------------------------------------
};
......
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