Skip to content

GitLab

Commit 890963b9 authored by peastman's avatar peastman
Browse files

Beginning of CUDA implementation of CustomManyParticleForce

parent 0d5e0b55
Hide whitespace changes
Inline Side-by-side
Showing with 1001 additions and 3 deletions
platforms/cuda/include/CudaKernels.h
View file @ 890963b9
...@@ -925,6 +925,54 @@ private: ...@@ -925,6 +925,54 @@ private:
const System& system; const System& system;
}; };
/**
* This kernel is invoked by CustomManyParticleForce to calculate the forces acting on the system.
*/
class CudaCalcCustomManyParticleForceKernel : public CalcCustomManyParticleForceKernel {
public:
CudaCalcCustomManyParticleForceKernel(std::string name, const Platform& platform, CudaContext& cu, const System& system) : CalcCustomManyParticleForceKernel(name, platform),
hasInitializedKernel(false), cu(cu), params(NULL), globals(NULL), particleTypes(NULL), system(system) {
}
~CudaCalcCustomManyParticleForceKernel();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param force the CustomManyParticleForce this kernel will be used for
*/
void initialize(const System& system, const CustomManyParticleForce& force);
/**
* Execute the kernel to calculate the forces and/or energy.
*
* @param context the context in which to execute this kernel
* @param includeForces true if forces should be calculated
* @param includeEnergy true if the energy should be calculated
* @return the potential energy due to the force
*/
double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
/**
* Copy changed parameters over to a context.
*
* @param context the context to copy parameters to
* @param force the CustomManyParticleForce to copy the parameters from
*/
void copyParametersToContext(ContextImpl& context, const CustomManyParticleForce& force);
private:
CudaContext& cu;
bool hasInitializedKernel;
NonbondedMethod nonbondedMethod;
CudaParameterSet* params;
CudaArray* globals;
CudaArray* particleTypes;
std::vector<std::string> globalParamNames;
std::vector<float> globalParamValues;
std::vector<CudaArray*> tabulatedFunctions;
std::vector<void*> forceArgs;
const System& system;
CUfunction forceKernel;
};
/** /**
* This kernel is invoked by VerletIntegrator to take one time step. * This kernel is invoked by VerletIntegrator to take one time step.
*/ */
......
platforms/cuda/src/CudaKernelFactory.cpp
View file @ 890963b9
...@@ -106,6 +106,8 @@ KernelImpl* CudaKernelFactory::createKernelImpl(std::string name, const Platform ...@@ -106,6 +106,8 @@ KernelImpl* CudaKernelFactory::createKernelImpl(std::string name, const Platform
return new CudaCalcCustomHbondForceKernel(name, platform, cu, context.getSystem()); return new CudaCalcCustomHbondForceKernel(name, platform, cu, context.getSystem());
if (name == CalcCustomCompoundBondForceKernel::Name()) if (name == CalcCustomCompoundBondForceKernel::Name())
return new CudaCalcCustomCompoundBondForceKernel(name, platform, cu, context.getSystem()); return new CudaCalcCustomCompoundBondForceKernel(name, platform, cu, context.getSystem());
if (name == CalcCustomManyParticleForceKernel::Name())
return new CudaCalcCustomManyParticleForceKernel(name, platform, cu, context.getSystem());
if (name == IntegrateVerletStepKernel::Name()) if (name == IntegrateVerletStepKernel::Name())
return new CudaIntegrateVerletStepKernel(name, platform, cu); return new CudaIntegrateVerletStepKernel(name, platform, cu);
if (name == IntegrateLangevinStepKernel::Name()) if (name == IntegrateLangevinStepKernel::Name())
......
platforms/cuda/src/CudaKernels.cpp
View file @ 890963b9
...@@ -33,6 +33,7 @@ ...@@ -33,6 +33,7 @@
#include "openmm/internal/ContextImpl.h" #include "openmm/internal/ContextImpl.h"
#include "openmm/internal/CustomCompoundBondForceImpl.h" #include "openmm/internal/CustomCompoundBondForceImpl.h"
#include "openmm/internal/CustomHbondForceImpl.h" #include "openmm/internal/CustomHbondForceImpl.h"
#include "openmm/internal/CustomManyParticleForceImpl.h"
#include "openmm/internal/CustomNonbondedForceImpl.h" #include "openmm/internal/CustomNonbondedForceImpl.h"
#include "openmm/internal/NonbondedForceImpl.h" #include "openmm/internal/NonbondedForceImpl.h"
#include "CudaBondedUtilities.h" #include "CudaBondedUtilities.h"
...@@ -1963,7 +1964,7 @@ void CudaCalcCustomNonbondedForceKernel::initialize(const System& system, const ...@@ -1963,7 +1964,7 @@ void CudaCalcCustomNonbondedForceKernel::initialize(const System& system, const
map<string, Lepton::CustomFunction*> functions; map<string, Lepton::CustomFunction*> functions;
vector<pair<string, string> > functionDefinitions; vector<pair<string, string> > functionDefinitions;
vector<const TabulatedFunction*> functionList; vector<const TabulatedFunction*> functionList;
for (int i = 0; i < force.getNumFunctions(); i++) { for (int i = 0; i < force.getNumTabulatedFunctions(); i++) {
functionList.push_back(&force.getTabulatedFunction(i)); functionList.push_back(&force.getTabulatedFunction(i));
string name = force.getTabulatedFunctionName(i); string name = force.getTabulatedFunctionName(i);
string arrayName = prefix+"table"+cu.intToString(i); string arrayName = prefix+"table"+cu.intToString(i);
...@@ -3766,7 +3767,7 @@ void CudaCalcCustomHbondForceKernel::initialize(const System& system, const Cust ...@@ -3766,7 +3767,7 @@ void CudaCalcCustomHbondForceKernel::initialize(const System& system, const Cust
vector<pair<string, string> > functionDefinitions; vector<pair<string, string> > functionDefinitions;
vector<const TabulatedFunction*> functionList; vector<const TabulatedFunction*> functionList;
stringstream tableArgs; stringstream tableArgs;
for (int i = 0; i < force.getNumFunctions(); i++) { for (int i = 0; i < force.getNumTabulatedFunctions(); i++) {
functionList.push_back(&force.getTabulatedFunction(i)); functionList.push_back(&force.getTabulatedFunction(i));
string name = force.getTabulatedFunctionName(i); string name = force.getTabulatedFunctionName(i);
string arrayName = "table"+cu.intToString(i); string arrayName = "table"+cu.intToString(i);
...@@ -4147,7 +4148,7 @@ void CudaCalcCustomCompoundBondForceKernel::initialize(const System& system, con ...@@ -4147,7 +4148,7 @@ void CudaCalcCustomCompoundBondForceKernel::initialize(const System& system, con
vector<pair<string, string> > functionDefinitions; vector<pair<string, string> > functionDefinitions;
vector<const TabulatedFunction*> functionList; vector<const TabulatedFunction*> functionList;
stringstream tableArgs; stringstream tableArgs;
for (int i = 0; i < force.getNumFunctions(); i++) { for (int i = 0; i < force.getNumTabulatedFunctions(); i++) {
functionList.push_back(&force.getTabulatedFunction(i)); functionList.push_back(&force.getTabulatedFunction(i));
string name = force.getTabulatedFunctionName(i); string name = force.getTabulatedFunctionName(i);
functions[name] = cu.getExpressionUtilities().getFunctionPlaceholder(force.getTabulatedFunction(i)); functions[name] = cu.getExpressionUtilities().getFunctionPlaceholder(force.getTabulatedFunction(i));
...@@ -4395,6 +4396,387 @@ void CudaCalcCustomCompoundBondForceKernel::copyParametersToContext(ContextImpl& ...@@ -4395,6 +4396,387 @@ void CudaCalcCustomCompoundBondForceKernel::copyParametersToContext(ContextImpl&
cu.invalidateMolecules(); cu.invalidateMolecules();
} }
class CudaCustomManyParticleForceInfo : public CudaForceInfo {
public:
CudaCustomManyParticleForceInfo(const CustomManyParticleForce& force) : force(force) {
}
bool areParticlesIdentical(int particle1, int particle2) {
vector<double> params1, params2;
int type1, type2;
force.getParticleParameters(particle1, params1, type1);
force.getParticleParameters(particle2, params2, type2);
if (type1 != type2)
return false;
for (int i = 0; i < (int) params1.size(); i++)
if (params1[i] != params2[i])
return false;
return true;
}
int getNumParticleGroups() {
return 0;
}
void getParticlesInGroup(int index, vector<int>& particles) {
}
bool areGroupsIdentical(int group1, int group2) {
return true;
}
private:
const CustomManyParticleForce& force;
};
CudaCalcCustomManyParticleForceKernel::~CudaCalcCustomManyParticleForceKernel() {
cu.setAsCurrent();
if (params != NULL)
delete params;
if (globals != NULL)
delete globals;
if (particleTypes != NULL)
delete particleTypes;
for (int i = 0; i < (int) tabulatedFunctions.size(); i++)
delete tabulatedFunctions[i];
}
void CudaCalcCustomManyParticleForceKernel::initialize(const System& system, const CustomManyParticleForce& force) {
cu.setAsCurrent();
int numParticles = force.getNumParticles();
int particlesPerSet = force.getNumParticlesPerSet();
nonbondedMethod = CalcCustomManyParticleForceKernel::NonbondedMethod(force.getNonbondedMethod());
// Record parameter values.
params = new CudaParameterSet(cu, force.getNumPerParticleParameters(), numParticles, "customManyParticleParameters");
if (force.getNumGlobalParameters() > 0)
globals = CudaArray::create<float>(cu, force.getNumGlobalParameters(), "customManyParticleGlobals");
vector<vector<float> > paramVector(numParticles);
for (int i = 0; i < numParticles; i++) {
vector<double> parameters;
int type;
force.getParticleParameters(i, parameters, type);
paramVector[i].resize(parameters.size());
for (int j = 0; j < (int) parameters.size(); j++)
paramVector[i][j] = (float) parameters[j];
}
params->setParameterValues(paramVector);
cu.addForce(new CudaCustomManyParticleForceInfo(force));
// Record the tabulated functions.
map<string, Lepton::CustomFunction*> functions;
vector<pair<string, string> > functionDefinitions;
vector<const TabulatedFunction*> functionList;
stringstream tableArgs;
for (int i = 0; i < force.getNumTabulatedFunctions(); i++) {
functionList.push_back(&force.getTabulatedFunction(i));
string name = force.getTabulatedFunctionName(i);
functions[name] = cu.getExpressionUtilities().getFunctionPlaceholder(force.getTabulatedFunction(i));
int width;
vector<float> f = cu.getExpressionUtilities().computeFunctionCoefficients(force.getTabulatedFunction(i), width);
CudaArray* array = CudaArray::create<float>(cu, f.size(), "TabulatedFunction");
tabulatedFunctions.push_back(array);
array->upload(f);
string arrayName = cu.getBondedUtilities().addArgument(array->getDevicePointer(), width == 1 ? "float" : "float"+cu.intToString(width));
functionDefinitions.push_back(make_pair(name, arrayName));
}
// Record information about parameters.
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);
}
map<string, string> variables;
for (int i = 0; i < particlesPerSet; i++) {
string index = cu.intToString(i+1);
variables["x"+index] = "pos"+index+".x";
variables["y"+index] = "pos"+index+".y";
variables["z"+index] = "pos"+index+".z";
}
for (int i = 0; i < force.getNumPerParticleParameters(); i++) {
const string& name = force.getPerParticleParameterName(i);
variables[name] = "params"+params->getParameterSuffix(i);
}
if (force.getNumGlobalParameters() > 0) {
globals = CudaArray::create<float>(cu, force.getNumGlobalParameters(), "customManyParticleGlobals");
globals->upload(globalParamValues);
for (int i = 0; i < force.getNumGlobalParameters(); i++) {
const string& name = force.getGlobalParameterName(i);
string value = "globals["+cu.intToString(i)+"]";
variables[name] = value;
}
}
// Now to generate the kernel. First, it needs to calculate all distances, angles,
// and dihedrals the expression depends on.
map<string, vector<int> > distances;
map<string, vector<int> > angles;
map<string, vector<int> > dihedrals;
Lepton::ParsedExpression energyExpression = CustomManyParticleForceImpl::prepareExpression(force, functions, distances, angles, dihedrals);
map<string, Lepton::ParsedExpression> forceExpressions;
set<string> computedDeltas;
vector<string> atomNames, posNames;
for (int i = 0; i < particlesPerSet; i++) {
string index = cu.intToString(i+1);
atomNames.push_back("P"+index);
posNames.push_back("pos"+index);
}
stringstream compute;
int index = 0;
for (map<string, vector<int> >::const_iterator iter = distances.begin(); iter != distances.end(); ++iter, ++index) {
const vector<int>& atoms = iter->second;
string deltaName = atomNames[atoms[0]]+atomNames[atoms[1]];
if (computedDeltas.count(deltaName) == 0) {
compute<<"real4 delta"<<deltaName<<" = delta("<<posNames[atoms[0]]<<", "<<posNames[atoms[1]]<<");\n";
computedDeltas.insert(deltaName);
}
compute<<"real r_"<<deltaName<<" = sqrt(delta"<<deltaName<<".w);\n";
variables[iter->first] = "r_"+deltaName;
forceExpressions["real dEdDistance"+cu.intToString(index)+" = "] = energyExpression.differentiate(iter->first).optimize();
}
index = 0;
for (map<string, vector<int> >::const_iterator iter = angles.begin(); iter != angles.end(); ++iter, ++index) {
const vector<int>& atoms = iter->second;
string deltaName1 = atomNames[atoms[1]]+atomNames[atoms[0]];
string deltaName2 = atomNames[atoms[1]]+atomNames[atoms[2]];
string angleName = "angle_"+atomNames[atoms[0]]+atomNames[atoms[1]]+atomNames[atoms[2]];
if (computedDeltas.count(deltaName1) == 0) {
compute<<"real4 delta"<<deltaName1<<" = delta("<<posNames[atoms[1]]<<", "<<posNames[atoms[0]]<<");\n";
computedDeltas.insert(deltaName1);
}
if (computedDeltas.count(deltaName2) == 0) {
compute<<"real4 delta"<<deltaName2<<" = delta("<<posNames[atoms[1]]<<", "<<posNames[atoms[2]]<<");\n";
computedDeltas.insert(deltaName2);
}
compute<<"real "<<angleName<<" = computeAngle(delta"<<deltaName1<<", delta"<<deltaName2<<");\n";
variables[iter->first] = angleName;
forceExpressions["real dEdAngle"+cu.intToString(index)+" = "] = energyExpression.differentiate(iter->first).optimize();
}
index = 0;
for (map<string, vector<int> >::const_iterator iter = dihedrals.begin(); iter != dihedrals.end(); ++iter, ++index) {
const vector<int>& atoms = iter->second;
string deltaName1 = atomNames[atoms[0]]+atomNames[atoms[1]];
string deltaName2 = atomNames[atoms[2]]+atomNames[atoms[1]];
string deltaName3 = atomNames[atoms[2]]+atomNames[atoms[3]];
string crossName1 = "cross_"+deltaName1+"_"+deltaName2;
string crossName2 = "cross_"+deltaName2+"_"+deltaName3;
string dihedralName = "dihedral_"+atomNames[atoms[0]]+atomNames[atoms[1]]+atomNames[atoms[2]]+atomNames[atoms[3]];
if (computedDeltas.count(deltaName1) == 0) {
compute<<"real4 delta"<<deltaName1<<" = delta("<<posNames[atoms[0]]<<", "<<posNames[atoms[1]]<<");\n";
computedDeltas.insert(deltaName1);
}
if (computedDeltas.count(deltaName2) == 0) {
compute<<"real4 delta"<<deltaName2<<" = delta("<<posNames[atoms[2]]<<", "<<posNames[atoms[1]]<<");\n";
computedDeltas.insert(deltaName2);
}
if (computedDeltas.count(deltaName3) == 0) {
compute<<"real4 delta"<<deltaName3<<" = delta("<<posNames[atoms[2]]<<", "<<posNames[atoms[3]]<<");\n";
computedDeltas.insert(deltaName3);
}
compute<<"real4 "<<crossName1<<" = computeCross(delta"<<deltaName1<<", delta"<<deltaName2<<");\n";
compute<<"real4 "<<crossName2<<" = computeCross(delta"<<deltaName2<<", delta"<<deltaName3<<");\n";
compute<<"real "<<dihedralName<<" = computeAngle("<<crossName1<<", "<<crossName2<<");\n";
compute<<dihedralName<<" *= (delta"<<deltaName1<<".x*"<<crossName2<<".x + delta"<<deltaName1<<".y*"<<crossName2<<".y + delta"<<deltaName1<<".z*"<<crossName2<<".z < 0 ? -1 : 1);\n";
variables[iter->first] = dihedralName;
forceExpressions["real dEdDihedral"+cu.intToString(index)+" = "] = energyExpression.differentiate(iter->first).optimize();
}
// Now evaluate the expressions.
for (int i = 0; i < (int) params->getBuffers().size(); i++) {
CudaNonbondedUtilities::ParameterInfo& buffer = params->getBuffers()[i];
compute<<buffer.getType()<<" params"<<(i+1)<<" = global_params"<<(i+1)<<"[index];\n";
}
forceExpressions["energy += "] = energyExpression;
compute << cu.getExpressionUtilities().createExpressions(forceExpressions, variables, functionList, functionDefinitions, "temp");
// Apply forces to atoms.
vector<string> forceNames;
for (int i = 0; i < particlesPerSet; i++) {
string istr = cu.intToString(i+1);
string forceName = "force"+istr;
forceNames.push_back(forceName);
compute<<"real3 "<<forceName<<" = make_real3(0);\n";
compute<<"{\n";
Lepton::ParsedExpression forceExpressionX = energyExpression.differentiate("x"+istr).optimize();
Lepton::ParsedExpression forceExpressionY = energyExpression.differentiate("y"+istr).optimize();
Lepton::ParsedExpression forceExpressionZ = energyExpression.differentiate("z"+istr).optimize();
map<string, Lepton::ParsedExpression> expressions;
if (!isZeroExpression(forceExpressionX))
expressions[forceName+".x -= "] = forceExpressionX;
if (!isZeroExpression(forceExpressionY))
expressions[forceName+".y -= "] = forceExpressionY;
if (!isZeroExpression(forceExpressionZ))
expressions[forceName+".z -= "] = forceExpressionZ;
if (expressions.size() > 0)
compute<<cu.getExpressionUtilities().createExpressions(expressions, variables, functionList, functionDefinitions, "coordtemp");
compute<<"}\n";
}
index = 0;
for (map<string, vector<int> >::const_iterator iter = distances.begin(); iter != distances.end(); ++iter, ++index) {
const vector<int>& atoms = iter->second;
string deltaName = atomNames[atoms[0]]+atomNames[atoms[1]];
string value = "(dEdDistance"+cu.intToString(index)+"/r_"+deltaName+")*trim(delta"+deltaName+")";
compute<<forceNames[atoms[0]]<<" += "<<"-"<<value<<";\n";
compute<<forceNames[atoms[1]]<<" += "<<value<<";\n";
}
index = 0;
for (map<string, vector<int> >::const_iterator iter = angles.begin(); iter != angles.end(); ++iter, ++index) {
const vector<int>& atoms = iter->second;
string deltaName1 = atomNames[atoms[1]]+atomNames[atoms[0]];
string deltaName2 = atomNames[atoms[1]]+atomNames[atoms[2]];
compute<<"{\n";
compute<<"real3 crossProd = cross(delta"<<deltaName2<<", delta"<<deltaName1<<");\n";
compute<<"real lengthCross = max(SQRT(dot(crossProd, crossProd)), 1e-6f);\n";
compute<<"real3 deltaCross0 = -cross(trim(delta"<<deltaName1<<"), crossProd)*dEdAngle"<<cu.intToString(index)<<"/(delta"<<deltaName1<<".w*lengthCross);\n";
compute<<"real3 deltaCross2 = cross(trim(delta"<<deltaName2<<"), crossProd)*dEdAngle"<<cu.intToString(index)<<"/(delta"<<deltaName2<<".w*lengthCross);\n";
compute<<"real3 deltaCross1 = -(deltaCross0+deltaCross2);\n";
compute<<forceNames[atoms[0]]<<" += deltaCross0;\n";
compute<<forceNames[atoms[1]]<<" += deltaCross1;\n";
compute<<forceNames[atoms[2]]<<" += deltaCross2;\n";
compute<<"}\n";
}
index = 0;
for (map<string, vector<int> >::const_iterator iter = dihedrals.begin(); iter != dihedrals.end(); ++iter, ++index) {
const vector<int>& atoms = iter->second;
string deltaName1 = atomNames[atoms[0]]+atomNames[atoms[1]];
string deltaName2 = atomNames[atoms[2]]+atomNames[atoms[1]];
string deltaName3 = atomNames[atoms[2]]+atomNames[atoms[3]];
string crossName1 = "cross_"+deltaName1+"_"+deltaName2;
string crossName2 = "cross_"+deltaName2+"_"+deltaName3;
compute<<"{\n";
compute<<"real r = sqrt(delta"<<deltaName2<<".w);\n";
compute<<"real4 ff;\n";
compute<<"ff.x = (-dEdDihedral"<<cu.intToString(index)<<"*r)/"<<crossName1<<".w;\n";
compute<<"ff.y = (delta"<<deltaName1<<".x*delta"<<deltaName2<<".x + delta"<<deltaName1<<".y*delta"<<deltaName2<<".y + delta"<<deltaName1<<".z*delta"<<deltaName2<<".z)/delta"<<deltaName2<<".w;\n";
compute<<"ff.z = (delta"<<deltaName3<<".x*delta"<<deltaName2<<".x + delta"<<deltaName3<<".y*delta"<<deltaName2<<".y + delta"<<deltaName3<<".z*delta"<<deltaName2<<".z)/delta"<<deltaName2<<".w;\n";
compute<<"ff.w = (dEdDihedral"<<cu.intToString(index)<<"*r)/"<<crossName2<<".w;\n";
compute<<"real3 internalF0 = ff.x*trim("<<crossName1<<");\n";
compute<<"real3 internalF3 = ff.w*trim("<<crossName2<<");\n";
compute<<"real3 s = ff.y*internalF0 - ff.z*internalF3;\n";
compute<<forceNames[atoms[0]]<<" += internalF0;\n";
compute<<forceNames[atoms[1]]<<" += s-internalF0;\n";
compute<<forceNames[atoms[2]]<<" += -s-internalF3;\n";
compute<<forceNames[atoms[3]]<<" += internalF3;\n";
compute<<"}\n";
}
// Store forces to global memory.
for (int i = 0; i < particlesPerSet; i++)
compute<<"storeForce(atom"<<(i+1)<<", "<<forceNames[i]<<", forceBuffers);\n";
// Create other replacements that depend on the number of particles per set.
stringstream numCombinations, atomsForCombination, isValidCombination, permute, loadData;
for (int i = 0; i < particlesPerSet; i++) {
permute<<"int atom"<<(i+1)<<" = p"<<(i+1)<<";\n";
loadData<<"real4 pos"<<(i+1)<<" = posq[atom"<<(i+1)<<"];\n";
}
for (int i = 2; i < particlesPerSet; i++) {
if (i > 2)
isValidCombination<<" && ";
isValidCombination<<"p"<<(i+1)<<">p"<<i;
}
atomsForCombination<<"int tempIndex = index;\n";
for (int i = 1; i < particlesPerSet; i++) {
if (i > 1)
numCombinations<<"*";
numCombinations<<"numNeighbors";
atomsForCombination<<"int p"<<(i+1)<<" = p1+1+tempIndex%numNeighbors;\n";
atomsForCombination<<"tempIndex /= numNeighbors;\n";
}
// Create replacements for extra arguments.
stringstream extraArgs;
if (force.getNumGlobalParameters() > 0)
extraArgs<<", const float* __restrict__ globals";
for (int i = 0; i < (int) params->getBuffers().size(); i++) {
CudaNonbondedUtilities::ParameterInfo& buffer = params->getBuffers()[i];
extraArgs<<", const "<<buffer.getType()<<"* __restrict__ global_params";
}
// Create the kernels.
map<string, string> replacements;
replacements["COMPUTE_INTERACTION"] = compute.str();
replacements["NUM_CANDIDATE_COMBINATIONS"] = numCombinations.str();
replacements["FIND_ATOMS_FOR_COMBINATION_INDEX"] = atomsForCombination.str();
replacements["IS_VALID_COMBINATION"] = isValidCombination.str();
replacements["PERMUTE_ATOMS"] = permute.str();
replacements["LOAD_PARTICLE_DATA"] = loadData.str();
replacements["PARAMETER_ARGUMENTS"] = extraArgs.str()+tableArgs.str();
map<string, string> defines;
defines["NUM_ATOMS"] = cu.intToString(cu.getNumAtoms());
defines["PADDED_NUM_ATOMS"] = cu.intToString(cu.getPaddedNumAtoms());
defines["M_PI"] = cu.doubleToString(M_PI);
std::cout << cu.replaceStrings(CudaKernelSources::vectorOps+CudaKernelSources::customManyParticle, replacements)<< std::endl;
CUmodule module = cu.createModule(cu.replaceStrings(CudaKernelSources::vectorOps+CudaKernelSources::customManyParticle, replacements), defines);
forceKernel = cu.getKernel(module, "computeInteraction");
}
double CudaCalcCustomManyParticleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
if (!hasInitializedKernel) {
hasInitializedKernel = true;
int index = 0;
forceArgs.push_back(&cu.getForce().getDevicePointer());
forceArgs.push_back(&cu.getEnergyBuffer().getDevicePointer());
forceArgs.push_back(&cu.getPosq().getDevicePointer());
if (nonbondedMethod != NoCutoff) {
forceArgs.push_back(cu.getPeriodicBoxSizePointer());
forceArgs.push_back(cu.getInvPeriodicBoxSizePointer());
}
if (globals != NULL)
forceArgs.push_back(&globals->getDevicePointer());
for (int i = 0; i < (int) params->getBuffers().size(); i++) {
CudaNonbondedUtilities::ParameterInfo& buffer = params->getBuffers()[i];
forceArgs.push_back(&buffer.getMemory());
}
for (int i = 0; i < (int) tabulatedFunctions.size(); i++)
forceArgs.push_back(&tabulatedFunctions[i]->getDevicePointer());
}
if (globals != NULL) {
bool changed = false;
for (int i = 0; i < (int) globalParamNames.size(); i++) {
float value = (float) context.getParameter(globalParamNames[i]);
if (value != globalParamValues[i])
changed = true;
globalParamValues[i] = value;
}
if (changed)
globals->upload(globalParamValues);
}
cu.executeKernel(forceKernel, &forceArgs[0], cu.getNumAtoms()*CudaContext::ThreadBlockSize, CudaContext::ThreadBlockSize);
return 0.0;
}
void CudaCalcCustomManyParticleForceKernel::copyParametersToContext(ContextImpl& context, const CustomManyParticleForce& force) {
cu.setAsCurrent();
int numParticles = force.getNumParticles();
if (numParticles != cu.getNumAtoms())
throw OpenMMException("updateParametersInContext: The number of particles has changed");
// Record the per-particle parameters.
vector<vector<float> > paramVector(numParticles);
vector<double> parameters;
int type;
for (int i = 0; i < numParticles; i++) {
force.getParticleParameters(i, parameters, type);
paramVector[i].resize(parameters.size());
for (int j = 0; j < (int) parameters.size(); j++)
paramVector[i][j] = (float) parameters[j];
}
params->setParameterValues(paramVector);
// Mark that the current reordering may be invalid.
cu.invalidateMolecules();
}
CudaIntegrateVerletStepKernel::~CudaIntegrateVerletStepKernel() { CudaIntegrateVerletStepKernel::~CudaIntegrateVerletStepKernel() {
} }
......
platforms/cuda/src/CudaPlatform.cpp
View file @ 890963b9
...@@ -81,6 +81,7 @@ CudaPlatform::CudaPlatform() { ...@@ -81,6 +81,7 @@ CudaPlatform::CudaPlatform() {
registerKernelFactory(CalcCustomExternalForceKernel::Name(), factory); registerKernelFactory(CalcCustomExternalForceKernel::Name(), factory);
registerKernelFactory(CalcCustomHbondForceKernel::Name(), factory); registerKernelFactory(CalcCustomHbondForceKernel::Name(), factory);
registerKernelFactory(CalcCustomCompoundBondForceKernel::Name(), factory); registerKernelFactory(CalcCustomCompoundBondForceKernel::Name(), factory);
registerKernelFactory(CalcCustomManyParticleForceKernel::Name(), factory);
registerKernelFactory(IntegrateVerletStepKernel::Name(), factory); registerKernelFactory(IntegrateVerletStepKernel::Name(), factory);
registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory); registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);
registerKernelFactory(IntegrateBrownianStepKernel::Name(), factory); registerKernelFactory(IntegrateBrownianStepKernel::Name(), factory);
......
platforms/cuda/src/kernels/customManyParticle.cu 0 → 100644
View file @ 890963b9
/**
* Record the force on an atom to global memory.
*/
inline __device__ void storeForce(int atom, real3 force, unsigned long long* __restrict__ forceBuffers) {
atomicAdd(&forceBuffers[atom], static_cast<unsigned long long>((long long) (force.x*0x100000000)));
atomicAdd(&forceBuffers[atom+PADDED_NUM_ATOMS], static_cast<unsigned long long>((long long) (force.y*0x100000000)));
atomicAdd(&forceBuffers[atom+2*PADDED_NUM_ATOMS], static_cast<unsigned long long>((long long) (force.z*0x100000000)));
}
/**
* Convert a real4 to a real3 by removing its last element.
*/
inline __device__ real3 trim(real4 v) {
return make_real3(v.x, v.y, v.z);
}
/**
* Compute the difference between two vectors, setting the fourth component to the squared magnitude.
*/
inline __device__ real4 delta(real4 vec1, real4 vec2) {
real4 result = make_real4(vec1.x-vec2.x, vec1.y-vec2.y, vec1.z-vec2.z, 0);
result.w = result.x*result.x + result.y*result.y + result.z*result.z;
return result;
}
/**
* Compute the angle between two vectors. The w component of each vector should contain the squared magnitude.
*/
__device__ real computeAngle(real4 vec1, real4 vec2) {
real dotProduct = vec1.x*vec2.x + vec1.y*vec2.y + vec1.z*vec2.z;
real cosine = dotProduct*RSQRT(vec1.w*vec2.w);
real angle;
if (cosine > 0.99f || cosine < -0.99f) {
// We're close to the singularity in acos(), so take the cross product and use asin() instead.
real3 crossProduct = cross(vec1, vec2);
real scale = vec1.w*vec2.w;
angle = ASIN(SQRT(dot(crossProduct, crossProduct)/scale));
if (cosine < 0.0f)
angle = M_PI-angle;
}
else
angle = ACOS(cosine);
return angle;
}
/**
* Compute the cross product of two vectors, setting the fourth component to the squared magnitude.
*/
inline __device__ real4 computeCross(real4 vec1, real4 vec2) {
real3 cp = cross(vec1, vec2);
return make_real4(cp.x, cp.y, cp.z, cp.x*cp.x+cp.y*cp.y+cp.z*cp.z);
}
/**
* Compute the interaction.
*/
extern "C" __global__ void computeInteraction(
unsigned long long* __restrict__ forceBuffers, real* __restrict__ energyBuffer, const real4* __restrict__ posq
#ifdef USE_CUTOFF
, real4 periodicBoxSize, real4 invPeriodicBoxSize
#endif
PARAMETER_ARGUMENTS) {
real energy = 0.0f;
// Loop over particles to be the first one in the set.
for (int p1 = blockIdx.x; p1 < NUM_ATOMS; p1 += gridDim.x) {
int numNeighbors = NUM_ATOMS-p1-1;
int numCombinations = NUM_CANDIDATE_COMBINATIONS;
for (int index = threadIdx.x; index < numCombinations; index += blockDim.x) {
FIND_ATOMS_FOR_COMBINATION_INDEX;
bool includeInteraction = IS_VALID_COMBINATION;
if (includeInteraction) {
PERMUTE_ATOMS;
LOAD_PARTICLE_DATA;
COMPUTE_INTERACTION;
}
}
}
energyBuffer[blockIdx.x*blockDim.x+threadIdx.x] += energy;
}
\ No newline at end of file
platforms/cuda/tests/TestCudaCustomManyParticleForce.cpp 0 → 100644
View file @ 890963b9
/* -------------------------------------------------------------------------- *
* 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) 2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the CPU implementation of CustomManyParticleForce.
*/
#ifdef WIN32
#define _USE_MATH_DEFINES // Needed to get M_PI
#endif
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "CudaPlatform.h"
#include "openmm/CustomCompoundBondForce.h"
#include "openmm/CustomManyParticleForce.h"
#include "openmm/System.h"
#include "openmm/TabulatedFunction.h"
#include "openmm/VerletIntegrator.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
CudaPlatform platform;
void validateAxilrodTeller(CustomManyParticleForce* force, const vector<Vec3>& positions, const vector<const int*>& expectedSets, double boxSize) {
// Create a System and Context.
int numParticles = force->getNumParticles();
CustomManyParticleForce::NonbondedMethod nonbondedMethod = force->getNonbondedMethod();
System system;
for (int i = 0; i < numParticles; i++)
system.addParticle(1.0);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
system.addForce(force);
VerletIntegrator integrator(0.001);
Context context(system, integrator, platform);
context.setPositions(positions);
State state1 = context.getState(State::Forces | State::Energy);
double c = context.getParameter("C");
// See if the energy matches the expected value.
double expectedEnergy = 0;
for (int i = 0; i < (int) expectedSets.size(); i++) {
int p1 = expectedSets[i][0];
int p2 = expectedSets[i][1];
int p3 = expectedSets[i][2];
Vec3 d12 = positions[p2]-positions[p1];
Vec3 d13 = positions[p3]-positions[p1];
Vec3 d23 = positions[p3]-positions[p2];
if (nonbondedMethod == CustomManyParticleForce::CutoffPeriodic) {
for (int j = 0; j < 3; j++) {
d12[j] -= floor(d12[j]/boxSize+0.5f)*boxSize;
d13[j] -= floor(d13[j]/boxSize+0.5f)*boxSize;
d23[j] -= floor(d23[j]/boxSize+0.5f)*boxSize;
}
}
double r12 = sqrt(d12.dot(d12));
double r13 = sqrt(d13.dot(d13));
double r23 = sqrt(d23.dot(d23));
double ctheta1 = d12.dot(d13)/(r12*r13);
double ctheta2 = -d12.dot(d23)/(r12*r23);
double ctheta3 = d13.dot(d23)/(r13*r23);
double rprod = r12*r13*r23;
expectedEnergy += c*(1+3*ctheta1*ctheta2*ctheta3)/(rprod*rprod*rprod);
}
ASSERT_EQUAL_TOL(expectedEnergy, state1.getPotentialEnergy(), 1e-5);
// Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.
const vector<Vec3>& forces = state1.getForces();
double norm = 0.0;
for (int i = 0; i < (int) forces.size(); ++i)
norm += forces[i].dot(forces[i]);
norm = std::sqrt(norm);
const double stepSize = 1e-3;
double step = 0.5*stepSize/norm;
vector<Vec3> positions2(numParticles), positions3(numParticles);
for (int i = 0; i < (int) positions.size(); ++i) {
Vec3 p = positions[i];
Vec3 f = forces[i];
positions2[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
positions3[i] = Vec3(p[0]+f[0]*step, p[1]+f[1]*step, p[2]+f[2]*step);
}
context.setPositions(positions2);
State state2 = context.getState(State::Energy);
context.setPositions(positions3);
State state3 = context.getState(State::Energy);
ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state3.getPotentialEnergy())/stepSize, 1e-4);
}
void testNoCutoff() {
CustomManyParticleForce* force = new CustomManyParticleForce(3,
"C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
"theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
"r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
force->addGlobalParameter("C", 1.5);
vector<double> params;
force->addParticle(params);
force->addParticle(params);
force->addParticle(params);
force->addParticle(params);
vector<Vec3> positions;
positions.push_back(Vec3(0, 0, 0));
positions.push_back(Vec3(1, 0, 0));
positions.push_back(Vec3(0, 1.1, 0.3));
positions.push_back(Vec3(0.4, 0, -0.8));
int sets[4][3] = {{0,1,2}, {1,2,3}, {2,3,0}, {3,0,1}};
vector<const int*> expectedSets(&sets[0], &sets[4]);
validateAxilrodTeller(force, positions, expectedSets, 2.0);
}
void testCutoff() {
CustomManyParticleForce* force = new CustomManyParticleForce(3,
"C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
"theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
"r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
force->addGlobalParameter("C", 1.5);
force->setNonbondedMethod(CustomManyParticleForce::CutoffNonPeriodic);
force->setCutoffDistance(1.55);
vector<double> params;
force->addParticle(params);
force->addParticle(params);
force->addParticle(params);
force->addParticle(params);
force->addParticle(params);
vector<Vec3> positions;
positions.push_back(Vec3(0, 0, 0));
positions.push_back(Vec3(1, 0, 0));
positions.push_back(Vec3(0, 1.1, 0.3));
positions.push_back(Vec3(0.4, 0, -0.8));
positions.push_back(Vec3(0.2, 0.5, -0.1));
int sets[7][3] = {{0,1,2}, {0,1,3}, {0,1,4}, {0,2,4}, {0,3,4}, {1,2,4}, {1,3,4}};
vector<const int*> expectedSets(&sets[0], &sets[7]);
validateAxilrodTeller(force, positions, expectedSets, 2.0);
}
void testPeriodic() {
CustomManyParticleForce* force = new CustomManyParticleForce(3,
"C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
"theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
"r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
force->addGlobalParameter("C", 1.5);
force->setNonbondedMethod(CustomManyParticleForce::CutoffPeriodic);
force->setCutoffDistance(1.05);
vector<double> params;
force->addParticle(params);
force->addParticle(params);
force->addParticle(params);
force->addParticle(params);
force->addParticle(params);
vector<Vec3> positions;
positions.push_back(Vec3(0, 0, 0));
positions.push_back(Vec3(1, 0, 0));
positions.push_back(Vec3(0, 1.1, 0.3));
positions.push_back(Vec3(0.4, 0, -0.8));
positions.push_back(Vec3(0.2, 0.5, -0.1));
double boxSize = 2.1;
int sets[5][3] = {{0,1,3}, {0,1,4}, {0,2,4}, {0,3,4}, {1,3,4}};
vector<const int*> expectedSets(&sets[0], &sets[5]);
validateAxilrodTeller(force, positions, expectedSets, boxSize);
}
void testExclusions() {
CustomManyParticleForce* force = new CustomManyParticleForce(3,
"C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
"theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
"r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
force->addGlobalParameter("C", 1.5);
vector<double> params;
force->addParticle(params);
force->addParticle(params);
force->addParticle(params);
force->addParticle(params);
force->addParticle(params);
vector<Vec3> positions;
positions.push_back(Vec3(0, 0, 0));
positions.push_back(Vec3(1, 0, 0));
positions.push_back(Vec3(0, 1.1, 0.3));
positions.push_back(Vec3(0.4, 0, -0.8));
positions.push_back(Vec3(0.2, 0.5, -0.1));
force->addExclusion(0, 2);
force->addExclusion(0, 3);
int sets[5][3] = {{0,1,4}, {1,2,3}, {1,2,4}, {1,3,4}, {2,3,4}};
vector<const int*> expectedSets(&sets[0], &sets[5]);
validateAxilrodTeller(force, positions, expectedSets, 2.0);
}
void testAllTerms() {
int numParticles = 4;
// Create a system with a CustomManyParticleForce.
System system1;
CustomManyParticleForce* force1 = new CustomManyParticleForce(4,
"distance(p1,p2)+angle(p1,p4,p3)+dihedral(p1,p3,p2,p4)+x1+y4+z3");
system1.addForce(force1);
vector<double> params;
for (int i = 0; i < numParticles; i++) {
system1.addParticle(1.0);
force1->addParticle(params, i);
}
set<int> filter;
filter.insert(0);
force1->setTypeFilter(0, filter);
filter.clear();
filter.insert(1);
force1->setTypeFilter(1, filter);
filter.clear();
filter.insert(3);
force1->setTypeFilter(2, filter);
filter.clear();
filter.insert(2);
force1->setTypeFilter(3, filter);
// Create a system that use a CustomCompoundBondForce to compute exactly the same interactions.
System system2;
CustomCompoundBondForce* force2 = new CustomCompoundBondForce(4,
"distance(p1,p2)+angle(p1,p3,p4)+dihedral(p1,p4,p2,p3)+x1+y3+z4");
system2.addForce(force2);
vector<int> particles;
particles.push_back(0);
particles.push_back(1);
particles.push_back(2);
particles.push_back(3);
force2->addBond(particles, params);
for (int i = 0; i < numParticles; i++)
system2.addParticle(1.0);
// Create contexts for both of them.
vector<Vec3> positions;
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; i++)
positions.push_back(Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt)));
VerletIntegrator integrator1(0.001);
VerletIntegrator integrator2(0.001);
Context context1(system1, integrator1, platform);
Context context2(system2, integrator2, platform);
context1.setPositions(positions);
context2.setPositions(positions);
// See if they produce identical forces and energies.
State state1 = context1.getState(State::Forces | State::Energy);
State state2 = context2.getState(State::Forces | State::Energy);
ASSERT_EQUAL_TOL(state2.getPotentialEnergy(), state1.getPotentialEnergy(), 1e-4);
for (int i = 0; i < numParticles; i++)
ASSERT_EQUAL_VEC(state2.getForces()[i], state1.getForces()[i], 1e-4);
}
void testParameters() {
// Create a system.
int numParticles = 5;
System system;
CustomManyParticleForce* force = new CustomManyParticleForce(3, "C*scale1*scale2*scale3*(distance(p1,p2)+distance(p2,p3)+distance(p1,p3))");
force->addGlobalParameter("C", 2.0);
force->addPerParticleParameter("scale");
vector<double> params(1);
vector<Vec3> positions;
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < numParticles; i++) {
params[0] = i+1;
force->addParticle(params);
positions.push_back(Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt)));
system.addParticle(1.0);
}
system.addForce(force);
VerletIntegrator integrator(0.001);
Context context(system, integrator, platform);
context.setPositions(positions);
// See if the energy is correct.
State state = context.getState(State::Energy);
double expectedEnergy = 0;
for (int i = 0; i < numParticles; i++)
for (int j = i+1; j < numParticles; j++)
for (int k = j+1; k < numParticles; k++) {
Vec3 d12 = positions[j]-positions[i];
Vec3 d13 = positions[k]-positions[i];
Vec3 d23 = positions[k]-positions[j];
double r12 = sqrt(d12.dot(d12));
double r13 = sqrt(d13.dot(d13));
double r23 = sqrt(d23.dot(d23));
expectedEnergy += 2.0*(i+1)*(j+1)*(k+1)*(r12+r13+r23);
}
ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), 1e-5);
// Modify the parameters.
context.setParameter("C", 3.5);
for (int i = 0; i < numParticles; i++) {
params[0] = 0.5*i-0.1;
force->setParticleParameters(i, params, 0);
}
force->updateParametersInContext(context);
// See if the energy is still correct.
state = context.getState(State::Energy);
expectedEnergy = 0;
for (int i = 0; i < numParticles; i++)
for (int j = i+1; j < numParticles; j++)
for (int k = j+1; k < numParticles; k++) {
Vec3 d12 = positions[j]-positions[i];
Vec3 d13 = positions[k]-positions[i];
Vec3 d23 = positions[k]-positions[j];
double r12 = sqrt(d12.dot(d12));
double r13 = sqrt(d13.dot(d13));
double r23 = sqrt(d23.dot(d23));
expectedEnergy += 3.5*(0.5*i-0.1)*(0.5*j-0.1)*(0.5*k-0.1)*(r12+r13+r23);
}
ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), 1e-5);
}
void testTabulatedFunctions() {
int numParticles = 5;
// Create two tabulated functions.
vector<double> values;
values.push_back(0.0);
values.push_back(50.0);
Continuous1DFunction* f1 = new Continuous1DFunction(values, 0, 100);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<double> c(numParticles);
for (int i = 0; i < numParticles; i++)
c[i] = genrand_real2(sfmt);
values.resize(numParticles*numParticles*numParticles);
for (int i = 0; i < numParticles; i++)
for (int j = 0; j < numParticles; j++)
for (int k = 0; k < numParticles; k++)
values[i+numParticles*j+numParticles*numParticles*k] = c[i]+c[j]+c[k];
Discrete3DFunction* f2 = new Discrete3DFunction(numParticles, numParticles, numParticles, values);
// Create a system.
System system;
CustomManyParticleForce* force = new CustomManyParticleForce(3, "f1(distance(p1,p2)+distance(p2,p3)+distance(p1,p3))*f2(atom1, atom2, atom3)");
force->addPerParticleParameter("atom");
force->addTabulatedFunction("f1", f1);
force->addTabulatedFunction("f2", f2);
vector<double> params(1);
vector<Vec3> positions;
for (int i = 0; i < numParticles; i++) {
params[0] = i;
force->addParticle(params);
positions.push_back(Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt)));
system.addParticle(1.0);
}
system.addForce(force);
VerletIntegrator integrator(0.001);
Context context(system, integrator, platform);
context.setPositions(positions);
// See if the energy is correct.
State state = context.getState(State::Energy);
double expectedEnergy = 0;
for (int i = 0; i < numParticles; i++)
for (int j = i+1; j < numParticles; j++)
for (int k = j+1; k < numParticles; k++) {
Vec3 d12 = positions[j]-positions[i];
Vec3 d13 = positions[k]-positions[i];
Vec3 d23 = positions[k]-positions[j];
double r12 = sqrt(d12.dot(d12));
double r13 = sqrt(d13.dot(d13));
double r23 = sqrt(d23.dot(d23));
expectedEnergy += 0.5*(r12+r13+r23)*(c[i]+c[j]+c[k]);
}
ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), 1e-5);
}
void testTypeFilters() {
// Create a system.
System system;
for (int i = 0; i < 5; i++)
system.addParticle(1.0);
CustomManyParticleForce* force = new CustomManyParticleForce(3, "c1*(distance(p1,p2)+distance(p1,p3))");
force->addPerParticleParameter("c");
double c[] = {1.0, 2.0, 1.3, 1.5, -2.1};
int type[] = {0, 1, 0, 1, 5};
vector<double> params(1);
for (int i = 0; i < 5; i++) {
params[0] = c[i];
force->addParticle(params, type[i]);
}
vector<Vec3> positions;
positions.push_back(Vec3(0, 0, 0));
positions.push_back(Vec3(1, 0, 0));
positions.push_back(Vec3(0, 1.1, 0.3));
positions.push_back(Vec3(0.4, 0, -0.8));
positions.push_back(Vec3(0.2, 0.5, -0.1));
set<int> f1, f2;
f1.insert(0);
f2.insert(1);
f2.insert(5);
force->setTypeFilter(0, f1);
force->setTypeFilter(1, f2);
force->setTypeFilter(2, f2);
system.addForce(force);
VerletIntegrator integrator(0.001);
Context context(system, integrator, platform);
context.setPositions(positions);
// See if the energy is correct.
State state = context.getState(State::Energy);
double expectedEnergy = 0;
int sets[6][3] = {{0,1,3}, {0,1,4}, {0,3,4}, {2,1,3}, {2,1,4}, {2,3,4}};
for (int i = 0; i < 6; i++) {
int p1 = sets[i][0];
int p2 = sets[i][1];
int p3 = sets[i][2];
Vec3 d12 = positions[p2]-positions[p1];
Vec3 d13 = positions[p3]-positions[p1];
double r12 = sqrt(d12.dot(d12));
double r13 = sqrt(d13.dot(d13));
expectedEnergy += c[p1]*(r12+r13);
}
ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), 1e-5);
}
int main(int argc, char* argv[]) {
try {
if (argc > 1)
platform.setPropertyDefaultValue("CudaPrecision", string(argv[1]));
testNoCutoff();
// testCutoff();
// testPeriodic();
// testExclusions();
// testAllTerms();
// testParameters();
// testTabulatedFunctions();
// testTypeFilters();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment