Merge pull request #2570 from peastman/periodicexceptions

NonbondedForce can optionally apply periodic boundary conditions to exceptions

openmmapi/include/openmm/NonbondedForce.h

@@ -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

@@ -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

@@ -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

@@ -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

@@ -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

 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

@@ -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

 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

@@ -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

@@ -32,7 +32,7 @@ namespace OpenMM {
 
 class OPENMM_EXPORT ReferenceLJCoulomb14 : public ReferenceBondIxn {
 
-   public:
+public:
 
     /**---------------------------------------------------------------------------------------
 
@@ -52,11 +52,21 @@ class OPENMM_EXPORT ReferenceLJCoulomb14 : public ReferenceBondIxn {
 
      /**---------------------------------------------------------------------------------------
 
-         Calculate Ryckaert-Bellemans bond ixn
+       Set the force to use periodic boundary conditions.
 
-         @param atomIndices      atom indices of 4 atoms in bond
+       @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       six RB parameters
+       @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
 
@@ -66,6 +76,9 @@ class OPENMM_EXPORT ReferenceLJCoulomb14 : public ReferenceBondIxn {
                           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

@@ -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

@@ -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
@@ -74,6 +81,9 @@ void ReferenceLJCoulomb14::calculateBondIxn(vector<int>& atomIndices, vector<Vec
 
     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]);

serialization/src/NonbondedForceProxy.cpp

@@ -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

@@ -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

@@ -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();