Throw exception if the box size is reduced to less than twice the cutoff

platforms/cuda/src/CudaKernels.cpp

@@ -50,6 +50,12 @@ void CudaCalcForcesAndEnergyKernel::beginComputation(ContextImpl& context, bool
     _gpuContext* gpu = data.gpu;
     if (data.nonbondedMethod != NO_CUTOFF && data.computeForceCount%100 == 0)
         gpuReorderAtoms(gpu);
+    if ((data.hasNonbonded && data.nonbondedMethod != NO_CUTOFF && data.nonbondedMethod != CUTOFF) ||
+        (data.hasCustomNonbonded && data.customNonbondedMethod != NO_CUTOFF && data.customNonbondedMethod != CUTOFF)) {
+        double minAllowedSize = 2*gpu->sim.nonbondedCutoff;
+        if (gpu->sim.periodicBoxSizeX < minAllowedSize || gpu->sim.periodicBoxSizeY < minAllowedSize || gpu->sim.periodicBoxSizeZ < minAllowedSize)
+            throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
+    }
     data.computeForceCount++;
     if (gpu->bIncludeGBSA || gpu->bIncludeGBVI)
         kClearBornSumAndForces(gpu);

platforms/cuda/tests/TestCudaMonteCarloBarostat.cpp

@@ -54,18 +54,35 @@ void testChangingBoxSize() {
     System system;
     system.setDefaultPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 5, 0), Vec3(0, 0, 6));
     system.addParticle(1.0);
+    NonbondedForce* nb = new NonbondedForce();
+    nb->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    nb->setCutoffDistance(2.0);
+    nb->addParticle(1, 0.5, 0.5);
+    system.addForce(nb);
     LangevinIntegrator integrator(300.0, 1.0, 0.01);
     Context context(system, integrator, platform);
     Vec3 x, y, z;
-    context.getState(0).getPeriodicBoxVectors(x, y, z);
+    context.getState(State::Forces).getPeriodicBoxVectors(x, y, z);
     ASSERT_EQUAL_VEC(Vec3(4, 0, 0), x, 0);
     ASSERT_EQUAL_VEC(Vec3(0, 5, 0), y, 0);
     ASSERT_EQUAL_VEC(Vec3(0, 0, 6), z, 0);
     context.setPeriodicBoxVectors(Vec3(7, 0, 0), Vec3(0, 8, 0), Vec3(0, 0, 9));
-    context.getState(0).getPeriodicBoxVectors(x, y, z);
+    context.getState(State::Forces).getPeriodicBoxVectors(x, y, z);
     ASSERT_EQUAL_VEC(Vec3(7, 0, 0), x, 0);
     ASSERT_EQUAL_VEC(Vec3(0, 8, 0), y, 0);
     ASSERT_EQUAL_VEC(Vec3(0, 0, 9), z, 0);
