Commit 8a1641ff authored by Peter Eastman's avatar Peter Eastman
Browse files

Custom forces in OpenCL platform support an unlimited number of parameters

parent b321a758
......@@ -288,10 +288,8 @@ OpenCLCalcCustomBondForceKernel::~OpenCLCalcCustomBondForceKernel() {
}
void OpenCLCalcCustomBondForceKernel::initialize(const System& system, const CustomBondForce& force) {
if (force.getNumPerBondParameters() > 4)
throw OpenMMException("OpenCLPlatform only supports four per-bond parameters for custom bonded forces");
numBonds = force.getNumBonds();
params = new OpenCLArray<mm_float4>(cl, numBonds, "customBondParams");
params = new OpenCLParameterSet(cl, force.getNumPerBondParameters(), numBonds, "customBondParams");
indices = new OpenCLArray<mm_int4>(cl, numBonds, "customBondIndices");
string extraArguments;
if (force.getNumGlobalParameters() > 0) {
......@@ -299,23 +297,18 @@ void OpenCLCalcCustomBondForceKernel::initialize(const System& system, const Cus
extraArguments += ", __constant float* globals";
}
vector<int> forceBufferCounter(system.getNumParticles(), 0);
vector<mm_float4> paramVector(numBonds);
vector<vector<cl_float> > paramVector(numBonds);
vector<mm_int4> indicesVector(numBonds);
for (int i = 0; i < numBonds; i++) {
int particle1, particle2;
vector<double> parameters;
force.getBondParameters(i, particle1, particle2, parameters);
if (parameters.size() > 0)
paramVector[i].x = (cl_float) parameters[0];
if (parameters.size() > 1)
paramVector[i].y = (cl_float) parameters[1];
if (parameters.size() > 2)
paramVector[i].z = (cl_float) parameters[2];
if (parameters.size() > 3)
paramVector[i].w = (cl_float) parameters[3];
paramVector[i].resize(parameters.size());
for (int j = 0; j < parameters.size(); j++)
paramVector[i][j] = (cl_float) parameters[j];
indicesVector[i] = (mm_int4) {particle1, particle2, forceBufferCounter[particle1]++, forceBufferCounter[particle2]++};
}
params->upload(paramVector);
params->setParameterValues(paramVector);
indices->upload(indicesVector);
int maxBuffers = 1;
for (int i = 0; i < forceBufferCounter.size(); i++)
......@@ -345,20 +338,24 @@ void OpenCLCalcCustomBondForceKernel::initialize(const System& system, const Cus
map<string, string> variables;
variables["r"] = "r";
string suffixes[] = {".x", ".y", ".z", ".w"};
for (int i = 0; i < force.getNumPerBondParameters(); i++) {
const string& name = force.getPerBondParameterName(i);
variables[name] = "bondParams"+suffixes[i];
variables[name] = "bondParams"+params->getParameterSuffix(i);
}
for (int i = 0; i < force.getNumGlobalParameters(); i++) {
const string& name = force.getGlobalParameterName(i);
string value = "globals["+intToString(i)+"]";
variables[name] = value;
}
map<string, string> replacements;
stringstream compute;
for (int i = 0; i < (int) params->getBuffers().size(); i++) {
const OpenCLNonbondedUtilities::ParameterInfo& buffer = params->getBuffers()[i];
extraArguments += ", __global "+buffer.getType()+"* "+buffer.getName();
compute<<buffer.getType()<<" bondParams"<<(i+1)<<" = "<<buffer.getName()<<"[index];\n";
}
vector<pair<string, string> > functions;
compute << OpenCLExpressionUtilities::createExpressions(expressions, variables, functions, "temp", "");
map<string, string> replacements;
replacements["COMPUTE_FORCE"] = compute.str();
replacements["EXTRA_ARGUMENTS"] = extraArguments;
cl::Program program = cl.createProgram(cl.loadSourceFromFile("customBondForce.cl", replacements));
......@@ -384,10 +381,14 @@ void OpenCLCalcCustomBondForceKernel::executeForces(ContextImpl& context) {
kernel.setArg<cl::Buffer>(2, cl.getForceBuffers().getDeviceBuffer());
kernel.setArg<cl::Buffer>(3, cl.getEnergyBuffer().getDeviceBuffer());
kernel.setArg<cl::Buffer>(4, cl.getPosq().getDeviceBuffer());
kernel.setArg<cl::Buffer>(5, params->getDeviceBuffer());
kernel.setArg<cl::Buffer>(6, indices->getDeviceBuffer());
kernel.setArg<cl::Buffer>(5, indices->getDeviceBuffer());
int nextIndex = 6;
if (globals != NULL)
kernel.setArg<cl::Buffer>(7, globals->getDeviceBuffer());
kernel.setArg<cl::Buffer>(nextIndex++, globals->getDeviceBuffer());
for (int i = 0; i < (int) params->getBuffers().size(); i++) {
const OpenCLNonbondedUtilities::ParameterInfo& buffer = params->getBuffers()[i];
kernel.setArg<cl::Buffer>(nextIndex++, buffer.getBuffer());
}
}
cl.executeKernel(kernel, numBonds);
}
......@@ -887,8 +888,6 @@ OpenCLCalcCustomNonbondedForceKernel::~OpenCLCalcCustomNonbondedForceKernel() {
}
void OpenCLCalcCustomNonbondedForceKernel::initialize(const System& system, const CustomNonbondedForce& force) {
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)
;
......@@ -898,24 +897,19 @@ void OpenCLCalcCustomNonbondedForceKernel::initialize(const System& system, cons
int numParticles = force.getNumParticles();
string extraArguments;
params = new OpenCLArray<mm_float4>(cl, numParticles, "customNonbondedParameters");
params = new OpenCLParameterSet(cl, force.getNumPerParticleParameters(), numParticles, "customNonbondedParameters");
if (force.getNumGlobalParameters() > 0) {
globals = new OpenCLArray<cl_float>(cl, force.getNumGlobalParameters(), "customNonbondedGlobals", false, CL_MEM_READ_ONLY);
extraArguments += ", __constant float* globals";
}
vector<mm_float4> paramVec(numParticles);
vector<vector<cl_float> > paramVector(numParticles);
vector<vector<int> > exclusionList(numParticles);
for (int i = 0; i < numParticles; i++) {
vector<double> parameters;
force.getParticleParameters(i, parameters);
if (parameters.size() > 0)
paramVec[i].x = (cl_float) parameters[0];
if (parameters.size() > 1)
paramVec[i].y = (cl_float) parameters[1];
if (parameters.size() > 2)
paramVec[i].z = (cl_float) parameters[2];
if (parameters.size() > 3)
paramVec[i].w = (cl_float) parameters[3];
paramVector[i].resize(parameters.size());
for (int j = 0; j < parameters.size(); j++)
paramVector[i][j] = (cl_float) parameters[j];
exclusionList[i].push_back(i);
}
for (int i = 0; i < force.getNumExclusions(); i++) {
......@@ -924,7 +918,7 @@ void OpenCLCalcCustomNonbondedForceKernel::initialize(const System& system, cons
exclusionList[particle1].push_back(particle2);
exclusionList[particle2].push_back(particle1);
}
params->upload(paramVec);
params->setParameterValues(paramVector);
// This class serves as a placeholder for custom functions in expressions.
......@@ -1021,11 +1015,10 @@ void OpenCLCalcCustomNonbondedForceKernel::initialize(const System& system, cons
map<string, string> variables;
variables["r"] = "r";
string suffixes[] = {".x", ".y", ".z", ".w"};
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];
variables[name+"1"] = prefix+"params"+params->getParameterSuffix(i, "1");
variables[name+"2"] = prefix+"params"+params->getParameterSuffix(i, "2");
}
for (int i = 0; i < force.getNumGlobalParameters(); i++) {
const string& name = force.getGlobalParameterName(i);
......@@ -1038,7 +1031,10 @@ void OpenCLCalcCustomNonbondedForceKernel::initialize(const System& system, cons
replacements["COMPUTE_FORCE"] = compute.str();
string source = cl.loadSourceFromFile("customNonbonded.cl", replacements);
cl.getNonbondedUtilities().addInteraction(useCutoff, usePeriodic, true, force.getCutoffDistance(), exclusionList, source);
cl.getNonbondedUtilities().addParameter(OpenCLNonbondedUtilities::ParameterInfo(prefix+"params", "float4", sizeof(cl_float4), params->getDeviceBuffer()));
for (int i = 0; i < (int) params->getBuffers().size(); i++) {
const OpenCLNonbondedUtilities::ParameterInfo& buffer = params->getBuffers()[i];
cl.getNonbondedUtilities().addParameter(OpenCLNonbondedUtilities::ParameterInfo(prefix+"params"+intToString(i+1), buffer.getType(), buffer.getSize(), buffer.getBuffer()));
}
if (globals != NULL) {
globals->upload(globalParamValues);
cl.getNonbondedUtilities().addArgument(OpenCLNonbondedUtilities::ParameterInfo(prefix+"globals", "float", sizeof(cl_float), globals->getDeviceBuffer()));
......@@ -1256,31 +1252,24 @@ OpenCLCalcCustomExternalForceKernel::~OpenCLCalcCustomExternalForceKernel() {
}
void OpenCLCalcCustomExternalForceKernel::initialize(const System& system, const CustomExternalForce& force) {
if (force.getNumPerParticleParameters() > 4)
throw OpenMMException("OpenCLPlatform only supports four per-particle parameters for custom external forces");
numParticles = force.getNumParticles();
params = new OpenCLArray<mm_float4>(cl, numParticles, "customExternalParams");
params = new OpenCLParameterSet(cl, force.getNumPerParticleParameters(), numParticles, "customExternalParams");
indices = new OpenCLArray<cl_int>(cl, numParticles, "customExternalIndices");
string extraArguments;
if (force.getNumGlobalParameters() > 0) {
globals = new OpenCLArray<cl_float>(cl, force.getNumGlobalParameters(), "customExternalGlobals", false, CL_MEM_READ_ONLY);
extraArguments += ", __constant float* globals";
}
vector<mm_float4> paramVector(numParticles);
vector<vector<cl_float> > paramVector(numParticles);
vector<cl_int> indicesVector(numParticles);
for (int i = 0; i < numParticles; i++) {
vector<double> parameters;
force.getParticleParameters(i, indicesVector[i], parameters);
if (parameters.size() > 0)
paramVector[i].x = (cl_float) parameters[0];
if (parameters.size() > 1)
paramVector[i].y = (cl_float) parameters[1];
if (parameters.size() > 2)
paramVector[i].z = (cl_float) parameters[2];
if (parameters.size() > 3)
paramVector[i].w = (cl_float) parameters[3];
paramVector[i].resize(parameters.size());
for (int j = 0; j < parameters.size(); j++)
paramVector[i][j] = (cl_float) parameters[j];
}
params->upload(paramVector);
params->setParameterValues(paramVector);
indices->upload(indicesVector);
cl.addForce(new OpenCLCustomExternalForceInfo(force, system.getNumParticles()));
......@@ -1313,20 +1302,24 @@ void OpenCLCalcCustomExternalForceKernel::initialize(const System& system, const
variables["x"] = "pos.x";
variables["y"] = "pos.y";
variables["z"] = "pos.z";
string suffixes[] = {".x", ".y", ".z", ".w"};
for (int i = 0; i < force.getNumPerParticleParameters(); i++) {
const string& name = force.getPerParticleParameterName(i);
variables[name] = "particleParams"+suffixes[i];
variables[name] = "particleParams"+params->getParameterSuffix(i);
}
for (int i = 0; i < force.getNumGlobalParameters(); i++) {
const string& name = force.getGlobalParameterName(i);
string value = "globals["+intToString(i)+"]";
variables[name] = value;
}
map<string, string> replacements;
stringstream compute;
for (int i = 0; i < (int) params->getBuffers().size(); i++) {
const OpenCLNonbondedUtilities::ParameterInfo& buffer = params->getBuffers()[i];
extraArguments += ", __global "+buffer.getType()+"* "+buffer.getName();
compute<<buffer.getType()<<" particleParams"<<(i+1)<<" = "<<buffer.getName()<<"[index];\n";
}
vector<pair<string, string> > functions;
compute << OpenCLExpressionUtilities::createExpressions(expressions, variables, functions, "temp", "");
map<string, string> replacements;
replacements["COMPUTE_FORCE"] = compute.str();
replacements["EXTRA_ARGUMENTS"] = extraArguments;
cl::Program program = cl.createProgram(cl.loadSourceFromFile("customExternalForce.cl", replacements));
......@@ -1351,10 +1344,14 @@ void OpenCLCalcCustomExternalForceKernel::executeForces(ContextImpl& context) {
kernel.setArg<cl::Buffer>(1, cl.getForceBuffers().getDeviceBuffer());
kernel.setArg<cl::Buffer>(2, cl.getEnergyBuffer().getDeviceBuffer());
kernel.setArg<cl::Buffer>(3, cl.getPosq().getDeviceBuffer());
kernel.setArg<cl::Buffer>(4, params->getDeviceBuffer());
kernel.setArg<cl::Buffer>(5, indices->getDeviceBuffer());
kernel.setArg<cl::Buffer>(4, indices->getDeviceBuffer());
int nextIndex = 5;
if (globals != NULL)
kernel.setArg<cl::Buffer>(6, globals->getDeviceBuffer());
kernel.setArg<cl::Buffer>(nextIndex++, globals->getDeviceBuffer());
for (int i = 0; i < (int) params->getBuffers().size(); i++) {
const OpenCLNonbondedUtilities::ParameterInfo& buffer = params->getBuffers()[i];
kernel.setArg<cl::Buffer>(nextIndex++, buffer.getBuffer());
}
}
cl.executeKernel(kernel, numParticles);
}
......
......@@ -30,6 +30,7 @@
#include "OpenCLPlatform.h"
#include "OpenCLArray.h"
#include "OpenCLContext.h"
#include "OpenCLParameterSet.h"
#include "openmm/kernels.h"
#include "openmm/System.h"
......@@ -218,7 +219,7 @@ private:
bool hasInitializedKernel;
OpenCLContext& cl;
System& system;
OpenCLArray<mm_float4>* params;
OpenCLParameterSet* params;
OpenCLArray<mm_int4>* indices;
OpenCLArray<cl_float>* globals;
std::vector<std::string> globalParamNames;
......@@ -417,7 +418,7 @@ public:
private:
bool hasInitializedKernel;
OpenCLContext& cl;
OpenCLArray<mm_float4>* params;
OpenCLParameterSet* params;
OpenCLArray<cl_float>* globals;
OpenCLArray<mm_float4>* tabulatedFunctionParams;
std::vector<std::string> globalParamNames;
......@@ -503,7 +504,7 @@ private:
bool hasInitializedKernel;
OpenCLContext& cl;
System& system;
OpenCLArray<mm_float4>* params;
OpenCLParameterSet* params;
OpenCLArray<cl_int>* indices;
OpenCLArray<cl_float>* globals;
std::vector<std::string> globalParamNames;
......
/* -------------------------------------------------------------------------- *
* 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) 2009 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU Lesser General Public License as published *
* by the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
* -------------------------------------------------------------------------- */
#include "OpenCLParameterSet.h"
#include "openmm/OpenMMException.h"
#include <cmath>
#include <sstream>
using namespace OpenMM;
using namespace std;
OpenCLParameterSet::OpenCLParameterSet(OpenCLContext& context, int numParameters, int numObjects, const string& name) :
context(context), numParameters(numParameters), numObjects(numObjects), name(name) {
int params = numParameters;
int bufferCount = 0;
try {
while (params > 2) {
cl::Buffer* buf = new cl::Buffer(context.getContext(), CL_MEM_READ_WRITE, numObjects*sizeof(mm_float4));
std::stringstream name;
name << "param" << (++bufferCount);
buffers.push_back(OpenCLNonbondedUtilities::ParameterInfo(name.str(), "float4", sizeof(mm_float4), *buf));
params -= 4;
}
if (params > 1) {
cl::Buffer* buf = new cl::Buffer(context.getContext(), CL_MEM_READ_WRITE, numObjects*sizeof(mm_float2));
std::stringstream name;
name << "param" << (++bufferCount);
buffers.push_back(OpenCLNonbondedUtilities::ParameterInfo(name.str(), "float2", sizeof(mm_float2), *buf));
params -= 2;
}
if (params > 0) {
cl::Buffer* buf = new cl::Buffer(context.getContext(), CL_MEM_READ_WRITE, numObjects*sizeof(cl_float));
std::stringstream name;
name << "param" << (++bufferCount);
buffers.push_back(OpenCLNonbondedUtilities::ParameterInfo(name.str(), "float", sizeof(cl_float), *buf));
}
}
catch (cl::Error err) {
stringstream str;
str<<"Error creating parameter set "<<name<<": "<<err.what()<<" ("<<err.err()<<")";
throw OpenMMException(str.str());
}
}
OpenCLParameterSet::~OpenCLParameterSet() {
for (int i = 0; i < (int) buffers.size(); i++)
delete &buffers[i].getBuffer();
}
void OpenCLParameterSet::setParameterValues(const vector<vector<cl_float> >& values) {
try {
int base = 0;
for (int i = 0; i < (int) buffers.size(); i++) {
if (buffers[i].getType() == "float4") {
vector<mm_float4> data(numObjects);
for (int j = 0; j < numObjects; j++)
data[j] = (mm_float4) {values[j][base], values[j][base+1], values[j][base+2], values[j][base+3]};
context.getQueue().enqueueWriteBuffer(buffers[i].getBuffer(), CL_TRUE, 0, numObjects*buffers[i].getSize(), &data[0]);
base += 4;
}
else if (buffers[i].getType() == "float2") {
vector<mm_float2> data(numObjects);
for (int j = 0; j < numObjects; j++)
data[j] = (mm_float2) {values[j][base], values[j][base+1]};
context.getQueue().enqueueWriteBuffer(buffers[i].getBuffer(), CL_TRUE, 0, numObjects*buffers[i].getSize(), &data[0]);
base += 2;
}
else if (buffers[i].getType() == "float") {
vector<cl_float> data(numObjects);
for (int j = 0; j < numObjects; j++)
data[j] = values[j][base];
context.getQueue().enqueueWriteBuffer(buffers[i].getBuffer(), CL_TRUE, 0, numObjects*buffers[i].getSize(), &data[0]);
}
else
throw OpenMMException("Internal error: Unknown buffer type in OpenCLParameterSet");
}
}
catch (cl::Error err) {
stringstream str;
str<<"Error uploading parameter set "<<name<<": "<<err.what()<<" ("<<err.err()<<")";
throw OpenMMException(str.str());
}
}
string OpenCLParameterSet::getParameterSuffix(int index, const std::string& extraSuffix) const {
const string suffixes[] = {".x", ".y", ".z", ".w"};
int buffer = -1;
for (int i = 0; buffer == -1 && i < (int) buffers.size(); i++) {
if (index*sizeof(cl_float) < buffers[i].getSize())
buffer = i;
else
index -= buffers[i].getSize()/sizeof(cl_float);
}
if (buffer == -1)
throw OpenMMException("Internal error: Illegal argument to OpenCLParameterSet::getParameterSuffix() ("+name+")");
stringstream suffix;
suffix << (buffer+1) << extraSuffix;
if (buffers[buffer].getType() != "float")
suffix << suffixes[index];
return suffix.str();
}
#ifndef OPENMM_OPENCLPARAMETERSET_H_
#define OPENMM_OPENCLPARAMETERSET_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) 2009 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU Lesser General Public License as published *
* by the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU Lesser General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
* -------------------------------------------------------------------------- */
#include "OpenCLContext.h"
#include "OpenCLNonbondedUtilities.h"
namespace OpenMM {
class OpenCLNonbondedUtilities;
/**
* This class represents a set of floating point parameter values for a set of objects (particles, bonds, etc.).
* It automatically creates an appropriate set of cl::Buffers to hold the parameter values, based
* on the number of parameters required.
*/
class OpenCLParameterSet {
public:
/**
* Create an OpenCLParameterSet.
*
* @param context the context for which to create the parameter set
* @param numParameters the number of parameters for each object
* @param numObjects the number of objects to store parameter values for
* @param name the name of the parameter set
*/
OpenCLParameterSet(OpenCLContext& context, int numParameters, int numObjects, const std::string& name);
~OpenCLParameterSet();
/**
* Get the number of parameters.
*/
int getNumParameters() const {
return numParameters;
}
/**
* Get the number of objects.
*/
int getNumObjects() const {
return numObjects;
}
/**
* Set the values of all parameters.
*
* @param values values[i][j] contains the value of parameter j for object i
*/
void setParameterValues(const std::vector<std::vector<cl_float> >& values);
/**
* Get a set of OpenCLNonbondedUtilities::ParameterInfo objects which describe the Buffers
* containing the data.
*/
const std::vector<OpenCLNonbondedUtilities::ParameterInfo>& getBuffers() const {
return buffers;
}
/**
* Get a suffix to add to variable names when accessing a certain parameter.
*
* @param index the index of the parameter
* @param extraSuffix an extra suffix to add to the variable name
* @return the suffix to append
*/
std::string getParameterSuffix(int index, const std::string& extraSuffix = "") const;
private:
OpenCLContext& context;
int numParameters;
int numObjects;
std::string name;
std::vector<OpenCLNonbondedUtilities::ParameterInfo> buffers;
};
} // namespace OpenMM
#endif /*OPENMM_OPENCLPARAMETERSET_H_*/
......@@ -3,14 +3,13 @@
*/
__kernel void computeCustomBondForces(int numAtoms, int numBonds, __global float4* forceBuffers, __global float* energyBuffer,
__global float4* posq, __global float4* params, __global int4* indices
__global float4* posq, __global int4* indices
EXTRA_ARGUMENTS) {
float energy = 0.0f;
for (int index = get_global_id(0); index < numBonds; index += get_global_size(0)) {
// Look up the data for this bond.
int4 atoms = indices[index];
float4 bondParams = params[index];
float4 delta = posq[atoms.y]-posq[atoms.x];
// Compute the force.
......
......@@ -3,14 +3,13 @@
*/
__kernel void computeCustomExternalForces(int numTerms, __global float4* forceBuffers, __global float* energyBuffer,
__global float4* posq, __global float4* params, __global int* indices
__global float4* posq, __global int* indices
EXTRA_ARGUMENTS) {
float energy = 0.0f;
for (int index = get_global_id(0); index < numTerms; index += get_global_size(0)) {
// Look up the data for this particle.
int atom = indices[index];
float4 particleParams = params[index];
float4 pos = posq[atom];
// Compute the force.
......
......@@ -81,9 +81,44 @@ void testBonds() {
ASSERT_EQUAL_TOL(0.5*0.8*0.5*0.5 + 0.5*0.7*0.2*0.2, state.getPotentialEnergy(), TOL);
}
void testManyParameters() {
OpenCLPlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomBondForce* forceField = new CustomBondForce("(a+b+c+d+e+f+g+h+i)*r");
forceField->addPerBondParameter("a");
forceField->addPerBondParameter("b");
forceField->addPerBondParameter("c");
forceField->addPerBondParameter("d");
forceField->addPerBondParameter("e");
forceField->addPerBondParameter("f");
forceField->addPerBondParameter("g");
forceField->addPerBondParameter("h");
forceField->addPerBondParameter("i");
vector<double> parameters(forceField->getNumPerBondParameters());
for (int i = 0; i < parameters.size(); i++)
parameters[i] = i;
forceField->addBond(0, 1, parameters);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(0, 2.5, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
double f = 1+2+3+4+5+6+7+8;
ASSERT_EQUAL_VEC(Vec3(0, f, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(0, -f, 0), forces[1], TOL);
ASSERT_EQUAL_TOL(f*2.5, state.getPotentialEnergy(), TOL);
}
int main() {
try {
testBonds();
testManyParameters();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
......
......@@ -81,9 +81,41 @@ void testForce() {
ASSERT_EQUAL_TOL(0.5*(1.0 + 2.0*1.5*1.5 + 3.0*1.5*1.5), state.getPotentialEnergy(), TOL);
}
void testManyParameters() {
OpenCLPlatform platform;
System system;
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomExternalForce* forceField = new CustomExternalForce("xscale*(x-x0)^2+yscale*(y-y0)^2+zscale*(z-z0)^2");
forceField->addPerParticleParameter("x0");
forceField->addPerParticleParameter("y0");
forceField->addPerParticleParameter("z0");
forceField->addPerParticleParameter("xscale");
forceField->addPerParticleParameter("yscale");
forceField->addPerParticleParameter("zscale");
vector<double> parameters(6);
parameters[0] = 1.0;
parameters[1] = 2.0;
parameters[2] = 3.0;
parameters[3] = 0.1;
parameters[4] = 0.2;
parameters[5] = 0.3;
forceField->addParticle(0, parameters);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(1);
positions[0] = Vec3(0, -1, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
ASSERT_EQUAL_VEC(Vec3(2*0.1*1.0, 2*0.2*3.0, 2*0.3*3.0), forces[0], TOL);
ASSERT_EQUAL_TOL(0.1*1*1 + 0.2*3*3 + 0.3*3*3, state.getPotentialEnergy(), TOL);
}
int main() {
try {
testForce();
testManyParameters();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
......
......@@ -115,6 +115,45 @@ void testParameters() {
ASSERT_EQUAL_TOL(1.5*3.0*(12*12*12), state.getPotentialEnergy(), TOL);
}
void testManyParameters() {
OpenCLPlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
CustomNonbondedForce* forceField = new CustomNonbondedForce("(a1*a2+b1*b2+c1*c2+d1*d2+e1*e2)*r");
forceField->addPerParticleParameter("a");
forceField->addPerParticleParameter("b");
forceField->addPerParticleParameter("c");
forceField->addPerParticleParameter("d");
forceField->addPerParticleParameter("e");
vector<double> params(5);
params[0] = 1.0;
params[1] = 2.0;
params[2] = 3.0;
params[3] = 4.0;
params[4] = 5.0;
forceField->addParticle(params);
params[0] = 1.1;
params[1] = 1.2;
params[2] = 1.3;
params[3] = 1.4;
params[4] = 1.5;
forceField->addParticle(params);
system.addForce(forceField);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(2, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
vector<Vec3> forces = state.getForces();
double force = 1*1.1 + 2*1.2 + 3*1.3 + 4*1.4 + 5*1.5;
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[1], TOL);
ASSERT_EQUAL_TOL(2*force, state.getPotentialEnergy(), TOL);
}
void testExclusions() {
OpenCLPlatform platform;
System system;
......@@ -312,6 +351,7 @@ int main() {
try {
testSimpleExpression();
testParameters();
testManyParameters();
testExclusions();
testCutoff();
testPeriodic();
......
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