Skip to content

GitLab

Commit f5acdd7a authored by peastman's avatar peastman
Browse files

CustomManyParticleForce offers two different "permutation modes". Implemented... 

CustomManyParticleForce offers two different "permutation modes".  Implemented it for Reference and CPU platforms.
parent bf464848
Hide whitespace changes
Inline Side-by-side
Showing with 409 additions and 38 deletions
openmmapi/include/openmm/CustomManyParticleForce.h
View file @ f5acdd7a
......@@ -75,8 +75,8 @@ namespace OpenMM {
* you can modify its parameters by calling setParticleParameters(). This will have no effect on Contexts that already exist
* unless you call updateParametersInContext().
*
* Multi-particle interactions can be very expensive to evaluate, so they are usually used with a cutoff distance. If two particles
* are further apart than the cutoff, <i>all</i> sets that include those two particles will be omitted.
* Multi-particle interactions can be very expensive to evaluate, so they are usually used with a cutoff distance. The exact
* interpretation of the cutoff depends on the permutation mode, as discussed below.
*
* CustomManyParticleForce also lets you specify "exclusions", particular pairs of particles whose interactions should be
* omitted from force and energy calculations. This is most often used for particles that are bonded to each other.
......@@ -89,7 +89,8 @@ namespace OpenMM {
* "C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
* "theta1=angle(p1,p2,p3); theta2=angle(p2,p3,p1); theta3=angle(p3,p1,p2);"
* "r12=distance(p1,p2); r13=distance(p1,p3); r23=distance(p2,p3)");
* <pre></tt>
* force->setPermutationMode(CustomManyParticleForce::SinglePermutation);
* </pre></tt>
*
* This force depends on one parameter, C. The following code defines it as a global parameter:
*
......@@ -97,16 +98,38 @@ namespace OpenMM {
* force->addGlobalParameter("C", 1.0);
* </pre></tt>
*
* The expression <i>must</i> be symmetric with respect to the particles. It typically will only be evaluated once for
* each set of particles, and no guarantee is made about which particle will be identified as "particle 1". In the above
* example, the energy only depends on the products cos(theta1)*cos(theta2)*cos(theta3) and r12*r13*r23, both of which are
* unchanged if the labels p1, p2, and p3 are permuted. If that were not true, the results would be undefined, because
* Notice that the expression is symmetric with respect to the particles. It only depends on the products
* cos(theta1)*cos(theta2)*cos(theta3) and r12*r13*r23, both of which are unchanged if the labels p1, p2, and p3 are permuted.
* This is required because we specified SinglePermutation as the permutation mode. (This is the default, so we did not
* really need to set it, but doing so makes the example clearer.) In this mode, the expression is only evaluated once for
* each set of particles. No guarantee is made about which particle will be identified as p1, p2, etc. Therefore, the
* energy <i>must</i> be symmetric with respect to exchange of particles. Otherwise, the results would be undefined because
* permuting the labels would change the energy.
*
* In some cases this requirement is overly restrictive. When some particles are fundamentally different from others,
* the expression may be inherently non-symmetric. An example would be a water model that involves three particles,
* two of which <i>must</i> be hydrogen and one of which <i>must</i> be oxygen. Cases like this can be implemented
* using particle types.
* Not all many-particle interactions work this way. Another common pattern is for the expression to describe an interaction
* between one central particle and other nearby particles. An example of this is the 3-particle piece of the Stillinger-Weber
* potential:
*
* <tt><pre>CustomManyParticleForce* force = new CustomManyParticleForce(3,
* "L*eps*(cos(theta1)+1/3)^2*exp(sigma*gamma/(r12-a*sigma))*exp(sigma*gamma/(r13-a*sigma));"
"r12 = distance(p1,p2); r13 = distance(p1,p3); theta1 = angle(p3,p1,p2)");
* force->setPermutationMode(CustomManyParticleForce::UniqueCentralParticle);
* </pre></tt>
*
* When the permutation mode is set to UniqueCentralParticle, particle p1 is treated as the central particle. For a set of
* N particles, the expression is evaluated N times, once with each particle as p1. The expression can therefore treat
* p1 differently from the other particles. Notice that it is still symmetric with respect to p2 and p3, however. There
* is no guarantee about how those labels will be assigned to particles.
*
* Distance cutoffs are applied in different ways depending on the permutation mode. In SinglePermutation mode, every particle
* in the set must be within the cutoff distance of every other particle. If <i>any</i> two particles are further apart than
* the cutoff distance, the interaction is skipped. In UniqueCentralParticle mode, each particle must be within the cutoff
* distance of the central particle, but not necessarily of all the other particles. The cutoff may therefore exclude a subset
* of the permutations of a set of particles.
*
* Another common situation is that some particles are fundamentally different from others, causing the expression to be
* inherently non-symmetric. An example would be a water model that involves three particles, two of which <i>must</i> be
* hydrogen and one of which <i>must</i> be oxygen. Cases like this can be implemented using particle types.
*
* A particle type is an integer that you specify when you call addParticle(). (If you omit the argument, it defaults
* to 0.) For the water model, you could specify 0 for all oxygen atoms and 1 for all hydrogen atoms. You can then
......@@ -157,6 +180,26 @@ public:
*/
CutoffPeriodic = 2,
};
/**
* This is an enumeration of the different modes for selecting which permutations of a set of particles to evaluate the
* interaction for.
*/
enum PermutationMode {
/**
* For any set of particles, the interaction is evaluated only once for a single permutation of the particles.
* There is no guarantee about which permutation will be used (aside from the requirement to satisfy type filters),
* so the expression must be symmetric. If cutoffs are used, then every particle in the set must be within the
* cutoff distance of every other particle.
*/
SinglePermutation = 0,
/**
* The interaction is treated as an interaction between one central particle (p1) and various other nearby particles
* (p2, p3, ...). For a set of N particles it will be evaluated N times, once with each particle as p1. The expression
* must be symmetric with respect to the other particles, but may treat p1 differently. If cutoffs are used, then
* every particle must be within the cutoff distance of p1.
*/
UniqueCentralParticle = 1
};
/**
* Create a CustomManyParticleForce.
*
......@@ -218,6 +261,14 @@ public:
* Set the method used for handling long range nonbonded interactions.
*/
void setNonbondedMethod(NonbondedMethod method);
/**
* Get the mode that selects which permutations of a set of particles to evaluate the interaction for.
*/
PermutationMode getPermutationMode() const;
/**
* Set the mode that selects which permutations of a set of particles to evaluate the interaction for.
*/
void setPermutationMode(PermutationMode mode);
/**
* Get the cutoff distance (in nm) being used for nonbonded interactions. If the NonbondedMethod in use
* is NoCutoff, this value will have no effect.
......@@ -420,6 +471,7 @@ private:
class FunctionInfo;
int particlesPerSet;
NonbondedMethod nonbondedMethod;
PermutationMode permutationMode;
double cutoffDistance;
std::string energyExpression;
std::vector<ParticleParameterInfo> particleParameters;
......
openmmapi/src/CustomManyParticleForce.cpp
View file @ f5acdd7a
......@@ -44,7 +44,7 @@ using namespace OpenMM;
using namespace std;
CustomManyParticleForce::CustomManyParticleForce(int particlesPerSet, const string& energy) :
particlesPerSet(particlesPerSet), energyExpression(energy), nonbondedMethod(NoCutoff), cutoffDistance(1.0), typeFilters(particlesPerSet) {
particlesPerSet(particlesPerSet), energyExpression(energy), nonbondedMethod(NoCutoff), permutationMode(SinglePermutation), cutoffDistance(1.0), typeFilters(particlesPerSet) {
}
CustomManyParticleForce::~CustomManyParticleForce() {
......@@ -68,6 +68,14 @@ void CustomManyParticleForce::setNonbondedMethod(NonbondedMethod method) {
nonbondedMethod = method;
}
CustomManyParticleForce::PermutationMode CustomManyParticleForce::getPermutationMode() const {
return permutationMode;
}
void CustomManyParticleForce::setPermutationMode(PermutationMode mode) {
permutationMode = mode;
}
double CustomManyParticleForce::getCutoffDistance() const {
return cutoffDistance;
}
......
platforms/cpu/include/CpuCustomManyParticleForce.h
View file @ f5acdd7a
......@@ -51,7 +51,7 @@ private:
class ComputeForceTask;
class ThreadData;
int numParticles, numParticlesPerSet, numPerParticleParameters, numTypes;
bool useCutoff, usePeriodic;
bool useCutoff, usePeriodic, centralParticleMode;
RealOpenMM cutoffDistance;
RealOpenMM periodicBoxSize[3];
CpuNeighborList* neighborList;
......@@ -60,6 +60,7 @@ private:
std::vector<int> particleTypes;
std::vector<int> orderIndex;
std::vector<std::vector<int> > particleOrder;
std::vector<std::vector<int> > particleNeighbors;
std::vector<ThreadData*> threadData;
// The following variables are used to make information accessible to the individual threads.
float* posq;
......
platforms/cpu/src/CpuCustomManyParticleForce.cpp
View file @ f5acdd7a
......@@ -54,6 +54,7 @@ CpuCustomManyParticleForce::CpuCustomManyParticleForce(const CustomManyParticleF
numParticles = force.getNumParticles();
numParticlesPerSet = force.getNumParticlesPerSet();
numPerParticleParameters = force.getNumPerParticleParameters();
centralParticleMode = (force.getPermutationMode() == CustomManyParticleForce::UniqueCentralParticle);
// Create custom functions for the tabulated functions.
......@@ -114,10 +115,31 @@ void CpuCustomManyParticleForce::calculateIxn(AlignedArray<float>& posq, RealOpe
gmx_atomic_set(&counter, 0);
this->atomicCounter = &counter;
if (useCutoff) {
// Construct a neighbor list.
// Construct a neighbor list. We use CpuNeighborList to do this, but then copy the result
// into a new data structure. This is needed because in UniqueCentralParticle mode, the
// the neighbor list needs to include symmetric pairs.
particleNeighbors.resize(numParticles);
for (int i = 0; i < numParticles; i++)
particleNeighbors[i].clear();
float boxSizeFloat[] = {(float) periodicBoxSize[0], (float) periodicBoxSize[1], (float) periodicBoxSize[2]};
neighborList->computeNeighborList(numParticles, posq, exclusions, boxSizeFloat, usePeriodic, cutoffDistance, threads);
for (int blockIndex = 0; blockIndex < neighborList->getNumBlocks(); blockIndex++) {
const vector<int>& neighbors = neighborList->getBlockNeighbors(blockIndex);
const vector<char>& exclusions = neighborList->getBlockExclusions(blockIndex);
int numNeighbors = neighbors.size();
for (int i = 0; i < 4; i++) {
int p1 = neighborList->getSortedAtoms()[4*blockIndex+i];
for (int j = 0; j < numNeighbors; j++) {
if ((exclusions[j] & (1<<i)) == 0) {
int p2 = neighbors[j];
particleNeighbors[p1].push_back(p2);
if (centralParticleMode)
particleNeighbors[p2].push_back(p1);
}
}
}
}
}
// Signal the threads to start running and wait for them to finish.
......@@ -147,27 +169,12 @@ void CpuCustomManyParticleForce::threadComputeForce(ThreadPool& threads, int thr
if (useCutoff) {
// Loop over interactions from the neighbor list.
vector<int> particles;
while (true) {
int blockIndex = gmx_atomic_fetch_add(reinterpret_cast<gmx_atomic_t*>(atomicCounter), 1);
if (blockIndex >= neighborList->getNumBlocks())
int i = gmx_atomic_fetch_add(reinterpret_cast<gmx_atomic_t*>(atomicCounter), 1);
if (i >= numParticles)
break;
const vector<int>& neighbors = neighborList->getBlockNeighbors(blockIndex);
const vector<char>& exclusions = neighborList->getBlockExclusions(blockIndex);
int numNeighbors = neighbors.size();
for (int i = 0; i < 4; i++) {
particleIndices[0] = neighborList->getSortedAtoms()[4*blockIndex+i];
// Build a filtered list of neighbors after removing exclusions. We'll check for actual exclusions
// again later, but the neighbor list also includes padding atoms that it marks as exclusions, so
// we need to remove those now.
particles.resize(0);
for (int j = 0; j < numNeighbors; j++)
if ((exclusions[j] & (1<<i)) == 0)
particles.push_back(neighbors[j]);
loopOverInteractions(particles, particleIndices, 1, 0, particleParameters, forces, data, boxSize, invBoxSize);
}
particleIndices[0] = i;
loopOverInteractions(particleNeighbors[i], particleIndices, 1, 0, particleParameters, forces, data, boxSize, invBoxSize);
}
}
else {
......@@ -181,7 +188,8 @@ void CpuCustomManyParticleForce::threadComputeForce(ThreadPool& threads, int thr
if (i >= numParticles)
break;
particleIndices[0] = i;
loopOverInteractions(particles, particleIndices, 1, i+1, particleParameters, forces, data, boxSize, invBoxSize);
int startIndex = (centralParticleMode ? 0 : i+1);
loopOverInteractions(particles, particleIndices, 1, startIndex, particleParameters, forces, data, boxSize, invBoxSize);
}
}
}
......@@ -208,6 +216,7 @@ void CpuCustomManyParticleForce::loopOverInteractions(vector<int>& availablePart
RealOpenMM** particleParameters, float* forces, ThreadData& data, const fvec4& boxSize, const fvec4& invBoxSize) {
int numParticles = availableParticles.size();
double cutoff2 = cutoffDistance*cutoffDistance;
int checkRange = (centralParticleMode ? 1 : loopIndex);
for (int i = startIndex; i < numParticles; i++) {
int particle = availableParticles[i];
......@@ -218,7 +227,7 @@ void CpuCustomManyParticleForce::loopOverInteractions(vector<int>& availablePart
fvec4 deltaR;
fvec4 pos1(posq+4*particle);
float r2;
for (int j = 0; j < loopIndex && include; j++) {
for (int j = 0; j < checkRange && include; j++) {
fvec4 pos2(posq+4*particleSet[j]);
computeDelta(pos1, pos2, deltaR, r2, boxSize, invBoxSize);
include &= (r2 < cutoff2);
......@@ -227,6 +236,8 @@ void CpuCustomManyParticleForce::loopOverInteractions(vector<int>& availablePart
for (int j = 0; j < loopIndex && include; j++)
include &= (exclusions[particle].find(particleSet[j]) == exclusions[particle].end());
if (include) {
if (loopIndex > 0 && availableParticles[i] == particleSet[0])
continue;
particleSet[loopIndex] = availableParticles[i];
if (loopIndex == numParticlesPerSet-1)
calculateOneIxn(particleSet, particleParameters, forces, data, boxSize, invBoxSize);
......
platforms/cpu/tests/TestCpuCustomManyParticleForce.cpp
View file @ f5acdd7a
......@@ -121,6 +121,77 @@ void validateAxilrodTeller(CustomManyParticleForce* force, const vector<Vec3>& p
ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state3.getPotentialEnergy())/stepSize, 1e-4);
}
void validateStillingerWeber(CustomManyParticleForce* force, const vector<Vec3>& positions, const vector<const int*>& expectedSets, double boxSize) {
// Create a System and Context.
int numParticles = force->getNumParticles();
CustomManyParticleForce::NonbondedMethod nonbondedMethod = force->getNonbondedMethod();
System system;
for (int i = 0; i < numParticles; i++)
system.addParticle(1.0);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
system.addForce(force);
VerletIntegrator integrator(0.001);
CpuPlatform platform;
Context context(system, integrator, platform);
context.setPositions(positions);
State state1 = context.getState(State::Forces | State::Energy);
double L = context.getParameter("L");
double eps = context.getParameter("eps");
double a = context.getParameter("a");
double gamma = context.getParameter("gamma");
double sigma = context.getParameter("sigma");
// See if the energy matches the expected value.
double expectedEnergy = 0;
for (int i = 0; i < (int) expectedSets.size(); i++) {
int p1 = expectedSets[i][0];
int p2 = expectedSets[i][1];
int p3 = expectedSets[i][2];
Vec3 d12 = positions[p2]-positions[p1];
Vec3 d13 = positions[p3]-positions[p1];
Vec3 d23 = positions[p3]-positions[p2];
if (nonbondedMethod == CustomManyParticleForce::CutoffPeriodic) {
for (int j = 0; j < 3; j++) {
d12[j] -= floor(d12[j]/boxSize+0.5f)*boxSize;
d13[j] -= floor(d13[j]/boxSize+0.5f)*boxSize;
d23[j] -= floor(d23[j]/boxSize+0.5f)*boxSize;
}
}
double r12 = sqrt(d12.dot(d12));
double r13 = sqrt(d13.dot(d13));
double r23 = sqrt(d23.dot(d23));
double ctheta1 = d12.dot(d13)/(r12*r13);
double ctheta2 = -d12.dot(d23)/(r12*r23);
double ctheta3 = d13.dot(d23)/(r13*r23);
expectedEnergy += L*eps*(ctheta1+1.0/3.0)*(ctheta1+1.0/3.0)*exp(sigma*gamma/(r12-a*sigma))*exp(sigma*gamma/(r13-a*sigma));
}
ASSERT_EQUAL_TOL(expectedEnergy, state1.getPotentialEnergy(), 1e-5);
// Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.
const vector<Vec3>& forces = state1.getForces();
double norm = 0.0;
for (int i = 0; i < (int) forces.size(); ++i)
norm += forces[i].dot(forces[i]);
norm = std::sqrt(norm);
const double stepSize = 1e-3;
double step = 0.5*stepSize/norm;
vector<Vec3> positions2(numParticles), positions3(numParticles);
for (int i = 0; i < (int) positions.size(); ++i) {
Vec3 p = positions[i];
Vec3 f = forces[i];
positions2[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
positions3[i] = Vec3(p[0]+f[0]*step, p[1]+f[1]*step, p[2]+f[2]*step);
}
context.setPositions(positions2);
State state2 = context.getState(State::Energy);
context.setPositions(positions3);
State state3 = context.getState(State::Energy);
ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state3.getPotentialEnergy())/stepSize, 1e-4);
}
void testNoCutoff() {
CustomManyParticleForce* force = new CustomManyParticleForce(3,
"C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
......@@ -504,6 +575,102 @@ void testLargeSystem() {
ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-4);
}
void testCentralParticleModeNoCutoff() {
CustomManyParticleForce* force = new CustomManyParticleForce(3,
"L*eps*(cos(theta1)+1/3)^2*exp(sigma*gamma/(r12-a*sigma))*exp(sigma*gamma/(r13-a*sigma));"
"r12 = distance(p1,p2); r13 = distance(p1,p3); theta1 = angle(p3,p1,p2)");
force->setPermutationMode(CustomManyParticleForce::UniqueCentralParticle);
force->addGlobalParameter("L", 23.13);
force->addGlobalParameter("eps", 25.894776);
force->addGlobalParameter("a", 1.8);
force->addGlobalParameter("sigma", 0.23925);
force->addGlobalParameter("gamma", 1.2);
vector<double> params;
force->addParticle(params);
force->addParticle(params);
force->addParticle(params);
force->addParticle(params);
vector<Vec3> positions;
positions.push_back(Vec3(0, 0, 0));
positions.push_back(Vec3(1, 0, 0));
positions.push_back(Vec3(0, 1.1, 0.3));
positions.push_back(Vec3(0.4, 0, -0.8));
int sets[12][3] = {{0,1,2}, {0,1,3}, {0,2,3}, {1,0,2}, {1,0,3}, {1, 2, 3}, {2,0,1}, {2,0,3}, {2, 1, 3}, {3,0,1}, {3,0,2}, {3,1,2}};
vector<const int*> expectedSets(&sets[0], &sets[12]);
validateStillingerWeber(force, positions, expectedSets, 2.0);
}
void testCentralParticleModeCutoff() {
CustomManyParticleForce* force = new CustomManyParticleForce(3,
"L*eps*(cos(theta1)+1/3)^2*exp(sigma*gamma/(r12-a*sigma))*exp(sigma*gamma/(r13-a*sigma));"
"r12 = distance(p1,p2); r13 = distance(p1,p3); theta1 = angle(p3,p1,p2)");
force->setPermutationMode(CustomManyParticleForce::UniqueCentralParticle);
force->addGlobalParameter("L", 23.13);
force->addGlobalParameter("eps", 25.894776);
force->addGlobalParameter("a", 1.8);
force->addGlobalParameter("sigma", 0.23925);
force->addGlobalParameter("gamma", 1.2);
force->setNonbondedMethod(CustomManyParticleForce::CutoffNonPeriodic);
force->setCutoffDistance(1.55);
vector<double> params;
force->addParticle(params);
force->addParticle(params);
force->addParticle(params);
force->addParticle(params);
vector<Vec3> positions;
positions.push_back(Vec3(0, 0, 0));
positions.push_back(Vec3(1, 0, 0));
positions.push_back(Vec3(0, 1.1, 0.3));
positions.push_back(Vec3(0.4, 0, -0.8));
int sets[8][3] = {{0,1,2}, {0,1,3}, {0,2,3}, {1,0,2}, {1,0,3}, {1, 2, 3}, {2,0,1}, {3,0,1}};
vector<const int*> expectedSets(&sets[0], &sets[8]);
validateStillingerWeber(force, positions, expectedSets, 2.0);
}
void testCentralParticleModeLargeSystem() {
int gridSize = 8;
int numParticles = gridSize*gridSize*gridSize;
double boxSize = 2.0;
double spacing = boxSize/gridSize;
CpuPlatform platform;
CustomManyParticleForce* force = new CustomManyParticleForce(3,
"L*eps*(cos(theta1)+1/3)^2*exp(sigma*gamma/(r12-a*sigma))*exp(sigma*gamma/(r13-a*sigma));"
"r12 = distance(p1,p2); r13 = distance(p1,p3); theta1 = angle(p3,p1,p2)");
force->setPermutationMode(CustomManyParticleForce::UniqueCentralParticle);
force->addGlobalParameter("L", 23.13);
force->addGlobalParameter("eps", 25.894776);
force->addGlobalParameter("a", 1.8);
force->addGlobalParameter("sigma", 0.23925);
force->addGlobalParameter("gamma", 1.2);
force->setNonbondedMethod(CustomManyParticleForce::CutoffPeriodic);
force->setCutoffDistance(1.8*0.23925);
vector<double> params;
vector<Vec3> positions;
System system;
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < gridSize; i++)
for (int j = 0; j < gridSize; j++)
for (int k = 0; k < gridSize; k++) {
force->addParticle(params);
positions.push_back(Vec3((i+0.4*genrand_real2(sfmt))*spacing, (j+0.4*genrand_real2(sfmt))*spacing, (k+0.4*genrand_real2(sfmt))*spacing));
system.addParticle(1.0);
}
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
system.addForce(force);
VerletIntegrator integrator1(0.001);
VerletIntegrator integrator2(0.001);
Context context1(system, integrator1, Platform::getPlatformByName("Reference"));
Context context2(system, integrator2, platform);
context1.setPositions(positions);
context2.setPositions(positions);
State state1 = context1.getState(State::Forces | State::Energy);
State state2 = context2.getState(State::Forces | State::Energy);
ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-4);
for (int i = 0; i < numParticles; i++)
ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-4);
}
int main() {
try {
testNoCutoff();
......@@ -515,6 +682,9 @@ int main() {
testTabulatedFunctions();
testTypeFilters();
testLargeSystem();
testCentralParticleModeNoCutoff();
testCentralParticleModeCutoff();
testCentralParticleModeLargeSystem();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
......
platforms/reference/include/ReferenceCustomManyParticleIxn.h
View file @ f5acdd7a
......@@ -44,7 +44,7 @@ class ReferenceCustomManyParticleIxn {
class AngleTermInfo;
class DihedralTermInfo;
int numParticlesPerSet, numPerParticleParameters, numTypes;
bool useCutoff, usePeriodic;
bool useCutoff, usePeriodic, centralParticleMode;
RealOpenMM cutoffDistance;
RealOpenMM periodicBoxSize[3];
Lepton::ExpressionProgram energyExpression;
......
platforms/reference/src/SimTKReference/ReferenceCustomManyParticleIxn.cpp
View file @ f5acdd7a
......@@ -41,6 +41,7 @@ using namespace std;
ReferenceCustomManyParticleIxn::ReferenceCustomManyParticleIxn(const CustomManyParticleForce& force) : useCutoff(false), usePeriodic(false) {
numParticlesPerSet = force.getNumParticlesPerSet();
numPerParticleParameters = force.getNumPerParticleParameters();
centralParticleMode = (force.getPermutationMode() == CustomManyParticleForce::UniqueCentralParticle);
// Create custom functions for the tabulated functions.
......@@ -134,8 +135,11 @@ void ReferenceCustomManyParticleIxn::loopOverInteractions(vector<int>& particles
RealOpenMM** particleParameters, map<string, double>& variables, vector<OpenMM::RealVec>& forces,
RealOpenMM* totalEnergy) const {
int numParticles = atomCoordinates.size();
int start = (loopIndex == 0 ? 0 : particles[loopIndex-1]+1);
int firstPartialLoop = (centralParticleMode ? 2 : 1);
int start = (loopIndex < firstPartialLoop ? 0 : particles[loopIndex-1]+1);
for (int i = start; i < numParticles; i++) {
if (loopIndex > 0 && i == particles[0])
continue;
particles[loopIndex] = i;
if (loopIndex == numParticlesPerSet-1)
calculateOneIxn(particles, atomCoordinates, particleParameters, variables, forces, totalEnergy);
......@@ -173,7 +177,7 @@ void ReferenceCustomManyParticleIxn::calculateOneIxn(const vector<int>& particle
int p2 = permutedParticles[j];
if (exclusions[p1].find(p2) != exclusions[p1].end())
return;
if (useCutoff) {
if (useCutoff && (i == 0 || !centralParticleMode)) {
RealOpenMM delta[ReferenceForce::LastDeltaRIndex];
computeDelta(p1, p2, delta, atomCoordinates);
if (delta[ReferenceForce::RIndex] >= cutoffDistance)
......
platforms/reference/tests/TestReferenceCustomManyParticleForce.cpp
View file @ f5acdd7a
......@@ -121,6 +121,77 @@ void validateAxilrodTeller(CustomManyParticleForce* force, const vector<Vec3>& p
ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state3.getPotentialEnergy())/stepSize, 1e-4);
}
void validateStillingerWeber(CustomManyParticleForce* force, const vector<Vec3>& positions, const vector<const int*>& expectedSets, double boxSize) {
// Create a System and Context.
int numParticles = force->getNumParticles();
CustomManyParticleForce::NonbondedMethod nonbondedMethod = force->getNonbondedMethod();
System system;
for (int i = 0; i < numParticles; i++)
system.addParticle(1.0);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
system.addForce(force);
VerletIntegrator integrator(0.001);
ReferencePlatform platform;
Context context(system, integrator, platform);
context.setPositions(positions);
State state1 = context.getState(State::Forces | State::Energy);
double L = context.getParameter("L");
double eps = context.getParameter("eps");
double a = context.getParameter("a");
double gamma = context.getParameter("gamma");
double sigma = context.getParameter("sigma");
// See if the energy matches the expected value.
double expectedEnergy = 0;
for (int i = 0; i < (int) expectedSets.size(); i++) {
int p1 = expectedSets[i][0];
int p2 = expectedSets[i][1];
int p3 = expectedSets[i][2];
Vec3 d12 = positions[p2]-positions[p1];
Vec3 d13 = positions[p3]-positions[p1];
Vec3 d23 = positions[p3]-positions[p2];
if (nonbondedMethod == CustomManyParticleForce::CutoffPeriodic) {
for (int j = 0; j < 3; j++) {
d12[j] -= floor(d12[j]/boxSize+0.5f)*boxSize;
d13[j] -= floor(d13[j]/boxSize+0.5f)*boxSize;
d23[j] -= floor(d23[j]/boxSize+0.5f)*boxSize;
}
}
double r12 = sqrt(d12.dot(d12));
double r13 = sqrt(d13.dot(d13));
double r23 = sqrt(d23.dot(d23));
double ctheta1 = d12.dot(d13)/(r12*r13);
double ctheta2 = -d12.dot(d23)/(r12*r23);
double ctheta3 = d13.dot(d23)/(r13*r23);
expectedEnergy += L*eps*(ctheta1+1.0/3.0)*(ctheta1+1.0/3.0)*exp(sigma*gamma/(r12-a*sigma))*exp(sigma*gamma/(r13-a*sigma));
}
ASSERT_EQUAL_TOL(expectedEnergy, state1.getPotentialEnergy(), 1e-5);
// Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.
const vector<Vec3>& forces = state1.getForces();
double norm = 0.0;
for (int i = 0; i < (int) forces.size(); ++i)
norm += forces[i].dot(forces[i]);
norm = std::sqrt(norm);
const double stepSize = 1e-3;
double step = 0.5*stepSize/norm;
vector<Vec3> positions2(numParticles), positions3(numParticles);
for (int i = 0; i < (int) positions.size(); ++i) {
Vec3 p = positions[i];
Vec3 f = forces[i];
positions2[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
positions3[i] = Vec3(p[0]+f[0]*step, p[1]+f[1]*step, p[2]+f[2]*step);
}
context.setPositions(positions2);
State state2 = context.getState(State::Energy);
context.setPositions(positions3);
State state3 = context.getState(State::Energy);
ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state3.getPotentialEnergy())/stepSize, 1e-4);
}
void testNoCutoff() {
CustomManyParticleForce* force = new CustomManyParticleForce(3,
"C*(1+3*cos(theta1)*cos(theta2)*cos(theta3))/(r12*r13*r23)^3;"
......@@ -464,6 +535,58 @@ void testTypeFilters() {
ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), 1e-5);
}
void testCentralParticleModeNoCutoff() {
CustomManyParticleForce* force = new CustomManyParticleForce(3,
"L*eps*(cos(theta1)+1/3)^2*exp(sigma*gamma/(r12-a*sigma))*exp(sigma*gamma/(r13-a*sigma));"
"r12 = distance(p1,p2); r13 = distance(p1,p3); theta1 = angle(p3,p1,p2)");
force->setPermutationMode(CustomManyParticleForce::UniqueCentralParticle);
force->addGlobalParameter("L", 23.13);
force->addGlobalParameter("eps", 25.894776);
force->addGlobalParameter("a", 1.8);
force->addGlobalParameter("sigma", 0.23925);
force->addGlobalParameter("gamma", 1.2);
vector<double> params;
force->addParticle(params);
force->addParticle(params);
force->addParticle(params);
force->addParticle(params);
vector<Vec3> positions;
positions.push_back(Vec3(0, 0, 0));
positions.push_back(Vec3(1, 0, 0));
positions.push_back(Vec3(0, 1.1, 0.3));
positions.push_back(Vec3(0.4, 0, -0.8));
int sets[12][3] = {{0,1,2}, {0,1,3}, {0,2,3}, {1,0,2}, {1,0,3}, {1, 2, 3}, {2,0,1}, {2,0,3}, {2, 1, 3}, {3,0,1}, {3,0,2}, {3,1,2}};
vector<const int*> expectedSets(&sets[0], &sets[12]);
validateStillingerWeber(force, positions, expectedSets, 2.0);
}
void testCentralParticleModeCutoff() {
CustomManyParticleForce* force = new CustomManyParticleForce(3,
"L*eps*(cos(theta1)+1/3)^2*exp(sigma*gamma/(r12-a*sigma))*exp(sigma*gamma/(r13-a*sigma));"
"r12 = distance(p1,p2); r13 = distance(p1,p3); theta1 = angle(p3,p1,p2)");
force->setPermutationMode(CustomManyParticleForce::UniqueCentralParticle);
force->addGlobalParameter("L", 23.13);
force->addGlobalParameter("eps", 25.894776);
force->addGlobalParameter("a", 1.8);
force->addGlobalParameter("sigma", 0.23925);
force->addGlobalParameter("gamma", 1.2);
force->setNonbondedMethod(CustomManyParticleForce::CutoffNonPeriodic);
force->setCutoffDistance(1.55);
vector<double> params;
force->addParticle(params);
force->addParticle(params);
force->addParticle(params);
force->addParticle(params);
vector<Vec3> positions;
positions.push_back(Vec3(0, 0, 0));
positions.push_back(Vec3(1, 0, 0));
positions.push_back(Vec3(0, 1.1, 0.3));
positions.push_back(Vec3(0.4, 0, -0.8));
int sets[8][3] = {{0,1,2}, {0,1,3}, {0,2,3}, {1,0,2}, {1,0,3}, {1, 2, 3}, {2,0,1}, {3,0,1}};
vector<const int*> expectedSets(&sets[0], &sets[8]);
validateStillingerWeber(force, positions, expectedSets, 2.0);
}
int main() {
try {
testNoCutoff();
......@@ -474,6 +597,8 @@ int main() {
testParameters();
testTabulatedFunctions();
testTypeFilters();
testCentralParticleModeNoCutoff();
testCentralParticleModeCutoff();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
......
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