Merge pull request #235 from peastman/master

Minor optimizations to CPU platform

platforms/cpu/src/CpuGBSAOBCForce.cpp

@@ -361,14 +361,13 @@ void CpuGBSAOBCForce::threadComputeForce(ThreadPool& threads, int threadIndex) {
             fvec4 logRatio = fastLog(u_ij/l_ij);
             fvec4 t3 = 0.125f*(1.0f + scaledRadiusJ2*r2Inverse)*(l_ij2 - u_ij2) + 0.25f*logRatio*r2Inverse;
             fvec4 de = bornForce*t3*rInverse;
+            de = blend(0.0f, de, include);
             fvec4 result[4] = {dx*de, dy*de, dz*de, 0.0f};
             transpose(result[0], result[1], result[2], result[3]);
             fvec4 atomForce(forces+4*atomJ);
             for (int j = 0; j < 4; j++) {
-                if (include[j]) {
-                    blockAtomForce[j] += result[j];
-                    atomForce -= result[j];
-                }
+                blockAtomForce[j] += result[j];
+                atomForce -= result[j];
             }
             atomForce.store(forces+4*atomJ);
         }

platforms/cpu/src/CpuKernels.cpp

@@ -84,13 +84,13 @@ void CpuCalcForcesAndEnergyKernel::beginComputation(ContextImpl& context, bool i
     AlignedArray<float>& posq = data.posq;
     vector<RealVec>& posData = extractPositions(context);
     RealVec boxSize = extractBoxSize(context);
-    float floatBoxSize[3] = {(float) boxSize[0], (float) boxSize[1], (float) boxSize[2]};
+    double invBoxSize[3] = {1/boxSize[0], 1/boxSize[1], 1/boxSize[2]};
     int numParticles = context.getSystem().getNumParticles();
     if (data.isPeriodic)
         for (int i = 0; i < numParticles; i++)
             for (int j = 0; j < 3; j++) {
                 RealOpenMM x = posData[i][j];
-                double base = floor(x/boxSize[j])*boxSize[j];
+                double base = floor(x*invBoxSize[j])*boxSize[j];
                 posq[4*i+j] = (float) (x-base);
             }
     else

platforms/cpu/src/CpuNeighborList.cpp

@@ -229,6 +229,8 @@ public:
                         maxz = max(maxz, atomPos[2]+dist);
                     }
                 }
+                if (minz == maxz)
+                    continue;
                 bool needPeriodic = (centerPos[0]-blockWidth[0] < maxDistance || centerPos[0]+blockWidth[0] > periodicBoxSize[0]-maxDistance ||
                                      centerPos[1]-blockWidth[1] < maxDistance || centerPos[1]+blockWidth[1] > periodicBoxSize[1]-maxDistance ||
                                      minz < 0.0f || maxz > periodicBoxSize[2]);

platforms/cpu/src/CpuNonbondedForce.cpp

@@ -472,6 +472,7 @@ void CpuNonbondedForce::calculateBlockIxn(int blockIndex, float* forces, double*
                     break;
                 }
     }
+    const float invSwitchingInterval = 1/(cutoffDistance-switchingDistance);
     
     // Loop over neighbors for this block.
     
@@ -512,9 +513,9 @@ void CpuNonbondedForce::calculateBlockIxn(int blockIndex, float* forces, double*
             dEdR = epsSig6*(12.0f*sig6 - 6.0f);
             energy = epsSig6*(sig6-1.0f);
             if (useSwitch) {
-                fvec4 t = (r>switchingDistance) & ((r-switchingDistance)/(cutoffDistance-switchingDistance));
+                fvec4 t = (r>switchingDistance) & ((r-switchingDistance)*invSwitchingInterval);
                 fvec4 switchValue = 1+t*t*t*(-10.0f+t*(15.0f-t*6.0f));
-                fvec4 switchDeriv = t*t*(-30.0f+t*(60.0f-t*30.0f))/(cutoffDistance-switchingDistance);
+                fvec4 switchDeriv = t*t*(-30.0f+t*(60.0f-t*30.0f))*invSwitchingInterval;
                 dEdR = switchValue*dEdR - energy*switchDeriv*r;
                 energy *= switchValue;
             }
@@ -584,6 +585,7 @@ void CpuNonbondedForce::calculateBlockEwaldIxn(int blockIndex, float* forces, do
                 needPeriodic = true;
                 break;
             }
+    const float invSwitchingInterval = 1/(cutoffDistance-switchingDistance);
     
     // Loop over neighbors for this block.
     
@@ -624,9 +626,9 @@ void CpuNonbondedForce::calculateBlockEwaldIxn(int blockIndex, float* forces, do
             dEdR = epsSig6*(12.0f*sig6 - 6.0f);
             energy = epsSig6*(sig6-1.0f);
             if (useSwitch) {
-                fvec4 t = (r>switchingDistance) & ((r-switchingDistance)/(cutoffDistance-switchingDistance));
+                fvec4 t = (r>switchingDistance) & ((r-switchingDistance)*invSwitchingInterval);
                 fvec4 switchValue = 1+t*t*t*(-10.0f+t*(15.0f-t*6.0f));
-                fvec4 switchDeriv = t*t*(-30.0f+t*(60.0f-t*30.0f))/(cutoffDistance-switchingDistance);
+                fvec4 switchDeriv = t*t*(-30.0f+t*(60.0f-t*30.0f))*invSwitchingInterval;
                 dEdR = switchValue*dEdR - energy*switchDeriv*r;
                 energy *= switchValue;
             }