diff --git a/src/input.cpp b/src/input.cpp
index fb1dd7db7f..399499a5c2 100644
--- a/src/input.cpp
+++ b/src/input.cpp
@@ -1350,10 +1350,7 @@ void Input::suffix()
 
 void Input::thermo()
 {
-  if (narg != 1) error->all(FLERR,"Illegal thermo command");
-  int n = atoi(arg[0]);
-  if (n < 0) error->all(FLERR,"Illegal thermo command");
-  output->thermo_every = n;
+  output->set_thermo(narg,arg);
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/min.cpp b/src/min.cpp
index bb084e7f1c..e5c2cb46af 100644
--- a/src/min.cpp
+++ b/src/min.cpp
@@ -278,7 +278,7 @@ void Min::setup()
       requestor[m]->min_xf_get(m);
 
   modify->setup(vflag);
-  output->setup(1,1);
+  output->setup();
 
   // stats for Finish to print
 
diff --git a/src/output.cpp b/src/output.cpp
index 3a0f832f82..b7a70222fe 100644
--- a/src/output.cpp
+++ b/src/output.cpp
@@ -122,8 +122,7 @@ Output::~Output()
 void Output::init()
 {
   thermo->init();
-  if (thermo_every) delete [] var_thermo;
-  else if (var_thermo) {
+  if (var_thermo) {
     ivar_thermo = input->variable->find(var_thermo);
     if (ivar_thermo < 0)
       error->all(FLERR,"Variable name for thermo every does not exist");
@@ -161,31 +160,30 @@ void Output::init()
    perform output for setup of run/min
    do dump first, so memory_usage will include dump allocation
    do thermo last, so will print after memory_usage
-   headflag = 0/1 for printing out header line
    memflag = 0/1 for printing out memory usage
 ------------------------------------------------------------------------- */
 
-void Output::setup(int headflag, int memflag)
+void Output::setup(int memflag)
 {
   bigint ntimestep = update->ntimestep;
 
-  // perform dump at start of run if current timestep is multiple of every
-  //   and last dump was not on this timestep
-  // perform dump on first step ever even if not multiple of every
-  //   if firstflag is set
-  // do not write on last step of run unless multiple of every
-  // NOTE: what if dump freq is set by variable, then what happens?
+  // perform dump at start of run only if:
+  //   current timestep is multiple of every and last dump not on this step
+  //   this is first run after dump created and firstflag is set
+  //   note that variable freq will not write unless triggered by firstflag
   // set next_dump to multiple of every or variable value
   // set next_dump_any to smallest next_dump
+  // wrap dumps that invoke computes and variable eval with clear/add
+  // if dump not written now, use addstep_compute_all() since don't know
+  //   what computes the dump write would invoke
   // if no dumps, set next_dump_any to last+1 so will not influence next
-  // wrap dumps that invoke computes with clear/add
-  // if dump not written now, add_all of next step since clear/add is no-op
 
   int writeflag;
 
   if (ndump && update->restrict_output == 0) {
     for (int idump = 0; idump < ndump; idump++) {
-      if (dump[idump]->clearstep) modify->clearstep_compute();
+      if (dump[idump]->clearstep || every_dump[idump] == 0)
+        modify->clearstep_compute();
       writeflag = 0;
       if (every_dump[idump] && ntimestep % every_dump[idump] == 0 &&
           last_dump[idump] != ntimestep) writeflag = 1;
@@ -205,7 +203,7 @@ void Output::setup(int headflag, int memflag)
           error->all(FLERR,"Dump every variable returned a bad timestep");
         next_dump[idump] = nextdump;
       }
-      if (dump[idump]->clearstep) {
+      if (dump[idump]->clearstep || every_dump[idump] == 0) {
         if (writeflag) modify->addstep_compute(next_dump[idump]);
         else modify->addstep_compute_all(next_dump[idump]);
       }
@@ -216,8 +214,8 @@ void Output::setup(int headflag, int memflag)
 
   // do not write restart files at start of run
   // set next_restart values to multiple of every or variable value
-  // do 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
+  // wrap variable eval with clear/add
+  // if no restarts, set next_restart to last+1 so will not influence next
 
   if (restart_flag && update->restrict_output == 0) {
     if (restart_flag_single) {
@@ -253,25 +251,24 @@ void Output::setup(int headflag, int memflag)
 
   if (memflag) memory_usage();
 
-  // always do thermo with header at start of run unless called via rerun
-  // set next_thermo to multiple of every or last step of run (if smaller)
-  // if every = 0, set next_thermo to last step of run
+  // set next_thermo to multiple of every or variable eval if var defined
+  // insure thermo output on last step of run
   // thermo may invoke computes so wrap with clear/add
 
   modify->clearstep_compute();
 
-  if (headflag) thermo->header();
+  thermo->header();
   thermo->compute(0);
   last_thermo = ntimestep;
 
-  if (thermo_every) {
-    next_thermo = (ntimestep/thermo_every)*thermo_every + thermo_every;
-    next_thermo = MIN(next_thermo,update->laststep);
-  } else if (var_thermo) {
+  if (var_thermo) {
     next_thermo = static_cast<bigint>
       (input->variable->compute_equal(ivar_thermo));
     if (next_thermo <= ntimestep)
       error->all(FLERR,"Thermo every variable returned a bad timestep");
+  } else if (thermo_every) {
+    next_thermo = (ntimestep/thermo_every)*thermo_every + thermo_every;
+    next_thermo = MIN(next_thermo,update->laststep);
   } else next_thermo = update->laststep;
 
   modify->addstep_compute(next_thermo);
@@ -284,14 +281,14 @@ void Output::setup(int headflag, int memflag)
 
 /* ----------------------------------------------------------------------
    perform all output for this timestep
-   only perform output if next matches current step and last doesn't
+   only perform output if next matches current step and last output doesn't
    do dump/restart before thermo so thermo CPU time will include them
 ------------------------------------------------------------------------- */
 
 void Output::write(bigint ntimestep)
 {
   // next_dump does not force output on last step of run
-  // wrap dumps that invoke computes with clear/add
+  // wrap dumps that invoke computes or eval of variable with clear/add
   // download data from GPU if necessary
 
   if (next_dump_any == ntimestep) {
@@ -299,7 +296,8 @@ void Output::write(bigint ntimestep)
 
     for (int idump = 0; idump < ndump; idump++) {
       if (next_dump[idump] == ntimestep && last_dump[idump] != ntimestep) {
-        if (dump[idump]->clearstep) modify->clearstep_compute();
+        if (dump[idump]->clearstep || every_dump[idump] == 0)
+          modify->clearstep_compute();
         dump[idump]->write();
         last_dump[idump] = ntimestep;
         if (every_dump[idump]) next_dump[idump] += every_dump[idump];
@@ -310,7 +308,8 @@ void Output::write(bigint ntimestep)
             error->all(FLERR,"Dump every variable returned a bad timestep");
           next_dump[idump] = nextdump;
         }
-        if (dump[idump]->clearstep) modify->addstep_compute(next_dump[idump]);
+        if (dump[idump]->clearstep || every_dump[idump] == 0)
+          modify->addstep_compute(next_dump[idump]);
       }
       if (idump) next_dump_any = MIN(next_dump_any,next_dump[idump]);
       else next_dump_any = next_dump[0];
@@ -320,6 +319,7 @@ void Output::write(bigint ntimestep)
   // next_restart does not force output on last step of run
   // for toggle = 0, replace "*" with current timestep in restart filename
   // download data from GPU if necessary
+  // eval of variable may invoke computes so wrap with clear/add
 
   if (next_restart == ntimestep && last_restart != ntimestep) {
     if (lmp->cuda && !lmp->cuda->oncpu) lmp->cuda->downloadAll();
@@ -334,11 +334,13 @@ void Output::write(bigint ntimestep)
       delete [] file;
       if (restart_every_single) next_restart_single += restart_every_single;
       else {
+        modify->clearstep_compute();
         bigint nextrestart = static_cast<bigint>
           (input->variable->compute_equal(ivar_restart_single));
         if (nextrestart <= ntimestep)
           error->all(FLERR,"Restart variable returned a bad timestep");
         next_restart_single = nextrestart;
+        modify->addstep_compute(next_restart_single);
       }
     }
     if (next_restart_double == ntimestep) {
@@ -351,11 +353,13 @@ void Output::write(bigint ntimestep)
       }
       if (restart_every_double) next_restart_double += restart_every_double;
       else {
+        modify->clearstep_compute();
         bigint nextrestart = static_cast<bigint>
           (input->variable->compute_equal(ivar_restart_double));
         if (nextrestart <= ntimestep)
           error->all(FLERR,"Restart variable returned a bad timestep");
         next_restart_double = nextrestart;
+        modify->addstep_compute(next_restart_double);
       }
     }
     last_restart = ntimestep;
@@ -369,13 +373,13 @@ void Output::write(bigint ntimestep)
     modify->clearstep_compute();
     thermo->compute(1);
     last_thermo = ntimestep;
-    if (thermo_every) next_thermo += thermo_every;
-    else if (var_thermo) {
+    if (var_thermo) {
       next_thermo = static_cast<bigint>
         (input->variable->compute_equal(ivar_thermo));
       if (next_thermo <= ntimestep)
         error->all(FLERR,"Thermo every variable returned a bad timestep");
-    } else next_thermo = update->laststep;
+    } else if (thermo_every) next_thermo += thermo_every;
+    else next_thermo = update->laststep;
     next_thermo = MIN(next_thermo,update->laststep);
     modify->addstep_compute(next_thermo);
   }
@@ -429,28 +433,31 @@ void Output::write_restart(bigint ntimestep)
 }
 
 /* ----------------------------------------------------------------------
-   timestep is being changed
+   timestep is being changed, called by update->reset_timestep()
    reset next timestep values for dumps, restart, thermo output
    reset to smallest value >= new timestep
-   called by update->reset_timestep()
+   if next timestep set by varaible evaluation,
+     eval for ntimestep-1, so current ntimestep can be returned
+     no guarantee that variable can be evaluated for ntimestep-1
+       if it depends on computes, but live with that rare case for now
 ------------------------------------------------------------------------- */
 
 void Output::reset_timestep(bigint ntimestep)
 {
-  // worry about clear/add of varialbe?
-
   for (int idump = 0; idump < ndump; idump++) {
     if (every_dump[idump]) {
       next_dump[idump] = (ntimestep/every_dump[idump])*every_dump[idump];
       if (next_dump[idump] < ntimestep) next_dump[idump] += every_dump[idump];
     } else {
-      /*
-        bigint nextdump = static_cast<bigint>
+      modify->clearstep_compute();
+      update->ntimestep--;
+      bigint nextdump = static_cast<bigint>
         (input->variable->compute_equal(ivar_dump[idump]));
-        if (nextdump <= ntimestep)
+      if (nextdump <= ntimestep)
         error->all(FLERR,"Dump every variable returned a bad timestep");
-        next_dump[idump] = nextdump;
-      */
+      update->ntimestep++;
+      next_dump[idump] = nextdump;
+      modify->addstep_compute(next_dump[idump]);
     }
     if (idump) next_dump_any = MIN(next_dump_any,next_dump[idump]);
     else next_dump_any = next_dump[0];
@@ -463,13 +470,15 @@ void Output::reset_timestep(bigint ntimestep)
       if (next_restart_single < ntimestep) 
         next_restart_single += restart_every_single;
     } else {
-      /*
-        bigint nextrestart = static_cast<bigint>
+      modify->clearstep_compute();
+      update->ntimestep--;
+      bigint nextrestart = static_cast<bigint>
         (input->variable->compute_equal(ivar_restart_single));
-        if (nextrestart <= ntimestep)
+      if (nextrestart <= ntimestep)
         error->all(FLERR,"Restart variable returned a bad timestep");
-        next_restart_single = nextrestart;
-      */
+      update->ntimestep++;
+      next_restart_single = nextrestart;
+      modify->addstep_compute(next_restart_single);
     }
   } else next_restart_single = update->laststep + 1;
 
@@ -480,22 +489,35 @@ void Output::reset_timestep(bigint ntimestep)
       if (next_restart_double < ntimestep) 
         next_restart_double += restart_every_double;
     } else {
-    /*
+      modify->clearstep_compute();
+      update->ntimestep--;
       bigint nextrestart = static_cast<bigint>
         (input->variable->compute_equal(ivar_restart_double));
-        if (nextrestart <= ntimestep)
+      if (nextrestart <= ntimestep)
         error->all(FLERR,"Restart variable returned a bad timestep");
-        next_restart_double = nextrestart;
-      */
+      update->ntimestep++;
+      next_restart_double = nextrestart;
+      modify->addstep_compute(next_restart_double);
     }
   } else next_restart_double = update->laststep + 1;
 
   next_restart = MIN(next_restart_single,next_restart_double);
-
-  next_thermo = (ntimestep/thermo_every)*thermo_every;
-  if (next_thermo < ntimestep) next_thermo += thermo_every;
-
-  // worry about thermo output on last step?
+  
+  if (var_thermo) {
+    modify->clearstep_compute();
+    update->ntimestep--;
+    next_thermo = static_cast<bigint>
+      (input->variable->compute_equal(ivar_thermo));
+    if (next_thermo < ntimestep)
+      error->all(FLERR,"Thermo every variable returned a bad timestep");
+    update->ntimestep++;
+    next_thermo = MIN(next_thermo,update->laststep);
+    modify->addstep_compute(next_thermo);
+  } else {
+    next_thermo = (ntimestep/thermo_every)*thermo_every;
+    if (next_thermo < ntimestep) next_thermo += thermo_every;
+    next_thermo = MIN(next_thermo,update->laststep);
+  }
 
   next = MIN(next_dump_any,next_restart);
   next = MIN(next,next_thermo);
@@ -598,6 +620,25 @@ void Output::delete_dump(char *id)
   ndump--;
 }
 
+/* ----------------------------------------------------------------------
+   set thermo output frequency from input script
+------------------------------------------------------------------------- */
+
+void Output::set_thermo(int narg, char **arg)
+{
+  if (narg != 1) error->all(FLERR,"Illegal thermo command");
+
+  if (strstr(arg[0],"v_") == arg[0]) {
+    delete [] var_thermo;
+    int n = strlen(&arg[0][2]) + 1;
+    var_thermo = new char[n];
+    strcpy(var_thermo,&arg[0][2]);
+  } else {
+    thermo_every = atoi(arg[0]);
+    if (thermo_every < 0) error->all(FLERR,"Illegal thermo command");
+  }
+}
+
 /* ----------------------------------------------------------------------
    new Thermo style
 ------------------------------------------------------------------------- */
diff --git a/src/output.h b/src/output.h
index cfcf3438ef..dd06b41e42 100644
--- a/src/output.h
+++ b/src/output.h
@@ -23,16 +23,16 @@ class Output : protected Pointers {
   bigint next;                 // next timestep for any kind of output
 
   bigint next_thermo;          // next timestep for thermo output
-  int thermo_every;            // thermo output every this many steps
+  int thermo_every;            // output freq for thermo, 0 if first/last only
   bigint last_thermo;          // last timestep thermo was output
-  char *var_thermo;            // variable name for thermo frequency
+  char *var_thermo;            // variable name for thermo freq, NULL if every
   int ivar_thermo;             // variable index for thermo frequency
   class Thermo *thermo;        // Thermodynamic computations
 
   int ndump;                   // # of Dumps defined
   int max_dump;                // max size of Dump list
   bigint next_dump_any;        // next timestep for any Dump
-  int *every_dump;             // output of each Dump every this many steps
+  int *every_dump;             // write freq for each Dump, 0 if var
   bigint *next_dump;           // next timestep to do each Dump
   bigint *last_dump;           // last timestep each snapshot was output
   char **var_dump;             // variable name for dump frequency
@@ -60,7 +60,7 @@ class Output : protected Pointers {
   Output(class LAMMPS *);
   ~Output();
   void init();
-  void setup(int,int);               // initial output before run/min
+  void setup(int memflag = 1);       // initial output before run/min
   void write(bigint);                // output for current timestep
   void write_dump(bigint);           // force output of dump snapshots
   void write_restart(bigint);        // force output of a restart file
@@ -70,6 +70,7 @@ class Output : protected Pointers {
   void modify_dump(int, char **);    // modify a Dump
   void delete_dump(char *);          // delete a Dump from Dump list
 
+  void set_thermo(int, char **);     // set thermo output freqquency
   void create_thermo(int, char **);  // create a thermo style
   void create_restart(int, char **); // create Restart and restart files
 
diff --git a/src/rerun.cpp b/src/rerun.cpp
index 3c3f744d2d..6cd5b2e5b4 100644
--- a/src/rerun.cpp
+++ b/src/rerun.cpp
@@ -123,7 +123,6 @@ void Rerun::command(int narg, char **arg)
   // read all relevant snapshots
   // uset setup_minimal() since atoms are already owned by correct procs
   // addstep_compute_all() insures energy/virial computed on every snapshot
-  // set update->nsteps to ndump for Finish stats to print
 
   update->whichflag = 1;
 
@@ -153,7 +152,7 @@ void Rerun::command(int narg, char **arg)
     modify->addstep_compute_all(ntimestep);
     update->integrate->setup_minimal(1);
     modify->end_of_step();
-    if (firstflag) output->setup(firstflag,firstflag);
+    if (firstflag) output->setup();
     else if (output->next) output->write(ntimestep);
 
     firstflag = 0;
@@ -161,10 +160,17 @@ void Rerun::command(int narg, char **arg)
     if (ntimestep < 0) break;
   }
 
+  // insure thermo output on last dump timestep
+
+  output->next_thermo = ntimestep;
+  output->write(ntimestep);
+
   timer->barrier_stop(TIME_LOOP);
 
   update->integrate->cleanup();
 
+  // set update->nsteps to ndump for Finish stats to print
+
   update->nsteps = ndump;
 
   Finish finish(lmp);
diff --git a/src/respa.cpp b/src/respa.cpp
index 192cb2696d..ba3e575933 100644
--- a/src/respa.cpp
+++ b/src/respa.cpp
@@ -379,7 +379,7 @@ void Respa::setup()
 
   modify->setup(vflag);
   sum_flevel_f();
-  output->setup(1,1);
+  output->setup();
   update->setupflag = 0;
 }
 
diff --git a/src/run.cpp b/src/run.cpp
index 87368b3a24..781001ae9f 100644
--- a/src/run.cpp
+++ b/src/run.cpp
@@ -167,7 +167,7 @@ void Run::command(int narg, char **arg)
     if (preflag || update->first_update == 0) {
       lmp->init();
       update->integrate->setup();
-    } else output->setup(1,0);
+    } else output->setup(0);
 
     timer->init();
     timer->barrier_start(TIME_LOOP);
@@ -205,7 +205,7 @@ void Run::command(int narg, char **arg)
       if (preflag || iter == 0) {
         lmp->init();
         update->integrate->setup();
-      } else output->setup(1,0);
+      } else output->setup(0);
 
       timer->init();
       timer->barrier_start(TIME_LOOP);
diff --git a/src/thermo.cpp b/src/thermo.cpp
index 1a92c7e8ee..52d87e19ad 100644
--- a/src/thermo.cpp
+++ b/src/thermo.cpp
@@ -412,17 +412,7 @@ void Thermo::modify_params(int narg, char **arg)
 
   int iarg = 0;
   while (iarg < narg) {
-    if (strcmp(arg[iarg],"every") == 0) {
-      if (iarg+2 > narg) error->all(FLERR,"Illegal thermo_modify command");
-      if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
-        delete [] output->var_thermo;
-        int n = strlen(&arg[iarg+1][2]) + 1;
-        output->var_thermo = new char[n];
-        strcpy(output->var_thermo,&arg[iarg+1][2]);
-      } else error->all(FLERR,"Illegal thermo_modify command");
-      output->thermo_every = 0;
-      iarg += 2;
-    } else if (strcmp(arg[iarg],"temp") == 0) {
+    if (strcmp(arg[iarg],"temp") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal thermo_modify command");
       if (index_temp < 0) error->all(FLERR,"Thermo style does not use temp");
       delete [] id_compute[index_temp];
diff --git a/src/update.cpp b/src/update.cpp
index 8f6d3f2129..1eda53689e 100644
--- a/src/update.cpp
+++ b/src/update.cpp
@@ -362,6 +362,7 @@ void Update::reset_timestep(int narg, char **arg)
    reset invoked flags of computes,
      so nothing will think they are current between runs
    clear timestep list of computes that store future invocation times
+   called from input script and rerun command
 ------------------------------------------------------------------------- */
 
 void Update::reset_timestep(bigint newstep)
@@ -376,7 +377,7 @@ void Update::reset_timestep(bigint newstep)
     if (modify->fix[i]->time_depend)
       error->all(FLERR,
                  "Cannot reset timestep with a time-dependent fix defined");
-    modify->fix[i]->reset_timestep(newstep);
+    modify->fix[i]->reset_timestep(ntimestep);
   }
 
   eflag_global = vflag_global = -1;
diff --git a/src/verlet.cpp b/src/verlet.cpp
index dac1aab7dd..095dccbd6a 100644
--- a/src/verlet.cpp
+++ b/src/verlet.cpp
@@ -139,7 +139,7 @@ void Verlet::setup()
   if (force->newton) comm->reverse_comm();
 
   modify->setup(vflag);
-  output->setup(1,1);
+  output->setup();
   update->setupflag = 0;
 }