Periodic box dimensions can be changed in the middle of a simulation

platforms/opencl/src/OpenCLKernels.cpp

@@ -70,6 +70,9 @@ void OpenCLCalcForcesAndEnergyKernel::initialize(const System& system) {
 }
 
 void OpenCLCalcForcesAndEnergyKernel::beginForceComputation(ContextImpl& context) {
+    Vec3 boxVectors[3];
+    context.getOwner().getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
+    cl.setPeriodicBoxSize(boxVectors[0][0], boxVectors[1][1], boxVectors[2][2]);
     if (cl.getNonbondedUtilities().getUseCutoff() && cl.getComputeForceCount()%100 == 0)
         cl.reorderAtoms();
     cl.setComputeForceCount(cl.getComputeForceCount()+1);
@@ -83,6 +86,9 @@ void OpenCLCalcForcesAndEnergyKernel::finishForceComputation(ContextImpl& contex
 }
 
 void OpenCLCalcForcesAndEnergyKernel::beginEnergyComputation(ContextImpl& context) {
+    Vec3 boxVectors[3];
+    context.getOwner().getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
+    cl.setPeriodicBoxSize(boxVectors[0][0], boxVectors[1][1], boxVectors[2][2]);
     if (cl.getNonbondedUtilities().getUseCutoff() && cl.getComputeForceCount()%100 == 0)
         cl.reorderAtoms();
     cl.setComputeForceCount(cl.getComputeForceCount()+1);
@@ -117,7 +123,7 @@ void OpenCLUpdateStateDataKernel::getPositions(ContextImpl& context, std::vector
     OpenCLArray<cl_int>& order = cl.getAtomIndex();
     int numParticles = context.getSystem().getNumParticles();
     positions.resize(numParticles);
-    mm_float4 periodicBoxSize = cl.getNonbondedUtilities().getPeriodicBoxSize();
+    mm_float4 periodicBoxSize = cl.getPeriodicBoxSize();
     for (int i = 0; i < numParticles; ++i) {
         mm_float4 pos = posq[i];
         mm_int4 offset = cl.getPosCellOffsets()[i];
@@ -1076,8 +1082,6 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
     sigmaEpsilon->upload(sigmaEpsilonVector);
     bool useCutoff = (force.getNonbondedMethod() != NonbondedForce::NoCutoff);
     bool usePeriodic = (force.getNonbondedMethod() != NonbondedForce::NoCutoff && force.getNonbondedMethod() != NonbondedForce::CutoffNonPeriodic);
-    Vec3 boxVectors[3];
-    system.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
     map<string, string> defines;
     defines["HAS_COULOMB"] = (hasCoulomb ? "1" : "0");
     defines["HAS_LENNARD_JONES"] = (hasLJ ? "1" : "0");
@@ -1107,10 +1111,6 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
         replacements["KMAX_X"] = intToString(kmaxx);
         replacements["KMAX_Y"] = intToString(kmaxy);
         replacements["KMAX_Z"] = intToString(kmaxz);
-        replacements["RECIPROCAL_BOX_SIZE_X"] = doubleToString(2.0*M_PI/boxVectors[0][0]);
-        replacements["RECIPROCAL_BOX_SIZE_Y"] = doubleToString(2.0*M_PI/boxVectors[1][1]);
-        replacements["RECIPROCAL_BOX_SIZE_Z"] = doubleToString(2.0*M_PI/boxVectors[2][2]);
-        replacements["RECIPROCAL_COEFFICIENT"] = doubleToString(ONE_4PI_EPS0*4*M_PI/(boxVectors[0][0]*boxVectors[1][1]*boxVectors[2][2]));
         replacements["EXP_COEFFICIENT"] = doubleToString(-1.0/(4.0*alpha*alpha));
         cl::Program program = cl.createProgram(OpenCLKernelSources::ewald, replacements);
         ewaldSumsKernel = cl::Kernel(program, "calculateEwaldCosSinSums");
@@ -1282,14 +1282,6 @@ void OpenCLCalcNonbondedForceKernel::executeForces(ContextImpl& context) {
             ewaldForcesKernel.setArg<cl::Buffer>(2, cosSinSums->getDeviceBuffer());
         }
         if (pmeGrid != NULL) {
-            mm_float4 boxSize = cl.getNonbondedUtilities().getPeriodicBoxSize();
-            pmeDefines["PERIODIC_BOX_SIZE_X"] = doubleToString(boxSize.x);
-            pmeDefines["PERIODIC_BOX_SIZE_Y"] = doubleToString(boxSize.y);
-            pmeDefines["PERIODIC_BOX_SIZE_Z"] = doubleToString(boxSize.z);
-            pmeDefines["INV_PERIODIC_BOX_SIZE_X"] = doubleToString(1.0/boxSize.x);
-            pmeDefines["INV_PERIODIC_BOX_SIZE_Y"] = doubleToString(1.0/boxSize.y);
-            pmeDefines["INV_PERIODIC_BOX_SIZE_Z"] = doubleToString(1.0/boxSize.z);
-            pmeDefines["RECIP_SCALE_FACTOR"] = doubleToString(1.0/(M_PI*boxSize.x*boxSize.y*boxSize.z));
             cl::Program program = cl.createProgram(OpenCLKernelSources::pme, pmeDefines);
             pmeGridIndexKernel = cl::Kernel(program, "updateGridIndexAndFraction");
             pmeAtomRangeKernel = cl::Kernel(program, "findAtomRangeForGrid");
@@ -1326,18 +1318,32 @@ void OpenCLCalcNonbondedForceKernel::executeForces(ContextImpl& context) {
     if (exceptionIndices != NULL)
         cl.executeKernel(exceptionsKernel, exceptionIndices->getSize());
     if (cosSinSums != NULL) {
+        mm_float4 boxSize = cl.getPeriodicBoxSize();
+        mm_float4 recipBoxSize = mm_float4((float) (2*M_PI/boxSize.x), (float) (2*M_PI/boxSize.y), (float) (2*M_PI/boxSize.z), 0);
+        float recipCoefficient = ONE_4PI_EPS0*4*M_PI/(boxSize.x*boxSize.y*boxSize.z);
+        ewaldSumsKernel.setArg<mm_float4>(3, recipBoxSize);
+        ewaldSumsKernel.setArg<cl_float>(4, recipCoefficient);
         cl.executeKernel(ewaldSumsKernel, cosSinSums->getSize());
+        ewaldForcesKernel.setArg<mm_float4>(3, recipBoxSize);
+        ewaldForcesKernel.setArg<cl_float>(4, recipCoefficient);
         cl.executeKernel(ewaldForcesKernel, cl.getNumAtoms());
     }
     if (pmeGrid != NULL) {
+        mm_float4 boxSize = cl.getPeriodicBoxSize();
+        mm_float4 invBoxSize = cl.getInvPeriodicBoxSize();
+        pmeGridIndexKernel.setArg<mm_float4>(2, invBoxSize);
         cl.executeKernel(pmeGridIndexKernel, cl.getNumAtoms());
         sort->sort(*pmeAtomGridIndex);
         cl.executeKernel(pmeAtomRangeKernel, cl.getNumAtoms());
+        pmeUpdateBsplinesKernel.setArg<mm_float4>(5, invBoxSize);
         cl.executeKernel(pmeUpdateBsplinesKernel, cl.getNumAtoms());
         cl.executeKernel(pmeSpreadChargeKernel, cl.getNumAtoms());
         fft->execFFT(*pmeGrid, true);
+        pmeConvolutionKernel.setArg<mm_float4>(5, invBoxSize);
+        pmeConvolutionKernel.setArg<cl_float>(6, (float) (1.0/(M_PI*boxSize.x*boxSize.y*boxSize.z)));
         cl.executeKernel(pmeConvolutionKernel, cl.getNumAtoms());
         fft->execFFT(*pmeGrid, false);
+        pmeInterpolateForceKernel.setArg<mm_float4>(5, invBoxSize);
         cl.executeKernel(pmeInterpolateForceKernel, cl.getNumAtoms());
     }
 }
@@ -1584,12 +1590,6 @@ void OpenCLCalcGBSAOBCForceKernel::executeForces(ContextImpl& context) {
             defines["USE_CUTOFF"] = "1";
         if (nb.getUsePeriodic())
             defines["USE_PERIODIC"] = "1";
-        defines["PERIODIC_BOX_SIZE_X"] = doubleToString(nb.getPeriodicBoxSize().x);
-        defines["PERIODIC_BOX_SIZE_Y"] = doubleToString(nb.getPeriodicBoxSize().y);
-        defines["PERIODIC_BOX_SIZE_Z"] = doubleToString(nb.getPeriodicBoxSize().z);
-        defines["INV_PERIODIC_BOX_SIZE_X"] = doubleToString(1.0/nb.getPeriodicBoxSize().x);
-        defines["INV_PERIODIC_BOX_SIZE_Y"] = doubleToString(1.0/nb.getPeriodicBoxSize().y);
-        defines["INV_PERIODIC_BOX_SIZE_Z"] = doubleToString(1.0/nb.getPeriodicBoxSize().z);
         defines["CUTOFF_SQUARED"] = doubleToString(nb.getCutoffDistance()*nb.getCutoffDistance());
         defines["PREFACTOR"] = doubleToString(prefactor);
         defines["NUM_ATOMS"] = intToString(cl.getNumAtoms());
@@ -1652,6 +1652,12 @@ void OpenCLCalcGBSAOBCForceKernel::executeForces(ContextImpl& context) {
     }
     cl.clearBuffer(*bornSum);
     cl.clearBuffer(*bornForce);
+    if (nb.getUseCutoff()) {
+        computeBornSumKernel.setArg<mm_float4>(8, cl.getPeriodicBoxSize());
+        computeBornSumKernel.setArg<mm_float4>(9, cl.getInvPeriodicBoxSize());
+        force1Kernel.setArg<mm_float4>(10, cl.getPeriodicBoxSize());
+        force1Kernel.setArg<mm_float4>(11, cl.getInvPeriodicBoxSize());
+    }
     cl.executeKernel(computeBornSumKernel, nb.getTiles().getSize()*OpenCLContext::TileSize);
     cl.executeKernel(reduceBornSumKernel, cl.getPaddedNumAtoms());
     cl.executeKernel(force1Kernel, nb.getTiles().getSize()*OpenCLContext::TileSize);
@@ -1883,14 +1889,6 @@ void OpenCLCalcCustomGBForceKernel::initialize(const System& system, const Custo
             defines["USE_PERIODIC"] = "1";
         if (useExclusionsForValue)
             defines["USE_EXCLUSIONS"] = "1";
-        Vec3 boxVectors[3];
-        system.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
-        defines["PERIODIC_BOX_SIZE_X"] = doubleToString(boxVectors[0][0]);
-        defines["PERIODIC_BOX_SIZE_Y"] = doubleToString(boxVectors[1][1]);
-        defines["PERIODIC_BOX_SIZE_Z"] = doubleToString(boxVectors[2][2]);
-        defines["INV_PERIODIC_BOX_SIZE_X"] = doubleToString(1.0/boxVectors[0][0]);
-        defines["INV_PERIODIC_BOX_SIZE_Y"] = doubleToString(1.0/boxVectors[1][1]);
-        defines["INV_PERIODIC_BOX_SIZE_Z"] = doubleToString(1.0/boxVectors[2][2]);
         defines["CUTOFF_SQUARED"] = doubleToString(force.getCutoffDistance()*force.getCutoffDistance());
         defines["NUM_ATOMS"] = intToString(cl.getNumAtoms());
         defines["PADDED_NUM_ATOMS"] = intToString(cl.getPaddedNumAtoms());
@@ -2037,14 +2035,6 @@ void OpenCLCalcCustomGBForceKernel::initialize(const System& system, const Custo
             defines["USE_PERIODIC"] = "1";
         if (anyExclusions)
             defines["USE_EXCLUSIONS"] = "1";
-        Vec3 boxVectors[3];
-        system.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
-        defines["PERIODIC_BOX_SIZE_X"] = doubleToString(boxVectors[0][0]);
-        defines["PERIODIC_BOX_SIZE_Y"] = doubleToString(boxVectors[1][1]);
-        defines["PERIODIC_BOX_SIZE_Z"] = doubleToString(boxVectors[2][2]);
-        defines["INV_PERIODIC_BOX_SIZE_X"] = doubleToString(1.0/boxVectors[0][0]);
-        defines["INV_PERIODIC_BOX_SIZE_Y"] = doubleToString(1.0/boxVectors[1][1]);
-        defines["INV_PERIODIC_BOX_SIZE_Z"] = doubleToString(1.0/boxVectors[2][2]);
         defines["CUTOFF_SQUARED"] = doubleToString(force.getCutoffDistance()*force.getCutoffDistance());
         defines["NUM_ATOMS"] = intToString(cl.getNumAtoms());
         defines["PADDED_NUM_ATOMS"] = intToString(cl.getPaddedNumAtoms());
@@ -2290,6 +2280,7 @@ void OpenCLCalcCustomGBForceKernel::executeForces(ContextImpl& context) {
             pairValueKernel.setArg<cl::Buffer>(index++, nb.getInteractingTiles().getDeviceBuffer());
             pairValueKernel.setArg<cl::Buffer>(index++, nb.getInteractionFlags().getDeviceBuffer());
             pairValueKernel.setArg<cl::Buffer>(index++, nb.getInteractionCount().getDeviceBuffer());
+            index += 2; // Periodic box size arguments are set when the kernel is executed.
         }
         else {
             pairValueKernel.setArg<cl::Buffer>(index++, nb.getTiles().getDeviceBuffer());
@@ -2336,6 +2327,7 @@ void OpenCLCalcCustomGBForceKernel::executeForces(ContextImpl& context) {
             pairEnergyKernel.setArg<cl::Buffer>(index++, nb.getInteractingTiles().getDeviceBuffer());
             pairEnergyKernel.setArg<cl::Buffer>(index++, nb.getInteractionFlags().getDeviceBuffer());
             pairEnergyKernel.setArg<cl::Buffer>(index++, nb.getInteractionCount().getDeviceBuffer());
+            index += 2; // Periodic box size arguments are set when the kernel is executed.
         }
         else {
             pairEnergyKernel.setArg<cl::Buffer>(index++, nb.getTiles().getDeviceBuffer());
@@ -2400,6 +2392,16 @@ void OpenCLCalcCustomGBForceKernel::executeForces(ContextImpl& context) {
         const OpenCLNonbondedUtilities::ParameterInfo& buffer = energyDerivs->getBuffers()[i];
         cl.clearBuffer(buffer.getMemory(), buffer.getSize()*energyDerivs->getNumObjects()/sizeof(cl_float));
     }
+    if (nb.getUseCutoff()) {
+        pairValueKernel.setArg<mm_float4>(10, cl.getPeriodicBoxSize());
+        pairValueKernel.setArg<mm_float4>(11, cl.getInvPeriodicBoxSize());
+        pairEnergyKernel.setArg<mm_float4>(11, cl.getPeriodicBoxSize());
+        pairEnergyKernel.setArg<mm_float4>(12, cl.getInvPeriodicBoxSize());
+        if (separateChainRuleKernel) {
+            chainRuleKernel.setArg<mm_float4>(10, cl.getPeriodicBoxSize());
+            chainRuleKernel.setArg<mm_float4>(11, cl.getInvPeriodicBoxSize());
+        }
+    }
     cl.executeKernel(pairValueKernel, nb.getTiles().getSize()*OpenCLContext::TileSize);
     cl.executeKernel(perParticleValueKernel, cl.getPaddedNumAtoms());
     cl.executeKernel(pairEnergyKernel, nb.getTiles().getSize()*OpenCLContext::TileSize);
@@ -3011,6 +3013,7 @@ void OpenCLCalcCustomHbondForceKernel::executeForces(ContextImpl& context) {
         donorKernel.setArg<cl::Buffer>(index++, acceptors->getDeviceBuffer());
         donorKernel.setArg<cl::Buffer>(index++, donorBufferIndices->getDeviceBuffer());
         donorKernel.setArg(index++, 3*OpenCLContext::ThreadBlockSize*sizeof(mm_float4), NULL);
+        index += 2; // Periodic box size arguments are set when the kernel is executed.
         if (globals != NULL)
             donorKernel.setArg<cl::Buffer>(index++, globals->getDeviceBuffer());
         for (int i = 0; i < (int) donorParams->getBuffers().size(); i++) {
@@ -3035,6 +3038,7 @@ void OpenCLCalcCustomHbondForceKernel::executeForces(ContextImpl& context) {
         acceptorKernel.setArg<cl::Buffer>(index++, acceptors->getDeviceBuffer());
         acceptorKernel.setArg<cl::Buffer>(index++, acceptorBufferIndices->getDeviceBuffer());
         acceptorKernel.setArg(index++, 3*OpenCLContext::ThreadBlockSize*sizeof(mm_float4), NULL);
+        index += 2; // Periodic box size arguments are set when the kernel is executed.
         if (globals != NULL)
             acceptorKernel.setArg<cl::Buffer>(index++, globals->getDeviceBuffer());
         for (int i = 0; i < (int) donorParams->getBuffers().size(); i++) {
@@ -3051,7 +3055,11 @@ void OpenCLCalcCustomHbondForceKernel::executeForces(ContextImpl& context) {
             acceptorKernel.setArg<cl::Buffer>(index++, tabulatedFunctionParams->getDeviceBuffer());
         }
     }
+    donorKernel.setArg<mm_float4>(8, cl.getPeriodicBoxSize());
+    donorKernel.setArg<mm_float4>(9, cl.getInvPeriodicBoxSize());
     cl.executeKernel(donorKernel, std::max(numDonors, numAcceptors));
+    acceptorKernel.setArg<mm_float4>(8, cl.getPeriodicBoxSize());
+    acceptorKernel.setArg<mm_float4>(9, cl.getInvPeriodicBoxSize());
     cl.executeKernel(acceptorKernel, std::max(numDonors, numAcceptors));
 }
 

platforms/opencl/src/OpenCLNonbondedUtilities.cpp

@@ -207,12 +207,6 @@ void OpenCLNonbondedUtilities::initialize(const System& system) {
     exclusions->upload(exclusionVec);
     exclusionIndex->upload(exclusionIndexVec);
 
-    // Record the periodic box size.
-
-    Vec3 boxVectors[3];
-    system.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
-    periodicBoxSize = mm_float4((float) boxVectors[0][0], (float) boxVectors[1][1], (float) boxVectors[2][2], 0.0f);
-
     // Create data structures for the neighbor list.
 
     if (useCutoff) {
@@ -236,28 +230,25 @@ void OpenCLNonbondedUtilities::initialize(const System& system) {
         cl::Program interactingBlocksProgram = context.createProgram(OpenCLKernelSources::findInteractingBlocks, defines);
         findBlockBoundsKernel = cl::Kernel(interactingBlocksProgram, "findBlockBounds");
         findBlockBoundsKernel.setArg<cl_int>(0, context.getNumAtoms());
-        findBlockBoundsKernel.setArg<mm_float4>(1, periodicBoxSize);
-        findBlockBoundsKernel.setArg<cl::Buffer>(2, context.getPosq().getDeviceBuffer());
-        findBlockBoundsKernel.setArg<cl::Buffer>(3, blockCenter->getDeviceBuffer());
-        findBlockBoundsKernel.setArg<cl::Buffer>(4, blockBoundingBox->getDeviceBuffer());
+        findBlockBoundsKernel.setArg<cl::Buffer>(3, context.getPosq().getDeviceBuffer());
+        findBlockBoundsKernel.setArg<cl::Buffer>(4, blockCenter->getDeviceBuffer());
+        findBlockBoundsKernel.setArg<cl::Buffer>(5, blockBoundingBox->getDeviceBuffer());
         findInteractingBlocksKernel = cl::Kernel(interactingBlocksProgram, "findBlocksWithInteractions");
         findInteractingBlocksKernel.setArg<cl_int>(0, tiles->getSize());
         findInteractingBlocksKernel.setArg<cl_float>(1, (cl_float) (cutoff*cutoff));
-        findInteractingBlocksKernel.setArg<mm_float4>(2, periodicBoxSize);
-        findInteractingBlocksKernel.setArg<cl::Buffer>(3, tiles->getDeviceBuffer());
-        findInteractingBlocksKernel.setArg<cl::Buffer>(4, blockCenter->getDeviceBuffer());
-        findInteractingBlocksKernel.setArg<cl::Buffer>(5, blockBoundingBox->getDeviceBuffer());
-        findInteractingBlocksKernel.setArg<cl::Buffer>(6, interactionFlags->getDeviceBuffer());
+        findInteractingBlocksKernel.setArg<cl::Buffer>(4, tiles->getDeviceBuffer());
+        findInteractingBlocksKernel.setArg<cl::Buffer>(5, blockCenter->getDeviceBuffer());
+        findInteractingBlocksKernel.setArg<cl::Buffer>(6, blockBoundingBox->getDeviceBuffer());
+        findInteractingBlocksKernel.setArg<cl::Buffer>(7, interactionFlags->getDeviceBuffer());
         findInteractionsWithinBlocksKernel = cl::Kernel(interactingBlocksProgram, "findInteractionsWithinBlocks");
         findInteractionsWithinBlocksKernel.setArg<cl_float>(0, (cl_float) (cutoff*cutoff));
-        findInteractionsWithinBlocksKernel.setArg<mm_float4>(1, periodicBoxSize);
-        findInteractionsWithinBlocksKernel.setArg<cl::Buffer>(2, context.getPosq().getDeviceBuffer());
-        findInteractionsWithinBlocksKernel.setArg<cl::Buffer>(3, interactingTiles->getDeviceBuffer());
-        findInteractionsWithinBlocksKernel.setArg<cl::Buffer>(4, blockCenter->getDeviceBuffer());
-        findInteractionsWithinBlocksKernel.setArg<cl::Buffer>(5, blockBoundingBox->getDeviceBuffer());
-        findInteractionsWithinBlocksKernel.setArg<cl::Buffer>(6, interactionFlags->getDeviceBuffer());
-        findInteractionsWithinBlocksKernel.setArg<cl::Buffer>(7, interactionCount->getDeviceBuffer());
-        findInteractionsWithinBlocksKernel.setArg(8, OpenCLContext::ThreadBlockSize*sizeof(cl_uint), NULL);
+        findInteractionsWithinBlocksKernel.setArg<cl::Buffer>(3, context.getPosq().getDeviceBuffer());
+        findInteractionsWithinBlocksKernel.setArg<cl::Buffer>(4, interactingTiles->getDeviceBuffer());
+        findInteractionsWithinBlocksKernel.setArg<cl::Buffer>(5, blockCenter->getDeviceBuffer());
+        findInteractionsWithinBlocksKernel.setArg<cl::Buffer>(6, blockBoundingBox->getDeviceBuffer());
+        findInteractionsWithinBlocksKernel.setArg<cl::Buffer>(7, interactionFlags->getDeviceBuffer());
+        findInteractionsWithinBlocksKernel.setArg<cl::Buffer>(8, interactionCount->getDeviceBuffer());
+        findInteractionsWithinBlocksKernel.setArg(9, OpenCLContext::ThreadBlockSize*sizeof(cl_uint), NULL);
     }
 }
 
@@ -267,16 +258,28 @@ void OpenCLNonbondedUtilities::prepareInteractions() {
 
     // Compute the neighbor list.
 
+    findBlockBoundsKernel.setArg<mm_float4>(1, context.getPeriodicBoxSize());
+    findBlockBoundsKernel.setArg<mm_float4>(2, context.getInvPeriodicBoxSize());
     context.executeKernel(findBlockBoundsKernel, context.getNumAtoms());
+    findInteractingBlocksKernel.setArg<mm_float4>(2, context.getPeriodicBoxSize());
+    findInteractingBlocksKernel.setArg<mm_float4>(3, context.getInvPeriodicBoxSize());
     context.executeKernel(findInteractingBlocksKernel, context.getNumAtoms());
     compact->compactStream(*interactingTiles, *tiles, *interactionFlags, *interactionCount);
-    if (context.getSIMDWidth() == 32)
+    if (context.getSIMDWidth() == 32) {
+        findInteractionsWithinBlocksKernel.setArg<mm_float4>(1, context.getPeriodicBoxSize());
+        findInteractionsWithinBlocksKernel.setArg<mm_float4>(2, context.getInvPeriodicBoxSize());
         context.executeKernel(findInteractionsWithinBlocksKernel, context.getNumAtoms());
+    }
 }
 
 void OpenCLNonbondedUtilities::computeInteractions() {
-    if (tiles != NULL)
+    if (tiles != NULL) {
+        if (useCutoff) {
+            forceKernel.setArg<mm_float4>(10, context.getPeriodicBoxSize());
+            forceKernel.setArg<mm_float4>(11, context.getInvPeriodicBoxSize());
+        }
         context.executeKernel(forceKernel, tiles->getSize()*OpenCLContext::TileSize);
+    }
 }
 
 cl::Kernel OpenCLNonbondedUtilities::createInteractionKernel(const string& source, const vector<ParameterInfo>& params, const vector<ParameterInfo>& arguments, bool useExclusions, bool isSymmetric) const {
@@ -382,12 +385,6 @@ cl::Kernel OpenCLNonbondedUtilities::createInteractionKernel(const string& sourc
         defines["USE_EXCLUSIONS"] = "1";
     if (isSymmetric)
         defines["USE_SYMMETRIC"] = "1";
-    defines["PERIODIC_BOX_SIZE_X"] = OpenCLExpressionUtilities::doubleToString(periodicBoxSize.x);
-    defines["PERIODIC_BOX_SIZE_Y"] = OpenCLExpressionUtilities::doubleToString(periodicBoxSize.y);
-    defines["PERIODIC_BOX_SIZE_Z"] = OpenCLExpressionUtilities::doubleToString(periodicBoxSize.z);
-    defines["INV_PERIODIC_BOX_SIZE_X"] = OpenCLExpressionUtilities::doubleToString(1.0/periodicBoxSize.x);
-    defines["INV_PERIODIC_BOX_SIZE_Y"] = OpenCLExpressionUtilities::doubleToString(1.0/periodicBoxSize.y);
-    defines["INV_PERIODIC_BOX_SIZE_Z"] = OpenCLExpressionUtilities::doubleToString(1.0/periodicBoxSize.z);
     defines["CUTOFF_SQUARED"] = OpenCLExpressionUtilities::doubleToString(cutoff*cutoff);
     defines["NUM_ATOMS"] = OpenCLExpressionUtilities::intToString(context.getNumAtoms());
     defines["PADDED_NUM_ATOMS"] = OpenCLExpressionUtilities::intToString(context.getPaddedNumAtoms());
@@ -409,6 +406,7 @@ cl::Kernel OpenCLNonbondedUtilities::createInteractionKernel(const string& sourc
         kernel.setArg<cl::Buffer>(index++, interactingTiles->getDeviceBuffer());
         kernel.setArg<cl::Buffer>(index++, interactionFlags->getDeviceBuffer());
         kernel.setArg<cl::Buffer>(index++, interactionCount->getDeviceBuffer());
+        index += 2; // The periodic box size arguments are set when the kernel is executed.
     }
     else {
         kernel.setArg<cl::Buffer>(index++, tiles->getDeviceBuffer());

platforms/opencl/src/OpenCLNonbondedUtilities.h

@@ -120,12 +120,6 @@ public:
     double getCutoffDistance() {
         return cutoff;
     }
-    /**
-     * Get the periodic box size.
-     */
-    mm_float4 getPeriodicBoxSize() {
-        return periodicBoxSize;
-    }
     /**
      * Prepare to compute interactions.  This updates the neighbor list.
      */
@@ -218,7 +212,6 @@ private:
     double cutoff;
     bool useCutoff, usePeriodic, forceBufferPerAtomBlock;
     int numForceBuffers;
-    mm_float4 periodicBoxSize;
 };
 
 /**

platforms/opencl/src/kernels/customGBEnergyN2_default.cl

@@ -11,7 +11,7 @@ void computeN2Energy(__global float4* forceBuffers, __global float* energyBuffer
 	__global float4* posq, __local float4* local_posq, __global unsigned int* exclusions, __global unsigned int* exclusionIndices,
         __local float4* tempForceBuffer, __global unsigned int* tiles,
 #ifdef USE_CUTOFF
-        __global unsigned int* interactionFlags, __global unsigned int* interactionCount
+        __global unsigned int* interactionFlags, __global unsigned int* interactionCount, float4 periodicBoxSize, float4 invPeriodicBoxSize
 #else
         unsigned int numTiles
 #endif
@@ -57,9 +57,9 @@ void computeN2Energy(__global float4* forceBuffers, __global float* energyBuffer
                 float4 posq2 = local_posq[atom2];
                 float4 delta = (float4) (posq2.xyz - posq1.xyz, 0.0f);
 #ifdef USE_PERIODIC
-                delta.x -= floor(delta.x*INV_PERIODIC_BOX_SIZE_X+0.5f)*PERIODIC_BOX_SIZE_X;
-                delta.y -= floor(delta.y*INV_PERIODIC_BOX_SIZE_Y+0.5f)*PERIODIC_BOX_SIZE_Y;
-                delta.z -= floor(delta.z*INV_PERIODIC_BOX_SIZE_Z+0.5f)*PERIODIC_BOX_SIZE_Z;
+                delta.x -= floor(delta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+                delta.y -= floor(delta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+                delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
 #endif
                 float r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
 #ifdef USE_CUTOFF
@@ -133,9 +133,9 @@ void computeN2Energy(__global float4* forceBuffers, __global float* energyBuffer
                 float4 posq2 = local_posq[atom2];
                 float4 delta = (float4) (posq2.xyz - posq1.xyz, 0.0f);
 #ifdef USE_PERIODIC
-                delta.x -= floor(delta.x*INV_PERIODIC_BOX_SIZE_X+0.5f)*PERIODIC_BOX_SIZE_X;
-                delta.y -= floor(delta.y*INV_PERIODIC_BOX_SIZE_Y+0.5f)*PERIODIC_BOX_SIZE_Y;
-                delta.z -= floor(delta.z*INV_PERIODIC_BOX_SIZE_Z+0.5f)*PERIODIC_BOX_SIZE_Z;
+                delta.x -= floor(delta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+                delta.y -= floor(delta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+                delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
 #endif
                 float r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
 #ifdef USE_CUTOFF

platforms/opencl/src/kernels/customGBEnergyN2_nvidia.cl

@@ -11,7 +11,7 @@ void computeN2Energy(__global float4* forceBuffers, __global float* energyBuffer
 	__global float4* posq, __local float4* local_posq, __global unsigned int* exclusions, __global unsigned int* exclusionIndices,
         __local float4* tempBuffer, __global unsigned int* tiles,
 #ifdef USE_CUTOFF
-        __global unsigned int* interactionFlags, __global unsigned int* interactionCount
+        __global unsigned int* interactionFlags, __global unsigned int* interactionCount, float4 periodicBoxSize, float4 invPeriodicBoxSize
 #else
         unsigned int numTiles
 #endif
@@ -56,9 +56,9 @@ void computeN2Energy(__global float4* forceBuffers, __global float* energyBuffer
                 float4 posq2 = local_posq[atom2];
                 float4 delta = (float4) (posq2.xyz - posq1.xyz, 0.0f);
 #ifdef USE_PERIODIC
-                delta.x -= floor(delta.x*INV_PERIODIC_BOX_SIZE_X+0.5f)*PERIODIC_BOX_SIZE_X;
-                delta.y -= floor(delta.y*INV_PERIODIC_BOX_SIZE_Y+0.5f)*PERIODIC_BOX_SIZE_Y;
-                delta.z -= floor(delta.z*INV_PERIODIC_BOX_SIZE_Z+0.5f)*PERIODIC_BOX_SIZE_Z;
+                delta.x -= floor(delta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+                delta.y -= floor(delta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+                delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
 #endif
                 float r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
 #ifdef USE_CUTOFF
@@ -129,9 +129,9 @@ void computeN2Energy(__global float4* forceBuffers, __global float* energyBuffer
                     float4 posq2 = local_posq[atom2];
                     float4 delta = (float4) (posq2.xyz - posq1.xyz, 0.0f);
 #ifdef USE_PERIODIC
-                    delta.x -= floor(delta.x*INV_PERIODIC_BOX_SIZE_X+0.5f)*PERIODIC_BOX_SIZE_X;
-                    delta.y -= floor(delta.y*INV_PERIODIC_BOX_SIZE_Y+0.5f)*PERIODIC_BOX_SIZE_Y;
-                    delta.z -= floor(delta.z*INV_PERIODIC_BOX_SIZE_Z+0.5f)*PERIODIC_BOX_SIZE_Z;
+                    delta.x -= floor(delta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+                    delta.y -= floor(delta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+                    delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
 #endif
                     float r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
 #ifdef USE_CUTOFF

platforms/opencl/src/kernels/customGBValueN2_default.cl

@@ -9,7 +9,7 @@ void computeN2Value(__global float4* posq, __local float4* local_posq, __global
         __global unsigned int* exclusionIndices, __global float* global_value, __local float* local_value,
         __local float* tempBuffer, __global unsigned int* tiles,
 #ifdef USE_CUTOFF
-        __global unsigned int* interactionFlags, __global unsigned int* interactionCount
+        __global unsigned int* interactionFlags, __global unsigned int* interactionCount, float4 periodicBoxSize, float4 invPeriodicBoxSize
 #else
         unsigned int numTiles
 #endif
@@ -54,9 +54,9 @@ void computeN2Value(__global float4* posq, __local float4* local_posq, __global
                 float4 posq2 = local_posq[atom2];
                 float4 delta = (float4) (posq2.xyz - posq1.xyz, 0.0f);
 #ifdef USE_PERIODIC
-                delta.x -= floor(delta.x*INV_PERIODIC_BOX_SIZE_X+0.5f)*PERIODIC_BOX_SIZE_X;
-                delta.y -= floor(delta.y*INV_PERIODIC_BOX_SIZE_Y+0.5f)*PERIODIC_BOX_SIZE_Y;
-                delta.z -= floor(delta.z*INV_PERIODIC_BOX_SIZE_Z+0.5f)*PERIODIC_BOX_SIZE_Z;
+                delta.x -= floor(delta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+                delta.y -= floor(delta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+                delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
 #endif
                 float r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
 #ifdef USE_CUTOFF
@@ -127,9 +127,9 @@ void computeN2Value(__global float4* posq, __local float4* local_posq, __global
                 float4 posq2 = local_posq[atom2];
                 float4 delta = (float4) (posq2.xyz - posq1.xyz, 0.0f);
 #ifdef USE_PERIODIC
-                delta.x -= floor(delta.x*INV_PERIODIC_BOX_SIZE_X+0.5f)*PERIODIC_BOX_SIZE_X;
-                delta.y -= floor(delta.y*INV_PERIODIC_BOX_SIZE_Y+0.5f)*PERIODIC_BOX_SIZE_Y;
-                delta.z -= floor(delta.z*INV_PERIODIC_BOX_SIZE_Z+0.5f)*PERIODIC_BOX_SIZE_Z;
+                delta.x -= floor(delta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+                delta.y -= floor(delta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+                delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
 #endif
                 float r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
 #ifdef USE_CUTOFF

platforms/opencl/src/kernels/customGBValueN2_nvidia.cl

@@ -9,7 +9,7 @@ void computeN2Value(__global float4* posq, __local float4* local_posq, __global
         __global unsigned int* exclusionIndices, __global float* global_value, __local float* local_value,
         __local float* tempBuffer, __global unsigned int* tiles,
 #ifdef USE_CUTOFF
-        __global unsigned int* interactionFlags, __global unsigned int* interactionCount
+        __global unsigned int* interactionFlags, __global unsigned int* interactionCount, float4 periodicBoxSize, float4 invPeriodicBoxSize
 #else
         unsigned int numTiles
 #endif
@@ -54,9 +54,9 @@ void computeN2Value(__global float4* posq, __local float4* local_posq, __global
                 float4 posq2 = local_posq[atom2];
                 float4 delta = (float4) (posq2.xyz - posq1.xyz, 0.0f);
 #ifdef USE_PERIODIC
-                delta.x -= floor(delta.x*INV_PERIODIC_BOX_SIZE_X+0.5f)*PERIODIC_BOX_SIZE_X;
-                delta.y -= floor(delta.y*INV_PERIODIC_BOX_SIZE_Y+0.5f)*PERIODIC_BOX_SIZE_Y;
-                delta.z -= floor(delta.z*INV_PERIODIC_BOX_SIZE_Z+0.5f)*PERIODIC_BOX_SIZE_Z;
+                delta.x -= floor(delta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+                delta.y -= floor(delta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+                delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
 #endif
                 float r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
 #ifdef USE_CUTOFF
@@ -115,9 +115,9 @@ void computeN2Value(__global float4* posq, __local float4* local_posq, __global
                             float4 posq2 = local_posq[atom2];
                             float4 delta = (float4) (posq2.xyz - posq1.xyz, 0.0f);
 #ifdef USE_PERIODIC
-                            delta.x -= floor(delta.x*INV_PERIODIC_BOX_SIZE_X+0.5f)*PERIODIC_BOX_SIZE_X;
-                            delta.y -= floor(delta.y*INV_PERIODIC_BOX_SIZE_Y+0.5f)*PERIODIC_BOX_SIZE_Y;
-                            delta.z -= floor(delta.z*INV_PERIODIC_BOX_SIZE_Z+0.5f)*PERIODIC_BOX_SIZE_Z;
+                            delta.x -= floor(delta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+                            delta.y -= floor(delta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+                            delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
 #endif
                             float r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
                             float tempValue1 = 0.0f;
@@ -170,9 +170,9 @@ void computeN2Value(__global float4* posq, __local float4* local_posq, __global
                     float4 posq2 = local_posq[atom2];
                     float4 delta = (float4) (posq2.xyz - posq1.xyz, 0.0f);
 #ifdef USE_PERIODIC
-                    delta.x -= floor(delta.x*INV_PERIODIC_BOX_SIZE_X+0.5f)*PERIODIC_BOX_SIZE_X;
-                    delta.y -= floor(delta.y*INV_PERIODIC_BOX_SIZE_Y+0.5f)*PERIODIC_BOX_SIZE_Y;
-                    delta.z -= floor(delta.z*INV_PERIODIC_BOX_SIZE_Z+0.5f)*PERIODIC_BOX_SIZE_Z;
+                    delta.x -= floor(delta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+                    delta.y -= floor(delta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+                    delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
 #endif
                     float r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
 #ifdef USE_CUTOFF

platforms/opencl/src/kernels/customHbondForce.cl

@@ -11,12 +11,12 @@ float4 delta(float4 vec1, float4 vec2) {
  * Compute the difference between two vectors, taking periodic boundary conditions into account
  * and setting the fourth component to the squared magnitude.
  */
-float4 deltaPeriodic(float4 vec1, float4 vec2) {
+float4 deltaPeriodic(float4 vec1, float4 vec2, float4 periodicBoxSize, float4 invPeriodicBoxSize) {
     float4 result = (float4) (vec1.x-vec2.x, vec1.y-vec2.y, vec1.z-vec2.z, 0.0f);
 #ifdef USE_PERIODIC
-    result.x -= floor(result.x*INV_PERIODIC_BOX_SIZE_X+0.5f)*PERIODIC_BOX_SIZE_X;
-    result.y -= floor(result.y*INV_PERIODIC_BOX_SIZE_Y+0.5f)*PERIODIC_BOX_SIZE_Y;
-    result.z -= floor(result.z*INV_PERIODIC_BOX_SIZE_Z+0.5f)*PERIODIC_BOX_SIZE_Z;
+    result.x -= floor(result.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+    result.y -= floor(result.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+    result.z -= floor(result.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
 #endif
     result.w = result.x*result.x + result.y*result.y + result.z*result.z;
     return result;
@@ -56,7 +56,7 @@ float4 computeCross(float4 vec1, float4 vec2) {
  * Compute forces on donors.
  */
 __kernel void computeDonorForces(__global float4* forceBuffers, __global float* energyBuffer, __global float4* posq, __global int4* exclusions,
-        __global int4* donorAtoms, __global int4* acceptorAtoms, __global int4* donorBufferIndices, __local float4* posBuffer
+        __global int4* donorAtoms, __global int4* acceptorAtoms, __global int4* donorBufferIndices, __local float4* posBuffer, float4 periodicBoxSize, float4 invPeriodicBoxSize
         PARAMETER_ARGUMENTS) {
     float energy = 0.0f;
     float4 f1 = (float4) 0;
@@ -102,7 +102,7 @@ __kernel void computeDonorForces(__global float4* forceBuffers, __global float*
                     float4 a1 = posBuffer[3*index];
                     float4 a2 = posBuffer[3*index+1];
                     float4 a3 = posBuffer[3*index+2];
-                    float4 deltaD1A1 = deltaPeriodic(d1, a1);
+                    float4 deltaD1A1 = deltaPeriodic(d1, a1, periodicBoxSize, invPeriodicBoxSize);
 #ifdef USE_CUTOFF
                     if (deltaD1A1.w < CUTOFF_SQUARED) {
 #endif
@@ -142,7 +142,7 @@ __kernel void computeDonorForces(__global float4* forceBuffers, __global float*
  * Compute forces on acceptors.
  */
 __kernel void computeAcceptorForces(__global float4* forceBuffers, __global float* energyBuffer, __global float4* posq, __global int4* exclusions,
-        __global int4* donorAtoms, __global int4* acceptorAtoms, __global int4* acceptorBufferIndices, __local float4* posBuffer
+        __global int4* donorAtoms, __global int4* acceptorAtoms, __global int4* acceptorBufferIndices, __local float4* posBuffer, float4 periodicBoxSize, float4 invPeriodicBoxSize
         PARAMETER_ARGUMENTS) {
     float4 f1 = (float4) 0;
     float4 f2 = (float4) 0;
@@ -187,7 +187,7 @@ __kernel void computeAcceptorForces(__global float4* forceBuffers, __global floa
                     float4 d1 = posBuffer[3*index];
                     float4 d2 = posBuffer[3*index+1];
                     float4 d3 = posBuffer[3*index+2];
-                    float4 deltaD1A1 = deltaPeriodic(d1, a1);
+                    float4 deltaD1A1 = deltaPeriodic(d1, a1, periodicBoxSize, invPeriodicBoxSize);
 #ifdef USE_CUTOFF
                     if (deltaD1A1.w < CUTOFF_SQUARED) {
 #endif

platforms/opencl/src/kernels/ewald.cl

@@ -7,7 +7,7 @@ float2 multofFloat2(float2 a, float2 b) {
  * Precompute the cosine and sine sums which appear in each force term.
  */
 
-__kernel void calculateEwaldCosSinSums(__global float* energyBuffer, __global float4* posq, __global float2* cosSinSum) {
+__kernel void calculateEwaldCosSinSums(__global float* energyBuffer, __global float4* posq, __global float2* cosSinSum, float4 reciprocalPeriodicBoxSize, float reciprocalCoefficient) {
     const unsigned int ksizex = 2*KMAX_X-1;
     const unsigned int ksizey = 2*KMAX_Y-1;
     const unsigned int ksizez = 2*KMAX_Z-1;
@@ -24,9 +24,9 @@ __kernel void calculateEwaldCosSinSums(__global float* energyBuffer, __global fl
         int ry = remainder/ksizez;
         int rz = remainder - ry*ksizez - KMAX_Z + 1;
         ry += -KMAX_Y + 1;
-        float kx = rx*RECIPROCAL_BOX_SIZE_X;
-        float ky = ry*RECIPROCAL_BOX_SIZE_Y;
-        float kz = rz*RECIPROCAL_BOX_SIZE_Z;
+        float kx = rx*reciprocalPeriodicBoxSize.x;
+        float ky = ry*reciprocalPeriodicBoxSize.y;
+        float kz = rz*reciprocalPeriodicBoxSize.z;
 
         // Compute the sum for this wave vector.
 
@@ -47,7 +47,7 @@ __kernel void calculateEwaldCosSinSums(__global float* energyBuffer, __global fl
 
         float k2 = kx*kx + ky*ky + kz*kz;
         float ak = exp(k2*EXP_COEFFICIENT) / k2;
-        energy += RECIPROCAL_COEFFICIENT*ak*(sum.x*sum.x + sum.y*sum.y);
+        energy += reciprocalCoefficient*ak*(sum.x*sum.x + sum.y*sum.y);
         index += get_global_size(0);
     }
     energyBuffer[get_global_id(0)] += energy;
@@ -58,7 +58,7 @@ __kernel void calculateEwaldCosSinSums(__global float* energyBuffer, __global fl
  * previous routine.
  */
 
-__kernel void calculateEwaldForces(__global float4* forceBuffers, __global float4* posq, __global float2* cosSinSum) {
+__kernel void calculateEwaldForces(__global float4* forceBuffers, __global float4* posq, __global float2* cosSinSum, float4 reciprocalPeriodicBoxSize, float reciprocalCoefficient) {
     unsigned int atom = get_global_id(0);
     while (atom < NUM_ATOMS) {
         float4 force = forceBuffers[atom];
@@ -69,15 +69,15 @@ __kernel void calculateEwaldForces(__global float4* forceBuffers, __global float
         int lowry = 0;
         int lowrz = 1;
         for (int rx = 0; rx < KMAX_X; rx++) {
-            float kx = rx*RECIPROCAL_BOX_SIZE_X;
+            float kx = rx*reciprocalPeriodicBoxSize.x;
             for (int ry = lowry; ry < KMAX_Y; ry++) {
-                float ky = ry*RECIPROCAL_BOX_SIZE_Y;
+                float ky = ry*reciprocalPeriodicBoxSize.y;
                 float phase = apos.x*kx;
                 float2 tab_xy = (float2) (cos(phase), sin(phase));
                 phase = apos.y*ky;
                 tab_xy = multofFloat2(tab_xy, (float2) (cos(phase), sin(phase)));
                 for (int rz = lowrz; rz < KMAX_Z; rz++) {
-                    float kz = rz*RECIPROCAL_BOX_SIZE_Z;
+                    float kz = rz*reciprocalPeriodicBoxSize.z;
 
                     // Compute the force contribution of this wave vector.
 
@@ -87,7 +87,7 @@ __kernel void calculateEwaldForces(__global float4* forceBuffers, __global float
                     phase = apos.z*kz;
                     float2 structureFactor = multofFloat2(tab_xy, (float2) (cos(phase), sin(phase)));
                     float2 sum = cosSinSum[index];
-                    float dEdR = 2*RECIPROCAL_COEFFICIENT*ak*apos.w*(sum.x*structureFactor.y - sum.y*structureFactor.x);
+                    float dEdR = 2*reciprocalCoefficient*ak*apos.w*(sum.x*structureFactor.y - sum.y*structureFactor.x);
                     force.x += dEdR*kx;
                     force.y += dEdR*ky;
                     force.z += dEdR*kz;

platforms/opencl/src/kernels/findInteractingBlocks.cl

@@ -3,15 +3,15 @@
 /**
  * Find a bounding box for the atoms in each block.
  */
-__kernel void findBlockBounds(int numAtoms, float4 periodicBoxSize, __global float4* posq, __global float4* blockCenter, __global float4* blockBoundingBox) {
+__kernel void findBlockBounds(int numAtoms, float4 periodicBoxSize, float4 invPeriodicBoxSize, __global float4* posq, __global float4* blockCenter, __global float4* blockBoundingBox) {
     int index = get_global_id(0);
     int base = index*TILE_SIZE;
     while (base < numAtoms) {
         float4 pos = posq[base];
 #ifdef USE_PERIODIC
-        pos.x -= floor(pos.x/periodicBoxSize.x)*periodicBoxSize.x;
-        pos.y -= floor(pos.y/periodicBoxSize.y)*periodicBoxSize.y;
-        pos.z -= floor(pos.z/periodicBoxSize.z)*periodicBoxSize.z;
+        pos.x -= floor(pos.x*invPeriodicBoxSize.x)*periodicBoxSize.x;
+        pos.y -= floor(pos.y*invPeriodicBoxSize.y)*periodicBoxSize.y;
+        pos.z -= floor(pos.z*invPeriodicBoxSize.z)*periodicBoxSize.z;
         float4 firstPoint = pos;
 #endif
         float4 minPos = pos;
@@ -20,9 +20,9 @@ __kernel void findBlockBounds(int numAtoms, float4 periodicBoxSize, __global flo
         for (int i = base+1; i < last; i++) {
             pos = posq[i];
 #ifdef USE_PERIODIC
-            pos.x -= floor((pos.x-firstPoint.x)/periodicBoxSize.x+0.5f)*periodicBoxSize.x;
-            pos.y -= floor((pos.y-firstPoint.y)/periodicBoxSize.y+0.5f)*periodicBoxSize.y;
-            pos.z -= floor((pos.z-firstPoint.z)/periodicBoxSize.z+0.5f)*periodicBoxSize.z;
+            pos.x -= floor((pos.x-firstPoint.x)*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+            pos.y -= floor((pos.y-firstPoint.y)*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+            pos.z -= floor((pos.z-firstPoint.z)*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
 #endif
             minPos = min(minPos, pos);
             maxPos = max(maxPos, pos);
@@ -38,7 +38,7 @@ __kernel void findBlockBounds(int numAtoms, float4 periodicBoxSize, __global flo
  * Compare the bounding boxes for each pair of blocks.  If they are sufficiently far apart,
  * mark them as non-interacting.
  */
-__kernel void findBlocksWithInteractions(int numTiles, float cutoffSquared, float4 periodicBoxSize, __global unsigned int* tiles, __global float4* blockCenter,
+__kernel void findBlocksWithInteractions(int numTiles, float cutoffSquared, float4 periodicBoxSize, float4 invPeriodicBoxSize, __global unsigned int* tiles, __global float4* blockCenter,
         __global float4* blockBoundingBox, __global unsigned int* interactionFlag) {
     int index = get_global_id(0);
     while (index < numTiles) {
@@ -52,9 +52,9 @@ __kernel void findBlocksWithInteractions(int numTiles, float cutoffSquared, floa
 
         float4 delta = blockCenter[x]-blockCenter[y];
 #ifdef USE_PERIODIC
-        delta.x -= floor(delta.x/periodicBoxSize.x+0.5f)*periodicBoxSize.x;
-        delta.y -= floor(delta.y/periodicBoxSize.y+0.5f)*periodicBoxSize.y;
-        delta.z -= floor(delta.z/periodicBoxSize.z+0.5f)*periodicBoxSize.z;
+        delta.x -= floor(delta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+        delta.y -= floor(delta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+        delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
 #endif
         float4 boxSizea = blockBoundingBox[x];
         float4 boxSizeb = blockBoundingBox[y];
@@ -70,7 +70,7 @@ __kernel void findBlocksWithInteractions(int numTiles, float cutoffSquared, floa
  * Compare each atom in one block to the bounding box of another block, and set
  * flags for which ones are interacting.
  */
-__kernel void findInteractionsWithinBlocks(float cutoffSquared, float4 periodicBoxSize, __global float4* posq, __global unsigned int* tiles, __global float4* blockCenter,
+__kernel void findInteractionsWithinBlocks(float cutoffSquared, float4 periodicBoxSize, float4 invPeriodicBoxSize, __global float4* posq, __global unsigned int* tiles, __global float4* blockCenter,
             __global float4* blockBoundingBox, __global unsigned int* interactionFlags, __global unsigned int* interactionCount, __local unsigned int* flags) {
     unsigned int totalWarps = get_global_size(0)/TILE_SIZE;
     unsigned int warp = get_global_id(0)/TILE_SIZE;
@@ -105,9 +105,9 @@ __kernel void findInteractionsWithinBlocks(float cutoffSquared, float4 periodicB
 
             float4 delta = apos-center;
 #ifdef USE_PERIODIC
-            delta.x -= floor(delta.x/periodicBoxSize.x+0.5f)*periodicBoxSize.x;
-            delta.y -= floor(delta.y/periodicBoxSize.y+0.5f)*periodicBoxSize.y;
-            delta.z -= floor(delta.z/periodicBoxSize.z+0.5f)*periodicBoxSize.z;
+            delta.x -= floor(delta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+            delta.y -= floor(delta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+            delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
 #endif
             delta = max((float4) 0.0f, fabs(delta)-boxSize);
             int thread = get_local_id(0);

platforms/opencl/src/kernels/gbsaObc_default.cl

@@ -17,7 +17,7 @@ typedef struct {
 __kernel __attribute__((reqd_work_group_size(WORK_GROUP_SIZE, 1, 1)))
 void computeBornSum(__global float* global_bornSum, __global float4* posq, __global float2* global_params, __local AtomData* localData, __local float* tempBuffer, __global unsigned int* tiles,
 #ifdef USE_CUTOFF
-        __global unsigned int* interactionFlags, __global unsigned int* interactionCount) {
+        __global unsigned int* interactionFlags, __global unsigned int* interactionCount, float4 periodicBoxSize, float4 invPeriodicBoxSize) {
 #else
         unsigned int numTiles) {
 #endif
@@ -56,9 +56,9 @@ void computeBornSum(__global float* global_bornSum, __global float4* posq, __glo
             for (unsigned int j = 0; j < TILE_SIZE/2; j++) {
                 float4 delta = (float4) (localData[baseLocalAtom+j].x-posq1.x, localData[baseLocalAtom+j].y-posq1.y, localData[baseLocalAtom+j].z-posq1.z, 0.0f);
 #ifdef USE_PERIODIC
-                delta.x -= floor(delta.x*INV_PERIODIC_BOX_SIZE_X+0.5f)*PERIODIC_BOX_SIZE_X;
-                delta.y -= floor(delta.y*INV_PERIODIC_BOX_SIZE_Y+0.5f)*PERIODIC_BOX_SIZE_Y;
-                delta.z -= floor(delta.z*INV_PERIODIC_BOX_SIZE_Z+0.5f)*PERIODIC_BOX_SIZE_Z;
+                delta.x -= floor(delta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+                delta.y -= floor(delta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+                delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
 #endif
                 float r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
 #ifdef USE_CUTOFF
@@ -124,9 +124,9 @@ void computeBornSum(__global float* global_bornSum, __global float4* posq, __glo
             for (unsigned int j = 0; j < TILE_SIZE/2; j++) {
                 float4 delta = (float4) (localData[baseLocalAtom+tj].x-posq1.x, localData[baseLocalAtom+tj].y-posq1.y, localData[baseLocalAtom+tj].z-posq1.z, 0.0f);
 #ifdef USE_PERIODIC
-                delta.x -= floor(delta.x*INV_PERIODIC_BOX_SIZE_X+0.5f)*PERIODIC_BOX_SIZE_X;
-                delta.y -= floor(delta.y*INV_PERIODIC_BOX_SIZE_Y+0.5f)*PERIODIC_BOX_SIZE_Y;
-                delta.z -= floor(delta.z*INV_PERIODIC_BOX_SIZE_Z+0.5f)*PERIODIC_BOX_SIZE_Z;
+                delta.x -= floor(delta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+                delta.y -= floor(delta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+                delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
 #endif
                 float r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
 #ifdef USE_CUTOFF
@@ -198,7 +198,7 @@ void computeGBSAForce1(__global float4* forceBuffers, __global float* energyBuff
         __global float4* posq, __global float* global_bornRadii,
         __global float* global_bornForce, __local AtomData* localData, __local float4* tempBuffer, __global unsigned int* tiles,
 #ifdef USE_CUTOFF
-        __global unsigned int* interactionFlags, __global unsigned int* interactionCount) {
+        __global unsigned int* interactionFlags, __global unsigned int* interactionCount, float4 periodicBoxSize, float4 invPeriodicBoxSize) {
 #else
         unsigned int numTiles) {
 #endif
@@ -238,9 +238,9 @@ void computeGBSAForce1(__global float4* forceBuffers, __global float* energyBuff
                     float4 posq2 = (float4) (localData[baseLocalAtom+j].x, localData[baseLocalAtom+j].y, localData[baseLocalAtom+j].z, localData[baseLocalAtom+j].q);
                     float4 delta = (float4) (posq2.xyz - posq1.xyz, 0.0f);
 #ifdef USE_PERIODIC
-                    delta.x -= floor(delta.x*INV_PERIODIC_BOX_SIZE_X+0.5f)*PERIODIC_BOX_SIZE_X;
-                    delta.y -= floor(delta.y*INV_PERIODIC_BOX_SIZE_Y+0.5f)*PERIODIC_BOX_SIZE_Y;
-                    delta.z -= floor(delta.z*INV_PERIODIC_BOX_SIZE_Z+0.5f)*PERIODIC_BOX_SIZE_Z;
+                    delta.x -= floor(delta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+                    delta.y -= floor(delta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+                    delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
 #endif
                     float r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
                     float invR = RSQRT(r2);
@@ -312,9 +312,9 @@ void computeGBSAForce1(__global float4* forceBuffers, __global float* energyBuff
                     float4 posq2 = (float4) (localData[baseLocalAtom+tj].x, localData[baseLocalAtom+tj].y, localData[baseLocalAtom+tj].z, localData[baseLocalAtom+tj].q);
                     float4 delta = (float4) (posq2.xyz - posq1.xyz, 0.0f);
 #ifdef USE_PERIODIC
-                    delta.x -= floor(delta.x*INV_PERIODIC_BOX_SIZE_X+0.5f)*PERIODIC_BOX_SIZE_X;
-                    delta.y -= floor(delta.y*INV_PERIODIC_BOX_SIZE_Y+0.5f)*PERIODIC_BOX_SIZE_Y;
-                    delta.z -= floor(delta.z*INV_PERIODIC_BOX_SIZE_Z+0.5f)*PERIODIC_BOX_SIZE_Z;
+                    delta.x -= floor(delta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+                    delta.y -= floor(delta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+                    delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
 #endif
                     float r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
                     float invR = RSQRT(r2);

platforms/opencl/src/kernels/gbsaObc_nvidia.cl

@@ -17,7 +17,7 @@ typedef struct {
 __kernel __attribute__((reqd_work_group_size(WORK_GROUP_SIZE, 1, 1)))
 void computeBornSum(__global float* global_bornSum, __global float4* posq, __global float2* global_params, __local AtomData* localData, __local float* tempBuffer, __global unsigned int* tiles,
 #ifdef USE_CUTOFF
-        __global unsigned int* interactionFlags, __global unsigned int* interactionCount) {
+        __global unsigned int* interactionFlags, __global unsigned int* interactionCount, float4 periodicBoxSize, float4 invPeriodicBoxSize) {
 #else
         unsigned int numTiles) {
 #endif
@@ -56,9 +56,9 @@ void computeBornSum(__global float* global_bornSum, __global float4* posq, __glo
             for (unsigned int j = 0; j < TILE_SIZE; j++) {
                 float4 delta = (float4) (localData[tbx+j].x-posq1.x, localData[tbx+j].y-posq1.y, localData[tbx+j].z-posq1.z, 0.0f);
 #ifdef USE_PERIODIC
-                delta.x -= floor(delta.x*INV_PERIODIC_BOX_SIZE_X+0.5f)*PERIODIC_BOX_SIZE_X;
-                delta.y -= floor(delta.y*INV_PERIODIC_BOX_SIZE_Y+0.5f)*PERIODIC_BOX_SIZE_Y;
-                delta.z -= floor(delta.z*INV_PERIODIC_BOX_SIZE_Z+0.5f)*PERIODIC_BOX_SIZE_Z;
+                delta.x -= floor(delta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+                delta.y -= floor(delta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+                delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
 #endif
                 float r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
 #ifdef USE_CUTOFF
@@ -120,9 +120,9 @@ void computeBornSum(__global float* global_bornSum, __global float4* posq, __glo
                         if ((flags&(1<<j)) != 0) {
                             float4 delta = (float4) (localData[tbx+j].x-posq1.x, localData[tbx+j].y-posq1.y, localData[tbx+j].z-posq1.z, 0.0f);
 #ifdef USE_PERIODIC
-                            delta.x -= floor(delta.x*INV_PERIODIC_BOX_SIZE_X+0.5f)*PERIODIC_BOX_SIZE_X;
-                            delta.y -= floor(delta.y*INV_PERIODIC_BOX_SIZE_Y+0.5f)*PERIODIC_BOX_SIZE_Y;
-                            delta.z -= floor(delta.z*INV_PERIODIC_BOX_SIZE_Z+0.5f)*PERIODIC_BOX_SIZE_Z;
+                            delta.x -= floor(delta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+                            delta.y -= floor(delta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+                            delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
 #endif
                             float r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
                             tempBuffer[get_local_id(0)] = 0.0f;
@@ -189,9 +189,9 @@ void computeBornSum(__global float* global_bornSum, __global float4* posq, __glo
                 for (unsigned int j = 0; j < TILE_SIZE; j++) {
                     float4 delta = (float4) (localData[tbx+tj].x-posq1.x, localData[tbx+tj].y-posq1.y, localData[tbx+tj].z-posq1.z, 0.0f);
 #ifdef USE_PERIODIC
-                    delta.x -= floor(delta.x*INV_PERIODIC_BOX_SIZE_X+0.5f)*PERIODIC_BOX_SIZE_X;
-                    delta.y -= floor(delta.y*INV_PERIODIC_BOX_SIZE_Y+0.5f)*PERIODIC_BOX_SIZE_Y;
-                    delta.z -= floor(delta.z*INV_PERIODIC_BOX_SIZE_Z+0.5f)*PERIODIC_BOX_SIZE_Z;
+                    delta.x -= floor(delta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+                    delta.y -= floor(delta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+                    delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
 #endif
                     float r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
 #ifdef USE_CUTOFF
@@ -258,7 +258,7 @@ void computeGBSAForce1(__global float4* forceBuffers, __global float* energyBuff
         __global float4* posq, __global float* global_bornRadii,
         __global float* global_bornForce, __local AtomData* localData, __local float4* tempBuffer, __global unsigned int* tiles,
 #ifdef USE_CUTOFF
-        __global unsigned int* interactionFlags, __global unsigned int* interactionCount) {
+        __global unsigned int* interactionFlags, __global unsigned int* interactionCount, float4 periodicBoxSize, float4 invPeriodicBoxSize) {
 #else
         unsigned int numTiles) {
 #endif
@@ -298,9 +298,9 @@ void computeGBSAForce1(__global float4* forceBuffers, __global float* energyBuff
                     float4 posq2 = (float4) (localData[tbx+j].x, localData[tbx+j].y, localData[tbx+j].z, localData[tbx+j].q);
                     float4 delta = (float4) (posq2.xyz - posq1.xyz, 0.0f);
 #ifdef USE_PERIODIC
-                    delta.x -= floor(delta.x*INV_PERIODIC_BOX_SIZE_X+0.5f)*PERIODIC_BOX_SIZE_X;
-                    delta.y -= floor(delta.y*INV_PERIODIC_BOX_SIZE_Y+0.5f)*PERIODIC_BOX_SIZE_Y;
-                    delta.z -= floor(delta.z*INV_PERIODIC_BOX_SIZE_Z+0.5f)*PERIODIC_BOX_SIZE_Z;
+                    delta.x -= floor(delta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+                    delta.y -= floor(delta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+                    delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
 #endif
                     float r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
                     float invR = RSQRT(r2);
@@ -367,9 +367,9 @@ void computeGBSAForce1(__global float4* forceBuffers, __global float* energyBuff
                             float4 posq2 = (float4) (localData[tbx+j].x, localData[tbx+j].y, localData[tbx+j].z, localData[tbx+j].q);
                             float4 delta = (float4) (posq2.xyz - posq1.xyz, 0.0f);
 #ifdef USE_PERIODIC
-                            delta.x -= floor(delta.x*INV_PERIODIC_BOX_SIZE_X+0.5f)*PERIODIC_BOX_SIZE_X;
-                            delta.y -= floor(delta.y*INV_PERIODIC_BOX_SIZE_Y+0.5f)*PERIODIC_BOX_SIZE_Y;
-                            delta.z -= floor(delta.z*INV_PERIODIC_BOX_SIZE_Z+0.5f)*PERIODIC_BOX_SIZE_Z;
+                            delta.x -= floor(delta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+                            delta.y -= floor(delta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+                            delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
 #endif
                             float r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
                             float invR = RSQRT(r2);
@@ -433,9 +433,9 @@ void computeGBSAForce1(__global float4* forceBuffers, __global float* energyBuff
                         float4 posq2 = (float4) (localData[tbx+tj].x, localData[tbx+tj].y, localData[tbx+tj].z, localData[tbx+tj].q);
                         float4 delta = (float4) (posq2.xyz - posq1.xyz, 0.0f);
 #ifdef USE_PERIODIC
-                        delta.x -= floor(delta.x*INV_PERIODIC_BOX_SIZE_X+0.5f)*PERIODIC_BOX_SIZE_X;
-                        delta.y -= floor(delta.y*INV_PERIODIC_BOX_SIZE_Y+0.5f)*PERIODIC_BOX_SIZE_Y;
-                        delta.z -= floor(delta.z*INV_PERIODIC_BOX_SIZE_Z+0.5f)*PERIODIC_BOX_SIZE_Z;
+                        delta.x -= floor(delta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+                        delta.y -= floor(delta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+                        delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
 #endif
                         float r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
                         float invR = RSQRT(r2);

platforms/opencl/src/kernels/nonbonded_default.cl

@@ -15,7 +15,7 @@ __kernel __attribute__((reqd_work_group_size(WORK_GROUP_SIZE, 1, 1)))
 void computeNonbonded(__global float4* forceBuffers, __global float* energyBuffer, __global float4* posq, __global unsigned int* exclusions,
         __global unsigned int* exclusionIndices, __local AtomData* localData, __local float4* tempBuffer, __global unsigned int* tiles,
 #ifdef USE_CUTOFF
-        __global unsigned int* interactionFlags, __global unsigned int* interactionCount
+        __global unsigned int* interactionFlags, __global unsigned int* interactionCount, float4 periodicBoxSize, float4 invPeriodicBoxSize
 #else
         unsigned int numTiles
 #endif
@@ -63,9 +63,9 @@ void computeNonbonded(__global float4* forceBuffers, __global float* energyBuffe
                 float4 posq2 = (float4) (localData[atom2].x, localData[atom2].y, localData[atom2].z, localData[atom2].q);
                 float4 delta = (float4) (posq2.xyz - posq1.xyz, 0.0f);
 #ifdef USE_PERIODIC
-                delta.x -= floor(delta.x*INV_PERIODIC_BOX_SIZE_X+0.5f)*PERIODIC_BOX_SIZE_X;
-                delta.y -= floor(delta.y*INV_PERIODIC_BOX_SIZE_Y+0.5f)*PERIODIC_BOX_SIZE_Y;
-                delta.z -= floor(delta.z*INV_PERIODIC_BOX_SIZE_Z+0.5f)*PERIODIC_BOX_SIZE_Z;
+                delta.x -= floor(delta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+                delta.y -= floor(delta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+                delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
 #endif
                 float r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
                 float invR = RSQRT(r2);
@@ -139,9 +139,9 @@ void computeNonbonded(__global float4* forceBuffers, __global float* energyBuffe
                 float4 posq2 = (float4) (localData[atom2].x, localData[atom2].y, localData[atom2].z, localData[atom2].q);
                 float4 delta = (float4) (posq2.xyz - posq1.xyz, 0.0f);
 #ifdef USE_PERIODIC
-                delta.x -= floor(delta.x*INV_PERIODIC_BOX_SIZE_X+0.5f)*PERIODIC_BOX_SIZE_X;
-                delta.y -= floor(delta.y*INV_PERIODIC_BOX_SIZE_Y+0.5f)*PERIODIC_BOX_SIZE_Y;
-                delta.z -= floor(delta.z*INV_PERIODIC_BOX_SIZE_Z+0.5f)*PERIODIC_BOX_SIZE_Z;
+                delta.x -= floor(delta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+                delta.y -= floor(delta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+                delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
 #endif
                 float r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
                 float invR = RSQRT(r2);

platforms/opencl/src/kernels/nonbonded_nvidia.cl

@@ -15,7 +15,7 @@ __kernel __attribute__((reqd_work_group_size(WORK_GROUP_SIZE, 1, 1)))
 void computeNonbonded(__global float4* forceBuffers, __global float* energyBuffer, __global float4* posq, __global unsigned int* exclusions,
         __global unsigned int* exclusionIndices, __local AtomData* localData, __local float4* tempBuffer, __global unsigned int* tiles,
 #ifdef USE_CUTOFF
-        __global unsigned int* interactionFlags, __global unsigned int* interactionCount
+        __global unsigned int* interactionFlags, __global unsigned int* interactionCount, float4 periodicBoxSize, float4 invPeriodicBoxSize
 #else
         unsigned int numTiles
 #endif
@@ -63,9 +63,9 @@ void computeNonbonded(__global float4* forceBuffers, __global float* energyBuffe
                 float4 posq2 = (float4) (localData[atom2].x, localData[atom2].y, localData[atom2].z, localData[atom2].q);
                 float4 delta = (float4) (posq2.xyz - posq1.xyz, 0.0f);
 #ifdef USE_PERIODIC
-                delta.x -= floor(delta.x*INV_PERIODIC_BOX_SIZE_X+0.5f)*PERIODIC_BOX_SIZE_X;
-                delta.y -= floor(delta.y*INV_PERIODIC_BOX_SIZE_Y+0.5f)*PERIODIC_BOX_SIZE_Y;
-                delta.z -= floor(delta.z*INV_PERIODIC_BOX_SIZE_Z+0.5f)*PERIODIC_BOX_SIZE_Z;
+                delta.x -= floor(delta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+                delta.y -= floor(delta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+                delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
 #endif
                 float r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
                 float r = sqrt(r2);
@@ -128,9 +128,9 @@ void computeNonbonded(__global float4* forceBuffers, __global float* energyBuffe
                             float4 posq2 = (float4) (localData[atom2].x, localData[atom2].y, localData[atom2].z, localData[atom2].q);
                             float4 delta = (float4) (posq2.xyz - posq1.xyz, 0.0f);
 #ifdef USE_PERIODIC
-                            delta.x -= floor(delta.x*INV_PERIODIC_BOX_SIZE_X+0.5f)*PERIODIC_BOX_SIZE_X;
-                            delta.y -= floor(delta.y*INV_PERIODIC_BOX_SIZE_Y+0.5f)*PERIODIC_BOX_SIZE_Y;
-                            delta.z -= floor(delta.z*INV_PERIODIC_BOX_SIZE_Z+0.5f)*PERIODIC_BOX_SIZE_Z;
+                            delta.x -= floor(delta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+                            delta.y -= floor(delta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+                            delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
 #endif
                             float r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
                             float invR = RSQRT(r2);
@@ -195,9 +195,9 @@ void computeNonbonded(__global float4* forceBuffers, __global float* energyBuffe
                     float4 posq2 = (float4) (localData[atom2].x, localData[atom2].y, localData[atom2].z, localData[atom2].q);
                     float4 delta = (float4) (posq2.xyz - posq1.xyz, 0.0f);
 #ifdef USE_PERIODIC
-                    delta.x -= floor(delta.x*INV_PERIODIC_BOX_SIZE_X+0.5f)*PERIODIC_BOX_SIZE_X;
-                    delta.y -= floor(delta.y*INV_PERIODIC_BOX_SIZE_Y+0.5f)*PERIODIC_BOX_SIZE_Y;
-                    delta.z -= floor(delta.z*INV_PERIODIC_BOX_SIZE_Z+0.5f)*PERIODIC_BOX_SIZE_Z;
+                    delta.x -= floor(delta.x*invPeriodicBoxSize.x+0.5f)*periodicBoxSize.x;
+                    delta.y -= floor(delta.y*invPeriodicBoxSize.y+0.5f)*periodicBoxSize.y;
+                    delta.z -= floor(delta.z*invPeriodicBoxSize.z+0.5f)*periodicBoxSize.z;
 #endif
                     float r2 = delta.x*delta.x + delta.y*delta.y + delta.z*delta.z;
                     float invR = RSQRT(r2);

platforms/opencl/src/kernels/pme.cl

-__kernel void updateGridIndexAndFraction(__global float4* posq, __global float2* pmeAtomGridIndex) {
+__kernel void updateGridIndexAndFraction(__global float4* posq, __global float2* pmeAtomGridIndex, float4 invPeriodicBoxSize) {
     for (int i = get_global_id(0); i < NUM_ATOMS; i += get_global_size(0)) {
         float4 pos = posq[i];
-        float4 t = (float4) ((pos.x*INV_PERIODIC_BOX_SIZE_X+1.0f)*GRID_SIZE_X,
-                             (pos.y*INV_PERIODIC_BOX_SIZE_Y+1.0f)*GRID_SIZE_Y,
-                             (pos.z*INV_PERIODIC_BOX_SIZE_Z+1.0f)*GRID_SIZE_Z, 0.0f);
+        float4 t = (float4) ((pos.x*invPeriodicBoxSize.x+1.0f)*GRID_SIZE_X,
+                             (pos.y*invPeriodicBoxSize.y+1.0f)*GRID_SIZE_Y,
+                             (pos.z*invPeriodicBoxSize.z+1.0f)*GRID_SIZE_Z, 0.0f);
         int4 gridIndex = (int4) (((int) t.x) % GRID_SIZE_X,
                                  ((int) t.y) % GRID_SIZE_Y,
                                  ((int) t.z) % GRID_SIZE_Z, 0);
@@ -38,7 +38,7 @@ __kernel void findAtomRangeForGrid(__global float4* posq, __global float2* pmeAt
     }
 }
 
-__kernel void updateBsplines(__global float4* posq, __global float4* pmeBsplineTheta, __global float4* pmeBsplineDTheta, __local float4* bsplinesCache, __global float2* pmeAtomGridIndex) {
+__kernel void updateBsplines(__global float4* posq, __global float4* pmeBsplineTheta, __global float4* pmeBsplineDTheta, __local float4* bsplinesCache, __global float2* pmeAtomGridIndex, float4 invPeriodicBoxSize) {
     const float4 scale = 1.0f/(PME_ORDER-1);
     for (int i = get_global_id(0); i < NUM_ATOMS; i += get_global_size(0)) {
         __local float4* data = &bsplinesCache[get_local_id(0)*PME_ORDER];
@@ -48,9 +48,9 @@ __kernel void updateBsplines(__global float4* posq, __global float4* pmeBsplineT
             ddata[j] = 0.0f;
         }
         float4 pos = posq[i];
-        float4 t = (float4) ((pos.x*INV_PERIODIC_BOX_SIZE_X+1.0f)*GRID_SIZE_X,
-                             (pos.y*INV_PERIODIC_BOX_SIZE_Y+1.0f)*GRID_SIZE_Y,
-                             (pos.z*INV_PERIODIC_BOX_SIZE_Z+1.0f)*GRID_SIZE_Z, 0.0f);
+        float4 t = (float4) ((pos.x*invPeriodicBoxSize.x+1.0f)*GRID_SIZE_X,
+                             (pos.y*invPeriodicBoxSize.y+1.0f)*GRID_SIZE_Y,
+                             (pos.z*invPeriodicBoxSize.z+1.0f)*GRID_SIZE_Z, 0.0f);
         float4 dr = (float4) (t.x-(int) t.x, t.y-(int) t.y, t.z-(int) t.z, 0.0f);
         data[PME_ORDER-1] = 0.0f;
         data[1] = dr;
@@ -119,7 +119,7 @@ __kernel void gridSpreadCharge(__global float2* pmeAtomGridIndex, __global int*
 }
 
 __kernel void reciprocalConvolution(__global float2* pmeGrid, __global float* energyBuffer, __global float* pmeBsplineModuliX,
-        __global float* pmeBsplineModuliY, __global float* pmeBsplineModuliZ) {
+        __global float* pmeBsplineModuliY, __global float* pmeBsplineModuliZ, float4 invPeriodicBoxSize, float recipScaleFactor) {
     const unsigned int gridSize = GRID_SIZE_X*GRID_SIZE_Y*GRID_SIZE_Z;
     float energy = 0.0f;
     for (int index = get_global_id(0); index < gridSize; index += get_global_size(0)) {
@@ -132,29 +132,29 @@ __kernel void reciprocalConvolution(__global float2* pmeGrid, __global float* en
         int mx = (kx < (GRID_SIZE_X+1)/2) ? kx : (kx-GRID_SIZE_X);
         int my = (ky < (GRID_SIZE_Y+1)/2) ? ky : (ky-GRID_SIZE_Y);
         int mz = (kz < (GRID_SIZE_Z+1)/2) ? kz : (kz-GRID_SIZE_Z);
-        float mhx = mx*INV_PERIODIC_BOX_SIZE_X;
-        float mhy = my*INV_PERIODIC_BOX_SIZE_Y;
-        float mhz = mz*INV_PERIODIC_BOX_SIZE_Z;
+        float mhx = mx*invPeriodicBoxSize.x;
+        float mhy = my*invPeriodicBoxSize.y;
+        float mhz = mz*invPeriodicBoxSize.z;
         float bx = pmeBsplineModuliX[kx];
         float by = pmeBsplineModuliY[ky];
         float bz = pmeBsplineModuliZ[kz];
         float2 grid = pmeGrid[index];
         float m2 = mhx*mhx+mhy*mhy+mhz*mhz;
         float denom = m2*bx*by*bz;
-        float eterm = RECIP_SCALE_FACTOR*exp(-RECIP_EXP_FACTOR*m2)/denom;
+        float eterm = recipScaleFactor*exp(-RECIP_EXP_FACTOR*m2)/denom;
         pmeGrid[index] = (float2) (grid.x*eterm, grid.y*eterm);
         energy += eterm*(grid.x*grid.x + grid.y*grid.y);
     }
     energyBuffer[get_global_id(0)] += 0.5f*energy;
 }
 
-__kernel void gridInterpolateForce(__global float4* posq, __global float4* forceBuffers, __global float4* pmeBsplineTheta, __global float4* pmeBsplineDTheta, __global float2* pmeGrid) {
+__kernel void gridInterpolateForce(__global float4* posq, __global float4* forceBuffers, __global float4* pmeBsplineTheta, __global float4* pmeBsplineDTheta, __global float2* pmeGrid, float4 invPeriodicBoxSize) {
     for (int atom = get_global_id(0); atom < NUM_ATOMS; atom += get_global_size(0)) {
         float4 force = 0.0f;
         float4 pos = posq[atom];
-        float4 t = (float4) ((pos.x*INV_PERIODIC_BOX_SIZE_X+1.0f)*GRID_SIZE_X,
-                             (pos.y*INV_PERIODIC_BOX_SIZE_Y+1.0f)*GRID_SIZE_Y,
-                             (pos.z*INV_PERIODIC_BOX_SIZE_Z+1.0f)*GRID_SIZE_Z, 0.0f);
+        float4 t = (float4) ((pos.x*invPeriodicBoxSize.x+1.0f)*GRID_SIZE_X,
+                             (pos.y*invPeriodicBoxSize.y+1.0f)*GRID_SIZE_Y,
+                             (pos.z*invPeriodicBoxSize.z+1.0f)*GRID_SIZE_Z, 0.0f);
         int4 gridIndex = (int4) (((int) t.x) % GRID_SIZE_X,
                                  ((int) t.y) % GRID_SIZE_Y,
                                  ((int) t.z) % GRID_SIZE_Z, 0);
@@ -180,9 +180,9 @@ __kernel void gridInterpolateForce(__global float4* posq, __global float4* force
         }
         float4 totalForce = forceBuffers[atom];
         float q = pos.w*EPSILON_FACTOR;
-        totalForce.x -= q*force.x*GRID_SIZE_X*INV_PERIODIC_BOX_SIZE_X;
-        totalForce.y -= q*force.y*GRID_SIZE_Y*INV_PERIODIC_BOX_SIZE_Y;
-        totalForce.z -= q*force.z*GRID_SIZE_Z*INV_PERIODIC_BOX_SIZE_Z;
+        totalForce.x -= q*force.x*GRID_SIZE_X*invPeriodicBoxSize.x;
+        totalForce.y -= q*force.y*GRID_SIZE_Y*invPeriodicBoxSize.y;
+        totalForce.z -= q*force.z*GRID_SIZE_Z*invPeriodicBoxSize.z;
         forceBuffers[atom] = totalForce;
     }
 }

platforms/opencl/tests/TestOpenCLCustomGBForce.cpp

@@ -64,8 +64,8 @@ void testOBC(GBSAOBCForce::NonbondedMethod obcMethod, CustomGBForce::NonbondedMe
         standardSystem.addParticle(1.0);
         customSystem.addParticle(1.0);
     }
-    standardSystem.setPeriodicBoxVectors(Vec3(boxSize, 0.0, 0.0), Vec3(0.0, boxSize, 0.0), Vec3(0.0, 0.0, boxSize));
-    customSystem.setPeriodicBoxVectors(Vec3(boxSize, 0.0, 0.0), Vec3(0.0, boxSize, 0.0), Vec3(0.0, 0.0, boxSize));
+    standardSystem.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0.0, 0.0), Vec3(0.0, boxSize, 0.0), Vec3(0.0, 0.0, boxSize));
+    customSystem.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0.0, 0.0), Vec3(0.0, boxSize, 0.0), Vec3(0.0, 0.0, boxSize));
     GBSAOBCForce* obc = new GBSAOBCForce();
     CustomGBForce* custom = new CustomGBForce();
     obc->setCutoffDistance(2.0);

platforms/opencl/tests/TestOpenCLCustomNonbondedForce.cpp

@@ -226,7 +226,7 @@ void testPeriodic() {
     forceField->addParticle(vector<double>());
     forceField->setNonbondedMethod(CustomNonbondedForce::CutoffPeriodic);
     forceField->setCutoffDistance(2.0);
-    system.setPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 4, 0), Vec3(0, 0, 4));
+    system.setDefaultPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 4, 0), Vec3(0, 0, 4));
     system.addForce(forceField);
     Context context(system, integrator, platform);
     vector<Vec3> positions(3);

platforms/opencl/tests/TestOpenCLEwald.cpp

@@ -78,7 +78,7 @@ void testEwaldPME(bool includeExceptions) {
         nonbonded->addParticle(1.0, 1.0,0.0);
     for (int i = 0; i < numParticles/2; i++)
         nonbonded->addParticle(-1.0, 1.0,0.0);
-    system.setPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
     system.addForce(nonbonded);
 
     vector<Vec3> positions(numParticles);
@@ -183,7 +183,7 @@ void testEwald2Ions() {
     const double cutoff = 2.0;
     nonbonded->setCutoffDistance(cutoff);
     nonbonded->setEwaldErrorTolerance(TOL);
-    system.setPeriodicBoxVectors(Vec3(6, 0, 0), Vec3(0, 6, 0), Vec3(0, 0, 6));
+    system.setDefaultPeriodicBoxVectors(Vec3(6, 0, 0), Vec3(0, 6, 0), Vec3(0, 0, 6));
     system.addForce(nonbonded);
     Context context(system, integrator, platform);
     vector<Vec3> positions(2);
@@ -204,7 +204,7 @@ void testErrorTolerance(NonbondedForce::NonbondedMethod method) {
     const int numParticles = 51;
     const double boxWidth = 5.0;
     System system;
-    system.setPeriodicBoxVectors(Vec3(boxWidth, 0, 0), Vec3(0, boxWidth, 0), Vec3(0, 0, boxWidth));
+    system.setDefaultPeriodicBoxVectors(Vec3(boxWidth, 0, 0), Vec3(0, boxWidth, 0), Vec3(0, 0, boxWidth));
     NonbondedForce* force = new NonbondedForce();
     system.addForce(force);
     vector<Vec3> positions(numParticles);

platforms/opencl/tests/TestOpenCLGBSAOBCForce.cpp

@@ -91,7 +91,7 @@ void testCutoffAndPeriodic() {
     const double boxSize = 10.0;
     nonbonded->setCutoffDistance(cutoffDistance);
     gbsa->setCutoffDistance(cutoffDistance);
-    system.setPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
     system.addForce(gbsa);
     system.addForce(nonbonded);
     vector<Vec3> positions(2);
@@ -152,7 +152,7 @@ void testForce(int numParticles, NonbondedForce::NonbondedMethod method, GBSAOBC
     int grid = (int) floor(0.5+pow(numParticles, 1.0/3.0));
     if (method == NonbondedForce::CutoffPeriodic) {
         double boxSize = (grid+1)*1.1;
-        system.setPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+        system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
     }
     system.addForce(gbsa);
     system.addForce(nonbonded);

platforms/opencl/tests/TestOpenCLNonbondedForce.cpp

@@ -338,7 +338,7 @@ void testPeriodic() {
     nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
     const double cutoff = 2.0;
     nonbonded->setCutoffDistance(cutoff);
-    system.setPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 4, 0), Vec3(0, 0, 4));
+    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(3);
@@ -418,7 +418,7 @@ void testLargeSystem() {
     // Now do the same thing with periodic boundary conditions.
 
     nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
-    system.setPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
     clContext.reinitialize();
     referenceContext.reinitialize();
     clContext.setPositions(positions);
@@ -456,7 +456,7 @@ void testBlockInteractions(bool periodic) {
     }
     nonbonded->setNonbondedMethod(periodic ? NonbondedForce::CutoffPeriodic : NonbondedForce::CutoffNonPeriodic);
     nonbonded->setCutoffDistance(cutoff);
-    system.setPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
     system.addForce(nonbonded);
     Context context(system, integrator, cl);
     context.setPositions(positions);