Do not rebuild the neighbor list every time step

platforms/cpu/include/CpuKernels.h

@@ -86,6 +86,7 @@ private:
     std::vector<std::pair<float, float> > particleParams;
     std::vector<float> posq;
     std::vector<float> forces;
+    std::vector<RealVec> lastPositions;
     NonbondedMethod nonbondedMethod;
     CpuNeighborList neighborList;
     CpuNonbondedForce nonbonded;

platforms/cpu/include/CpuNonbondedForce.h

@@ -122,7 +122,6 @@ class CpuNonbondedForce {
          @param atomParameters   atom parameters (sigma/2, 2*sqrt(epsilon))
          @param exclusions       atom exclusion indices
                                  exclusions[atomIndex] contains the list of exclusions for that atom
-         @param fixedParameters  non atom parameters (not currently used)
          @param forces           force array (forces added)
          @param totalEnergy      total energy
             
@@ -130,7 +129,7 @@ class CpuNonbondedForce {
           
       void calculateReciprocalIxn(int numberOfAtoms, float* posq, std::vector<OpenMM::RealVec>& atomCoordinates,
                             const std::vector<std::pair<float, float> >& atomParameters, const std::vector<std::set<int> >& exclusions,
-                            float* fixedParameters, std::vector<OpenMM::RealVec>& forces, float* totalEnergy) const;
+                            std::vector<OpenMM::RealVec>& forces, float* totalEnergy) const;
       
       /**---------------------------------------------------------------------------------------
       
@@ -141,15 +140,13 @@ class CpuNonbondedForce {
          @param atomParameters   atom parameters (sigma/2, 2*sqrt(epsilon))
          @param exclusions       atom exclusion indices
                                  exclusions[atomIndex] contains the list of exclusions for that atom
-         @param fixedParameters  non atom parameters (not currently used)
          @param forces           force array (forces added)
          @param totalEnergy      total energy
       
          --------------------------------------------------------------------------------------- */
           
-      void calculateDirectIxn(int numberOfAtoms, float* posq,
-                            const std::vector<std::pair<float, float> >& atomParameters, const std::vector<std::set<int> >& exclusions,
-                            float* fixedParameters, float* forces, float* totalEnergy);
+      void calculateDirectIxn(int numberOfAtoms, float* posq, const std::vector<std::pair<float, float> >& atomParameters,
+            const std::vector<std::set<int> >& exclusions, float* forces, float* totalEnergy);
 
     /**
      * This routine contains the code executed by each thread.

platforms/cpu/src/CpuKernels.cpp

@@ -172,6 +172,7 @@ void CpuCalcNonbondedForceKernel::initialize(const System& system, const Nonbond
         dispersionCoefficient = NonbondedForceImpl::calcDispersionCorrection(system, force);
     else
         dispersionCoefficient = 0.0;
+    lastPositions.resize(numParticles, Vec3(1e10, 1e10, 1e10));
 }
 
 double CpuCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy, bool includeDirect, bool includeReciprocal) {
@@ -218,7 +219,21 @@ double CpuCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeFo
     for (int i = 0; i < 4*numParticles; i++)
         forces[i] = 0.0f;
     if (nonbondedMethod != NoCutoff) {
-        neighborList.computeNeighborList(numParticles, posq, exclusions, floatBoxSize, periodic || ewald || pme, nonbondedCutoff);
+        // Determine whether we need to recompute the neighbor list.
+        
+        double padding = 0.1*nonbondedCutoff;
+        bool needRecompute = false;
+        for (int i = 0; i < numParticles; i++) {
+            RealVec delta = posData[i]-lastPositions[i];
+            if (delta.dot(delta) > 0.25*padding*padding) {
+                needRecompute = true;
+                break;
+            }
+        }
+        if (needRecompute) {
+            neighborList.computeNeighborList(numParticles, posq, exclusions, floatBoxSize, periodic || ewald || pme, nonbondedCutoff+padding);
+            lastPositions = posData;
+        }
         nonbonded.setUseCutoff(nonbondedCutoff, neighborList.getNeighbors(), rfDielectric);
     }
     if (periodic || ewald || pme) {
@@ -235,7 +250,7 @@ double CpuCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeFo
         nonbonded.setUseSwitchingFunction(switchingDistance);
     float nonbondedEnergy = 0;
     if (includeDirect)
-        nonbonded.calculateDirectIxn(numParticles, &posq[0], particleParams, exclusions, 0, &forces[0], includeEnergy ? &nonbondedEnergy : NULL);
+        nonbonded.calculateDirectIxn(numParticles, &posq[0], particleParams, exclusions, &forces[0], includeEnergy ? &nonbondedEnergy : NULL);
     if (includeReciprocal) {
         if (useOptimizedPme) {
             PmeIO io(&posq[0], &forces[0], numParticles);
@@ -244,7 +259,7 @@ double CpuCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeFo
             optimizedPme.getAs<CalcPmeReciprocalForceKernel>().finishComputation(io);
         }
         else
-            nonbonded.calculateReciprocalIxn(numParticles, &posq[0], posData, particleParams, exclusions, 0, forceData, includeEnergy ? &nonbondedEnergy : NULL);
+            nonbonded.calculateReciprocalIxn(numParticles, &posq[0], posData, particleParams, exclusions, forceData, includeEnergy ? &nonbondedEnergy : NULL);
     }
     energy += nonbondedEnergy;
     for (int i = 0; i < numParticles; i++) {

platforms/cpu/src/CpuNonbondedForce.cpp

@@ -191,7 +191,7 @@ void CpuNonbondedForce::setUseSwitchingFunction(float distance) {
 
 void CpuNonbondedForce::calculateReciprocalIxn(int numberOfAtoms, float* posq, vector<OpenMM::RealVec>& atomCoordinates,
                                              const vector<pair<float, float> >& atomParameters, const vector<set<int> >& exclusions,
-                                             float* fixedParameters, vector<OpenMM::RealVec>& forces, float* totalEnergy) const {
+                                             vector<OpenMM::RealVec>& forces, float* totalEnergy) const {
     typedef std::complex<float> d_complex;
 
     static const float epsilon     =  1.0;
@@ -303,9 +303,8 @@ void CpuNonbondedForce::calculateReciprocalIxn(int numberOfAtoms, float* posq, v
 }
 
 
-void CpuNonbondedForce::calculateDirectIxn(int numberOfAtoms, float* posq,
-                                             const vector<pair<float, float> >& atomParameters, const vector<set<int> >& exclusions,
-                                             float* fixedParameters, float* forces, float* totalEnergy) {
+void CpuNonbondedForce::calculateDirectIxn(int numberOfAtoms, float* posq, const vector<pair<float, float> >& atomParameters,
+                const vector<set<int> >& exclusions, float* forces, float* totalEnergy) {
     // Record the parameters for the threads.
     
     this->posq = posq;
@@ -425,6 +424,8 @@ void CpuNonbondedForce::calculateOneIxn(int ii, int jj, float* forces, double* t
     __m128 posJ = _mm_loadu_ps(posq+4*jj);
     float r2;
     getDeltaR(posJ, posI, deltaR, r2, periodic);
+    if (cutoff && r2 >= cutoffDistance*cutoffDistance)
+        return;
     float r = sqrtf(r2);
     float inverseR = 1/r;
     float switchValue = 1, switchDeriv = 0;
@@ -475,6 +476,8 @@ void CpuNonbondedForce::calculateOneEwaldIxn(int ii, int jj, float* forces, doub
     __m128 posJ = _mm_loadu_ps(posq+4*jj);
     float r2;
     getDeltaR(posJ, posI, deltaR, r2, true);
+    if (r2 >= cutoffDistance*cutoffDistance)
+        return;
     float r         = sqrtf(r2);
     float inverseR  = 1/r;
     float switchValue = 1, switchDeriv = 0;