Skip to content

GitLab

Commit 81d60b1e authored by Lee-Ping Wang's avatar Lee-Ping Wang
Browse files

Merge branch 'master' of github.com:SimTk/openmm

parents d58cc041 68632d78
Hide whitespace changes
Inline Side-by-side
Showing with 1448 additions and 208 deletions
openmmapi/include/openmm/CustomNonbondedForce.h
View file @ 81d60b1e
......@@ -80,6 +80,25 @@ namespace OpenMM {
* if the labels 1 and 2 are reversed. In contrast, if it depended on the difference sigma1-sigma2, the results would
* be undefined, because reversing the labels 1 and 2 would change the energy.
*
* CustomNonbondedForce can operate in two modes. By default, it computes the interaction of every particle in the System
* with every other particle. Alternatively, you can restrict it to only a subset of particle pairs. To do this, specify
* one or more "interaction groups". An interaction group consists of two sets of particles that should interact with
* each other. Every particle in the first set interacts with every particle in the second set. For example, you might use
* this feature to compute a solute-solvent interaction energy, while omitting all interactions between two solute atoms
* or two solvent atoms.
*
* To create an interaction group, call addInteractionGroup(). You may add as many interaction groups as you want.
* Be aware of the following:
*
* <ul>
* <li>Exclusions are still taken into account, so the interactions between excluded pairs are omitted.</li>
* <li>Likewise, a particle will never interact with itself, even if it appears in both sets of an interaction group.</li>
* <li>If a particle pair appears in two different interaction groups, its interaction will be computed twice. This is
* sometimes useful, but be aware of it so you do not accidentally create unwanted duplicate interactions.</li>
* <li>If you do not add any interaction groups to a CustomNonbondedForce, it operates in the default mode where every
* particle interacts with every other particle.</li>
* </ul>
*
* When using a cutoff, by default the interaction is sharply truncated at the cutoff distance.
* Optionally you can instead use a switching function to make the interaction smoothly go to zero over a finite
* distance range. To enable this, call setUseSwitchingFunction(). You must also call setSwitchingDistance()
......@@ -167,6 +186,12 @@ public:
int getNumFunctions() const {
return functions.size();
}
/**
* Get the number of interaction groups that have been defined.
*/
int getNumInteractionGroups() const {
return interactionGroups.size();
}
/**
* Get the algebraic expression that gives the interaction energy between two particles
*/
......@@ -363,10 +388,34 @@ public:
* @param max the value of the independent variable corresponding to the last element of values
*/
void setFunctionParameters(int index, const std::string& name, const std::vector<double>& values, double min, double max);
/**
* Add an interaction group. An interaction will be computed between every particle in set1 and every particle in set2.
*
* @param set1 the first set of particles forming the interaction group
* @param set2 the second set of particles forming the interaction group
* @return the index of the interaction group that was added
*/
int addInteractionGroup(const std::set<int>& set1, const std::set<int>& set2);
/**
* Get the parameters for an interaction group.
*
* @param index the index of the interaction group for which to get parameters
* @param set1 the first set of particles forming the interaction group
* @param set2 the second set of particles forming the interaction group
*/
void getInteractionGroupParameters(int index, std::set<int>& set1, std::set<int>& set2) const;
/**
* Set the parameters for an interaction group.
*
* @param index the index of the interaction group for which to set parameters
* @param set1 the first set of particles forming the interaction group
* @param set2 the second set of particles forming the interaction group
*/
void setInteractionGroupParameters(int index, const std::set<int>& set1, const std::set<int>& set2);
/**
* Update the per-particle parameters in a Context to match those stored in this Force object. This method provides
* an efficient method to update certain parameters in an existing Context without needing to reinitialize it.
* Simply call setParticleParameters() to modify this object's parameters, then call updateParametersInState()
* Simply call setParticleParameters() to modify this object's parameters, then call updateParametersInContext()
* to copy them over to the Context.
*
* This method has several limitations. The only information it updates is the values of per-particle parameters.
......@@ -383,6 +432,7 @@ private:
class GlobalParameterInfo;
class ExclusionInfo;
class FunctionInfo;
class InteractionGroupInfo;
NonbondedMethod nonbondedMethod;
double cutoffDistance, switchingDistance;
bool useSwitchingFunction, useLongRangeCorrection;
......@@ -392,6 +442,7 @@ private:
std::vector<ParticleInfo> particles;
std::vector<ExclusionInfo> exclusions;
std::vector<FunctionInfo> functions;
std::vector<InteractionGroupInfo> interactionGroups;
};
/**
......@@ -465,6 +516,20 @@ public:
}
};
/**
* This is an internal class used to record information about an interaction group.
* @private
*/
class CustomNonbondedForce::InteractionGroupInfo {
public:
std::set<int> set1, set2;
InteractionGroupInfo() {
}
InteractionGroupInfo(const std::set<int>& set1, const std::set<int>& set2) :
set1(set1), set2(set2) {
}
};
} // namespace OpenMM
#endif /*OPENMM_CUSTOMNONBONDEDFORCE_H_*/
openmmapi/src/CustomNonbondedForce.cpp
View file @ 81d60b1e
......@@ -199,6 +199,23 @@ void CustomNonbondedForce::setFunctionParameters(int index, const std::string& n
functions[index].max = max;
}
int CustomNonbondedForce::addInteractionGroup(const std::set<int>& set1, const std::set<int>& set2) {
interactionGroups.push_back(InteractionGroupInfo(set1, set2));
return interactionGroups.size()-1;
}
void CustomNonbondedForce::getInteractionGroupParameters(int index, std::set<int>& set1, std::set<int>& set2) const {
ASSERT_VALID_INDEX(index, interactionGroups);
set1 = interactionGroups[index].set1;
set2 = interactionGroups[index].set2;
}
void CustomNonbondedForce::setInteractionGroupParameters(int index, const std::set<int>& set1, const std::set<int>& set2) {
ASSERT_VALID_INDEX(index, interactionGroups);
interactionGroups[index].set1 = set1;
interactionGroups[index].set2 = set2;
}
ForceImpl* CustomNonbondedForce::createImpl() const {
return new CustomNonbondedForceImpl(*this);
}
......
openmmapi/src/CustomNonbondedForceImpl.cpp
View file @ 81d60b1e
......@@ -42,14 +42,10 @@
#include "lepton/Parser.h"
#include <cmath>
#include <sstream>
#include <utility>
using namespace OpenMM;
using std::map;
using std::pair;
using std::vector;
using std::set;
using std::string;
using std::stringstream;
using namespace std;
CustomNonbondedForceImpl::CustomNonbondedForceImpl(const CustomNonbondedForce& owner) : owner(owner) {
}
......@@ -176,20 +172,6 @@ double CustomNonbondedForceImpl::calcLongRangeCorrection(const CustomNonbondedFo
if (force.getNonbondedMethod() == CustomNonbondedForce::NoCutoff || force.getNonbondedMethod() == CustomNonbondedForce::CutoffNonPeriodic)
return 0.0;
// Identify all particle classes (defined by parameters), and count the number of
// particles in each class.
map<vector<double>, int> classCounts;
for (int i = 0; i < force.getNumParticles(); i++) {
vector<double> parameters;
force.getParticleParameters(i, parameters);
map<vector<double>, int>::iterator entry = classCounts.find(parameters);
if (entry == classCounts.end())
classCounts[parameters] = 1;
else
entry->second++;
}
// Parse the energy expression.
map<string, Lepton::CustomFunction*> functions;
......@@ -201,20 +183,69 @@ double CustomNonbondedForceImpl::calcLongRangeCorrection(const CustomNonbondedFo
functions[name] = new TabulatedFunction(min, max, values);
}
Lepton::ExpressionProgram expression = Lepton::Parser::parse(force.getEnergyFunction(), functions).optimize().createProgram();
// Identify all particle classes (defined by parameters), and record the class of each particle.
int numParticles = force.getNumParticles();
vector<vector<double> > classes;
map<vector<double>, int> classIndex;
vector<int> atomClass(numParticles);
for (int i = 0; i < numParticles; i++) {
vector<double> parameters;
force.getParticleParameters(i, parameters);
if (classIndex.find(parameters) == classIndex.end()) {
classIndex[parameters] = classes.size();
classes.push_back(parameters);
}
atomClass[i] = classIndex[parameters];
}
int numClasses = classes.size();
// Count the total number of particle pairs for each pair of classes.
map<pair<int, int>, int> interactionCount;
if (force.getNumInteractionGroups() == 0) {
// Count the particles of each class.
vector<int> classCounts(numClasses, 0);
for (int i = 0; i < numParticles; i++)
classCounts[atomClass[i]]++;
for (int i = 0; i < numClasses; i++) {
interactionCount[make_pair(i, i)] = (classCounts[i]*(classCounts[i]+1))/2;
for (int j = i+1; j < numClasses; j++)
interactionCount[make_pair(i, j)] = classCounts[i]*classCounts[j];
}
}
else {
// Initialize the counts to 0.
for (int i = 0; i < numClasses; i++) {
for (int j = i; j < numClasses; j++)
interactionCount[make_pair(i, j)] = 0;
}
// Loop over interaction groups and count the interactions in each one.
for (int group = 0; group < force.getNumInteractionGroups(); group++) {
set<int> set1, set2;
force.getInteractionGroupParameters(group, set1, set2);
for (set<int>::const_iterator a1 = set1.begin(); a1 != set1.end(); ++a1)
for (set<int>::const_iterator a2 = set2.begin(); a2 != set2.end(); ++a2) {
if (*a1 >= *a2 && set1.find(*a2) != set1.end() && set2.find(*a1) != set2.end())
continue;
int class1 = atomClass[*a1];
int class2 = atomClass[*a2];
interactionCount[make_pair(min(class1, class2), max(class1, class2))]++;
}
}
}
// Loop over all pairs of classes to compute the coefficient.
double sum = 0;
for (map<vector<double>, int>::const_iterator entry = classCounts.begin(); entry != classCounts.end(); ++entry) {
int count = (entry->second*(entry->second+1))/2;
sum += count*integrateInteraction(expression, entry->first, entry->first, force, context);
}
for (map<vector<double>, int>::const_iterator class1 = classCounts.begin(); class1 != classCounts.end(); ++class1)
for (map<vector<double>, int>::const_iterator class2 = classCounts.begin(); class2 != class1; ++class2) {
int count = class1->second*class2->second;
sum += count*integrateInteraction(expression, class1->first, class2->first, force, context);
}
int numParticles = force.getNumParticles();
for (int i = 0; i < numClasses; i++)
for (int j = i; j < numClasses; j++)
sum += interactionCount[make_pair(i, j)]*integrateInteraction(expression, classes[i], classes[j], force, context);
int numInteractions = (numParticles*(numParticles+1))/2;
sum /= numInteractions;
return 2*M_PI*numParticles*numParticles*sum;
......
platforms/cuda/include/CudaKernels.h
View file @ 81d60b1e
......@@ -638,7 +638,7 @@ private:
class CudaCalcCustomNonbondedForceKernel : public CalcCustomNonbondedForceKernel {
public:
CudaCalcCustomNonbondedForceKernel(std::string name, const Platform& platform, CudaContext& cu, const System& system) : CalcCustomNonbondedForceKernel(name, platform),
cu(cu), params(NULL), globals(NULL), tabulatedFunctionParams(NULL), forceCopy(NULL), system(system) {
cu(cu), params(NULL), globals(NULL), tabulatedFunctionParams(NULL), interactionGroupData(NULL), forceCopy(NULL), system(system), hasInitializedKernel(false) {
}
~CudaCalcCustomNonbondedForceKernel();
/**
......@@ -665,15 +665,20 @@ public:
*/
void copyParametersToContext(ContextImpl& context, const CustomNonbondedForce& force);
private:
void initInteractionGroups(const CustomNonbondedForce& force, const std::string& interactionSource);
CudaContext& cu;
CudaParameterSet* params;
CudaArray* globals;
CudaArray* tabulatedFunctionParams;
CudaArray* interactionGroupData;
CUfunction interactionGroupKernel;
std::vector<void*> interactionGroupArgs;
std::vector<std::string> globalParamNames;
std::vector<float> globalParamValues;
std::vector<CudaArray*> tabulatedFunctions;
double longRangeCoefficient;
bool hasInitializedLongRangeCorrection;
bool hasInitializedLongRangeCorrection, hasInitializedKernel;
int numGroupThreadBlocks;
CustomNonbondedForce* forceCopy;
const System& system;
};
......
platforms/cuda/src/CudaKernels.cpp
View file @ 81d60b1e
......@@ -46,6 +46,7 @@
#include "lepton/ParsedExpression.h"
#include "SimTKOpenMMRealType.h"
#include "SimTKOpenMMUtilities.h"
#include <algorithm>
#include <cmath>
#include <set>
......@@ -1862,6 +1863,17 @@ void CudaCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& context,
class CudaCustomNonbondedForceInfo : public CudaForceInfo {
public:
CudaCustomNonbondedForceInfo(const CustomNonbondedForce& force) : force(force) {
if (force.getNumInteractionGroups() > 0) {
groupsForParticle.resize(force.getNumParticles());
for (int i = 0; i < force.getNumInteractionGroups(); i++) {
set<int> set1, set2;
force.getInteractionGroupParameters(i, set1, set2);
for (set<int>::const_iterator iter = set1.begin(); iter != set1.end(); ++iter)
groupsForParticle[*iter].insert(2*i);
for (set<int>::const_iterator iter = set2.begin(); iter != set2.end(); ++iter)
groupsForParticle[*iter].insert(2*i+1);
}
}
}
bool areParticlesIdentical(int particle1, int particle2) {
vector<double> params1;
......@@ -1871,6 +1883,8 @@ public:
for (int i = 0; i < (int) params1.size(); i++)
if (params1[i] != params2[i])
return false;
if (groupsForParticle.size() > 0 && groupsForParticle[particle1] != groupsForParticle[particle2])
return false;
return true;
}
int getNumParticleGroups() {
......@@ -1888,6 +1902,7 @@ public:
}
private:
const CustomNonbondedForce& force;
vector<set<int> > groupsForParticle;
};
CudaCalcCustomNonbondedForceKernel::~CudaCalcCustomNonbondedForceKernel() {
......@@ -1898,6 +1913,8 @@ CudaCalcCustomNonbondedForceKernel::~CudaCalcCustomNonbondedForceKernel() {
delete globals;
if (tabulatedFunctionParams != NULL)
delete tabulatedFunctionParams;
if (interactionGroupData != NULL)
delete interactionGroupData;
for (int i = 0; i < (int) tabulatedFunctions.size(); i++)
delete tabulatedFunctions[i];
if (forceCopy != NULL)
......@@ -1909,7 +1926,7 @@ void CudaCalcCustomNonbondedForceKernel::initialize(const System& system, const
int forceIndex;
for (forceIndex = 0; forceIndex < system.getNumForces() && &system.getForce(forceIndex) != &force; ++forceIndex)
;
string prefix = "custom"+cu.intToString(forceIndex)+"_";
string prefix = (force.getNumInteractionGroups() == 0 ? "custom"+cu.intToString(forceIndex)+"_" : "");
// Record parameters and exclusions.
......@@ -2010,14 +2027,18 @@ void CudaCalcCustomNonbondedForceKernel::initialize(const System& system, const
replacements["SWITCH_C5"] = cu.doubleToString(6/pow(force.getSwitchingDistance()-force.getCutoffDistance(), 5.0));
}
string source = cu.replaceStrings(CudaKernelSources::customNonbonded, replacements);
cu.getNonbondedUtilities().addInteraction(useCutoff, usePeriodic, true, force.getCutoffDistance(), exclusionList, source, force.getForceGroup());
for (int i = 0; i < (int) params->getBuffers().size(); i++) {
CudaNonbondedUtilities::ParameterInfo& buffer = params->getBuffers()[i];
cu.getNonbondedUtilities().addParameter(CudaNonbondedUtilities::ParameterInfo(prefix+"params"+cu.intToString(i+1), buffer.getComponentType(), buffer.getNumComponents(), buffer.getSize(), buffer.getMemory()));
}
if (globals != NULL) {
globals->upload(globalParamValues);
cu.getNonbondedUtilities().addArgument(CudaNonbondedUtilities::ParameterInfo(prefix+"globals", "float", 1, sizeof(float), globals->getDevicePointer()));
if (force.getNumInteractionGroups() > 0)
initInteractionGroups(force, source);
else {
cu.getNonbondedUtilities().addInteraction(useCutoff, usePeriodic, true, force.getCutoffDistance(), exclusionList, source, force.getForceGroup());
for (int i = 0; i < (int) params->getBuffers().size(); i++) {
CudaNonbondedUtilities::ParameterInfo& buffer = params->getBuffers()[i];
cu.getNonbondedUtilities().addParameter(CudaNonbondedUtilities::ParameterInfo(prefix+"params"+cu.intToString(i+1), buffer.getComponentType(), buffer.getNumComponents(), buffer.getSize(), buffer.getMemory()));
}
if (globals != NULL) {
globals->upload(globalParamValues);
cu.getNonbondedUtilities().addArgument(CudaNonbondedUtilities::ParameterInfo(prefix+"globals", "float", 1, sizeof(float), globals->getDevicePointer()));
}
}
cu.addForce(new CudaCustomNonbondedForceInfo(force));
......@@ -2033,6 +2054,242 @@ void CudaCalcCustomNonbondedForceKernel::initialize(const System& system, const
}
}
void CudaCalcCustomNonbondedForceKernel::initInteractionGroups(const CustomNonbondedForce& force, const string& interactionSource) {
// Process groups to form tiles.
vector<vector<int> > atomLists;
vector<pair<int, int> > tiles;
map<pair<int, int>, int> duplicateInteractions;
for (int group = 0; group < force.getNumInteractionGroups(); group++) {
// Get the list of atoms in this group and sort them.
set<int> set1, set2;
force.getInteractionGroupParameters(group, set1, set2);
vector<int> atoms1, atoms2;
atoms1.insert(atoms1.begin(), set1.begin(), set1.end());
atoms2.insert(atoms2.begin(), set2.begin(), set2.end());
sort(atoms1.begin(), atoms1.end());
sort(atoms2.begin(), atoms2.end());
// Find how many tiles we will create for this group.
int tileWidth = min(min(32, (int) atoms1.size()), (int) atoms2.size());
int numBlocks1 = (atoms1.size()+tileWidth-1)/tileWidth;
int numBlocks2 = (atoms2.size()+tileWidth-1)/tileWidth;
// Add the tiles.
for (int i = 0; i < numBlocks1; i++)
for (int j = 0; j < numBlocks2; j++)
tiles.push_back(make_pair(atomLists.size()+i, atomLists.size()+numBlocks1+j));
// Add the atom lists.
for (int i = 0; i < numBlocks1; i++) {
vector<int> atoms;
int first = i*tileWidth;
int last = min((i+1)*tileWidth, (int) atoms1.size());
for (int j = first; j < last; j++)
atoms.push_back(atoms1[j]);
atomLists.push_back(atoms);
}
for (int i = 0; i < numBlocks2; i++) {
vector<int> atoms;
int first = i*tileWidth;
int last = min((i+1)*tileWidth, (int) atoms2.size());
for (int j = first; j < last; j++)
atoms.push_back(atoms2[j]);
atomLists.push_back(atoms);
}
// If this group contains duplicate interactions, record that we need to skip them once.
for (int i = 0; i < (int) atoms1.size(); i++) {
int a1 = atoms1[i];
if (set2.find(a1) == set2.end())
continue;
for (int j = 0; j < (int) atoms2.size() && atoms2[j] < a1; j++) {
int a2 = atoms2[j];
if (set1.find(a2) != set1.end()) {
pair<int, int> key = make_pair(a2, a1);
if (duplicateInteractions.find(key) == duplicateInteractions.end())
duplicateInteractions[key] = 0;
duplicateInteractions[key]++;
}
}
}
}
// Build a lookup table for quickly identifying excluded interactions.
set<pair<int, int> > exclusions;
for (int i = 0; i < force.getNumExclusions(); i++) {
int p1, p2;
force.getExclusionParticles(i, p1, p2);
exclusions.insert(make_pair(min(p1, p2), max(p1, p2)));
}
// Build the exclusion flags for each tile. While we're at it, filter out tiles
// where all interactions are excluded, and sort the tiles by size.
vector<vector<int> > exclusionFlags(tiles.size());
vector<pair<int, int> > tileOrder;
for (int tile = 0; tile < tiles.size(); tile++) {
if (atomLists[tiles[tile].first].size() < atomLists[tiles[tile].second].size()) {
// For efficiency, we want the first axis to be the larger one.
int swap = tiles[tile].first;
tiles[tile].first = tiles[tile].second;
tiles[tile].second = swap;
}
vector<int>& atoms1 = atomLists[tiles[tile].first];
vector<int>& atoms2 = atomLists[tiles[tile].second];
vector<int> flags(atoms1.size(), (int) (1LL<<atoms2.size())-1);
int numExcluded = 0;
for (int i = 0; i < (int) atoms1.size(); i++)
for (int j = 0; j < (int) atoms2.size(); j++) {
int a1 = atoms1[i];
int a2 = atoms2[j];
bool isExcluded = false;
pair<int, int> key = make_pair(min(a1, a2), max(a1, a2));
if (a1 == a2 || exclusions.find(key) != exclusions.end())
isExcluded = true; // This is an excluded interaction.
else if (duplicateInteractions.find(key) != duplicateInteractions.end() && duplicateInteractions[key] > 0) {
// Both atoms are in both sets, so skip duplicate interactions.
isExcluded = true;
duplicateInteractions[key]--;
}
if (isExcluded) {
flags[i] &= -1-(1<<j);
numExcluded++;
}
}
if (numExcluded == atoms1.size()*atoms2.size())
continue; // All interactions are excluded.
tileOrder.push_back(make_pair((int) -atoms2.size(), tile));
exclusionFlags[tile] = flags;
}
sort(tileOrder.begin(), tileOrder.end());
// Merge tiles to get as close as possible to 32 along the first axis of each one.
vector<int> tileSetStart;
tileSetStart.push_back(0);
int tileSetSize = 0;
for (int i = 0; i < tileOrder.size(); i++) {
int tile = tileOrder[i].second;
int size = atomLists[tiles[tile].first].size();
if (tileSetSize+size > 32) {
tileSetStart.push_back(i);
tileSetSize = 0;
}
tileSetSize += size;
}
tileSetStart.push_back(tileOrder.size());
// Build the data structures.
int numTileSets = tileSetStart.size()-1;
vector<int4> groupData;
for (int tileSet = 0; tileSet < numTileSets; tileSet++) {
int indexInTileSet = 0;
for (int i = tileSetStart[tileSet]; i < tileSetStart[tileSet+1]; i++) {
int tile = tileOrder[i].second;
vector<int>& atoms1 = atomLists[tiles[tile].first];
vector<int>& atoms2 = atomLists[tiles[tile].second];
int range = indexInTileSet + ((indexInTileSet+atoms1.size())<<16);
int allFlags = (1<<atoms2.size())-1;
for (int j = 0; j < (int) atoms1.size(); j++) {
int a1 = atoms1[j];
int a2 = (j < atoms2.size() ? atoms2[j] : 0);
int flags = (exclusionFlags[tile].size() > 0 ? exclusionFlags[tile][j] : allFlags);
groupData.push_back(make_int4(a1, a2, range, flags<<indexInTileSet));
}
indexInTileSet += atoms1.size();
}
for (; indexInTileSet < 32; indexInTileSet++)
groupData.push_back(make_int4(0, 0, 0, 0));
}
interactionGroupData = CudaArray::create<int4>(cu, groupData.size(), "interactionGroupData");
interactionGroupData->upload(groupData);
// Create the kernel.
map<string, string> replacements;
replacements["COMPUTE_INTERACTION"] = interactionSource;
const string suffixes[] = {"x", "y", "z", "w"};
stringstream localData;
int localDataSize = 0;
vector<CudaNonbondedUtilities::ParameterInfo>& buffers = params->getBuffers();
for (int i = 0; i < (int) buffers.size(); i++) {
if (buffers[i].getNumComponents() == 1)
localData<<buffers[i].getComponentType()<<" params"<<(i+1)<<";\n";
else {
for (int j = 0; j < buffers[i].getNumComponents(); ++j)
localData<<buffers[i].getComponentType()<<" params"<<(i+1)<<"_"<<suffixes[j]<<";\n";
}
localDataSize += buffers[i].getSize();
}
replacements["ATOM_PARAMETER_DATA"] = localData.str();
stringstream args;
for (int i = 0; i < (int) buffers.size(); i++)
args<<", const "<<buffers[i].getType()<<"* __restrict__ global_params"<<(i+1);
if (globals != NULL)
args<<", const float* __restrict__ globals";
replacements["PARAMETER_ARGUMENTS"] = args.str();
stringstream load1;
for (int i = 0; i < (int) buffers.size(); i++)
load1<<buffers[i].getType()<<" params"<<(i+1)<<"1 = global_params"<<(i+1)<<"[atom1];\n";
replacements["LOAD_ATOM1_PARAMETERS"] = load1.str();
stringstream loadLocal2;
for (int i = 0; i < (int) buffers.size(); i++) {
if (buffers[i].getNumComponents() == 1)
loadLocal2<<"localData[threadIdx.x].params"<<(i+1)<<" = global_params"<<(i+1)<<"[atom2];\n";
else {
loadLocal2<<buffers[i].getType()<<" temp_params"<<(i+1)<<" = global_params"<<(i+1)<<"[atom2];\n";
for (int j = 0; j < buffers[i].getNumComponents(); ++j)
loadLocal2<<"localData[threadIdx.x].params"<<(i+1)<<"_"<<suffixes[j]<<" = temp_params"<<(i+1)<<"."<<suffixes[j]<<";\n";
}
}
replacements["LOAD_LOCAL_PARAMETERS"] = loadLocal2.str();
stringstream load2;
for (int i = 0; i < (int) buffers.size(); i++) {
if (buffers[i].getNumComponents() == 1)
load2<<buffers[i].getType()<<" params"<<(i+1)<<"2 = localData[localIndex].params"<<(i+1)<<";\n";
else {
load2<<buffers[i].getType()<<" params"<<(i+1)<<"2 = make_"<<buffers[i].getType()<<"(";
for (int j = 0; j < buffers[i].getNumComponents(); ++j) {
if (j > 0)
load2<<", ";
load2<<"localData[localIndex].params"<<(i+1)<<"_"<<suffixes[j];
}
load2<<");\n";
}
}
replacements["LOAD_ATOM2_PARAMETERS"] = load2.str();
map<string, string> defines;
if (force.getNonbondedMethod() != CustomNonbondedForce::NoCutoff)
defines["USE_CUTOFF"] = "1";
if (force.getNonbondedMethod() == CustomNonbondedForce::CutoffPeriodic)
defines["USE_PERIODIC"] = "1";
defines["THREAD_BLOCK_SIZE"] = cu.intToString(cu.getNonbondedUtilities().getForceThreadBlockSize());
double cutoff = force.getCutoffDistance();
defines["CUTOFF_SQUARED"] = cu.doubleToString(cutoff*cutoff);
defines["PADDED_NUM_ATOMS"] = cu.intToString(cu.getPaddedNumAtoms());
defines["TILE_SIZE"] = "32";
int numContexts = cu.getPlatformData().contexts.size();
int startIndex = cu.getContextIndex()*numTileSets/numContexts;
int endIndex = (cu.getContextIndex()+1)*numTileSets/numContexts;
defines["FIRST_TILE"] = cu.intToString(startIndex);
defines["LAST_TILE"] = cu.intToString(endIndex);
if ((localDataSize/4)%2 == 0 && !cu.getUseDoublePrecision())
defines["PARAMETER_SIZE_IS_EVEN"] = "1";
CUmodule program = cu.createModule(CudaKernelSources::vectorOps+cu.replaceStrings(CudaKernelSources::customNonbondedGroups, replacements), defines);
interactionGroupKernel = cu.getKernel(program, "computeInteractionGroups");
numGroupThreadBlocks = cu.getNonbondedUtilities().getNumForceThreadBlocks();
}
double CudaCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
if (globals != NULL) {
bool changed = false;
......@@ -2054,6 +2311,23 @@ double CudaCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bool in
longRangeCoefficient = CustomNonbondedForceImpl::calcLongRangeCorrection(*forceCopy, context.getOwner());
hasInitializedLongRangeCorrection = true;
}
if (interactionGroupData != NULL) {
if (!hasInitializedKernel) {
hasInitializedKernel = true;
interactionGroupArgs.push_back(&cu.getForce().getDevicePointer());
interactionGroupArgs.push_back(&cu.getEnergyBuffer().getDevicePointer());
interactionGroupArgs.push_back(&cu.getPosq().getDevicePointer());
interactionGroupArgs.push_back(&interactionGroupData->getDevicePointer());
interactionGroupArgs.push_back(cu.getPeriodicBoxSizePointer());
interactionGroupArgs.push_back(cu.getInvPeriodicBoxSizePointer());
for (int i = 0; i < (int) params->getBuffers().size(); i++)
interactionGroupArgs.push_back(&params->getBuffers()[i].getMemory());
if (globals != NULL)
interactionGroupArgs.push_back(&globals->getDevicePointer());
}
int forceThreadBlockSize = cu.getNonbondedUtilities().getForceThreadBlockSize();
cu.executeKernel(interactionGroupKernel, &interactionGroupArgs[0], numGroupThreadBlocks*forceThreadBlockSize, forceThreadBlockSize);
}
double4 boxSize = cu.getPeriodicBoxSize();
return longRangeCoefficient/(boxSize.x*boxSize.y*boxSize.z);
}
......
platforms/cuda/src/kernels/customNonbondedGroups.cu 0 → 100644
View file @ 81d60b1e
#define WARPS_PER_GROUP (THREAD_BLOCK_SIZE/TILE_SIZE)
typedef struct {
real x, y, z;
real q;
real fx, fy, fz;
ATOM_PARAMETER_DATA
#ifndef PARAMETER_SIZE_IS_EVEN
real padding;
#endif
} AtomData;
extern "C" __global__ void computeInteractionGroups(
unsigned long long* __restrict__ forceBuffers, real* __restrict__ energyBuffer, const real4* __restrict__ posq, const int4* __restrict__ groupData,
real4 periodicBoxSize, real4 invPeriodicBoxSize
PARAMETER_ARGUMENTS) {
const unsigned int totalWarps = (blockDim.x*gridDim.x)/TILE_SIZE;
const unsigned int warp = (blockIdx.x*blockDim.x+threadIdx.x)/TILE_SIZE; // global warpIndex
const unsigned int tgx = threadIdx.x & (TILE_SIZE-1); // index within the warp
const unsigned int tbx = threadIdx.x - tgx; // block warpIndex
real energy = 0.0f;
__shared__ AtomData localData[THREAD_BLOCK_SIZE];
const unsigned int startTile = FIRST_TILE+warp*(LAST_TILE-FIRST_TILE)/totalWarps;
const unsigned int endTile = FIRST_TILE+(warp+1)*(LAST_TILE-FIRST_TILE)/totalWarps;
for (int tile = startTile; tile < endTile; tile++) {
const int4 atomData = groupData[TILE_SIZE*tile+tgx];
const int atom1 = atomData.x;
const int atom2 = atomData.y;
const int rangeStart = atomData.z&0xFFFF;
const int rangeEnd = (atomData.z>>16)&0xFFFF;
const int exclusions = atomData.w;
real4 posq1 = posq[atom1];
LOAD_ATOM1_PARAMETERS
real3 force = make_real3(0);
real4 posq2 = posq[atom2];
localData[threadIdx.x].x = posq2.x;
localData[threadIdx.x].y = posq2.y;
localData[threadIdx.x].z = posq2.z;
localData[threadIdx.x].q = posq2.w;
LOAD_LOCAL_PARAMETERS
localData[threadIdx.x].fx = 0.0f;
localData[threadIdx.x].fy = 0.0f;
localData[threadIdx.x].fz = 0.0f;
int tj = tgx;
for (int j = rangeStart; j < rangeEnd; j++) {
bool isExcluded = (((exclusions>>tj)&1) == 0);
int localIndex = tbx+tj;
posq2 = make_real4(localData[localIndex].x, localData[localIndex].y, localData[localIndex].z, localData[localIndex].q);
real3 delta = make_real3(posq2.x-posq1.x, posq2.y-posq1.y, posq2.z-posq1.z);
#ifdef USE_PERIODIC
delta.x -= floor(delta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
delta.y -= floor(delta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
#endif
real r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
#ifdef USE_CUTOFF
if (!isExcluded && r2 < CUTOFF_SQUARED) {
#endif
real invR = RSQRT(r2);
real r = RECIP(invR);
LOAD_ATOM2_PARAMETERS
real dEdR = 0.0f;
real tempEnergy = 0.0f;
COMPUTE_INTERACTION
energy += tempEnergy;
delta *= dEdR;
force.x -= delta.x;
force.y -= delta.y;
force.z -= delta.z;
localData[localIndex].fx += delta.x;
localData[localIndex].fy += delta.y;
localData[localIndex].fz += delta.z;
#ifdef USE_CUTOFF
}
#endif
tj = (tj == rangeEnd-1 ? rangeStart : tj+1);
}
if (exclusions != 0) {
atomicAdd(&forceBuffers[atom1], static_cast<unsigned long long>((long long) (force.x*0x100000000)));
atomicAdd(&forceBuffers[atom1+PADDED_NUM_ATOMS], static_cast<unsigned long long>((long long) (force.y*0x100000000)));
atomicAdd(&forceBuffers[atom1+2*PADDED_NUM_ATOMS], static_cast<unsigned long long>((long long) (force.z*0x100000000)));
}
atomicAdd(&forceBuffers[atom2], static_cast<unsigned long long>((long long) (localData[threadIdx.x].fx*0x100000000)));
atomicAdd(&forceBuffers[atom2+PADDED_NUM_ATOMS], static_cast<unsigned long long>((long long) (localData[threadIdx.x].fy*0x100000000)));
atomicAdd(&forceBuffers[atom2+2*PADDED_NUM_ATOMS], static_cast<unsigned long long>((long long) (localData[threadIdx.x].fz*0x100000000)));
}
energyBuffer[blockIdx.x*blockDim.x+threadIdx.x] += energy;
}
\ No newline at end of file
platforms/cuda/tests/TestCudaCustomNonbondedForce.cpp
View file @ 81d60b1e
......@@ -34,6 +34,9 @@
* This tests all the different force terms in the CUDA implementation of CustomNonbondedForce.
*/
#ifdef WIN32
#define _USE_MATH_DEFINES // Needed to get M_PI
#endif
#include "openmm/internal/AssertionUtilities.h"
#include "sfmt/SFMT.h"
#include "openmm/Context.h"
......@@ -42,6 +45,7 @@
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include <cmath>
#include <iostream>
#include <vector>
......@@ -538,6 +542,179 @@ void testLongRangeCorrection() {
ASSERT_EQUAL_TOL(standardEnergy1-standardEnergy2, customEnergy1-customEnergy2, 1e-4);
}
void testInteractionGroups() {
const int numParticles = 6;
System system;
VerletIntegrator integrator(0.01);
CustomNonbondedForce* nonbonded = new CustomNonbondedForce("v1+v2");
nonbonded->addPerParticleParameter("v");
vector<double> params(1, 0.001);
for (int i = 0; i < numParticles; i++) {
system.addParticle(1.0);
nonbonded->addParticle(params);
params[0] *= 10;
}
set<int> set1, set2, set3, set4;
set1.insert(2);
set2.insert(0);
set2.insert(1);
set2.insert(2);
set2.insert(3);
set2.insert(4);
set2.insert(5);
nonbonded->addInteractionGroup(set1, set2); // Particle 2 interacts with every other particle.
set3.insert(0);
set3.insert(1);
set4.insert(4);
set4.insert(5);
nonbonded->addInteractionGroup(set3, set4); // Particles 0 and 1 interact with 4 and 5.
nonbonded->addExclusion(1, 2); // Add an exclusion to make sure it gets skipped.
system.addForce(nonbonded);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
context.setPositions(positions);
State state = context.getState(State::Energy);
double expectedEnergy = 331.423; // Each digit is the number of interactions a particle particle is involved in.
ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), TOL);
}
void testLargeInteractionGroup() {
const int numMolecules = 300;
const int numParticles = numMolecules*2;
const double boxSize = 20.0;
// Create a large system.
System system;
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
for (int i = 0; i < numParticles; i++)
system.addParticle(1.0);
CustomNonbondedForce* nonbonded = new CustomNonbondedForce("4*eps*((sigma/r)^12-(sigma/r)^6)+138.935456*q/r; q=q1*q2; sigma=0.5*(sigma1+sigma2); eps=sqrt(eps1*eps2)");
nonbonded->addPerParticleParameter("q");
nonbonded->addPerParticleParameter("sigma");
nonbonded->addPerParticleParameter("eps");
vector<Vec3> positions(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<double> params(3);
for (int i = 0; i < numMolecules; i++) {
if (i < numMolecules/2) {
params[0] = 1.0;
params[1] = 0.2;
params[2] = 0.1;
nonbonded->addParticle(params);
params[0] = -1.0;
params[1] = 0.1;
nonbonded->addParticle(params);
}
else {
params[0] = 1.0;
params[1] = 0.2;
params[2] = 0.2;
nonbonded->addParticle(params);
params[0] = -1.0;
params[1] = 0.1;
nonbonded->addParticle(params);
}
positions[2*i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
positions[2*i+1] = Vec3(positions[2*i][0]+1.0, positions[2*i][1], positions[2*i][2]);
nonbonded->addExclusion(2*i, 2*i+1);
}
nonbonded->setNonbondedMethod(CustomNonbondedForce::CutoffPeriodic);
system.addForce(nonbonded);
// Compute the forces.
VerletIntegrator integrator(0.01);
Context context(system, integrator, platform);
context.setPositions(positions);
State state1 = context.getState(State::Forces);
// Modify the force so only one particle interacts with everything else.
set<int> set1, set2;
set1.insert(151);
for (int i = 0; i < numParticles; i++)
set2.insert(i);
nonbonded->addInteractionGroup(set1, set2);
context.reinitialize();
context.setPositions(positions);
State state2 = context.getState(State::Forces);
// The force on that one particle should be the same.
ASSERT_EQUAL_VEC(state1.getForces()[151], state2.getForces()[151], 1e-4);
// Modify the interaction group so it includes all interactions. This should now reproduce the original forces
// on all atoms.
for (int i = 0; i < numParticles; i++)
set1.insert(i);
nonbonded->setInteractionGroupParameters(0, set1, set2);
context.reinitialize();
context.setPositions(positions);
State state3 = context.getState(State::Forces);
for (int i = 0; i < numParticles; i++)
ASSERT_EQUAL_VEC(state1.getForces()[i], state3.getForces()[i], 1e-4);
}
void testInteractionGroupLongRangeCorrection() {
const int numParticles = 10;
const double boxSize = 10.0;
const double cutoff = 0.5;
System system;
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
CustomNonbondedForce* nonbonded = new CustomNonbondedForce("c1*c2*r^-4");
nonbonded->addPerParticleParameter("c");
vector<Vec3> positions(numParticles);
vector<double> params(1);
for (int i = 0; i < numParticles; i++) {
system.addParticle(1.0);
params[0] = (i%2 == 0 ? 1.1 : 2.0);
nonbonded->addParticle(params);
positions[i] = Vec3(0.5*i, 0, 0);
}
nonbonded->setNonbondedMethod(CustomNonbondedForce::CutoffPeriodic);
nonbonded->setCutoffDistance(cutoff);
system.addForce(nonbonded);
// Setup nonbonded groups. They involve 1 interaction of type AA,
// 2 of type BB, and 5 of type AB.
set<int> set1, set2, set3, set4, set5;
set1.insert(0);
set1.insert(1);
set1.insert(2);
nonbonded->addInteractionGroup(set1, set1);
set2.insert(3);
set3.insert(4);
set3.insert(6);
set3.insert(8);
nonbonded->addInteractionGroup(set2, set3);
set4.insert(5);
set5.insert(7);
set5.insert(9);
nonbonded->addInteractionGroup(set4, set5);
// Compute energy with and without the correction.
VerletIntegrator integrator(0.01);
Context context(system, integrator, platform);
context.setPositions(positions);
double energy1 = context.getState(State::Energy).getPotentialEnergy();
nonbonded->setUseLongRangeCorrection(true);
context.reinitialize();
context.setPositions(positions);
double energy2 = context.getState(State::Energy).getPotentialEnergy();
// Check the result.
double sum = (1.1*1.1 + 2*2.0*2.0 + 5*1.1*2.0)*2.0;
int numPairs = (numParticles*(numParticles+1))/2;
double expected = 2*M_PI*numParticles*numParticles*sum/(numPairs*boxSize*boxSize*boxSize);
ASSERT_EQUAL_TOL(expected, energy2-energy1, 1e-4);
}
int main(int argc, char* argv[]) {
try {
if (argc > 1)
......@@ -553,6 +730,9 @@ int main(int argc, char* argv[]) {
testParallelComputation();
testSwitchingFunction();
testLongRangeCorrection();
testInteractionGroups();
testLargeInteractionGroup();
testInteractionGroupLongRangeCorrection();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
......
platforms/opencl/include/OpenCLKernels.h
View file @ 81d60b1e
......@@ -639,7 +639,7 @@ private:
class OpenCLCalcCustomNonbondedForceKernel : public CalcCustomNonbondedForceKernel {
public:
OpenCLCalcCustomNonbondedForceKernel(std::string name, const Platform& platform, OpenCLContext& cl, const System& system) : CalcCustomNonbondedForceKernel(name, platform),
cl(cl), params(NULL), globals(NULL), tabulatedFunctionParams(NULL), forceCopy(NULL), system(system) {
cl(cl), params(NULL), globals(NULL), tabulatedFunctionParams(NULL), interactionGroupData(NULL), forceCopy(NULL), system(system), hasInitializedKernel(false) {
}
~OpenCLCalcCustomNonbondedForceKernel();
/**
......@@ -666,15 +666,20 @@ public:
*/
void copyParametersToContext(ContextImpl& context, const CustomNonbondedForce& force);
private:
void initInteractionGroups(const CustomNonbondedForce& force, const std::string& interactionSource);
OpenCLContext& cl;
OpenCLParameterSet* params;
OpenCLArray* globals;
OpenCLArray* tabulatedFunctionParams;
OpenCLArray* interactionGroupData;
cl::Kernel interactionGroupKernel;
std::vector<void*> interactionGroupArgs;
std::vector<std::string> globalParamNames;
std::vector<cl_float> globalParamValues;
std::vector<OpenCLArray*> tabulatedFunctions;
double longRangeCoefficient;
bool hasInitializedLongRangeCorrection;
bool hasInitializedLongRangeCorrection, hasInitializedKernel;
int numGroupThreadBlocks;
CustomNonbondedForce* forceCopy;
const System& system;
};
......
platforms/opencl/include/OpenCLParameterSet.h
View file @ 81d60b1e
......@@ -81,6 +81,13 @@ public:
*/
template <class T>
void setParameterValues(const std::vector<std::vector<T> >& values);
/**
* Get a set of OpenCLNonbondedUtilities::ParameterInfo objects which describe the Buffers
* containing the data.
*/
std::vector<OpenCLNonbondedUtilities::ParameterInfo>& getBuffers() {
return buffers;
}
/**
* Get a set of OpenCLNonbondedUtilities::ParameterInfo objects which describe the Buffers
* containing the data.
......
platforms/opencl/src/OpenCLKernels.cpp
View file @ 81d60b1e
......@@ -46,6 +46,7 @@
#include "lepton/ParsedExpression.h"
#include "SimTKOpenMMRealType.h"
#include "SimTKOpenMMUtilities.h"
#include <algorithm>
#include <cmath>
#include <set>
......@@ -1875,6 +1876,17 @@ void OpenCLCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& contex
class OpenCLCustomNonbondedForceInfo : public OpenCLForceInfo {
public:
OpenCLCustomNonbondedForceInfo(int requiredBuffers, const CustomNonbondedForce& force) : OpenCLForceInfo(requiredBuffers), force(force) {
if (force.getNumInteractionGroups() > 0) {
groupsForParticle.resize(force.getNumParticles());
for (int i = 0; i < force.getNumInteractionGroups(); i++) {
set<int> set1, set2;
force.getInteractionGroupParameters(i, set1, set2);
for (set<int>::const_iterator iter = set1.begin(); iter != set1.end(); ++iter)
groupsForParticle[*iter].insert(2*i);
for (set<int>::const_iterator iter = set2.begin(); iter != set2.end(); ++iter)
groupsForParticle[*iter].insert(2*i+1);
}
}
}
bool areParticlesIdentical(int particle1, int particle2) {
vector<double> params1;
......@@ -1884,6 +1896,8 @@ public:
for (int i = 0; i < (int) params1.size(); i++)
if (params1[i] != params2[i])
return false;
if (groupsForParticle.size() > 0 && groupsForParticle[particle1] != groupsForParticle[particle2])
return false;
return true;
}
int getNumParticleGroups() {
......@@ -1901,6 +1915,7 @@ public:
}
private:
const CustomNonbondedForce& force;
vector<set<int> > groupsForParticle;
};
OpenCLCalcCustomNonbondedForceKernel::~OpenCLCalcCustomNonbondedForceKernel() {
......@@ -1910,6 +1925,8 @@ OpenCLCalcCustomNonbondedForceKernel::~OpenCLCalcCustomNonbondedForceKernel() {
delete globals;
if (tabulatedFunctionParams != NULL)
delete tabulatedFunctionParams;
if (interactionGroupData != NULL)
delete interactionGroupData;
for (int i = 0; i < (int) tabulatedFunctions.size(); i++)
delete tabulatedFunctions[i];
if (forceCopy != NULL)
......@@ -1920,7 +1937,7 @@ void OpenCLCalcCustomNonbondedForceKernel::initialize(const System& system, cons
int forceIndex;
for (forceIndex = 0; forceIndex < system.getNumForces() && &system.getForce(forceIndex) != &force; ++forceIndex)
;
string prefix = "custom"+cl.intToString(forceIndex)+"_";
string prefix = (force.getNumInteractionGroups() == 0 ? "custom"+cl.intToString(forceIndex)+"_" : "");
// Record parameters and exclusions.
......@@ -2021,14 +2038,18 @@ void OpenCLCalcCustomNonbondedForceKernel::initialize(const System& system, cons
replacements["SWITCH_C5"] = cl.doubleToString(6/pow(force.getSwitchingDistance()-force.getCutoffDistance(), 5.0));
}
string source = cl.replaceStrings(OpenCLKernelSources::customNonbonded, replacements);
cl.getNonbondedUtilities().addInteraction(useCutoff, usePeriodic, true, force.getCutoffDistance(), exclusionList, source, force.getForceGroup());
for (int i = 0; i < (int) params->getBuffers().size(); i++) {
const OpenCLNonbondedUtilities::ParameterInfo& buffer = params->getBuffers()[i];
cl.getNonbondedUtilities().addParameter(OpenCLNonbondedUtilities::ParameterInfo(prefix+"params"+cl.intToString(i+1), buffer.getComponentType(), buffer.getNumComponents(), buffer.getSize(), buffer.getMemory()));
}
if (globals != NULL) {
globals->upload(globalParamValues);
cl.getNonbondedUtilities().addArgument(OpenCLNonbondedUtilities::ParameterInfo(prefix+"globals", "float", 1, sizeof(cl_float), globals->getDeviceBuffer()));
if (force.getNumInteractionGroups() > 0)
initInteractionGroups(force, source);
else {
cl.getNonbondedUtilities().addInteraction(useCutoff, usePeriodic, true, force.getCutoffDistance(), exclusionList, source, force.getForceGroup());
for (int i = 0; i < (int) params->getBuffers().size(); i++) {
const OpenCLNonbondedUtilities::ParameterInfo& buffer = params->getBuffers()[i];
cl.getNonbondedUtilities().addParameter(OpenCLNonbondedUtilities::ParameterInfo(prefix+"params"+cl.intToString(i+1), buffer.getComponentType(), buffer.getNumComponents(), buffer.getSize(), buffer.getMemory()));
}
if (globals != NULL) {
globals->upload(globalParamValues);
cl.getNonbondedUtilities().addArgument(OpenCLNonbondedUtilities::ParameterInfo(prefix+"globals", "float", 1, sizeof(cl_float), globals->getDeviceBuffer()));
}
}
cl.addForce(new OpenCLCustomNonbondedForceInfo(cl.getNonbondedUtilities().getNumForceBuffers(), force));
......@@ -2044,6 +2065,250 @@ void OpenCLCalcCustomNonbondedForceKernel::initialize(const System& system, cons
}
}
void OpenCLCalcCustomNonbondedForceKernel::initInteractionGroups(const CustomNonbondedForce& force, const string& interactionSource) {
// Process groups to form tiles.
vector<vector<int> > atomLists;
vector<pair<int, int> > tiles;
map<pair<int, int>, int> duplicateInteractions;
for (int group = 0; group < force.getNumInteractionGroups(); group++) {
// Get the list of atoms in this group and sort them.
set<int> set1, set2;
force.getInteractionGroupParameters(group, set1, set2);
vector<int> atoms1, atoms2;
atoms1.insert(atoms1.begin(), set1.begin(), set1.end());
atoms2.insert(atoms2.begin(), set2.begin(), set2.end());
sort(atoms1.begin(), atoms1.end());
sort(atoms2.begin(), atoms2.end());
// Find how many tiles we will create for this group.
int tileWidth = min(min(32, (int) atoms1.size()), (int) atoms2.size());
int numBlocks1 = (atoms1.size()+tileWidth-1)/tileWidth;
int numBlocks2 = (atoms2.size()+tileWidth-1)/tileWidth;
// Add the tiles.
for (int i = 0; i < numBlocks1; i++)
for (int j = 0; j < numBlocks2; j++)
tiles.push_back(make_pair(atomLists.size()+i, atomLists.size()+numBlocks1+j));
// Add the atom lists.
for (int i = 0; i < numBlocks1; i++) {
vector<int> atoms;
int first = i*tileWidth;
int last = min((i+1)*tileWidth, (int) atoms1.size());
for (int j = first; j < last; j++)
atoms.push_back(atoms1[j]);
atomLists.push_back(atoms);
}
for (int i = 0; i < numBlocks2; i++) {
vector<int> atoms;
int first = i*tileWidth;
int last = min((i+1)*tileWidth, (int) atoms2.size());
for (int j = first; j < last; j++)
atoms.push_back(atoms2[j]);
atomLists.push_back(atoms);
}
// If this group contains duplicate interactions, record that we need to skip them once.
for (int i = 0; i < (int) atoms1.size(); i++) {
int a1 = atoms1[i];
if (set2.find(a1) == set2.end())
continue;
for (int j = 0; j < (int) atoms2.size() && atoms2[j] < a1; j++) {
int a2 = atoms2[j];
if (set1.find(a2) != set1.end()) {
pair<int, int> key = make_pair(a2, a1);
if (duplicateInteractions.find(key) == duplicateInteractions.end())
duplicateInteractions[key] = 0;
duplicateInteractions[key]++;
}
}
}
}
// Build a lookup table for quickly identifying excluded interactions.
set<pair<int, int> > exclusions;
for (int i = 0; i < force.getNumExclusions(); i++) {
int p1, p2;
force.getExclusionParticles(i, p1, p2);
exclusions.insert(make_pair(min(p1, p2), max(p1, p2)));
}
// Build the exclusion flags for each tile. While we're at it, filter out tiles
// where all interactions are excluded, and sort the tiles by size.
vector<vector<int> > exclusionFlags(tiles.size());
vector<pair<int, int> > tileOrder;
for (int tile = 0; tile < tiles.size(); tile++) {
if (atomLists[tiles[tile].first].size() < atomLists[tiles[tile].second].size()) {
// For efficiency, we want the first axis to be the larger one.
int swap = tiles[tile].first;
tiles[tile].first = tiles[tile].second;
tiles[tile].second = swap;
}
vector<int>& atoms1 = atomLists[tiles[tile].first];
vector<int>& atoms2 = atomLists[tiles[tile].second];
vector<int> flags(atoms1.size(), (int) (1LL<<atoms2.size())-1);
int numExcluded = 0;
for (int i = 0; i < (int) atoms1.size(); i++)
for (int j = 0; j < (int) atoms2.size(); j++) {
int a1 = atoms1[i];
int a2 = atoms2[j];
bool isExcluded = false;
pair<int, int> key = make_pair(min(a1, a2), max(a1, a2));
if (a1 == a2 || exclusions.find(key) != exclusions.end())
isExcluded = true; // This is an excluded interaction.
else if (duplicateInteractions.find(key) != duplicateInteractions.end() && duplicateInteractions[key] > 0) {
// Both atoms are in both sets, so skip duplicate interactions.
isExcluded = true;
duplicateInteractions[key]--;
}
if (isExcluded) {
flags[i] &= -1-(1<<j);
numExcluded++;
}
}
if (numExcluded == atoms1.size()*atoms2.size())
continue; // All interactions are excluded.
tileOrder.push_back(make_pair((int) -atoms2.size(), tile));
exclusionFlags[tile] = flags;
}
sort(tileOrder.begin(), tileOrder.end());
// Merge tiles to get as close as possible to 32 along the first axis of each one.
vector<int> tileSetStart;
tileSetStart.push_back(0);
int tileSetSize = 0;
for (int i = 0; i < tileOrder.size(); i++) {
int tile = tileOrder[i].second;
int size = atomLists[tiles[tile].first].size();
if (tileSetSize+size > 32) {
tileSetStart.push_back(i);
tileSetSize = 0;
}
tileSetSize += size;
}
tileSetStart.push_back(tileOrder.size());
// Build the data structures.
int numTileSets = tileSetStart.size()-1;
vector<mm_int4> groupData;
for (int tileSet = 0; tileSet < numTileSets; tileSet++) {
int indexInTileSet = 0;
int minSize = 0;
if (cl.getSIMDWidth() < 32) {
// We need to include a barrier inside the inner loop, so ensure that all
// threads will loop the same number of times.
for (int i = tileSetStart[tileSet]; i < tileSetStart[tileSet+1]; i++)
minSize = max(minSize, (int) atomLists[tiles[tileOrder[i].second].first].size());
}
for (int i = tileSetStart[tileSet]; i < tileSetStart[tileSet+1]; i++) {
int tile = tileOrder[i].second;
vector<int>& atoms1 = atomLists[tiles[tile].first];
vector<int>& atoms2 = atomLists[tiles[tile].second];
int range = indexInTileSet + ((indexInTileSet+max(minSize, (int) atoms1.size()))<<16);
int allFlags = (1<<atoms2.size())-1;
for (int j = 0; j < (int) atoms1.size(); j++) {
int a1 = atoms1[j];
int a2 = (j < atoms2.size() ? atoms2[j] : 0);
int flags = (exclusionFlags[tile].size() > 0 ? exclusionFlags[tile][j] : allFlags);
groupData.push_back(mm_int4(a1, a2, range, flags<<indexInTileSet));
}
indexInTileSet += atoms1.size();
}
for (; indexInTileSet < 32; indexInTileSet++)
groupData.push_back(mm_int4(0, 0, minSize<<16, 0));
}
interactionGroupData = OpenCLArray::create<mm_int4>(cl, groupData.size(), "interactionGroupData");
interactionGroupData->upload(groupData);
// Create the kernel.
map<string, string> replacements;
replacements["COMPUTE_INTERACTION"] = interactionSource;
const string suffixes[] = {"x", "y", "z", "w"};
stringstream localData;
int localDataSize = 0;
vector<OpenCLNonbondedUtilities::ParameterInfo>& buffers = params->getBuffers();
for (int i = 0; i < (int) buffers.size(); i++) {
if (buffers[i].getNumComponents() == 1)
localData<<buffers[i].getComponentType()<<" params"<<(i+1)<<";\n";
else {
for (int j = 0; j < buffers[i].getNumComponents(); ++j)
localData<<buffers[i].getComponentType()<<" params"<<(i+1)<<"_"<<suffixes[j]<<";\n";
}
localDataSize += buffers[i].getSize();
}
replacements["ATOM_PARAMETER_DATA"] = localData.str();
stringstream args;
for (int i = 0; i < (int) buffers.size(); i++)
args<<", __global const "<<buffers[i].getType()<<"* restrict global_params"<<(i+1);
if (globals != NULL)
args<<", __global const float* restrict globals";
replacements["PARAMETER_ARGUMENTS"] = args.str();
stringstream load1;
for (int i = 0; i < (int) buffers.size(); i++)
load1<<buffers[i].getType()<<" params"<<(i+1)<<"1 = global_params"<<(i+1)<<"[atom1];\n";
replacements["LOAD_ATOM1_PARAMETERS"] = load1.str();
stringstream loadLocal2;
for (int i = 0; i < (int) buffers.size(); i++) {
if (buffers[i].getNumComponents() == 1)
loadLocal2<<"localData[get_local_id(0)].params"<<(i+1)<<" = global_params"<<(i+1)<<"[atom2];\n";
else {
loadLocal2<<buffers[i].getType()<<" temp_params"<<(i+1)<<" = global_params"<<(i+1)<<"[atom2];\n";
for (int j = 0; j < buffers[i].getNumComponents(); ++j)
loadLocal2<<"localData[get_local_id(0)].params"<<(i+1)<<"_"<<suffixes[j]<<" = temp_params"<<(i+1)<<"."<<suffixes[j]<<";\n";
}
}
replacements["LOAD_LOCAL_PARAMETERS"] = loadLocal2.str();
stringstream load2;
for (int i = 0; i < (int) buffers.size(); i++) {
if (buffers[i].getNumComponents() == 1)
load2<<buffers[i].getType()<<" params"<<(i+1)<<"2 = localData[localIndex].params"<<(i+1)<<";\n";
else {
load2<<buffers[i].getType()<<" params"<<(i+1)<<"2 = ("<<buffers[i].getType()<<") (";
for (int j = 0; j < buffers[i].getNumComponents(); ++j) {
if (j > 0)
load2<<", ";
load2<<"localData[localIndex].params"<<(i+1)<<"_"<<suffixes[j];
}
load2<<");\n";
}
}
replacements["LOAD_ATOM2_PARAMETERS"] = load2.str();
map<string, string> defines;
if (force.getNonbondedMethod() != CustomNonbondedForce::NoCutoff)
defines["USE_CUTOFF"] = "1";
if (force.getNonbondedMethod() == CustomNonbondedForce::CutoffPeriodic)
defines["USE_PERIODIC"] = "1";
defines["THREAD_BLOCK_SIZE"] = cl.intToString(cl.getNonbondedUtilities().getForceThreadBlockSize());
double cutoff = force.getCutoffDistance();
defines["CUTOFF_SQUARED"] = cl.doubleToString(cutoff*cutoff);
defines["PADDED_NUM_ATOMS"] = cl.intToString(cl.getPaddedNumAtoms());
defines["TILE_SIZE"] = "32";
int numContexts = cl.getPlatformData().contexts.size();
int startIndex = cl.getContextIndex()*numTileSets/numContexts;
int endIndex = (cl.getContextIndex()+1)*numTileSets/numContexts;
defines["FIRST_TILE"] = cl.intToString(startIndex);
defines["LAST_TILE"] = cl.intToString(endIndex);
if ((localDataSize/4)%2 == 0 && !cl.getUseDoublePrecision())
defines["PARAMETER_SIZE_IS_EVEN"] = "1";
cl::Program program = cl.createProgram(cl.replaceStrings(OpenCLKernelSources::customNonbondedGroups, replacements), defines);
interactionGroupKernel = cl::Kernel(program, "computeInteractionGroups");
numGroupThreadBlocks = cl.getNonbondedUtilities().getNumForceThreadBlocks();
}
double OpenCLCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
if (globals != NULL) {
bool changed = false;
......@@ -2065,6 +2330,25 @@ double OpenCLCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bool
longRangeCoefficient = CustomNonbondedForceImpl::calcLongRangeCorrection(*forceCopy, context.getOwner());
hasInitializedLongRangeCorrection = true;
}
if (interactionGroupData != NULL) {
if (!hasInitializedKernel) {
hasInitializedKernel = true;
int index = 0;
bool useLong = cl.getSupports64BitGlobalAtomics();
interactionGroupKernel.setArg<cl::Buffer>(index++, (useLong ? cl.getLongForceBuffer() : cl.getForceBuffers()).getDeviceBuffer());
interactionGroupKernel.setArg<cl::Buffer>(index++, cl.getEnergyBuffer().getDeviceBuffer());
interactionGroupKernel.setArg<cl::Buffer>(index++, cl.getPosq().getDeviceBuffer());
interactionGroupKernel.setArg<cl::Buffer>(index++, interactionGroupData->getDeviceBuffer());
setPeriodicBoxSizeArg(cl, interactionGroupKernel, index++);
setInvPeriodicBoxSizeArg(cl, interactionGroupKernel, index++);
for (int i = 0; i < (int) params->getBuffers().size(); i++)
interactionGroupKernel.setArg<cl::Memory>(index++, params->getBuffers()[i].getMemory());
if (globals != NULL)
interactionGroupKernel.setArg<cl::Buffer>(index++, globals->getDeviceBuffer());
}
int forceThreadBlockSize = max(32, cl.getNonbondedUtilities().getForceThreadBlockSize());
cl.executeKernel(interactionGroupKernel, numGroupThreadBlocks*forceThreadBlockSize, forceThreadBlockSize);
}
mm_double4 boxSize = cl.getPeriodicBoxSizeDouble();
return longRangeCoefficient/(boxSize.x*boxSize.y*boxSize.z);
}
......
platforms/opencl/src/kernels/customNonbondedGroups.cl 0 → 100644
View file @ 81d60b1e
#ifdef SUPPORTS_64_BIT_ATOMICS
#pragma OPENCL EXTENSION cl_khr_int64_base_atomics : enable
#endif
#define WARPS_PER_GROUP (THREAD_BLOCK_SIZE/TILE_SIZE)
typedef struct {
real x, y, z;
real q;
real fx, fy, fz;
ATOM_PARAMETER_DATA
#ifndef PARAMETER_SIZE_IS_EVEN
real padding;
#endif
} AtomData;
/**
* This function is used on devices that don't support 64 bit atomics. Multiple threads within
* a single tile might have computed forces on the same atom. This loops over them and makes sure
* that only one thread updates the force on any given atom.
*/
void writeForces(__global real4* forceBuffers,__local AtomData* localData, int atomIndex) {
localData[get_local_id(0)].x = atomIndex;
SYNC_WARPS;
real4 forceSum = (real4) 0;
int start = (get_local_id(0)/TILE_SIZE)*TILE_SIZE;
int end = start+32;
bool isFirst = true;
for (int i = start; i < end; i++)
if (localData[i].x == atomIndex) {
forceSum += (real4) (localData[i].fx, localData[i].fy, localData[i].fz, 0);
isFirst &= (i >= get_local_id(0));
}
const unsigned int warp = get_global_id(0)/TILE_SIZE;
unsigned int offset = atomIndex + warp*PADDED_NUM_ATOMS;
if (isFirst)
forceBuffers[offset] += forceSum;
SYNC_WARPS;
}
__kernel void computeInteractionGroups(
#ifdef SUPPORTS_64_BIT_ATOMICS
__global long* restrict forceBuffers,
#else
__global real4* restrict forceBuffers,
#endif
__global real* restrict energyBuffer, __global const real4* restrict posq,
__global const int4* restrict groupData, real4 periodicBoxSize, real4 invPeriodicBoxSize
PARAMETER_ARGUMENTS) {
const unsigned int totalWarps = get_global_size(0)/TILE_SIZE;
const unsigned int warp = get_global_id(0)/TILE_SIZE; // global warpIndex
const unsigned int tgx = get_local_id(0) & (TILE_SIZE-1); // index within the warp
const unsigned int tbx = get_local_id(0) - tgx; // block warpIndex
real energy = 0.0f;
__local AtomData localData[THREAD_BLOCK_SIZE];
const unsigned int startTile = FIRST_TILE+warp*(LAST_TILE-FIRST_TILE)/totalWarps;
const unsigned int endTile = FIRST_TILE+(warp+1)*(LAST_TILE-FIRST_TILE)/totalWarps;
for (int tile = startTile; tile < endTile; tile++) {
const int4 atomData = groupData[TILE_SIZE*tile+tgx];
const int atom1 = atomData.x;
const int atom2 = atomData.y;
const int rangeStart = atomData.z&0xFFFF;
const int rangeEnd = (atomData.z>>16)&0xFFFF;
const int exclusions = atomData.w;
real4 posq1 = posq[atom1];
LOAD_ATOM1_PARAMETERS
real4 force = (real4) (0);
real4 posq2 = posq[atom2];
localData[get_local_id(0)].x = posq2.x;
localData[get_local_id(0)].y = posq2.y;
localData[get_local_id(0)].z = posq2.z;
localData[get_local_id(0)].q = posq2.w;
LOAD_LOCAL_PARAMETERS
localData[get_local_id(0)].fx = 0.0f;
localData[get_local_id(0)].fy = 0.0f;
localData[get_local_id(0)].fz = 0.0f;
int tj = tgx;
SYNC_WARPS;
for (int j = rangeStart; j < rangeEnd; j++) {
if (tj < rangeEnd) {
bool isExcluded = (((exclusions>>tj)&1) == 0);
int localIndex = tbx+tj;
posq2 = (real4) (localData[localIndex].x, localData[localIndex].y, localData[localIndex].z, localData[localIndex].q);
real4 delta = (real4) (posq2.xyz - posq1.xyz, 0);
#ifdef USE_PERIODIC
delta.xyz -= floor(delta.xyz*invPeriodicBoxSize.xyz+0.5f)*periodicBoxSize.xyz;
#endif
real r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
#ifdef USE_CUTOFF
if (!isExcluded && r2 < CUTOFF_SQUARED) {
#endif
real invR = RSQRT(r2);
real r = RECIP(invR);
LOAD_ATOM2_PARAMETERS
real dEdR = 0.0f;
real tempEnergy = 0.0f;
COMPUTE_INTERACTION
energy += tempEnergy;
delta *= dEdR;
force.xyz -= delta.xyz;
localData[localIndex].fx += delta.x;
localData[localIndex].fy += delta.y;
localData[localIndex].fz += delta.z;
#ifdef USE_CUTOFF
}
#endif
}
tj = (tj == rangeEnd-1 ? rangeStart : tj+1);
SYNC_WARPS;
}
#ifdef SUPPORTS_64_BIT_ATOMICS
if (exclusions != 0) {
atom_add(&forceBuffers[atom1], (long) (force.x*0x100000000));
atom_add(&forceBuffers[atom1+PADDED_NUM_ATOMS], (long) (force.y*0x100000000));
atom_add(&forceBuffers[atom1+2*PADDED_NUM_ATOMS], (long) (force.z*0x100000000));
}
atom_add(&forceBuffers[atom2], (long) (localData[get_local_id(0)].fx*0x100000000));
atom_add(&forceBuffers[atom2+PADDED_NUM_ATOMS], (long) (localData[get_local_id(0)].fy*0x100000000));
atom_add(&forceBuffers[atom2+2*PADDED_NUM_ATOMS], (long) (localData[get_local_id(0)].fz*0x100000000));
#else
writeForces(forceBuffers, localData, atom2);
localData[get_local_id(0)].fx = force.x;
localData[get_local_id(0)].fy = force.y;
localData[get_local_id(0)].fz = force.z;
writeForces(forceBuffers, localData, atom1);
#endif
}
energyBuffer[get_global_id(0)] += energy;
}
\ No newline at end of file
platforms/opencl/tests/TestOpenCLCustomNonbondedForce.cpp
View file @ 81d60b1e
......@@ -34,6 +34,9 @@
* This tests all the different force terms in the OpenCL implementation of CustomNonbondedForce.
*/
#ifdef WIN32
#define _USE_MATH_DEFINES // Needed to get M_PI
#endif
#include "openmm/internal/AssertionUtilities.h"
#include "sfmt/SFMT.h"
#include "openmm/Context.h"
......@@ -42,6 +45,7 @@
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include <cmath>
#include <iostream>
#include <vector>
......@@ -538,6 +542,179 @@ void testLongRangeCorrection() {
ASSERT_EQUAL_TOL(standardEnergy1-standardEnergy2, customEnergy1-customEnergy2, 1e-4);
}
void testInteractionGroups() {
const int numParticles = 6;
System system;
VerletIntegrator integrator(0.01);
CustomNonbondedForce* nonbonded = new CustomNonbondedForce("v1+v2");
nonbonded->addPerParticleParameter("v");
vector<double> params(1, 0.001);
for (int i = 0; i < numParticles; i++) {
system.addParticle(1.0);
nonbonded->addParticle(params);
params[0] *= 10;
}
set<int> set1, set2, set3, set4;
set1.insert(2);
set2.insert(0);
set2.insert(1);
set2.insert(2);
set2.insert(3);
set2.insert(4);
set2.insert(5);
nonbonded->addInteractionGroup(set1, set2); // Particle 2 interacts with every other particle.
set3.insert(0);
set3.insert(1);
set4.insert(4);
set4.insert(5);
nonbonded->addInteractionGroup(set3, set4); // Particles 0 and 1 interact with 4 and 5.
nonbonded->addExclusion(1, 2); // Add an exclusion to make sure it gets skipped.
system.addForce(nonbonded);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
context.setPositions(positions);
State state = context.getState(State::Energy);
double expectedEnergy = 331.423; // Each digit is the number of interactions a particle particle is involved in.
ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), TOL);
}
void testLargeInteractionGroup() {
const int numMolecules = 300;
const int numParticles = numMolecules*2;
const double boxSize = 20.0;
// Create a large system.
System system;
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
for (int i = 0; i < numParticles; i++)
system.addParticle(1.0);
CustomNonbondedForce* nonbonded = new CustomNonbondedForce("4*eps*((sigma/r)^12-(sigma/r)^6)+138.935456*q/r; q=q1*q2; sigma=0.5*(sigma1+sigma2); eps=sqrt(eps1*eps2)");
nonbonded->addPerParticleParameter("q");
nonbonded->addPerParticleParameter("sigma");
nonbonded->addPerParticleParameter("eps");
vector<Vec3> positions(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<double> params(3);
for (int i = 0; i < numMolecules; i++) {
if (i < numMolecules/2) {
params[0] = 1.0;
params[1] = 0.2;
params[2] = 0.1;
nonbonded->addParticle(params);
params[0] = -1.0;
params[1] = 0.1;
nonbonded->addParticle(params);
}
else {
params[0] = 1.0;
params[1] = 0.2;
params[2] = 0.2;
nonbonded->addParticle(params);
params[0] = -1.0;
params[1] = 0.1;
nonbonded->addParticle(params);
}
positions[2*i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
positions[2*i+1] = Vec3(positions[2*i][0]+1.0, positions[2*i][1], positions[2*i][2]);
nonbonded->addExclusion(2*i, 2*i+1);
}
nonbonded->setNonbondedMethod(CustomNonbondedForce::CutoffPeriodic);
system.addForce(nonbonded);
// Compute the forces.
VerletIntegrator integrator(0.01);
Context context(system, integrator, platform);
context.setPositions(positions);
State state1 = context.getState(State::Forces);
// Modify the force so only one particle interacts with everything else.
set<int> set1, set2;
set1.insert(151);
for (int i = 0; i < numParticles; i++)
set2.insert(i);
nonbonded->addInteractionGroup(set1, set2);
context.reinitialize();
context.setPositions(positions);
State state2 = context.getState(State::Forces);
// The force on that one particle should be the same.
ASSERT_EQUAL_VEC(state1.getForces()[151], state2.getForces()[151], 1e-4);
// Modify the interaction group so it includes all interactions. This should now reproduce the original forces
// on all atoms.
for (int i = 0; i < numParticles; i++)
set1.insert(i);
nonbonded->setInteractionGroupParameters(0, set1, set2);
context.reinitialize();
context.setPositions(positions);
State state3 = context.getState(State::Forces);
for (int i = 0; i < numParticles; i++)
ASSERT_EQUAL_VEC(state1.getForces()[i], state3.getForces()[i], 1e-4);
}
void testInteractionGroupLongRangeCorrection() {
const int numParticles = 10;
const double boxSize = 10.0;
const double cutoff = 0.5;
System system;
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
CustomNonbondedForce* nonbonded = new CustomNonbondedForce("c1*c2*r^-4");
nonbonded->addPerParticleParameter("c");
vector<Vec3> positions(numParticles);
vector<double> params(1);
for (int i = 0; i < numParticles; i++) {
system.addParticle(1.0);
params[0] = (i%2 == 0 ? 1.1 : 2.0);
nonbonded->addParticle(params);
positions[i] = Vec3(0.5*i, 0, 0);
}
nonbonded->setNonbondedMethod(CustomNonbondedForce::CutoffPeriodic);
nonbonded->setCutoffDistance(cutoff);
system.addForce(nonbonded);
// Setup nonbonded groups. They involve 1 interaction of type AA,
// 2 of type BB, and 5 of type AB.
set<int> set1, set2, set3, set4, set5;
set1.insert(0);
set1.insert(1);
set1.insert(2);
nonbonded->addInteractionGroup(set1, set1);
set2.insert(3);
set3.insert(4);
set3.insert(6);
set3.insert(8);
nonbonded->addInteractionGroup(set2, set3);
set4.insert(5);
set5.insert(7);
set5.insert(9);
nonbonded->addInteractionGroup(set4, set5);
// Compute energy with and without the correction.
VerletIntegrator integrator(0.01);
Context context(system, integrator, platform);
context.setPositions(positions);
double energy1 = context.getState(State::Energy).getPotentialEnergy();
nonbonded->setUseLongRangeCorrection(true);
context.reinitialize();
context.setPositions(positions);
double energy2 = context.getState(State::Energy).getPotentialEnergy();
// Check the result.
double sum = (1.1*1.1 + 2*2.0*2.0 + 5*1.1*2.0)*2.0;
int numPairs = (numParticles*(numParticles+1))/2;
double expected = 2*M_PI*numParticles*numParticles*sum/(numPairs*boxSize*boxSize*boxSize);
ASSERT_EQUAL_TOL(expected, energy2-energy1, 1e-4);
}
int main(int argc, char* argv[]) {
try {
if (argc > 1)
......@@ -553,6 +730,9 @@ int main(int argc, char* argv[]) {
testParallelComputation();
testSwitchingFunction();
testLongRangeCorrection();
testInteractionGroups();
testLargeInteractionGroup();
testInteractionGroupLongRangeCorrection();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
......
platforms/reference/include/ReferenceCustomHbondIxn.h
View file @ 81d60b1e
......@@ -29,6 +29,7 @@
#include "lepton/ExpressionProgram.h"
#include "lepton/ParsedExpression.h"
#include <map>
#include <set>
#include <vector>
// ---------------------------------------------------------------------------------------
......@@ -143,10 +144,8 @@ class ReferenceCustomHbondIxn : public ReferenceBondIxn {
@param atomCoordinates atom coordinates
@param donorParameters donor parameters values donorParameters[donorIndex][parameterIndex]
@param acceptorParameters acceptor parameters values acceptorParameters[acceptorIndex][parameterIndex]
@param exclusions exclusion indices exclusions[donorIndex][acceptorToExcludeIndex]
exclusions[donorIndex][0] = number of exclusions
exclusions[donorIndex][no.-1] = indices of acceptors to excluded from
interacting w/ donor donorIndex
@param exclusions exclusion indices
exclusions[donorIndex] contains the list of excluded acceptors for that donor
@param globalParameters the values of global parameters
@param forces force array (forces added)
@param totalEnergy total energy
......@@ -154,7 +153,7 @@ class ReferenceCustomHbondIxn : public ReferenceBondIxn {
--------------------------------------------------------------------------------------- */
void calculatePairIxn(std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMM** donorParameters, RealOpenMM** acceptorParameters,
int** exclusions, const std::map<std::string, double>& globalParameters,
std::vector<std::set<int> >& exclusions, const std::map<std::string, double>& globalParameters,
std::vector<OpenMM::RealVec>& forces, RealOpenMM* totalEnergy) const;
// ---------------------------------------------------------------------------------------
......
platforms/reference/include/ReferenceCustomNonbondedIxn.h
View file @ 81d60b1e
......@@ -29,6 +29,8 @@
#include "ReferenceNeighborList.h"
#include "lepton/ExpressionProgram.h"
#include <map>
#include <set>
#include <utility>
#include <vector>
// ---------------------------------------------------------------------------------------
......@@ -47,6 +49,7 @@ class ReferenceCustomNonbondedIxn {
Lepton::ExpressionProgram forceExpression;
std::vector<std::string> paramNames;
std::vector<std::string> particleParamNames;
std::vector<std::pair<std::set<int>, std::set<int> > > interactionGroups;
/**---------------------------------------------------------------------------------------
......@@ -97,6 +100,17 @@ class ReferenceCustomNonbondedIxn {
void setUseCutoff( RealOpenMM distance, const OpenMM::NeighborList& neighbors );
/**---------------------------------------------------------------------------------------
Restrict the force to a list of interaction groups.
@param distance the cutoff distance
@param neighbors the neighbor list to use
--------------------------------------------------------------------------------------- */
void setInteractionGroups(const std::vector<std::pair<std::set<int>, std::set<int> > >& groups);
/**---------------------------------------------------------------------------------------
Set the force to use a switching function.
......@@ -126,10 +140,8 @@ class ReferenceCustomNonbondedIxn {
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param atomParameters atom parameters (charges, c6, c12, ...) atomParameters[atomIndex][paramterIndex]
@param exclusions atom exclusion indices exclusions[atomIndex][atomToExcludeIndex]
exclusions[atomIndex][0] = number of exclusions
exclusions[atomIndex][1-no.] = atom indices of atoms to excluded from
interacting w/ atom atomIndex
@param exclusions atom exclusion indices
exclusions[atomIndex] contains the list of exclusions for that atom
@param fixedParameters non atom parameters (not currently used)
@param globalParameters the values of global parameters
@param forces force array (forces added)
......@@ -139,7 +151,7 @@ class ReferenceCustomNonbondedIxn {
--------------------------------------------------------------------------------------- */
void calculatePairIxn( int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
RealOpenMM** atomParameters, int** exclusions,
RealOpenMM** atomParameters, std::vector<std::set<int> >& exclusions,
RealOpenMM* fixedParameters, const std::map<std::string, double>& globalParameters,
std::vector<OpenMM::RealVec>& forces, RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const;
......
platforms/reference/include/ReferenceKernels.h
View file @ 81d60b1e
......@@ -572,7 +572,7 @@ public:
void copyParametersToContext(ContextImpl& context, const NonbondedForce& force);
private:
int numParticles, num14;
int **exclusionArray, **bonded14IndexArray;
int **bonded14IndexArray;
RealOpenMM **particleParamArray, **bonded14ParamArray;
RealOpenMM nonbondedCutoff, switchingDistance, rfDielectric, ewaldAlpha, dispersionCoefficient;
int kmax[3], gridSize[3];
......@@ -615,7 +615,6 @@ public:
void copyParametersToContext(ContextImpl& context, const CustomNonbondedForce& force);
private:
int numParticles;
int **exclusionArray;
RealOpenMM **particleParamArray;
RealOpenMM nonbondedCutoff, switchingDistance, periodicBoxSize[3], longRangeCoefficient;
bool useSwitchingFunction, hasInitializedLongRangeCorrection;
......@@ -624,6 +623,7 @@ private:
std::vector<std::set<int> > exclusions;
Lepton::ExpressionProgram energyExpression, forceExpression;
std::vector<std::string> parameterNames, globalParameterNames;
std::vector<std::pair<std::set<int>, std::set<int> > > interactionGroups;
NonbondedMethod nonbondedMethod;
NeighborList* neighborList;
};
......@@ -819,7 +819,6 @@ public:
private:
int numDonors, numAcceptors, numParticles;
bool isPeriodic;
int **exclusionArray;
RealOpenMM **donorParamArray, **acceptorParamArray;
RealOpenMM nonbondedCutoff;
ReferenceCustomHbondIxn* ixn;
......
platforms/reference/include/ReferenceLJCoulombIxn.h
View file @ 81d60b1e
......@@ -156,10 +156,8 @@ class ReferenceLJCoulombIxn {
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param atomParameters atom parameters (charges, c6, c12, ...) atomParameters[atomIndex][paramterIndex]
@param exclusions atom exclusion indices exclusions[atomIndex][atomToExcludeIndex]
exclusions[atomIndex][0] = number of exclusions
exclusions[atomIndex][1-no.] = atom indices of atoms to excluded from
interacting w/ atom atomIndex
@param exclusions atom exclusion indices
exclusions[atomIndex] contains the list of exclusions for that atom
@param fixedParameters non atom parameters (not currently used)
@param forces force array (forces added)
@param energyByAtom atom energy
......@@ -170,7 +168,7 @@ class ReferenceLJCoulombIxn {
--------------------------------------------------------------------------------------- */
void calculatePairIxn(int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
RealOpenMM** atomParameters, int** exclusions,
RealOpenMM** atomParameters, std::vector<std::set<int> >& exclusions,
RealOpenMM* fixedParameters, std::vector<OpenMM::RealVec>& forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy, bool includeDirect, bool includeReciprocal) const;
......@@ -182,10 +180,8 @@ private:
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param atomParameters atom parameters (charges, c6, c12, ...) atomParameters[atomIndex][paramterIndex]
@param exclusions atom exclusion indices exclusions[atomIndex][atomToExcludeIndex]
exclusions[atomIndex][0] = number of exclusions
exclusions[atomIndex][1-no.] = atom indices of atoms to excluded from
interacting w/ atom atomIndex
@param exclusions atom exclusion indices
exclusions[atomIndex] contains the list of exclusions for that atom
@param fixedParameters non atom parameters (not currently used)
@param forces force array (forces added)
@param energyByAtom atom energy
......@@ -196,7 +192,7 @@ private:
--------------------------------------------------------------------------------------- */
void calculateEwaldIxn(int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
RealOpenMM** atomParameters, int** exclusions,
RealOpenMM** atomParameters, std::vector<std::set<int> >& exclusions,
RealOpenMM* fixedParameters, std::vector<OpenMM::RealVec>& forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy, bool includeDirect, bool includeReciprocal) const;
};
......
platforms/reference/src/ReferenceKernels.cpp
View file @ 81d60b1e
......@@ -805,7 +805,6 @@ void ReferenceCalcCustomTorsionForceKernel::copyParametersToContext(ContextImpl&
ReferenceCalcNonbondedForceKernel::~ReferenceCalcNonbondedForceKernel() {
disposeRealArray(particleParamArray, numParticles);
disposeIntArray(exclusionArray, numParticles);
disposeIntArray(bonded14IndexArray, num14);
disposeRealArray(bonded14ParamArray, num14);
if (neighborList != NULL)
......@@ -843,14 +842,6 @@ void ReferenceCalcNonbondedForceKernel::initialize(const System& system, const N
particleParamArray[i][2] = static_cast<RealOpenMM>(charge);
}
this->exclusions = exclusions;
exclusionArray = new int*[numParticles];
for (int i = 0; i < numParticles; ++i) {
exclusionArray[i] = new int[exclusions[i].size()+1];
exclusionArray[i][0] = exclusions[i].size();
int index = 0;
for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter)
exclusionArray[i][++index] = *iter;
}
for (int i = 0; i < num14; ++i) {
int particle1, particle2;
double charge, radius, depth;
......@@ -914,7 +905,7 @@ double ReferenceCalcNonbondedForceKernel::execute(ContextImpl& context, bool inc
clj.setUsePME(ewaldAlpha, gridSize);
if (useSwitchingFunction)
clj.setUseSwitchingFunction(switchingDistance);
clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, includeEnergy ? &energy : NULL, includeDirect, includeReciprocal);
clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusions, 0, forceData, 0, includeEnergy ? &energy : NULL, includeDirect, includeReciprocal);
if (includeDirect) {
ReferenceBondForce refBondForce;
ReferenceLJCoulomb14 nonbonded14;
......@@ -1002,7 +993,6 @@ public:
ReferenceCalcCustomNonbondedForceKernel::~ReferenceCalcCustomNonbondedForceKernel() {
disposeRealArray(particleParamArray, numParticles);
disposeIntArray(exclusionArray, numParticles);
if (neighborList != NULL)
delete neighborList;
if (forceCopy != NULL)
......@@ -1032,14 +1022,6 @@ void ReferenceCalcCustomNonbondedForceKernel::initialize(const System& system, c
for (int j = 0; j < numParameters; j++)
particleParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
}
exclusionArray = new int*[numParticles];
for (int i = 0; i < numParticles; ++i) {
exclusionArray[i] = new int[exclusions[i].size()+1];
exclusionArray[i][0] = exclusions[i].size();
int index = 0;
for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter)
exclusionArray[i][++index] = *iter;
}
nonbondedMethod = CalcCustomNonbondedForceKernel::NonbondedMethod(force.getNonbondedMethod());
nonbondedCutoff = (RealOpenMM) force.getCutoffDistance();
if (nonbondedMethod == NoCutoff) {
......@@ -1090,6 +1072,14 @@ void ReferenceCalcCustomNonbondedForceKernel::initialize(const System& system, c
longRangeCoefficient = 0.0;
hasInitializedLongRangeCorrection = true;
}
// Record the interaction groups.
for (int i = 0; i < force.getNumInteractionGroups(); i++) {
set<int> set1, set2;
force.getInteractionGroupParameters(i, set1, set2);
interactionGroups.push_back(make_pair(set1, set2));
}
}
double ReferenceCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
......@@ -1109,6 +1099,8 @@ double ReferenceCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bo
throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
ixn.setPeriodic(box);
}
if (interactionGroups.size() > 0)
ixn.setInteractionGroups(interactionGroups);
bool globalParamsChanged = false;
for (int i = 0; i < (int) globalParameterNames.size(); i++) {
double value = context.getParameter(globalParameterNames[i]);
......@@ -1118,7 +1110,7 @@ double ReferenceCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bo
}
if (useSwitchingFunction)
ixn.setUseSwitchingFunction(switchingDistance);
ixn.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, globalParamValues, forceData, 0, includeEnergy ? &energy : NULL);
ixn.calculatePairIxn(numParticles, posData, particleParamArray, exclusions, 0, globalParamValues, forceData, 0, includeEnergy ? &energy : NULL);
// Add in the long range correction.
......@@ -1499,7 +1491,6 @@ void ReferenceCalcCustomExternalForceKernel::copyParametersToContext(ContextImpl
ReferenceCalcCustomHbondForceKernel::~ReferenceCalcCustomHbondForceKernel() {
disposeRealArray(donorParamArray, numDonors);
disposeRealArray(acceptorParamArray, numAcceptors);
disposeIntArray(exclusionArray, numDonors);
if (ixn != NULL)
delete ixn;
}
......@@ -1546,14 +1537,6 @@ void ReferenceCalcCustomHbondForceKernel::initialize(const System& system, const
for (int j = 0; j < numAcceptorParameters; j++)
acceptorParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
}
exclusionArray = new int*[numDonors];
for (int i = 0; i < numDonors; ++i) {
exclusionArray[i] = new int[exclusions[i].size()+1];
exclusionArray[i][0] = exclusions[i].size();
int index = 0;
for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter)
exclusionArray[i][++index] = *iter;
}
NonbondedMethod nonbondedMethod = CalcCustomHbondForceKernel::NonbondedMethod(force.getNonbondedMethod());
nonbondedCutoff = (RealOpenMM) force.getCutoffDistance();
......@@ -1602,7 +1585,7 @@ double ReferenceCalcCustomHbondForceKernel::execute(ContextImpl& context, bool i
map<string, double> globalParameters;
for (int i = 0; i < (int) globalParameterNames.size(); i++)
globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
ixn->calculatePairIxn(posData, donorParamArray, acceptorParamArray, exclusionArray, globalParameters, forceData, includeEnergy ? &energy : NULL);
ixn->calculatePairIxn(posData, donorParamArray, acceptorParamArray, exclusions, globalParameters, forceData, includeEnergy ? &energy : NULL);
return energy;
}
......
platforms/reference/src/SimTKReference/ReferenceCustomHbondIxn.cpp
View file @ 81d60b1e
......@@ -34,6 +34,7 @@
using std::map;
using std::pair;
using std::set;
using std::string;
using std::stringstream;
using std::vector;
......@@ -109,10 +110,8 @@ void ReferenceCustomHbondIxn::setPeriodic(RealVec& boxSize) {
@param atomCoordinates atom coordinates
@param donorParameters donor parameters values donorParameters[donorIndex][parameterIndex]
@param acceptorParameters acceptor parameters values acceptorParameters[acceptorIndex][parameterIndex]
@param exclusions exclusion indices exclusions[donorIndex][acceptorToExcludeIndex]
exclusions[donorIndex][0] = number of exclusions
exclusions[donorIndex][no.-1] = indices of acceptors to excluded from
interacting w/ donor donorIndex
@param exclusions exclusion indices
exclusions[donorIndex] contains the list of excluded acceptors for that donor
@param globalParameters the values of global parameters
@param forces force array (forces added)
@param totalEnergy total energy
......@@ -120,7 +119,7 @@ void ReferenceCustomHbondIxn::setPeriodic(RealVec& boxSize) {
--------------------------------------------------------------------------------------- */
void ReferenceCustomHbondIxn::calculatePairIxn(vector<RealVec>& atomCoordinates, RealOpenMM** donorParameters, RealOpenMM** acceptorParameters,
int** exclusions, const map<string, double>& globalParameters, vector<RealVec>& forces,
vector<set<int> >& exclusions, const map<string, double>& globalParameters, vector<RealVec>& forces,
RealOpenMM* totalEnergy) const {
map<string, double> variables = globalParameters;
......@@ -129,18 +128,8 @@ void ReferenceCustomHbondIxn::calculatePairIxn(vector<RealVec>& atomCoordinates,
int numDonors = donorAtoms.size();
int numAcceptors = acceptorAtoms.size();
int* exclusionIndices = new int[numAcceptors];
for( int ii = 0; ii < numAcceptors; ii++ ){
exclusionIndices[ii] = -1;
}
for( int donor = 0; donor < numDonors; donor++ ){
// set exclusions
for (int j = 1; j <= exclusions[donor][0]; j++)
exclusionIndices[exclusions[donor][j]] = donor;
// Initialize per-donor parameters.
for (int j = 0; j < (int) donorParamNames.size(); j++)
......@@ -149,16 +138,13 @@ void ReferenceCustomHbondIxn::calculatePairIxn(vector<RealVec>& atomCoordinates,
// loop over atom pairs
for( int acceptor = 0; acceptor < numAcceptors; acceptor++ ){
if( exclusionIndices[acceptor] != donor ){
if (exclusions[donor].find(acceptor) == exclusions[donor].end()) {
for (int j = 0; j < (int) acceptorParamNames.size(); j++)
variables[acceptorParamNames[j]] = acceptorParameters[acceptor][j];
calculateOneIxn(donor, acceptor, atomCoordinates, variables, forces, totalEnergy);
}
}
}
delete[] exclusionIndices;
}
/**---------------------------------------------------------------------------------------
......
platforms/reference/src/SimTKReference/ReferenceCustomNonbondedIxn.cpp
View file @ 81d60b1e
......@@ -32,8 +32,10 @@
#include "ReferenceCustomNonbondedIxn.h"
using std::map;
using std::pair;
using std::string;
using std::stringstream;
using std::set;
using std::vector;
using OpenMM::RealVec;
......@@ -94,6 +96,19 @@ ReferenceCustomNonbondedIxn::~ReferenceCustomNonbondedIxn( ){
neighborList = &neighbors;
}
/**---------------------------------------------------------------------------------------
Restrict the force to a list of interaction groups.
@param distance the cutoff distance
@param neighbors the neighbor list to use
--------------------------------------------------------------------------------------- */
void ReferenceCustomNonbondedIxn::setInteractionGroups(const vector<pair<set<int>, set<int> > >& groups) {
interactionGroups = groups;
}
/**---------------------------------------------------------------------------------------
Set the force to use a switching function.
......@@ -138,10 +153,8 @@ void ReferenceCustomNonbondedIxn::setUseSwitchingFunction( RealOpenMM distance )
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param atomParameters atom parameters atomParameters[atomIndex][paramterIndex]
@param exclusions atom exclusion indices exclusions[atomIndex][atomToExcludeIndex]
exclusions[atomIndex][0] = number of exclusions
exclusions[atomIndex][1-no.] = atom indices of atoms to excluded from
interacting w/ atom atomIndex
@param exclusions atom exclusion indices
exclusions[atomIndex] contains the list of exclusions for that atom
@param fixedParameters non atom parameters (not currently used)
@param globalParameters the values of global parameters
@param forces force array (forces added)
......@@ -151,53 +164,59 @@ void ReferenceCustomNonbondedIxn::setUseSwitchingFunction( RealOpenMM distance )
--------------------------------------------------------------------------------------- */
void ReferenceCustomNonbondedIxn::calculatePairIxn( int numberOfAtoms, vector<RealVec>& atomCoordinates,
RealOpenMM** atomParameters, int** exclusions,
RealOpenMM** atomParameters, vector<set<int> >& exclusions,
RealOpenMM* fixedParameters, const map<string, double>& globalParameters, vector<RealVec>& forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const {
map<string, double> variables = globalParameters;
if (cutoff) {
for (int i = 0; i < (int) neighborList->size(); i++) {
OpenMM::AtomPair pair = (*neighborList)[i];
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];
}
calculateOneIxn(pair.first, pair.second, atomCoordinates, variables, forces, energyByAtom, totalEnergy);
}
}
else {
// allocate and initialize exclusion array
int* exclusionIndices = new int[numberOfAtoms];
for( int ii = 0; ii < numberOfAtoms; ii++ ){
exclusionIndices[ii] = -1;
}
for( int ii = 0; ii < numberOfAtoms; ii++ ){
// set exclusions
for( int jj = 1; jj <= exclusions[ii][0]; jj++ ){
exclusionIndices[exclusions[ii][jj]] = ii;
}
// loop over atom pairs
for( int jj = ii+1; jj < numberOfAtoms; jj++ ){
if( exclusionIndices[jj] != ii ){
for (int j = 0; j < (int) paramNames.size(); j++) {
variables[particleParamNames[j*2]] = atomParameters[ii][j];
variables[particleParamNames[j*2+1]] = atomParameters[jj][j];
}
calculateOneIxn(ii, jj, atomCoordinates, variables, forces, energyByAtom, totalEnergy);
}
}
}
delete[] exclusionIndices;
}
map<string, double> variables = globalParameters;
if (interactionGroups.size() > 0) {
// The user has specified interaction groups, so compute only the requested interactions.
for (int group = 0; group < (int) interactionGroups.size(); group++) {
const set<int>& set1 = interactionGroups[group].first;
const set<int>& set2 = interactionGroups[group].second;
for (set<int>::const_iterator atom1 = set1.begin(); atom1 != set1.end(); ++atom1) {
for (set<int>::const_iterator atom2 = set2.begin(); atom2 != set2.end(); ++atom2) {
if (*atom1 == *atom2 || exclusions[*atom1].find(*atom2) != exclusions[*atom1].end())
continue; // This is an excluded interaction.
if (*atom1 > *atom2 && set1.find(*atom2) != set1.end() && set2.find(*atom1) != set2.end())
continue; // Both atoms are in both sets, so skip duplicate interactions.
for (int j = 0; j < (int) paramNames.size(); j++) {
variables[particleParamNames[j*2]] = atomParameters[*atom1][j];
variables[particleParamNames[j*2+1]] = atomParameters[*atom2][j];
}
calculateOneIxn(*atom1, *atom2, atomCoordinates, variables, forces, energyByAtom, totalEnergy);
}
}
}
}
else if (cutoff) {
// We are using a cutoff, so get the interactions from the neighbor list.
for (int i = 0; i < (int) neighborList->size(); i++) {
OpenMM::AtomPair pair = (*neighborList)[i];
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];
}
calculateOneIxn(pair.first, pair.second, atomCoordinates, variables, forces, energyByAtom, totalEnergy);
}
}
else {
// Every particle interacts with every other one.
for (int ii = 0; ii < numberOfAtoms; ii++) {
for (int jj = ii+1; jj < numberOfAtoms; jj++) {
if (exclusions[jj].find(ii) == exclusions[jj].end()) {
for (int j = 0; j < (int) paramNames.size(); j++) {
variables[particleParamNames[j*2]] = atomParameters[ii][j];
variables[particleParamNames[j*2+1]] = atomParameters[jj][j];
}
calculateOneIxn(ii, jj, atomCoordinates, variables, forces, energyByAtom, totalEnergy);
}
}
}
}
}
/**---------------------------------------------------------------------------------------
......
platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.cpp
View file @ 81d60b1e
......@@ -37,6 +37,7 @@
// make sure that erf() and erfc() are defined.
#include "openmm/internal/MSVC_erfc.h"
using std::set;
using std::vector;
using OpenMM::RealVec;
......@@ -169,10 +170,8 @@ void ReferenceLJCoulombIxn::setUseSwitchingFunction( RealOpenMM distance ) {
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param atomParameters atom parameters atomParameters[atomIndex][paramterIndex]
@param exclusions atom exclusion indices exclusions[atomIndex][atomToExcludeIndex]
exclusions[atomIndex][0] = number of exclusions
exclusions[atomIndex][1-no.] = atom indices of atoms to excluded from
interacting w/ atom atomIndex
@param exclusions atom exclusion indices
exclusions[atomIndex] contains the list of exclusions for that atom
@param fixedParameters non atom parameters (not currently used)
@param forces force array (forces added)
@param energyByAtom atom energy
......@@ -183,7 +182,7 @@ void ReferenceLJCoulombIxn::setUseSwitchingFunction( RealOpenMM distance ) {
--------------------------------------------------------------------------------------- */
void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>& atomCoordinates,
RealOpenMM** atomParameters, int** exclusions,
RealOpenMM** atomParameters, vector<set<int> >& exclusions,
RealOpenMM* fixedParameters, vector<RealVec>& forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy, bool includeDirect, bool includeReciprocal) const {
typedef std::complex<RealOpenMM> d_complex;
......@@ -423,10 +422,10 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
RealOpenMM totalExclusionEnergy = 0.0f;
for (int i = 0; i < numberOfAtoms; i++)
for (int j = 1; j <= exclusions[i][0]; j++)
if (exclusions[i][j] > i) {
for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter) {
if (*iter > i) {
int ii = i;
int jj = exclusions[i][j];
int jj = *iter;
RealOpenMM deltaR[2][ReferenceForce::LastDeltaRIndex];
ReferenceForce::getDeltaR( atomCoordinates[jj], atomCoordinates[ii], deltaR[0] );
......@@ -454,6 +453,7 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
energyByAtom[jj] -= realSpaceEwaldEnergy;
}
}
}
if( totalEnergy )
*totalEnergy -= totalExclusionEnergy;
......@@ -467,10 +467,8 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param atomParameters atom parameters atomParameters[atomIndex][paramterIndex]
@param exclusions atom exclusion indices exclusions[atomIndex][atomToExcludeIndex]
exclusions[atomIndex][0] = number of exclusions
exclusions[atomIndex][1-no.] = atom indices of atoms to excluded from
interacting w/ atom atomIndex
@param exclusions atom exclusion indices
exclusions[atomIndex] contains the list of exclusions for that atom
@param fixedParameters non atom parameters (not currently used)
@param forces force array (forces added)
@param energyByAtom atom energy
......@@ -481,7 +479,7 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
--------------------------------------------------------------------------------------- */
void ReferenceLJCoulombIxn::calculatePairIxn(int numberOfAtoms, vector<RealVec>& atomCoordinates,
RealOpenMM** atomParameters, int** exclusions,
RealOpenMM** atomParameters, vector<set<int> >& exclusions,
RealOpenMM* fixedParameters, vector<RealVec>& forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy, bool includeDirect, bool includeReciprocal) const {
......@@ -499,32 +497,13 @@ void ReferenceLJCoulombIxn::calculatePairIxn(int numberOfAtoms, vector<RealVec>&
}
}
else {
// allocate and initialize exclusion array
int* exclusionIndices = new int[numberOfAtoms];
for( int ii = 0; ii < numberOfAtoms; ii++ ){
exclusionIndices[ii] = -1;
}
for( int ii = 0; ii < numberOfAtoms; ii++ ){
// set exclusions
for( int jj = 1; jj <= exclusions[ii][0]; jj++ ){
exclusionIndices[exclusions[ii][jj]] = ii;
}
// loop over atom pairs
for( int jj = ii+1; jj < numberOfAtoms; jj++ ){
if( exclusionIndices[jj] != ii ){
calculateOneIxn(ii, jj, atomCoordinates, atomParameters, forces, energyByAtom, totalEnergy);
}
}
for( int jj = ii+1; jj < numberOfAtoms; jj++ )
if (exclusions[jj].find(ii) == exclusions[jj].end())
calculateOneIxn(ii, jj, atomCoordinates, atomParameters, forces, energyByAtom, totalEnergy);
}
delete[] exclusionIndices;
}
}
......
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