diff --git a/src/finish.cpp b/src/finish.cpp
index 2fe3dea480..d9708edb5b 100644
--- a/src/finish.cpp
+++ b/src/finish.cpp
@@ -38,53 +38,76 @@ Finish::Finish(LAMMPS *lmp) : Pointers(lmp) {}
 
 void Finish::end(int flag)
 {
-  int i,m;
+  int i,m,nneigh;
   int histo[10];
-  double time,tmp,ave,max,min;
-
+  int loopflag,minflag,prdflag,timeflag,fftflag,histoflag,neighflag;
+  double time,tmp,ave,max,min,natoms;
+  double time_loop,time_other;
+  
   int me,nprocs;
   MPI_Comm_rank(world,&me);
   MPI_Comm_size(world,&nprocs);
 
-  // deduce time_other
+  // choose flavors of statistical output
+  // flag determines caller
+  // flag = 0 = just loop summary
+  // flag = 1 = dynamics or minimization
+  // flag = 2 = PRD
+  
+  loopflag = 1;
+  minflag = prdflag = timeflag = fftflag = histoflag = neighflag = 0;
 
-  double time_other = timer->array[TIME_LOOP] -
-    (timer->array[TIME_PAIR] + timer->array[TIME_BOND] + 
-     timer->array[TIME_KSPACE] + timer->array[TIME_NEIGHBOR] +
-     timer->array[TIME_COMM] + timer->array[TIME_OUTPUT]);
-
-  double time_loop = timer->array[TIME_LOOP];
-  MPI_Allreduce(&time_loop,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
-  time_loop = tmp/nprocs;
-
-  // overall loop time
-  // use actual natoms, in case atoms were lost
-
-  double natoms;
-  double rlocal = atom->nlocal;
-  MPI_Allreduce(&rlocal,&natoms,1,MPI_DOUBLE,MPI_SUM,world);
-
-  if (me == 0) {
-    if (screen) 
-      fprintf(screen,
-	      "Loop time of %g on %d procs for %d steps with %.15g atoms\n",
-	      time_loop,nprocs,update->nsteps,natoms);
-    if (logfile)
-      fprintf(logfile,
-	      "Loop time of %g on %d procs for %d steps with %.15g atoms\n",
-	      time_loop,nprocs,update->nsteps,natoms);
+  if (flag == 1) {
+    if (update->whichflag == 2) minflag = 1;
+    timeflag = histoflag = neighflag = 1;
+    if (strstr(force->kspace_style,"pppm")) fftflag = 1;
+  }
+  if (flag == 2) {
+    prdflag = histoflag = neighflag = 1;
+    
   }
 
-  if (flag == 0) return;
+  // loop stats
 
-  if (me == 0) {
-    if (screen) fprintf(screen,"\n");
-    if (logfile) fprintf(logfile,"\n");
+  if (loopflag) {
+    time_other = timer->array[TIME_LOOP] -
+      (timer->array[TIME_PAIR] + timer->array[TIME_BOND] + 
+       timer->array[TIME_KSPACE] + timer->array[TIME_NEIGHBOR] +
+       timer->array[TIME_COMM] + timer->array[TIME_OUTPUT]);
+    
+    time_loop = timer->array[TIME_LOOP];
+    MPI_Allreduce(&time_loop,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
+    time_loop = tmp/nprocs;
+
+    // overall loop time
+    // use actual natoms, in case atoms were lost
+
+    natoms;
+    double rlocal = atom->nlocal;
+    MPI_Allreduce(&rlocal,&natoms,1,MPI_DOUBLE,MPI_SUM,world);
+    
+    if (me == 0) {
+      if (screen) 
+	fprintf(screen,
+		"Loop time of %g on %d procs for %d steps with %.15g atoms\n",
+		time_loop,nprocs,update->nsteps,natoms);
+      if (logfile)
+	fprintf(logfile,
+		"Loop time of %g on %d procs for %d steps with %.15g atoms\n",
+		time_loop,nprocs,update->nsteps,natoms);
+    }
+
+    if (time_loop == 0.0) time_loop = 1.0;
   }
 
   // minimization stats
 
-  if (update->whichflag == 2) {
+  if (minflag) {
+    if (me == 0) {
+      if (screen) fprintf(screen,"\n");
+      if (logfile) fprintf(logfile,"\n");
+    }
+
     if (me == 0) {
       if (screen) {
 	fprintf(screen,"Minimization stats:\n");
@@ -127,98 +150,170 @@ void Finish::end(int flag)
 		update->minimize->niter,update->minimize->neval);
       }
     }
+  }
+
+
+  // PRD stats using PAIR,BOND,KSPACE for dephase,dynamics,quench
+
+  if (prdflag) {
     if (me == 0) {
       if (screen) fprintf(screen,"\n");
       if (logfile) fprintf(logfile,"\n");
     }
-  }
 
-  // timing breakdowns
+    if (screen) fprintf(screen,"PRD stats:\n");
+    if (logfile) fprintf(logfile,"PRD stats:\n");
 
-  if (time_loop == 0.0) time_loop = 1.0;
+    time = timer->array[TIME_PAIR];
+    MPI_Allreduce(&time,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
+    time = tmp/nprocs;
+    if (me == 0) {
+      if (screen) 
+	fprintf(screen,"  Dephase  time (%%) = %g (%g)\n",
+		time,time/time_loop*100.0);
+      if (logfile) 
+	fprintf(logfile,"  Dephase  time (%%) = %g (%g)\n",
+		time,time/time_loop*100.0);
+    }
 
-  time = timer->array[TIME_PAIR];
-  MPI_Allreduce(&time,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
-  time = tmp/nprocs;
-  if (me == 0) {
-    if (screen) 
-      fprintf(screen,"Pair  time (%%) = %g (%g)\n",time,time/time_loop*100.0);
-    if (logfile) 
-      fprintf(logfile,"Pair  time (%%) = %g (%g)\n",time,time/time_loop*100.0);
-  }
-
-  if (atom->molecular) {
     time = timer->array[TIME_BOND];
     MPI_Allreduce(&time,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
     time = tmp/nprocs;
     if (me == 0) {
       if (screen) 
-	fprintf(screen,"Bond  time (%%) = %g (%g)\n",
+	fprintf(screen,"  Dynamics time (%%) = %g (%g)\n",
 		time,time/time_loop*100.0);
-      if (logfile)
-	fprintf(logfile,"Bond  time (%%) = %g (%g)\n",
+      if (logfile) 
+	fprintf(logfile,"  Dynamics time (%%) = %g (%g)\n",
 		time,time/time_loop*100.0);
     }
-  }
 
-  if (force->kspace) {
     time = timer->array[TIME_KSPACE];
     MPI_Allreduce(&time,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
     time = tmp/nprocs;
     if (me == 0) {
       if (screen) 
-	fprintf(screen,"Kspce time (%%) = %g (%g)\n",
+	fprintf(screen,"  Quench   time (%%) = %g (%g)\n",
 		time,time/time_loop*100.0);
-      if (logfile)
-	fprintf(logfile,"Kspce time (%%) = %g (%g)\n",
+      if (logfile) 
+	fprintf(logfile,"  Quench   time (%%) = %g (%g)\n",
+		time,time/time_loop*100.0);
+    }
+
+    time = time_other;
+    MPI_Allreduce(&time,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
+    time = tmp/nprocs;
+    if (me == 0) {
+      if (screen) 
+	fprintf(screen,"  Other    time (%%) = %g (%g)\n",
+		time,time/time_loop*100.0);
+      if (logfile) 
+	fprintf(logfile,"  Other    time (%%) = %g (%g)\n",
 		time,time/time_loop*100.0);
     }
   }
 
-  time = timer->array[TIME_NEIGHBOR];
-  MPI_Allreduce(&time,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
-  time = tmp/nprocs;
-  if (me == 0) {
-    if (screen) 
-      fprintf(screen,"Neigh time (%%) = %g (%g)\n",time,time/time_loop*100.0);
-    if (logfile) 
-      fprintf(logfile,"Neigh time (%%) = %g (%g)\n",time,time/time_loop*100.0);
-  }
+  // timing breakdowns
 
-  time = timer->array[TIME_COMM];
-  MPI_Allreduce(&time,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
-  time = tmp/nprocs;
-  if (me == 0) {
-    if (screen) 
-      fprintf(screen,"Comm  time (%%) = %g (%g)\n",time,time/time_loop*100.0);
-    if (logfile) 
-      fprintf(logfile,"Comm  time (%%) = %g (%g)\n",time,time/time_loop*100.0);
-  }
+  if (timeflag) {
+    if (me == 0) {
+      if (screen) fprintf(screen,"\n");
+      if (logfile) fprintf(logfile,"\n");
+    }
 
-  time = timer->array[TIME_OUTPUT];
-  MPI_Allreduce(&time,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
-  time = tmp/nprocs;
-  if (me == 0) {
-    if (screen) 
-      fprintf(screen,"Outpt time (%%) = %g (%g)\n",time,time/time_loop*100.0);
-    if (logfile) 
-      fprintf(logfile,"Outpt time (%%) = %g (%g)\n",time,time/time_loop*100.0);
-  }
+    time = timer->array[TIME_PAIR];
+    MPI_Allreduce(&time,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
+    time = tmp/nprocs;
+    if (me == 0) {
+      if (screen) 
+	fprintf(screen,"Pair  time (%%) = %g (%g)\n",
+		time,time/time_loop*100.0);
+      if (logfile) 
+	fprintf(logfile,"Pair  time (%%) = %g (%g)\n",
+		time,time/time_loop*100.0);
+    }
 
-  time = time_other;
-  MPI_Allreduce(&time,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
-  time = tmp/nprocs;
-  if (me == 0) {
-    if (screen) 
-      fprintf(screen,"Other time (%%) = %g (%g)\n",time,time/time_loop*100.0);
-    if (logfile) 
-      fprintf(logfile,"Other time (%%) = %g (%g)\n",time,time/time_loop*100.0);
+    if (atom->molecular) {
+      time = timer->array[TIME_BOND];
+      MPI_Allreduce(&time,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
+      time = tmp/nprocs;
+      if (me == 0) {
+	if (screen) 
+	  fprintf(screen,"Bond  time (%%) = %g (%g)\n",
+		  time,time/time_loop*100.0);
+	if (logfile)
+	  fprintf(logfile,"Bond  time (%%) = %g (%g)\n",
+		  time,time/time_loop*100.0);
+      }
+    }
+    
+    if (force->kspace) {
+      time = timer->array[TIME_KSPACE];
+      MPI_Allreduce(&time,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
+      time = tmp/nprocs;
+      if (me == 0) {
+	if (screen) 
+	  fprintf(screen,"Kspce time (%%) = %g (%g)\n",
+		  time,time/time_loop*100.0);
+	if (logfile)
+	  fprintf(logfile,"Kspce time (%%) = %g (%g)\n",
+		  time,time/time_loop*100.0);
+      }
+    }
+    
+    time = timer->array[TIME_NEIGHBOR];
+    MPI_Allreduce(&time,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
+    time = tmp/nprocs;
+    if (me == 0) {
+      if (screen) 
+	fprintf(screen,"Neigh time (%%) = %g (%g)\n",
+		time,time/time_loop*100.0);
+      if (logfile) 
+	fprintf(logfile,"Neigh time (%%) = %g (%g)\n",
+		time,time/time_loop*100.0);
+    }
+    
+    time = timer->array[TIME_COMM];
+    MPI_Allreduce(&time,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
+    time = tmp/nprocs;
+    if (me == 0) {
+      if (screen) 
+	fprintf(screen,"Comm  time (%%) = %g (%g)\n",
+		time,time/time_loop*100.0);
+      if (logfile) 
+	fprintf(logfile,"Comm  time (%%) = %g (%g)\n",
+		time,time/time_loop*100.0);
+    }
+    
+    time = timer->array[TIME_OUTPUT];
+    MPI_Allreduce(&time,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
+    time = tmp/nprocs;
+    if (me == 0) {
+      if (screen) 
+	fprintf(screen,"Outpt time (%%) = %g (%g)\n",
+		time,time/time_loop*100.0);
+      if (logfile) 
+	fprintf(logfile,"Outpt time (%%) = %g (%g)\n",
+		time,time/time_loop*100.0);
+    }
+    
+    time = time_other;
+    MPI_Allreduce(&time,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
+    time = tmp/nprocs;
+    if (me == 0) {
+      if (screen) 
+	fprintf(screen,"Other time (%%) = %g (%g)\n",
+		time,time/time_loop*100.0);
+      if (logfile) 
+	fprintf(logfile,"Other time (%%) = %g (%g)\n",
+		time,time/time_loop*100.0);
+    }
   }
-
+  
   // FFT timing statistics
   // time3d,time1d = total time during run for 3d and 1d FFTs
 
-  if (strstr(force->kspace_style,"pppm")) {
+  if (fftflag) {
     if (me == 0) {
       if (screen) fprintf(screen,"\n");
       if (logfile) fprintf(logfile,"\n");
@@ -265,151 +360,157 @@ void Finish::end(int flag)
     }
   }
 
-  if (me == 0) {
-    if (screen) fprintf(screen,"\n");
-    if (logfile) fprintf(logfile,"\n");
-  }
-
-  tmp = atom->nlocal;
-  stats(1,&tmp,&ave,&max,&min,10,histo);
-  if (me == 0) {
-    if (screen) {
-      fprintf(screen,"Nlocal:    %g ave %g max %g min\n",ave,max,min);
-      fprintf(screen,"Histogram:");
-      for (i = 0; i < 10; i++) fprintf(screen," %d",histo[i]);
-      fprintf(screen,"\n");
+  if (histoflag) {
+    if (me == 0) {
+      if (screen) fprintf(screen,"\n");
+      if (logfile) fprintf(logfile,"\n");
     }
-    if (logfile) {
-      fprintf(logfile,"Nlocal:    %g ave %g max %g min\n",ave,max,min);
-      fprintf(logfile,"Histogram:");
-      for (i = 0; i < 10; i++) fprintf(logfile," %d",histo[i]);
-      fprintf(logfile,"\n");
-    }
-  }
 
-  tmp = atom->nghost;
-  stats(1,&tmp,&ave,&max,&min,10,histo);
-  if (me == 0) {
-    if (screen) {
-      fprintf(screen,"Nghost:    %g ave %g max %g min\n",ave,max,min);
-      fprintf(screen,"Histogram:");
-      for (i = 0; i < 10; i++) fprintf(screen," %d",histo[i]);
-      fprintf(screen,"\n");
-    }
-    if (logfile) {
-      fprintf(logfile,"Nghost:    %g ave %g max %g min\n",ave,max,min);
-      fprintf(logfile,"Histogram:");
-      for (i = 0; i < 10; i++) fprintf(logfile," %d",histo[i]);
-      fprintf(logfile,"\n");
-    }
-  }
-
-  // find a non-skip neighbor list containing half the pairwise interactions
-  // count neighbors in that list for stats purposes
-
-  for (m = 0; m < neighbor->old_nrequest; m++)
-    if ((neighbor->old_requests[m]->half || neighbor->old_requests[m]->gran ||
-	 neighbor->old_requests[m]->respaouter ||
-	 neighbor->old_requests[m]->half_from_full) &&
-	neighbor->old_requests[m]->skip == 0) break;
-
-  int nneigh = 0;
-  if (m < neighbor->old_nrequest) {
-    int inum = neighbor->lists[m]->inum;
-    int *ilist = neighbor->lists[m]->ilist;
-    int *numneigh = neighbor->lists[m]->numneigh;
-    for (int ii = 0; ii < inum; ii++)
-      nneigh += numneigh[ilist[ii]];
-  }
-
-  tmp = nneigh;
-  stats(1,&tmp,&ave,&max,&min,10,histo);
-  if (me == 0) {
-    if (screen) {
-      fprintf(screen,"Neighs:    %g ave %g max %g min\n",ave,max,min);
-      fprintf(screen,"Histogram:");
-      for (i = 0; i < 10; i++) fprintf(screen," %d",histo[i]);
-      fprintf(screen,"\n");
-    }
-    if (logfile) {
-      fprintf(logfile,"Neighs:    %g ave %g max %g min\n",ave,max,min);
-      fprintf(logfile,"Histogram:");
-      for (i = 0; i < 10; i++) fprintf(logfile," %d",histo[i]);
-      fprintf(logfile,"\n");
-    }
-  }
-
-  // find a non-skip neighbor list containing full pairwise interactions
-
-  for (m = 0; m < neighbor->old_nrequest; m++)
-    if (neighbor->old_requests[m]->full &&
-	neighbor->old_requests[m]->skip == 0) break;
-
-  if (m < neighbor->old_nrequest) {
-
-    nneigh = 0;
-    for (i = 0; i < atom->nlocal; i++)
-      nneigh += neighbor->lists[m]->numneigh[i];
-
-    tmp = nneigh;
+    tmp = atom->nlocal;
     stats(1,&tmp,&ave,&max,&min,10,histo);
     if (me == 0) {
       if (screen) {
-	fprintf(screen,"FullNghs:  %g ave %g max %g min\n",ave,max,min);
+	fprintf(screen,"Nlocal:    %g ave %g max %g min\n",ave,max,min);
 	fprintf(screen,"Histogram:");
 	for (i = 0; i < 10; i++) fprintf(screen," %d",histo[i]);
 	fprintf(screen,"\n");
       }
       if (logfile) {
-	fprintf(logfile,"FullNghs: %g ave %g max %g min\n",ave,max,min);
+	fprintf(logfile,"Nlocal:    %g ave %g max %g min\n",ave,max,min);
 	fprintf(logfile,"Histogram:");
 	for (i = 0; i < 10; i++) fprintf(logfile," %d",histo[i]);
 	fprintf(logfile,"\n");
       }
     }
-  }
-
-  if (me == 0) {
-    if (screen) fprintf(screen,"\n");
-    if (logfile) fprintf(logfile,"\n");
-  }
-
-  tmp = nneigh;
-  double nall;
-  MPI_Allreduce(&tmp,&nall,1,MPI_DOUBLE,MPI_SUM,world);
-
-  int nspec;
-  double nspec_all;
-  if (atom->molecular) {
-    nspec = 0;
-    for (i = 0; i < atom->nlocal; i++) nspec += atom->nspecial[i][2];
-    tmp = nspec;
-    MPI_Allreduce(&tmp,&nspec_all,1,MPI_DOUBLE,MPI_SUM,world);
-  }
-
-  if (me == 0) {
-    if (screen) {
-      if (nall < 2.0e9) 
-	fprintf(screen,"Total # of neighbors = %d\n",static_cast<int> (nall));
-      else fprintf(screen,"Total # of neighbors = %g\n",nall);
-      if (natoms > 0) fprintf(screen,"Ave neighs/atom = %g\n",nall/natoms);
-      if (atom->molecular && natoms > 0) 
-	fprintf(screen,"Ave special neighs/atom = %g\n",nspec_all/natoms);
-      fprintf(screen,"Neighbor list builds = %d\n",neighbor->ncalls);
-      fprintf(screen,"Dangerous builds = %d\n",neighbor->ndanger);
+    
+    tmp = atom->nghost;
+    stats(1,&tmp,&ave,&max,&min,10,histo);
+    if (me == 0) {
+      if (screen) {
+	fprintf(screen,"Nghost:    %g ave %g max %g min\n",ave,max,min);
+	fprintf(screen,"Histogram:");
+	for (i = 0; i < 10; i++) fprintf(screen," %d",histo[i]);
+	fprintf(screen,"\n");
+      }
+      if (logfile) {
+	fprintf(logfile,"Nghost:    %g ave %g max %g min\n",ave,max,min);
+	fprintf(logfile,"Histogram:");
+	for (i = 0; i < 10; i++) fprintf(logfile," %d",histo[i]);
+	fprintf(logfile,"\n");
+      }
     }
-    if (logfile) {
-      if (nall < 2.0e9) 
-	fprintf(logfile,"Total # of neighbors = %d\n",static_cast<int> (nall));
-      else fprintf(logfile,"Total # of neighbors = %g\n",nall);
-      if (natoms > 0) fprintf(logfile,"Ave neighs/atom = %g\n",nall/natoms);
-      if (atom->molecular && natoms > 0) 
-	fprintf(logfile,"Ave special neighs/atom = %g\n",nspec_all/natoms);
-      fprintf(logfile,"Neighbor list builds = %d\n",neighbor->ncalls);
-      fprintf(logfile,"Dangerous builds = %d\n",neighbor->ndanger);
+    
+    // find a non-skip neighbor list containing half the pairwise interactions
+    // count neighbors in that list for stats purposes
+    
+    for (m = 0; m < neighbor->old_nrequest; m++)
+      if ((neighbor->old_requests[m]->half || 
+	   neighbor->old_requests[m]->gran ||
+	   neighbor->old_requests[m]->respaouter ||
+	   neighbor->old_requests[m]->half_from_full) &&
+	  neighbor->old_requests[m]->skip == 0) break;
+    
+    nneigh = 0;
+    if (m < neighbor->old_nrequest) {
+      int inum = neighbor->lists[m]->inum;
+      int *ilist = neighbor->lists[m]->ilist;
+      int *numneigh = neighbor->lists[m]->numneigh;
+      for (int ii = 0; ii < inum; ii++)
+	nneigh += numneigh[ilist[ii]];
+    }
+    
+    tmp = nneigh;
+    stats(1,&tmp,&ave,&max,&min,10,histo);
+    if (me == 0) {
+      if (screen) {
+	fprintf(screen,"Neighs:    %g ave %g max %g min\n",ave,max,min);
+	fprintf(screen,"Histogram:");
+	for (i = 0; i < 10; i++) fprintf(screen," %d",histo[i]);
+	fprintf(screen,"\n");
+      }
+      if (logfile) {
+	fprintf(logfile,"Neighs:    %g ave %g max %g min\n",ave,max,min);
+	fprintf(logfile,"Histogram:");
+	for (i = 0; i < 10; i++) fprintf(logfile," %d",histo[i]);
+	fprintf(logfile,"\n");
+      }
+    }
+    
+    // find a non-skip neighbor list containing full pairwise interactions
+    
+    for (m = 0; m < neighbor->old_nrequest; m++)
+      if (neighbor->old_requests[m]->full &&
+	  neighbor->old_requests[m]->skip == 0) break;
+    
+    if (m < neighbor->old_nrequest) {
+      nneigh = 0;
+      for (i = 0; i < atom->nlocal; i++)
+	nneigh += neighbor->lists[m]->numneigh[i];
+      
+      tmp = nneigh;
+      stats(1,&tmp,&ave,&max,&min,10,histo);
+      if (me == 0) {
+	if (screen) {
+	  fprintf(screen,"FullNghs:  %g ave %g max %g min\n",ave,max,min);
+	  fprintf(screen,"Histogram:");
+	  for (i = 0; i < 10; i++) fprintf(screen," %d",histo[i]);
+	  fprintf(screen,"\n");
+	}
+	if (logfile) {
+	  fprintf(logfile,"FullNghs: %g ave %g max %g min\n",ave,max,min);
+	  fprintf(logfile,"Histogram:");
+	  for (i = 0; i < 10; i++) fprintf(logfile," %d",histo[i]);
+	  fprintf(logfile,"\n");
+	}
+      }
     }
   }
 
+  if (neighflag) {
+    if (me == 0) {
+      if (screen) fprintf(screen,"\n");
+      if (logfile) fprintf(logfile,"\n");
+    }
+    
+    tmp = nneigh;
+    double nall;
+    MPI_Allreduce(&tmp,&nall,1,MPI_DOUBLE,MPI_SUM,world);
+    
+    int nspec;
+    double nspec_all;
+    if (atom->molecular) {
+      nspec = 0;
+      for (i = 0; i < atom->nlocal; i++) nspec += atom->nspecial[i][2];
+      tmp = nspec;
+      MPI_Allreduce(&tmp,&nspec_all,1,MPI_DOUBLE,MPI_SUM,world);
+    }
+    
+    if (me == 0) {
+      if (screen) {
+	if (nall < 2.0e9) 
+	  fprintf(screen,
+		  "Total # of neighbors = %d\n",static_cast<int> (nall));
+	else fprintf(screen,"Total # of neighbors = %g\n",nall);
+	if (natoms > 0) fprintf(screen,"Ave neighs/atom = %g\n",nall/natoms);
+	if (atom->molecular && natoms > 0) 
+	  fprintf(screen,"Ave special neighs/atom = %g\n",nspec_all/natoms);
+	fprintf(screen,"Neighbor list builds = %d\n",neighbor->ncalls);
+	fprintf(screen,"Dangerous builds = %d\n",neighbor->ndanger);
+      }
+      if (logfile) {
+	if (nall < 2.0e9) 
+	  fprintf(logfile,
+		  "Total # of neighbors = %d\n",static_cast<int> (nall));
+	else fprintf(logfile,"Total # of neighbors = %g\n",nall);
+	if (natoms > 0) fprintf(logfile,"Ave neighs/atom = %g\n",nall/natoms);
+	if (atom->molecular && natoms > 0) 
+	  fprintf(logfile,"Ave special neighs/atom = %g\n",nspec_all/natoms);
+	fprintf(logfile,"Neighbor list builds = %d\n",neighbor->ncalls);
+	fprintf(logfile,"Dangerous builds = %d\n",neighbor->ndanger);
+      }
+    }
+  }
+  
   if (logfile) fflush(logfile);
 }
 
diff --git a/src/fix_minimize.cpp b/src/fix_minimize.cpp
index b3ea54c2f5..0719b9857f 100644
--- a/src/fix_minimize.cpp
+++ b/src/fix_minimize.cpp
@@ -181,8 +181,9 @@ double FixMinimize::memory_usage()
 void FixMinimize::grow_arrays(int nmax)
 {
   for (int m = 0; m < nvector; m++)
-    vectors[m] = (double *) memory->srealloc(vectors[m],peratom[m]*nmax,
-					     "minimize:vector");
+    vectors[m] = (double *) 
+      memory->srealloc(vectors[m],peratom[m]*nmax*sizeof(double),
+		       "minimize:vector");
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/integrate.cpp b/src/integrate.cpp
index eaabb24c8e..f3749a37b9 100644
--- a/src/integrate.cpp
+++ b/src/integrate.cpp
@@ -79,7 +79,8 @@ void Integrate::ev_setup()
 }
 
 /* ----------------------------------------------------------------------
-   set eflag,vflag for current iteration with ntimestep
+   set eflag,vflag for current iteration
+   based on computes that need info on this ntimestep
    eflag = 0 = no energy computation
    eflag = 1 = global energy only
    eflag = 2 = per-atom energy only
diff --git a/src/integrate.h b/src/integrate.h
index ff23907bf7..88ede97adc 100644
--- a/src/integrate.h
+++ b/src/integrate.h
@@ -24,7 +24,8 @@ class Integrate : protected Pointers {
   virtual ~Integrate();
   virtual void init() = 0;
   virtual void setup() = 0;
-  virtual void iterate(int) = 0;
+  virtual void setup_minimal(int) = 0;
+  virtual void run(int) = 0;
   virtual void cleanup() {}
   virtual void reset_dt() {}
   virtual double memory_usage() {return 0.0;}
diff --git a/src/min.cpp b/src/min.cpp
index cc31f65bbe..a1a0264a2d 100644
--- a/src/min.cpp
+++ b/src/min.cpp
@@ -149,8 +149,6 @@ void Min::init()
   neigh_delay = neighbor->delay;
   neigh_dist_check = neighbor->dist_check;
   
-  // reset reneighboring criteria if necessary
-
   if (neigh_every != 1 || neigh_delay != 0 || neigh_dist_check != 1) {
     if (comm->me == 0) 
       error->warning("Resetting reneighboring criteria during minimization");
@@ -179,6 +177,12 @@ void Min::setup()
   nextra_global = modify->min_dof();
   if (nextra_global) fextra = new double[nextra_global];
 
+  // compute for potential energy
+
+  int id = modify->find_compute("thermo_pe");
+  if (id < 0) error->all("Minimization could not find thermo_pe compute");
+  pe_compute = modify->compute[id];
+
   // style-specific setup does two tasks
   // setup extra global dof vectors
   // setup extra per-atom dof vectors due to requests from Pair classes
@@ -198,7 +202,6 @@ void Min::setup()
   // setup domain, communication and neighboring
   // acquire ghosts
   // build neighbor lists
-  // reset gradient vector ptrs
 
   if (triclinic) domain->x2lamda(atom->nlocal);
   domain->pbc();
@@ -241,10 +244,63 @@ void Min::setup()
 }
 
 /* ----------------------------------------------------------------------
-   perform minimization, with setup first
+   setup without output or one-time post-init setup
+   flag = 0 = just force calculation
+   flag = 1 = reneighbor and force calculation
 ------------------------------------------------------------------------- */
 
-void Min::run()
+void Min::setup_minimal(int flag)
+{
+  // setup domain, communication and neighboring
+  // acquire ghosts
+  // build neighbor lists
+
+  if (flag) {
+    if (triclinic) domain->x2lamda(atom->nlocal);
+    domain->pbc();
+    domain->reset_box();
+    comm->setup();
+    if (neighbor->style) neighbor->setup_bins();
+    comm->exchange();
+    comm->borders();
+    if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
+    neighbor->build();
+    neighbor->ncalls = 0;
+  }
+
+  // compute all forces
+
+  ev_set(update->ntimestep);
+  force_clear();
+  
+  if (force->pair) force->pair->compute(eflag,vflag);
+
+  if (atom->molecular) {
+    if (force->bond) force->bond->compute(eflag,vflag);
+    if (force->angle) force->angle->compute(eflag,vflag);
+    if (force->dihedral) force->dihedral->compute(eflag,vflag);
+    if (force->improper) force->improper->compute(eflag,vflag);
+  }
+
+  if (force->kspace) {
+    force->kspace->setup();
+    force->kspace->compute(eflag,vflag);
+  }
+
+  if (force->newton) comm->reverse_communicate();
+
+  modify->setup(vflag);
+
+  // atoms may have migrated in comm->exchange()
+
+  reset_vectors();
+}
+
+/* ----------------------------------------------------------------------
+   perform minimization, calling iterate() for nsteps
+------------------------------------------------------------------------- */
+
+void Min::run(int nsteps)
 {
   // possible stop conditions
 
@@ -254,12 +310,6 @@ void Min::run()
 			 "linesearch alpha is zero",
 			 "forces are zero","quadratic factors are zero"};
 
-  // compute for potential energy
-
-  int id = modify->find_compute("thermo_pe");
-  if (id < 0) error->all("Minimization could not find thermo_pe compute");
-  pe_compute = modify->compute[id];
-
   // stats for Finish to print
 
   ecurrent = pe_compute->compute_scalar();
@@ -272,8 +322,7 @@ void Min::run()
 
   // minimizer iterations
 
-  timer->barrier_start(TIME_LOOP);
-  int stop_condition = iterate(update->nsteps);
+  int stop_condition = iterate(nsteps);
   stopstr = stopstrings[stop_condition];
 
   // account for early exit from iterate loop due to convergence
@@ -284,14 +333,17 @@ void Min::run()
   // add ntimestep to all computes that store invocation times
   //   since are hardwireing call to thermo/dumps and computes may not be ready
 
-  if (niter < update->nsteps) {
+  if (niter < nsteps) {
     niter++;
     update->nsteps = niter;
 
-    for (int idump = 0; idump < output->ndump; idump++)
-      output->next_dump[idump] = update->ntimestep;
-    output->next_dump_any = update->ntimestep;
-    if (output->restart_every) output->next_restart = update->ntimestep;
+    if (update->restrict_output == 0) {
+      for (int idump = 0; idump < output->ndump; idump++)
+	output->next_dump[idump] = update->ntimestep;
+      output->next_dump_any = update->ntimestep;
+      if (output->restart_every == 0) 
+	output->next_restart = update->ntimestep;
+    }
     output->next_thermo = update->ntimestep;
 
     modify->addstep_compute_all(update->ntimestep);
@@ -299,14 +351,6 @@ void Min::run()
     output->write(update->ntimestep);
   }
 
-  timer->barrier_stop(TIME_LOOP);
-
-  // reset reneighboring criteria
-
-  neighbor->every = neigh_every;
-  neighbor->delay = neigh_delay;
-  neighbor->dist_check = neigh_dist_check;
-
   // stats for Finish to print
 	
   efinal = ecurrent;
@@ -314,12 +358,18 @@ void Min::run()
   fnorminf_final = fnorm_inf();
 }
 
-/* ----------------------------------------------------------------------
-   delete fix at end of run, so its atom arrays won't persist
-------------------------------------------------------------------------- */
+/* ---------------------------------------------------------------------- */
 
 void Min::cleanup()
 {
+  // reset reneighboring criteria
+
+  neighbor->every = neigh_every;
+  neighbor->delay = neigh_delay;
+  neighbor->dist_check = neigh_dist_check;
+
+  // delete fix at end of run, so its atom arrays won't persist
+
   modify->delete_fix("MINIMIZE");
 }
 
@@ -527,7 +577,8 @@ void Min::ev_setup()
 }
 
 /* ----------------------------------------------------------------------
-   set eflag,vflag for current iteration with ntimestep
+   set eflag,vflag for current iteration
+   based on computes that need info on this ntimestep
    always set eflag_global = 1, since need energy every iteration
    eflag = 0 = no energy computation
    eflag = 1 = global energy only
diff --git a/src/min.h b/src/min.h
index 30f86fbefd..6091d2c3ad 100644
--- a/src/min.h
+++ b/src/min.h
@@ -30,7 +30,8 @@ class Min : protected Pointers {
   virtual ~Min();
   void init();
   void setup();
-  void run();
+  void setup_minimal(int);
+  void run(int);
   void cleanup();
   void request(class Pair *, int, double);
   double memory_usage() {return 0.0;}
diff --git a/src/minimize.cpp b/src/minimize.cpp
index a324e1f7bc..94d748f9dd 100644
--- a/src/minimize.cpp
+++ b/src/minimize.cpp
@@ -17,6 +17,7 @@
 #include "update.h"
 #include "min.h"
 #include "finish.h"
+#include "timer.h"
 #include "error.h"
 
 using namespace LAMMPS_NS;
@@ -48,7 +49,11 @@ void Minimize::command(int narg, char **arg)
 
   lmp->init();
   update->minimize->setup();
-  update->minimize->run();
+
+  timer->barrier_start(TIME_LOOP);
+  update->minimize->run(update->nsteps);
+  timer->barrier_stop(TIME_LOOP);
+
   update->minimize->cleanup();
 
   Finish finish(lmp);
diff --git a/src/output.cpp b/src/output.cpp
index a27518a5fb..1b896a9299 100644
--- a/src/output.cpp
+++ b/src/output.cpp
@@ -135,7 +135,7 @@ void Output::setup(int flag)
   // if no dumps, set next_dump_any to last+1 so will not influence next
   // dump custom may invoke computes so wrap with clear/add
 
-  if (ndump) {
+  if (ndump && update->restrict_output == 0) {
     for (int idump = 0; idump < ndump; idump++) {
       if (strcmp(dump[idump]->style,"custom") == 0)
 	modify->clearstep_compute();
@@ -158,7 +158,7 @@ void Output::setup(int flag)
   // will not write on last step of run unless multiple of every
   // if every = 0, set next_restart to last+1 so will not influence next
 
-  if (restart_every)
+  if (restart_every && update->restrict_output == 0)
     next_restart = (ntimestep/restart_every)*restart_every + restart_every;
   else next_restart = update->laststep + 1;
 
@@ -258,6 +258,43 @@ void Output::write(int ntimestep)
   next = MYMIN(next,next_thermo);
 }
 
+/* ----------------------------------------------------------------------
+   force a snapshot to be written for all dumps
+------------------------------------------------------------------------- */
+
+void Output::write_dump(int ntimestep)
+{
+  for (int idump = 0; idump < ndump; idump++) {
+    dump[idump]->write();
+    last_dump[idump] = ntimestep;
+  }
+}
+
+/* ----------------------------------------------------------------------
+   force a restart file to be written
+------------------------------------------------------------------------- */
+
+void Output::write_restart(int ntimestep)
+{
+  if (restart_toggle == 0) {
+    char *file = new char[strlen(restart1) + 16];
+    char *ptr = strchr(restart1,'*');
+    *ptr = '\0';
+    sprintf(file,"%s%d%s",restart1,ntimestep,ptr+1);
+    *ptr = '*';
+    restart->write(file);
+    delete [] file;
+  } else if (restart_toggle == 1) {
+    restart->write(restart1);
+    restart_toggle = 2;
+  } else if (restart_toggle == 2) {
+    restart->write(restart2);
+    restart_toggle = 1;
+  }
+
+  last_restart = ntimestep;
+}
+
 /* ----------------------------------------------------------------------
    add a Dump to list of Dumps 
 ------------------------------------------------------------------------- */
diff --git a/src/output.h b/src/output.h
index c6295471d9..2dedbf634e 100644
--- a/src/output.h
+++ b/src/output.h
@@ -48,6 +48,8 @@ class Output : protected Pointers {
   void init();
   void setup(int);                   // initial output before run/min
   void write(int);                   // output for current timestep
+  void write_dump(int);              // force output of dump snapshots
+  void write_restart(int);           // force output of a restart file
 
   void add_dump(int, char **);       // add a Dump to Dump list
   void modify_dump(int, char **);    // modify a Dump
diff --git a/src/respa.cpp b/src/respa.cpp
index 1e1350fdc6..7a7fd0eba3 100644
--- a/src/respa.cpp
+++ b/src/respa.cpp
@@ -363,10 +363,69 @@ void Respa::setup()
 }
 
 /* ----------------------------------------------------------------------
-   iterate for n steps
+   setup without output
+   flag = 0 = just force calculation
+   flag = 1 = reneighbor and force calculation
 ------------------------------------------------------------------------- */
 
-void Respa::iterate(int n)
+void Respa::setup_minimal(int flag)
+{
+  // setup domain, communication and neighboring
+  // acquire ghosts
+  // build neighbor lists
+
+  if (flag) {
+    if (triclinic) domain->x2lamda(atom->nlocal);
+    domain->pbc();
+    domain->reset_box();
+    comm->setup();
+    if (neighbor->style) neighbor->setup_bins();
+    comm->exchange();
+    comm->borders();
+    if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
+    neighbor->build();
+    neighbor->ncalls = 0;
+  }
+
+  // compute all forces
+
+  ev_set(update->ntimestep);
+
+  for (int ilevel = 0; ilevel < nlevels; ilevel++) {
+    force_clear(newton[ilevel]);
+    if (level_bond == ilevel && force->bond)
+      force->bond->compute(eflag,vflag);
+    if (level_angle == ilevel && force->angle) 
+      force->angle->compute(eflag,vflag);
+    if (level_dihedral == ilevel && force->dihedral) 
+      force->dihedral->compute(eflag,vflag);
+    if (level_improper == ilevel && force->improper) 
+      force->improper->compute(eflag,vflag);
+    if (level_pair == ilevel && force->pair)
+      force->pair->compute(eflag,vflag);
+    if (level_inner == ilevel && force->pair)
+      force->pair->compute_inner();
+    if (level_middle == ilevel && force->pair)
+      force->pair->compute_middle();
+    if (level_outer == ilevel && force->pair)
+      force->pair->compute_outer(eflag,vflag);
+    if (level_kspace == ilevel && force->kspace) {
+      force->kspace->setup();
+      force->kspace->compute(eflag,vflag);
+    }
+    if (newton[ilevel]) comm->reverse_communicate();
+    copy_f_flevel(ilevel);
+  }
+  
+  modify->setup(vflag);
+  sum_flevel_f();
+}
+
+/* ----------------------------------------------------------------------
+   run for N steps
+------------------------------------------------------------------------- */
+
+void Respa::run(int n)
 {
   int ntimestep;
 
diff --git a/src/respa.h b/src/respa.h
index 52c52b00e1..07bc2f47cb 100644
--- a/src/respa.h
+++ b/src/respa.h
@@ -37,7 +37,8 @@ class Respa : public Integrate {
   ~Respa();
   void init();
   void setup();
-  void iterate(int);
+  void setup_minimal(int);
+  void run(int);
   void cleanup();
   void reset_dt();
 
diff --git a/src/run.cpp b/src/run.cpp
index bb3ddc8f12..ef38951736 100644
--- a/src/run.cpp
+++ b/src/run.cpp
@@ -148,7 +148,7 @@ void Run::command(int narg, char **arg)
     }
 
     timer->barrier_start(TIME_LOOP);
-    update->integrate->iterate(nsteps);
+    update->integrate->run(nsteps);
     timer->barrier_stop(TIME_LOOP);
 
     update->integrate->cleanup();
@@ -186,7 +186,7 @@ void Run::command(int narg, char **arg)
       }
 
       timer->barrier_start(TIME_LOOP);
-      update->integrate->iterate(nsteps);
+      update->integrate->run(nsteps);
       timer->barrier_stop(TIME_LOOP);
 
       update->integrate->cleanup();
diff --git a/src/temper.cpp b/src/temper.cpp
index 0f3991a358..fb7431b763 100644
--- a/src/temper.cpp
+++ b/src/temper.cpp
@@ -196,7 +196,7 @@ void Temper::command(int narg, char **arg)
 	
     // run for nevery timesteps
 
-    update->integrate->iterate(nevery);
+    update->integrate->run(nevery);
 	
     // compute PE
     // notify compute it will be called at next swap
diff --git a/src/update.cpp b/src/update.cpp
index 1a52ca37e0..afb43199d1 100644
--- a/src/update.cpp
+++ b/src/update.cpp
@@ -45,6 +45,8 @@ Update::Update(LAMMPS *lmp) : Pointers(lmp)
   firststep = laststep = 0;
   beginstep = endstep = 0;
 
+  restrict_output = 0;
+
   eflag_global = vflag_global = -1;
 
   unit_style = NULL;
diff --git a/src/update.h b/src/update.h
index 41d3fab4be..d0d55c5cbe 100644
--- a/src/update.h
+++ b/src/update.h
@@ -29,8 +29,9 @@ class Update : protected Pointers {
   int beginstep,endstep;          // 1st and last step of multiple runs
   int first_update;               // 0 before initial update, 1 after
   int max_eval;                   // max force evaluations for minimizer
+  int restrict_output;            // 1 if output should not write dump/restart
 
-  int eflag_global,eflag_atom;    // timestep global/peratom eng is tallied
+  int eflag_global,eflag_atom;    // timestep global/peratom eng is tallied on
   int vflag_global,vflag_atom;    // ditto for virial
 
   char *unit_style;
diff --git a/src/verlet.cpp b/src/verlet.cpp
index 92beb1ce6e..954af1a6ac 100644
--- a/src/verlet.cpp
+++ b/src/verlet.cpp
@@ -122,10 +122,59 @@ void Verlet::setup()
 }
 
 /* ----------------------------------------------------------------------
-   iterate for n steps
+   setup without output
+   flag = 0 = just force calculation
+   flag = 1 = reneighbor and force calculation
 ------------------------------------------------------------------------- */
 
-void Verlet::iterate(int n)
+void Verlet::setup_minimal(int flag)
+{
+  // setup domain, communication and neighboring
+  // acquire ghosts
+  // build neighbor lists
+
+  if (flag) {
+    if (triclinic) domain->x2lamda(atom->nlocal);
+    domain->pbc();
+    domain->reset_box();
+    comm->setup();
+    if (neighbor->style) neighbor->setup_bins();
+    comm->exchange();
+    comm->borders();
+    if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
+    neighbor->build();
+    neighbor->ncalls = 0;
+  }
+
+  // compute all forces
+
+  ev_set(update->ntimestep);
+  force_clear();
+
+  if (force->pair) force->pair->compute(eflag,vflag);
+
+  if (atom->molecular) {
+    if (force->bond) force->bond->compute(eflag,vflag);
+    if (force->angle) force->angle->compute(eflag,vflag);
+    if (force->dihedral) force->dihedral->compute(eflag,vflag);
+    if (force->improper) force->improper->compute(eflag,vflag);
+  }
+
+  if (force->kspace) {
+    force->kspace->setup();
+    force->kspace->compute(eflag,vflag);
+  }
+
+  if (force->newton) comm->reverse_communicate();
+
+  modify->setup(vflag);
+}
+
+/* ----------------------------------------------------------------------
+   run for N steps
+------------------------------------------------------------------------- */
+
+void Verlet::run(int n)
 {
   int nflag,ntimestep;
 
diff --git a/src/verlet.h b/src/verlet.h
index 41b94dcb08..fc857d16e9 100644
--- a/src/verlet.h
+++ b/src/verlet.h
@@ -24,7 +24,8 @@ class Verlet : public Integrate {
   ~Verlet() {}
   void init();
   void setup();
-  void iterate(int);
+  void setup_minimal(int);
+  void run(int);
 
  private:
   int triclinic;                    // 0 if domain is orthog, 1 if triclinic