Fixed error using periodic exceptions with PME (#2935)

platforms/cpu/include/CpuNonbondedForce.h

 
-/* Portions copyright (c) 2006-2018 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2020 Stanford University and Simbios.
  * Contributors: Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -125,6 +125,14 @@ class CpuNonbondedForce {
 
       void setUseLJPME(float alpha, int meshSize[3]);
 
+      /**---------------------------------------------------------------------------------------
+
+         Set whether exceptions use periodic boundary conditions.
+
+         --------------------------------------------------------------------------------------- */
+
+      void setPeriodicExceptions(bool periodic);
+
       /**---------------------------------------------------------------------------------------
       
          Calculate Ewald ixn
@@ -172,7 +180,7 @@ class CpuNonbondedForce {
 protected:
         bool cutoff;
         bool useSwitch;
-        bool periodic;
+        bool periodic, periodicExceptions;
         bool triclinic;
         bool ewald;
         bool ljpme, pme;

platforms/cpu/src/CpuKernels.cpp

@@ -666,6 +666,7 @@ double CpuCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeFo
         if (boxVectors[0][0] < minAllowedSize || boxVectors[1][1] < minAllowedSize || boxVectors[2][2] < minAllowedSize)
             throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
         nonbonded->setPeriodic(boxVectors);
+        nonbonded->setPeriodicExceptions(exceptionsArePeriodic);
     }
     if (ewald)
         nonbonded->setUseEwald(ewaldAlpha, kmax[0], kmax[1], kmax[2]);

platforms/cpu/src/CpuNonbondedForce.cpp

 
-/* Portions copyright (c) 2006-2018 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2020 Stanford University and Simbios.
  * Contributors: Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -47,7 +47,7 @@ const int CpuNonbondedForce::NUM_TABLE_POINTS = 2048;
 
    --------------------------------------------------------------------------------------- */
 
-CpuNonbondedForce::CpuNonbondedForce() : cutoff(false), useSwitch(false), periodic(false), ewald(false), pme(false), ljpme(false), tableIsValid(false), expTableIsValid(false),
+CpuNonbondedForce::CpuNonbondedForce() : cutoff(false), useSwitch(false), periodic(false), periodicExceptions(false), ewald(false), pme(false), ljpme(false), tableIsValid(false), expTableIsValid(false),
     cutoffDistance(0.0f), alphaDispersionEwald(0.0f), alphaEwald(0.0f) {
 }
 
@@ -202,6 +202,9 @@ void CpuNonbondedForce::setUseLJPME(float alpha, int meshSize[3]) {
     }
 }
 