+    
+    // Shrinking the box too small should produce an exception.
+    
+    context.setPeriodicBoxVectors(Vec3(7, 0, 0), Vec3(0, 3.9, 0), Vec3(0, 0, 9));
+    bool ok = true;
+    try {
+        context.getState(State::Forces).getPeriodicBoxVectors(x, y, z);
+        ok = false;
+    }
+    catch (exception& ex) {
+    }
+    ASSERT(ok);
 }
 
 void testIdealGas() {

platforms/opencl/src/OpenCLNonbondedUtilities.cpp

@@ -24,6 +24,7 @@
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
  * -------------------------------------------------------------------------- */
 
+#include "openmm/OpenMMException.h"
 #include "OpenCLNonbondedUtilities.h"
 #include "OpenCLArray.h"
 #include "OpenCLKernelSources.h"
@@ -310,6 +311,12 @@ int OpenCLNonbondedUtilities::findExclusionIndex(int x, int y, const vector<cl_u
 void OpenCLNonbondedUtilities::prepareInteractions() {
     if (!useCutoff)
         return;
+    if (usePeriodic) {
+        mm_float4 box = context.getPeriodicBoxSize();
+        double minAllowedSize = 2*cutoff;
+        if (box.x < minAllowedSize || box.y < minAllowedSize || box.z < minAllowedSize)
+            throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
+    }
 
     // Compute the neighbor list.
 

platforms/opencl/tests/TestOpenCLMonteCarloBarostat.cpp

@@ -54,18 +54,35 @@ void testChangingBoxSize() {
     System system;
     system.setDefaultPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 5, 0), Vec3(0, 0, 6));
     system.addParticle(1.0);
+    NonbondedForce* nb = new NonbondedForce();
+    nb->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    nb->setCutoffDistance(2.0);
+    nb->addParticle(1, 0.5, 0.5);
+    system.addForce(nb);
     LangevinIntegrator integrator(300.0, 1.0, 0.01);
     Context context(system, integrator, platform);
     Vec3 x, y, z;
-    context.getState(0).getPeriodicBoxVectors(x, y, z);
+    context.getState(State::Forces).getPeriodicBoxVectors(x, y, z);
     ASSERT_EQUAL_VEC(Vec3(4, 0, 0), x, 0);
     ASSERT_EQUAL_VEC(Vec3(0, 5, 0), y, 0);
     ASSERT_EQUAL_VEC(Vec3(0, 0, 6), z, 0);
     context.setPeriodicBoxVectors(Vec3(7, 0, 0), Vec3(0, 8, 0), Vec3(0, 0, 9));
-    context.getState(0).getPeriodicBoxVectors(x, y, z);
+    context.getState(State::Forces).getPeriodicBoxVectors(x, y, z);
     ASSERT_EQUAL_VEC(Vec3(7, 0, 0), x, 0);
     ASSERT_EQUAL_VEC(Vec3(0, 8, 0), y, 0);
     ASSERT_EQUAL_VEC(Vec3(0, 0, 9), z, 0);
+    
+    // Shrinking the box too small should produce an exception.
+    
+    context.setPeriodicBoxVectors(Vec3(7, 0, 0), Vec3(0, 3.9, 0), Vec3(0, 0, 9));
+    bool ok = true;
+    try {
+        context.getState(State::Forces).getPeriodicBoxVectors(x, y, z);
+        ok = false;
+    }
+    catch (exception& ex) {
+    }
+    ASSERT(ok);
 }
 
 void testIdealGas() {

platforms/reference/src/ReferenceKernels.cpp

@@ -688,8 +688,13 @@ double ReferenceCalcNonbondedForceKernel::execute(ContextImpl& context, bool inc
         computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, extractBoxSize(context), periodic || ewald || pme, nonbondedCutoff, 0.0);
         clj.setUseCutoff(nonbondedCutoff, *neighborList, rfDielectric);
     }
-    if (periodic || ewald || pme)
-        clj.setPeriodic(extractBoxSize(context));
+    if (periodic || ewald || pme) {
+        RealVec& box = extractBoxSize(context);
+        double minAllowedSize = 2*nonbondedCutoff;
+        if (box[0] < minAllowedSize || box[1] < minAllowedSize || box[2] < minAllowedSize)
+            throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
+        clj.setPeriodic(box);
+    }
     if (ewald)
         clj.setUseEwald(ewaldAlpha, kmax[0], kmax[1], kmax[2]);
     if (pme)
@@ -821,8 +826,13 @@ double ReferenceCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bo
         computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, extractBoxSize(context), periodic, nonbondedCutoff, 0.0);
         ixn.setUseCutoff(nonbondedCutoff, *neighborList);
     }
-    if (periodic)
-        ixn.setPeriodic(extractBoxSize(context));
+    if (periodic) {
+        RealVec& box = extractBoxSize(context);
+        double minAllowedSize = 2*nonbondedCutoff;
+        if (box[0] < minAllowedSize || box[1] < minAllowedSize || box[2] < minAllowedSize)
+            throw OpenMMException("The periodic box size has decreased to less than twice the nonbonded cutoff.");
+        ixn.setPeriodic(box);
+    }
     map<string, double> globalParameters;
     for (int i = 0; i < (int) globalParameterNames.size(); i++)
         globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);

platforms/reference/tests/TestReferenceMonteCarloBarostat.cpp

@@ -54,18 +54,35 @@ void testChangingBoxSize() {
     System system;
     system.setDefaultPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 5, 0), Vec3(0, 0, 6));
     system.addParticle(1.0);
+    NonbondedForce* nb = new NonbondedForce();
+    nb->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    nb->setCutoffDistance(2.0);
+    nb->addParticle(1, 0.5, 0.5);
+    system.addForce(nb);
     LangevinIntegrator integrator(300.0, 1.0, 0.01);
     Context context(system, integrator, platform);
     Vec3 x, y, z;
-    context.getState(0).getPeriodicBoxVectors(x, y, z);
+    context.getState(State::Forces).getPeriodicBoxVectors(x, y, z);
     ASSERT_EQUAL_VEC(Vec3(4, 0, 0), x, 0);
     ASSERT_EQUAL_VEC(Vec3(0, 5, 0), y, 0);
     ASSERT_EQUAL_VEC(Vec3(0, 0, 6), z, 0);
     context.setPeriodicBoxVectors(Vec3(7, 0, 0), Vec3(0, 8, 0), Vec3(0, 0, 9));
-    context.getState(0).getPeriodicBoxVectors(x, y, z);
+    context.getState(State::Forces).getPeriodicBoxVectors(x, y, z);
     ASSERT_EQUAL_VEC(Vec3(7, 0, 0), x, 0);
     ASSERT_EQUAL_VEC(Vec3(0, 8, 0), y, 0);
     ASSERT_EQUAL_VEC(Vec3(0, 0, 9), z, 0);
+    
+    // Shrinking the box too small should produce an exception.
+    
+    context.setPeriodicBoxVectors(Vec3(7, 0, 0), Vec3(0, 3.9, 0), Vec3(0, 0, 9));
+    bool ok = true;
+    try {
+        context.getState(State::Forces).getPeriodicBoxVectors(x, y, z);
+        ok = false;
+    }
+    catch (exception& ex) {
+    }
+    ASSERT(ok);
 }
 
 void testIdealGas() {