Skip to content

GitLab

Commit ae07d525 authored by peastman's avatar peastman
Browse files

NonbondedForce can optionally apply periodic boundary conditions to exceptions

parent 40c04f74
Hide whitespace changes
Inline Side-by-side
Showing with 190 additions and 65 deletions
openmmapi/include/openmm/NonbondedForce.h
View file @ ae07d525
......@@ -9,7 +9,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2018 Stanford University and the Authors. *
* Portions copyright (c) 2008-2020 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -568,6 +568,32 @@ public:
nonbondedMethod == NonbondedForce::PME ||
nonbondedMethod == NonbondedForce::LJPME;
}
/**
* Get whether periodic boundary conditions should be applied to exceptions. Usually this is not
* appropriate, because exceptions are normally used to represent bonded interactions (1-2, 1-3, and
* 1-4 pairs), but there are situations when it does make sense. For example, you may want to simulate
* an infinite chain where one end of a molecule is bonded to the opposite end of the next periodic
* copy.
*
* Regardless of this value, periodic boundary conditions are only applied to exceptions if they also
* are applied to other interactions. If the nonbonded method is NoCutoff or CutoffNonPeriodic, this
* value is ignored. Also note that cutoffs are never applied to exceptions, again because they are
* normally used to represent bonded interactions.
*/
bool getExceptionsUsePeriodicBoundaryConditions() const;
/**
* Set whether periodic boundary conditions should be applied to exceptions. Usually this is not
* appropriate, because exceptions are normally used to represent bonded interactions (1-2, 1-3, and
* 1-4 pairs), but there are situations when it does make sense. For example, you may want to simulate
* an infinite chain where one end of a molecule is bonded to the opposite end of the next periodic
* copy.
*
* Regardless of this value, periodic boundary conditions are only applied to exceptions if they also
* get applied to other interactions. If the nonbonded method is NoCutoff or CutoffNonPeriodic, this
* value is ignored. Also note that cutoffs are never applied to exceptions, again because they are
* normally used to represent bonded interactions.
*/
void setExceptionsUsePeriodicBoundaryConditions(bool periodic);
protected:
ForceImpl* createImpl() const;
private:
......@@ -578,7 +604,7 @@ private:
class ExceptionOffsetInfo;
NonbondedMethod nonbondedMethod;
double cutoffDistance, switchingDistance, rfDielectric, ewaldErrorTol, alpha, dalpha;
bool useSwitchingFunction, useDispersionCorrection;
bool useSwitchingFunction, useDispersionCorrection, exceptionsUsePeriodic;
int recipForceGroup, nx, ny, nz, dnx, dny, dnz;
void addExclusionsToSet(const std::vector<std::set<int> >& bonded12, std::set<int>& exclusions, int baseParticle, int fromParticle, int currentLevel) const;
int getGlobalParameterIndex(const std::string& parameter) const;
......
openmmapi/src/NonbondedForce.cpp
View file @ ae07d525
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2018 Stanford University and the Authors. *
* Portions copyright (c) 2008-2020 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -48,7 +48,7 @@ using std::stringstream;
using std::vector;
NonbondedForce::NonbondedForce() : nonbondedMethod(NoCutoff), cutoffDistance(1.0), switchingDistance(-1.0), rfDielectric(78.3),
ewaldErrorTol(5e-4), alpha(0.0), dalpha(0.0), useSwitchingFunction(false), useDispersionCorrection(true), recipForceGroup(-1),
ewaldErrorTol(5e-4), alpha(0.0), dalpha(0.0), useSwitchingFunction(false), useDispersionCorrection(true), exceptionsUsePeriodic(false), recipForceGroup(-1),
nx(0), ny(0), nz(0), dnx(0), dny(0), dnz(0) {
}
......@@ -347,3 +347,11 @@ void NonbondedForce::setReciprocalSpaceForceGroup(int group) {
void NonbondedForce::updateParametersInContext(Context& context) {
dynamic_cast<NonbondedForceImpl&>(getImplInContext(context)).updateParametersInContext(getContextImpl(context));
}
bool NonbondedForce::getExceptionsUsePeriodicBoundaryConditions() const {
return exceptionsUsePeriodic;
}
void NonbondedForce::setExceptionsUsePeriodicBoundaryConditions(bool periodic) {
exceptionsUsePeriodic = periodic;
}
platforms/cpu/include/CpuKernels.h
View file @ ae07d525
......@@ -274,7 +274,7 @@ private:
std::vector<std::vector<double> > bonded14ParamArray;
double nonbondedCutoff, switchingDistance, rfDielectric, ewaldAlpha, ewaldDispersionAlpha, ewaldSelfEnergy, dispersionCoefficient;
int kmax[3], gridSize[3], dispersionGridSize[3];
bool useSwitchingFunction, useOptimizedPme, hasInitializedPme, hasInitializedDispersionPme, hasParticleOffsets, hasExceptionOffsets;
bool useSwitchingFunction, exceptionsArePeriodic, useOptimizedPme, hasInitializedPme, hasInitializedDispersionPme, hasParticleOffsets, hasExceptionOffsets;
std::vector<std::set<int> > exclusions;
std::vector<std::pair<float, float> > particleParams;
std::vector<float> C6params;
......
platforms/cpu/src/CpuKernels.cpp
View file @ ae07d525
......@@ -609,6 +609,10 @@ void CpuCalcNonbondedForceKernel::initialize(const System& system, const Nonbond
ewaldDispersionAlpha = alpha;
useSwitchingFunction = false;
}
if (nonbondedMethod == NoCutoff || nonbondedMethod == CutoffNonPeriodic)
exceptionsArePeriodic = false;
else
exceptionsArePeriodic = force.getExceptionsUsePeriodicBoundaryConditions();
rfDielectric = force.getReactionFieldDielectric();
if (force.getUseDispersionCorrection())
dispersionCoefficient = NonbondedForceImpl::calcDispersionCorrection(system, force);
......@@ -699,6 +703,10 @@ double CpuCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeFo
energy += nonbondedEnergy;
if (includeDirect) {
ReferenceLJCoulomb14 nonbonded14;
if (exceptionsArePeriodic) {
Vec3* boxVectors = extractBoxVectors(context);
nonbonded14.setPeriodic(boxVectors);
}
bondForce.calculateForce(posData, bonded14ParamArray, forceData, includeEnergy ? &energy : NULL, nonbonded14);
if (data.isPeriodic && nonbondedMethod != LJPME)
energy += dispersionCoefficient/(boxVectors[0][0]*boxVectors[1][1]*boxVectors[2][2]);
......
platforms/cuda/src/CudaKernels.cpp
View file @ ae07d525
......@@ -1104,6 +1104,7 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
}
baseExceptionParams.upload(baseExceptionParamsVec);
map<string, string> replacements;
replacements["APPLY_PERIODIC"] = (usePeriodic && force.getExceptionsUsePeriodicBoundaryConditions() ? "1" : "0");
replacements["PARAMS"] = cu.getBondedUtilities().addArgument(exceptionParams.getDevicePointer(), "float4");
cu.getBondedUtilities().addInteraction(atoms, cu.replaceStrings(CudaKernelSources::nonbondedExceptions, replacements), force.getForceGroup());
}
......
platforms/cuda/src/kernels/nonbondedExceptions.cu
View file @ ae07d525
float4 exceptionParams = PARAMS[index];
real3 delta = make_real3(pos2.x-pos1.x, pos2.y-pos1.y, pos2.z-pos1.z);
#if APPLY_PERIODIC
APPLY_PERIODIC_TO_DELTA(delta)
#endif
real r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
real invR = RSQRT(r2);
real sig2 = invR*exceptionParams.y;
......
platforms/opencl/src/OpenCLKernels.cpp
View file @ ae07d525
......@@ -1054,6 +1054,7 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
}
baseExceptionParams.upload(baseExceptionParamsVec);
map<string, string> replacements;
replacements["APPLY_PERIODIC"] = (usePeriodic && force.getExceptionsUsePeriodicBoundaryConditions() ? "1" : "0");
replacements["PARAMS"] = cl.getBondedUtilities().addArgument(exceptionParams.getDeviceBuffer(), "float4");
cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::nonbondedExceptions, replacements), force.getForceGroup());
}
......
platforms/opencl/src/kernels/nonbondedExceptions.cl
View file @ ae07d525
float4 exceptionParams = PARAMS[index];
real4 delta = pos2-pos1;
#if APPLY_PERIODIC
APPLY_PERIODIC_TO_DELTA(delta)
#endif
real r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
real invR = RSQRT(r2);
real sig2 = invR*exceptionParams.y;
......
platforms/reference/include/ReferenceKernels.h
View file @ ae07d525
......@@ -631,7 +631,7 @@ private:
std::map<std::pair<std::string, int>, std::array<double, 3> > particleParamOffsets, exceptionParamOffsets;
double nonbondedCutoff, switchingDistance, rfDielectric, ewaldAlpha, ewaldDispersionAlpha, dispersionCoefficient;
int kmax[3], gridSize[3], dispersionGridSize[3];
bool useSwitchingFunction;
bool useSwitchingFunction, exceptionsArePeriodic;
std::vector<std::set<int> > exclusions;
NonbondedMethod nonbondedMethod;
NeighborList* neighborList;
......
platforms/reference/include/ReferenceLJCoulomb14.h
View file @ ae07d525
......@@ -32,40 +32,53 @@ namespace OpenMM {
class OPENMM_EXPORT ReferenceLJCoulomb14 : public ReferenceBondIxn {
public:
/**---------------------------------------------------------------------------------------
Constructor
--------------------------------------------------------------------------------------- */
ReferenceLJCoulomb14();
/**---------------------------------------------------------------------------------------
Destructor
--------------------------------------------------------------------------------------- */
~ReferenceLJCoulomb14();
/**---------------------------------------------------------------------------------------
Calculate Ryckaert-Bellemans bond ixn
@param atomIndices atom indices of 4 atoms in bond
@param atomCoordinates atom coordinates
@param parameters six RB parameters
@param forces force array (forces added to current values)
@param totalEnergy if not null, the energy will be added to this
--------------------------------------------------------------------------------------- */
void calculateBondIxn(std::vector<int>& atomIndices, std::vector<OpenMM::Vec3>& atomCoordinates,
std::vector<double>& parameters, std::vector<OpenMM::Vec3>& forces,
double* totalEnergy, double* energyParamDerivs);
public:
/**---------------------------------------------------------------------------------------
Constructor
--------------------------------------------------------------------------------------- */
ReferenceLJCoulomb14();
/**---------------------------------------------------------------------------------------
Destructor
--------------------------------------------------------------------------------------- */
~ReferenceLJCoulomb14();
/**---------------------------------------------------------------------------------------
Set the force to use periodic boundary conditions.
@param vectors the vectors defining the periodic box
--------------------------------------------------------------------------------------- */
void setPeriodic(OpenMM::Vec3* vectors);
/**---------------------------------------------------------------------------------------
Calculate nonbonded 1-4 interactinos
@param atomIndices atom indices of the atoms in each pair
@param atomCoordinates atom coordinates
@param parameters (sigma, 4*epsilon, charge product) for each pair
@param forces force array (forces added to current values)
@param totalEnergy if not null, the energy will be added to this
--------------------------------------------------------------------------------------- */
void calculateBondIxn(std::vector<int>& atomIndices, std::vector<OpenMM::Vec3>& atomCoordinates,
std::vector<double>& parameters, std::vector<OpenMM::Vec3>& forces,
double* totalEnergy, double* energyParamDerivs);
private:
bool periodic;
OpenMM::Vec3 periodicBoxVectors[3];
};
} // namespace OpenMM
......
platforms/reference/src/ReferenceKernels.cpp
View file @ ae07d525
......@@ -991,6 +991,10 @@ void ReferenceCalcNonbondedForceKernel::initialize(const System& system, const N
ewaldDispersionAlpha = alpha;
useSwitchingFunction = false;
}
if (nonbondedMethod == NoCutoff || nonbondedMethod == CutoffNonPeriodic)
exceptionsArePeriodic = false;
else
exceptionsArePeriodic = force.getExceptionsUsePeriodicBoundaryConditions();
rfDielectric = force.getReactionFieldDielectric();
if (force.getUseDispersionCorrection())
dispersionCoefficient = NonbondedForceImpl::calcDispersionCorrection(system, force);
......@@ -1033,6 +1037,10 @@ double ReferenceCalcNonbondedForceKernel::execute(ContextImpl& context, bool inc
if (includeDirect) {
ReferenceBondForce refBondForce;
ReferenceLJCoulomb14 nonbonded14;
if (exceptionsArePeriodic) {
Vec3* boxVectors = extractBoxVectors(context);
nonbonded14.setPeriodic(boxVectors);
}
refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, includeEnergy ? &energy : NULL, nonbonded14);
if (periodic || ewald || pme) {
Vec3* boxVectors = extractBoxVectors(context);
......
platforms/reference/src/SimTKReference/ReferenceLJCoulomb14.cpp
View file @ ae07d525
......@@ -38,7 +38,7 @@ using namespace OpenMM;
--------------------------------------------------------------------------------------- */
ReferenceLJCoulomb14::ReferenceLJCoulomb14() {
ReferenceLJCoulomb14::ReferenceLJCoulomb14() : periodic(false) {
}
/**---------------------------------------------------------------------------------------
......@@ -50,6 +50,13 @@ ReferenceLJCoulomb14::ReferenceLJCoulomb14() {
ReferenceLJCoulomb14::~ReferenceLJCoulomb14() {
}
void ReferenceLJCoulomb14::setPeriodic(OpenMM::Vec3* vectors) {
periodic = true;
periodicBoxVectors[0] = vectors[0];
periodicBoxVectors[1] = vectors[1];
periodicBoxVectors[2] = vectors[2];
}
/**---------------------------------------------------------------------------------------
Calculate LJ 1-4 ixn
......@@ -68,33 +75,36 @@ ReferenceLJCoulomb14::~ReferenceLJCoulomb14() {
void ReferenceLJCoulomb14::calculateBondIxn(vector<int>& atomIndices, vector<Vec3>& atomCoordinates,
vector<double>& parameters, vector<Vec3>& forces,
double* totalEnergy, double* energyParamDerivs) {
double deltaR[2][ReferenceForce::LastDeltaRIndex];
double deltaR[2][ReferenceForce::LastDeltaRIndex];
// get deltaR, R2, and R between 2 atoms
// get deltaR, R2, and R between 2 atoms
int atomAIndex = atomIndices[0];
int atomBIndex = atomIndices[1];
ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], deltaR[0]);
int atomAIndex = atomIndices[0];
int atomBIndex = atomIndices[1];
if (periodic)
ReferenceForce::getDeltaRPeriodic(atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], periodicBoxVectors, deltaR[0]);
else
ReferenceForce::getDeltaR(atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], deltaR[0]);
double inverseR = 1.0/(deltaR[0][ReferenceForce::RIndex]);
double sig2 = inverseR*parameters[0];
sig2 *= sig2;
double sig6 = sig2*sig2*sig2;
double inverseR = 1.0/(deltaR[0][ReferenceForce::RIndex]);
double sig2 = inverseR*parameters[0];
sig2 *= sig2;
double sig6 = sig2*sig2*sig2;
double dEdR = parameters[1]*(12.0*sig6 - 6.0)*sig6;
dEdR += ONE_4PI_EPS0*parameters[2]*inverseR;
dEdR *= inverseR*inverseR;
double dEdR = parameters[1]*(12.0*sig6 - 6.0)*sig6;
dEdR += ONE_4PI_EPS0*parameters[2]*inverseR;
dEdR *= inverseR*inverseR;
// accumulate forces
// accumulate forces
for (int ii = 0; ii < 3; ii++) {
double force = dEdR*deltaR[0][ii];
forces[atomAIndex][ii] += force;
forces[atomBIndex][ii] -= force;
}
for (int ii = 0; ii < 3; ii++) {
double force = dEdR*deltaR[0][ii];
forces[atomAIndex][ii] += force;
forces[atomBIndex][ii] -= force;
}
// accumulate energies
// accumulate energies
if (totalEnergy != NULL)
*totalEnergy += parameters[1]*(sig6 - 1.0)*sig6 + (ONE_4PI_EPS0*parameters[2]*inverseR);
if (totalEnergy != NULL)
*totalEnergy += parameters[1]*(sig6 - 1.0)*sig6 + (ONE_4PI_EPS0*parameters[2]*inverseR);
}
serialization/src/NonbondedForceProxy.cpp
View file @ ae07d525
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2010-2018 Stanford University and the Authors. *
* Portions copyright (c) 2010-2020 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -42,7 +42,7 @@ NonbondedForceProxy::NonbondedForceProxy() : SerializationProxy("NonbondedForce"
}
void NonbondedForceProxy::serialize(const void* object, SerializationNode& node) const {
node.setIntProperty("version", 3);
node.setIntProperty("version", 4);
const NonbondedForce& force = *reinterpret_cast<const NonbondedForce*>(object);
node.setIntProperty("forceGroup", force.getForceGroup());
node.setIntProperty("method", (int) force.getNonbondedMethod());
......@@ -52,6 +52,7 @@ void NonbondedForceProxy::serialize(const void* object, SerializationNode& node)
node.setDoubleProperty("ewaldTolerance", force.getEwaldErrorTolerance());
node.setDoubleProperty("rfDielectric", force.getReactionFieldDielectric());
node.setIntProperty("dispersionCorrection", force.getUseDispersionCorrection());
node.setIntProperty("exceptionsUsePeriodic", force.getExceptionsUsePeriodicBoundaryConditions());
double alpha;
int nx, ny, nz;
force.getPMEParameters(alpha, nx, ny, nz);
......@@ -101,7 +102,7 @@ void NonbondedForceProxy::serialize(const void* object, SerializationNode& node)
void* NonbondedForceProxy::deserialize(const SerializationNode& node) const {
int version = node.getIntProperty("version");
if (version < 1 || version > 3)
if (version < 1 || version > 4)
throw OpenMMException("Unsupported version number");
NonbondedForce* force = new NonbondedForce();
try {
......@@ -137,6 +138,8 @@ void* NonbondedForceProxy::deserialize(const SerializationNode& node) const {
for (auto& offset : exceptionOffsets.getChildren())
force->addExceptionParameterOffset(offset.getStringProperty("parameter"), offset.getIntProperty("exception"), offset.getDoubleProperty("q"), offset.getDoubleProperty("sig"), offset.getDoubleProperty("eps"));
}
if (version >= 4)
force->setExceptionsUsePeriodicBoundaryConditions(node.getIntProperty("exceptionsUsePeriodic"));
const SerializationNode& particles = node.getChildNode("Particles");
for (auto& particle : particles.getChildren())
force->addParticle(particle.getDoubleProperty("q"), particle.getDoubleProperty("sig"), particle.getDoubleProperty("eps"));
......
serialization/tests/TestSerializeNonbondedForce.cpp
View file @ ae07d525
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2010-2018 Stanford University and the Authors. *
* Portions copyright (c) 2010-2020 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -50,6 +50,7 @@ void testSerialization() {
force.setEwaldErrorTolerance(1e-3);
force.setReactionFieldDielectric(50.0);
force.setUseDispersionCorrection(false);
force.setExceptionsUsePeriodicBoundaryConditions(true);
double alpha = 0.5;
int nx = 3, ny = 5, nz = 7;
force.setPMEParameters(alpha, nx, ny, nz);
......@@ -83,6 +84,7 @@ void testSerialization() {
ASSERT_EQUAL(force.getEwaldErrorTolerance(), force2.getEwaldErrorTolerance());
ASSERT_EQUAL(force.getReactionFieldDielectric(), force2.getReactionFieldDielectric());
ASSERT_EQUAL(force.getUseDispersionCorrection(), force2.getUseDispersionCorrection());
ASSERT_EQUAL(force.getExceptionsUsePeriodicBoundaryConditions(), force2.getExceptionsUsePeriodicBoundaryConditions());
ASSERT_EQUAL(force.getNumParticles(), force2.getNumParticles());
ASSERT_EQUAL(force.getNumExceptions(), force2.getNumExceptions());
ASSERT_EQUAL(force.getNumGlobalParameters(), force2.getNumGlobalParameters());
......
tests/TestNonbondedForce.h
View file @ ae07d525
......@@ -354,6 +354,44 @@ void testPeriodic() {
ASSERT_EQUAL_TOL(2*ONE_4PI_EPS0*(1.0)*(1.0+krf*1.0-crf), state.getPotentialEnergy(), TOL);
}
void testPeriodicExceptions() {
System system;
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->addParticle(1.0, 1, 0);
nonbonded->addParticle(1.0, 1, 0);
nonbonded->addException(0, 1, 1.0, 1.0, 0.0);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
const double cutoff = 2.0;
nonbonded->setCutoffDistance(cutoff);
system.setDefaultPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 4, 0), Vec3(0, 0, 4));
system.addForce(nonbonded);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(0, 0, 0);
positions[1] = Vec3(3, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
vector<Vec3> forces = state.getForces();
double force = ONE_4PI_EPS0/(3*3);
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], TOL);
ASSERT_EQUAL_TOL(ONE_4PI_EPS0/3, state.getPotentialEnergy(), TOL);
// Now make exceptions periodic and see if it changes correctly.
nonbonded->setExceptionsUsePeriodicBoundaryConditions(true);
context.reinitialize(true);
state = context.getState(State::Forces | State::Energy);
forces = state.getForces();
force = ONE_4PI_EPS0/(1*1);
ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[0], TOL);
ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[1], TOL);
ASSERT_EQUAL_TOL(ONE_4PI_EPS0/1, state.getPotentialEnergy(), TOL);
}
void testTriclinic() {
System system;
system.addParticle(1.0);
......@@ -809,6 +847,7 @@ int main(int argc, char* argv[]) {
testCutoff();
testCutoff14();
testPeriodic();
testPeriodicExceptions();
testTriclinic();
testLargeSystem();
testDispersionCorrection();
......
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