git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@3117 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/finish.cpp

@@ -38,29 +38,51 @@ 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
   
-  double time_other = timer->array[TIME_LOOP] -
+  loopflag = 1;
+  minflag = prdflag = timeflag = fftflag = histoflag = neighflag = 0;
+
+  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;
+    
+  }
+
+  // loop stats
+
+  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]);
     
-  double time_loop = timer->array[TIME_LOOP];
+    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;
+    natoms;
     double rlocal = atom->nlocal;
     MPI_Allreduce(&rlocal,&natoms,1,MPI_DOUBLE,MPI_SUM,world);
     
@@ -75,16 +97,17 @@ void Finish::end(int flag)
 		time_loop,nprocs,update->nsteps,natoms);
     }
 
-  if (flag == 0) return;
+    if (time_loop == 0.0) time_loop = 1.0;
+  }
 
+  // minimization stats
+
+  if (minflag) {
     if (me == 0) {
       if (screen) fprintf(screen,"\n");
       if (logfile) fprintf(logfile,"\n");
     }
 
-  // minimization stats
-
-  if (update->whichflag == 2) {
     if (me == 0) {
       if (screen) {
 	fprintf(screen,"Minimization stats:\n");
@@ -127,24 +150,87 @@ 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,"Pair  time (%%) = %g (%g)\n",time,time/time_loop*100.0);
+	fprintf(screen,"  Dephase  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);
+	fprintf(logfile,"  Dephase  time (%%) = %g (%g)\n",
+		time,time/time_loop*100.0);
+    }
+
+    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,"  Dynamics time (%%) = %g (%g)\n",
+		time,time/time_loop*100.0);
+      if (logfile) 
+	fprintf(logfile,"  Dynamics time (%%) = %g (%g)\n",
+		time,time/time_loop*100.0);
+    }
+
+    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,"  Quench   time (%%) = %g (%g)\n",
+		time,time/time_loop*100.0);
+      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);
+    }
+  }
+
+  // timing breakdowns
+
+  if (timeflag) {
+    if (me == 0) {
+      if (screen) fprintf(screen,"\n");
+      if (logfile) fprintf(logfile,"\n");
+    }
+
+    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) {
@@ -180,9 +266,11 @@ void Finish::end(int flag)
     time = tmp/nprocs;
     if (me == 0) {
       if (screen) 
-      fprintf(screen,"Neigh time (%%) = %g (%g)\n",time,time/time_loop*100.0);
+	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);
+	fprintf(logfile,"Neigh time (%%) = %g (%g)\n",
+		time,time/time_loop*100.0);
     }
     
     time = timer->array[TIME_COMM];
@@ -190,9 +278,11 @@ void Finish::end(int flag)
     time = tmp/nprocs;
     if (me == 0) {
       if (screen) 
-      fprintf(screen,"Comm  time (%%) = %g (%g)\n",time,time/time_loop*100.0);
+	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);
+	fprintf(logfile,"Comm  time (%%) = %g (%g)\n",
+		time,time/time_loop*100.0);
     }
     
     time = timer->array[TIME_OUTPUT];
@@ -200,9 +290,11 @@ void Finish::end(int flag)
     time = tmp/nprocs;
     if (me == 0) {
       if (screen) 
-      fprintf(screen,"Outpt time (%%) = %g (%g)\n",time,time/time_loop*100.0);
+	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);
+	fprintf(logfile,"Outpt time (%%) = %g (%g)\n",
+		time,time/time_loop*100.0);
     }
     
     time = time_other;
@@ -210,15 +302,18 @@ void Finish::end(int flag)
     time = tmp/nprocs;
     if (me == 0) {
       if (screen) 
-      fprintf(screen,"Other time (%%) = %g (%g)\n",time,time/time_loop*100.0);
+	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);
+	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,6 +360,7 @@ void Finish::end(int flag)
     }
   }
 
+  if (histoflag) {
     if (me == 0) {
       if (screen) fprintf(screen,"\n");
       if (logfile) fprintf(logfile,"\n");
@@ -308,12 +404,13 @@ void Finish::end(int flag)
     // 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 ||
+      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;
+    nneigh = 0;
     if (m < neighbor->old_nrequest) {
       int inum = neighbor->lists[m]->inum;
       int *ilist = neighbor->lists[m]->ilist;
@@ -346,7 +443,6 @@ void Finish::end(int flag)
 	  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];
@@ -368,7 +464,9 @@ void Finish::end(int flag)
 	}
       }
     }
+  }
 
+  if (neighflag) {
     if (me == 0) {
       if (screen) fprintf(screen,"\n");
       if (logfile) fprintf(logfile,"\n");
@@ -390,7 +488,8 @@ void Finish::end(int flag)
     if (me == 0) {
       if (screen) {
 	if (nall < 2.0e9) 
-	fprintf(screen,"Total # of neighbors = %d\n",static_cast<int> (nall));
+	  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) 
@@ -400,7 +499,8 @@ void Finish::end(int flag)
       }
       if (logfile) {
 	if (nall < 2.0e9) 
-	fprintf(logfile,"Total # of neighbors = %d\n",static_cast<int> (nall));
+	  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) 
@@ -409,6 +509,7 @@ void Finish::end(int flag)
 	fprintf(logfile,"Dangerous builds = %d\n",neighbor->ndanger);
       }
     }
+  }
   
   if (logfile) fflush(logfile);
 }

src/fix_minimize.cpp

@@ -181,7 +181,8 @@ 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,
+    vectors[m] = (double *) 
+      memory->srealloc(vectors[m],peratom[m]*nmax*sizeof(double),
 		       "minimize:vector");
 }
 

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

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;}

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(nsteps);
-  int stop_condition = iterate(update->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;
 
+    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) output->next_restart = 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

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;}

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);

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 
 ------------------------------------------------------------------------- */

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

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;
 

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();
 

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();

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

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;

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;

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;
 

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