+void CpuNonbondedForce::setPeriodicExceptions(bool periodic) {
+    periodicExceptions = periodic;
+}
 
 void CpuNonbondedForce::tabulateEwaldScaleFactor() {
     if (tableIsValid)
@@ -450,7 +453,7 @@ void CpuNonbondedForce::threadComputeDirect(ThreadPool& threads, int threadIndex
                         fvec4 deltaR;
                         fvec4 posJ((float) atomCoordinates[j][0], (float) atomCoordinates[j][1], (float) atomCoordinates[j][2], 0.0f);
                         float r2;
-                        getDeltaR(posJ, posI, deltaR, r2, false, boxSize, invBoxSize);
+                        getDeltaR(posJ, posI, deltaR, r2, periodicExceptions, boxSize, invBoxSize);
                         float r = sqrtf(r2);
                         float alphaR = alphaEwald*r;
                         float erfAlphaR = erf(alphaR);

platforms/cuda/src/CudaKernels.cpp

@@ -1052,6 +1052,7 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
             replacements["EWALD_ALPHA"] = cu.doubleToString(alpha);
             replacements["TWO_OVER_SQRT_PI"] = cu.doubleToString(2.0/sqrt(M_PI));
             replacements["DO_LJPME"] = doLJPME ? "1" : "0";
+            replacements["USE_PERIODIC"] = force.getExceptionsUsePeriodicBoundaryConditions() ? "1" : "0";
             if (doLJPME)
                 replacements["EWALD_DISPERSION_ALPHA"] = cu.doubleToString(dispersionAlpha);
             cu.getBondedUtilities().addInteraction(atoms, cu.replaceStrings(CudaKernelSources::pmeExclusions, replacements), force.getForceGroup());

platforms/cuda/src/kernels/pmeExclusions.cu

 const float4 exclusionParams = PARAMS[index];
 real3 delta = make_real3(pos2.x-pos1.x, pos2.y-pos1.y, pos2.z-pos1.z);
+#if USE_PERIODIC
+    APPLY_PERIODIC_TO_DELTA(delta)
+#endif
 const real r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
 const real r = SQRT(r2);
 const real invR = RECIP(r);

platforms/opencl/src/OpenCLKernels.cpp

@@ -1002,6 +1002,7 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
             replacements["EWALD_ALPHA"] = cl.doubleToString(alpha);
             replacements["TWO_OVER_SQRT_PI"] = cl.doubleToString(2.0/sqrt(M_PI));
             replacements["DO_LJPME"] = doLJPME ? "1" : "0";
+            replacements["USE_PERIODIC"] = force.getExceptionsUsePeriodicBoundaryConditions() ? "1" : "0";
             if (doLJPME)
                 replacements["EWALD_DISPERSION_ALPHA"] = cl.doubleToString(dispersionAlpha);
             cl.getBondedUtilities().addInteraction(atoms, cl.replaceStrings(OpenCLKernelSources::pmeExclusions, replacements), force.getForceGroup());

platforms/opencl/src/kernels/pmeExclusions.cl

 const float4 exclusionParams = PARAMS[index];
 real3 delta = (real3) (pos2.x-pos1.x, pos2.y-pos1.y, pos2.z-pos1.z);
+#if USE_PERIODIC
+    APPLY_PERIODIC_TO_DELTA(delta)
+#endif
 const real r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
 const real r = SQRT(r2);
 const real invR = RECIP(r);

platforms/reference/include/ReferenceLJCoulombIxn.h

 
-/* Portions copyright (c) 2006-2018 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2020 Stanford University and Simbios.
  * Contributors: Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -36,7 +36,7 @@ class ReferenceLJCoulombIxn {
        
       bool cutoff;
       bool useSwitch;
-      bool periodic;
+      bool periodic, periodicExceptions;
       bool ewald;
       bool pme, ljpme;
       const OpenMM::NeighborList* neighborList;
@@ -159,6 +159,14 @@ class ReferenceLJCoulombIxn {
 
       void setUseLJPME(double dalpha, int dmeshSize[3]);
       
+      /**---------------------------------------------------------------------------------------
+
+         Set whether exceptions use periodic boundary conditions.
+
+         --------------------------------------------------------------------------------------- */
+
+      void setPeriodicExceptions(bool periodic);
+
       /**---------------------------------------------------------------------------------------
       
          Calculate LJ Coulomb pair ixn

platforms/reference/src/ReferenceKernels.cpp

@@ -984,6 +984,7 @@ double ReferenceCalcNonbondedForceKernel::execute(ContextImpl& context, bool inc
         if (boxVectors[0][0] < minAllowedSize || boxVectors[1][1] < minAllowedSize || boxVectors[2][2] < minAllowedSize)
             throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
         clj.setPeriodic(boxVectors);
+        clj.setPeriodicExceptions(exceptionsArePeriodic);
     }
     if (ewald)
         clj.setUseEwald(ewaldAlpha, kmax[0], kmax[1], kmax[2]);

platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.cpp

 
-/* Portions copyright (c) 2006-2018 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2020 Stanford University and Simbios.
  * Contributors: Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -48,7 +48,7 @@ using namespace OpenMM;
 
    --------------------------------------------------------------------------------------- */
 
-ReferenceLJCoulombIxn::ReferenceLJCoulombIxn() : cutoff(false), useSwitch(false), periodic(false), ewald(false), pme(false), ljpme(false) {
+ReferenceLJCoulombIxn::ReferenceLJCoulombIxn() : cutoff(false), useSwitch(false), periodic(false), periodicExceptions(false), ewald(false), pme(false), ljpme(false) {
 }
 
 /**---------------------------------------------------------------------------------------
@@ -167,6 +167,10 @@ void ReferenceLJCoulombIxn::setUseLJPME(double alpha, int meshSize[3]) {
     ljpme = true;
 }
 
+void ReferenceLJCoulombIxn::setPeriodicExceptions(bool periodic) {
+    periodicExceptions = periodic;
+}
+
 /**---------------------------------------------------------------------------------------
 
    Calculate Ewald ixn
@@ -466,6 +470,9 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<Vec3>& a
                 int jj = exclusion;
 
                 double deltaR[2][ReferenceForce::LastDeltaRIndex];
+                if (periodicExceptions)
+                    ReferenceForce::getDeltaRPeriodic(atomCoordinates[jj], atomCoordinates[ii], periodicBoxVectors, deltaR[0]);
+                else
                     ReferenceForce::getDeltaR(atomCoordinates[jj], atomCoordinates[ii], deltaR[0]);
                 double r         = deltaR[0][ReferenceForce::RIndex];
                 double inverseR  = 1.0/(deltaR[0][ReferenceForce::RIndex]);

tests/TestNonbondedForce.h

@@ -893,6 +893,46 @@ void testParameterOffsets() {
     ASSERT_EQUAL_TOL(energy, context.getState(State::Energy).getPotentialEnergy(), 1e-5);
 }
 
+void testEwaldExceptions() {
+    // Create a minimal system using LJPME.
+
+    System system;
+    for (int i = 0; i < 4; i++)
+        system.addParticle(1.0);
+    system.setDefaultPeriodicBoxVectors(Vec3(2, 0, 0), Vec3(0, 2, 0), Vec3(0, 0, 2));
+    NonbondedForce* force = new NonbondedForce();
+    system.addForce(force);
+    force->setNonbondedMethod(NonbondedForce::LJPME);
+    force->setCutoffDistance(1.0);
+    force->addParticle(1.0, 0.5, 1.0);
+    force->addParticle(1.0, 0.5, 1.0);
+    force->addParticle(-1.0, 0.5, 1.0);
+    force->addParticle(-1.0, 0.5, 1.0);
+    vector<Vec3> positions = {
+        Vec3(0, 0, 0),
+        Vec3(1.5, 0, 0),
+        Vec3(0, 0.5, 0.5),
+        Vec3(0.2, 1.3, 0)
+    };
+    VerletIntegrator integrator(0.001);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    
+    // Compute the energy.
+    
+    double e1 = context.getState(State::Energy).getPotentialEnergy();
+
+    // Add a periodic exception and see if the energy changes by the correct amount.
+
+    force->addException(0, 1, 0.2, 0.8, 2.0);
+    force->setExceptionsUsePeriodicBoundaryConditions(true);
+    context.reinitialize(true);
+    double e2 = context.getState(State::Energy).getPotentialEnergy();
+    double r = 0.5;
+    double expectedChange = ONE_4PI_EPS0*(0.2-1.0)/r + 4*2.0*(pow(0.8/r, 12)-pow(0.8/r, 6)) - 4*1.0*(pow(0.5/r, 12)-pow(0.5/r, 6));
+    ASSERT_EQUAL_TOL(expectedChange, e2-e1, 1e-5);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
@@ -913,6 +953,7 @@ int main(int argc, char* argv[]) {
         testSwitchingFunction(NonbondedForce::PME);
         testTwoForces();
         testParameterOffsets();
+        testEwaldExceptions();
         runPlatformTests();
     }
     catch(const exception& e) {