diff --git a/doc/src/PDF/colvars-refman-lammps.pdf b/doc/src/PDF/colvars-refman-lammps.pdf
index 4afeedccc8..07d8254475 100644
Binary files a/doc/src/PDF/colvars-refman-lammps.pdf and b/doc/src/PDF/colvars-refman-lammps.pdf differ
diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp
index 22befdfd43..8b28eaa0df 100644
--- a/lib/colvars/colvar.cpp
+++ b/lib/colvars/colvar.cpp
@@ -21,10 +21,14 @@
 
 
 colvar::colvar()
-  : prev_timestep(-1)
 {
-  // Initialize static array once and for all
   runave_os = NULL;
+
+  prev_timestep = -1;
+  after_restart = false;
+  kinetic_energy = 0.0;
+  potential_energy = 0.0;
+
   init_cv_requires();
 }
 
@@ -38,6 +42,7 @@ namespace {
   }
 }
 
+
 int colvar::init(std::string const &conf)
 {
   cvm::log("Initializing a new collective variable.\n");
@@ -48,7 +53,7 @@ int colvar::init(std::string const &conf)
   colvarmodule *cv = cvm::main();
 
   get_keyval(conf, "name", this->name,
-             (std::string("colvar")+cvm::to_str(cv->variables()->size()+1)));
+             (std::string("colvar")+cvm::to_str(cv->variables()->size())));
 
   if ((cvm::colvar_by_name(this->name) != NULL) &&
       (cvm::colvar_by_name(this->name) != this)) {
@@ -63,9 +68,6 @@ int colvar::init(std::string const &conf)
 
   this->description = "colvar " + this->name;
 
-  kinetic_energy = 0.0;
-  potential_energy = 0.0;
-
   error_code |= init_components(conf);
   if (error_code != COLVARS_OK) {
     return cvm::get_error();
@@ -260,7 +262,6 @@ int colvar::init(std::string const &conf)
   f_old.reset();
 
   x_restart.type(value());
-  after_restart = false;
 
   reset_bias_force();
 
@@ -282,8 +283,7 @@ int colvar::init(std::string const &conf)
   // Now that the children are defined we can solve dependencies
   enable(f_cv_active);
 
-  if (cvm::b_analysis)
-    parse_analysis(conf);
+  error_code |= parse_analysis(conf);
 
   if (cvm::debug())
     cvm::log("Done initializing collective variable \""+this->name+"\".\n");
@@ -881,19 +881,7 @@ int colvar::parse_analysis(std::string const &conf)
       cvm::error("Error: runAveStride must be commensurate with the restart frequency.\n", INPUT_ERROR);
     }
 
-    get_keyval(conf, "runAveOutputFile", runave_outfile,
-                std::string(cvm::output_prefix()+"."+
-                             this->name+".runave.traj"));
-
-    size_t const this_cv_width = x.output_width(cvm::cv_width);
-    cvm::proxy->backup_file(runave_outfile);
-    runave_os = cvm::proxy->output_stream(runave_outfile);
-    *runave_os << "# " << cvm::wrap_string("step", cvm::it_width-2)
-               << "  "
-               << cvm::wrap_string("running average", this_cv_width)
-               << " "
-               << cvm::wrap_string("running stddev", this_cv_width)
-               << "\n";
+    get_keyval(conf, "runAveOutputFile", runave_outfile, runave_outfile);
   }
 
   acf_length = 0;
@@ -902,7 +890,6 @@ int colvar::parse_analysis(std::string const &conf)
 
     enable(f_cv_corrfunc);
 
-    std::string acf_colvar_name;
     get_keyval(conf, "corrFuncWithColvar", acf_colvar_name, this->name);
     if (acf_colvar_name == this->name) {
       cvm::log("Calculating auto-correlation function.\n");
@@ -918,8 +905,12 @@ int colvar::parse_analysis(std::string const &conf)
     } else if (acf_type_str == to_lower_cppstr(std::string("velocity"))) {
       acf_type = acf_vel;
       enable(f_cv_fdiff_velocity);
-      if (acf_colvar_name.size())
-        (cvm::colvar_by_name(acf_colvar_name))->enable(f_cv_fdiff_velocity);
+      colvar *cv2 = cvm::colvar_by_name(acf_colvar_name);
+      if (cv2 == NULL) {
+        return cvm::error("Error: collective variable \""+acf_colvar_name+
+                          "\" is not defined at this time.\n", INPUT_ERROR);
+      }
+      cv2->enable(f_cv_fdiff_velocity);
     } else if (acf_type_str == to_lower_cppstr(std::string("coordinate_p2"))) {
       acf_type = acf_p2coor;
     } else {
@@ -937,9 +928,7 @@ int colvar::parse_analysis(std::string const &conf)
     }
 
     get_keyval(conf, "corrFuncNormalize", acf_normalize, true);
-    get_keyval(conf, "corrFuncOutputFile", acf_outfile,
-                std::string(cvm::output_prefix()+"."+this->name+
-                             ".corrfunc.dat"));
+    get_keyval(conf, "corrFuncOutputFile", acf_outfile, acf_outfile);
   }
   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
 }
@@ -1389,9 +1378,9 @@ int colvar::calc_colvar_properties()
 {
   if (is_enabled(f_cv_fdiff_velocity)) {
     // calculate the velocity by finite differences
-    if (cvm::step_relative() == 0)
+    if (cvm::step_relative() == 0) {
       x_old = x;
-    else {
+    } else {
       v_fdiff = fdiff_velocity(x_old, x);
       v_reported = v_fdiff;
     }
@@ -1486,7 +1475,6 @@ cvm::real colvar::update_forces_energy()
         return 0.;
       }
     }
-    prev_timestep = cvm::step_relative();
 
     // Integrate with slow timestep (if time_step_factor != 1)
     cvm::real dt = cvm::dt() * cvm::real(time_step_factor);
@@ -1547,8 +1535,18 @@ cvm::real colvar::update_forces_energy()
   // bypass the extended Lagrangian mass)
   f += fb_actual;
 
+  if (cvm::debug())
+    cvm::log("Done updating colvar \""+this->name+"\".\n");
+  return (potential_energy + kinetic_energy);
+}
+
+
+int colvar::end_of_step()
+{
+  if (cvm::debug())
+    cvm::log("End of step for colvar \""+this->name+"\".\n");
+
   if (is_enabled(f_cv_fdiff_velocity)) {
-    // set it for the next step
     x_old = x;
   }
 
@@ -1556,9 +1554,9 @@ cvm::real colvar::update_forces_energy()
     f_old = f;
   }
 
-  if (cvm::debug())
-    cvm::log("Done updating colvar \""+this->name+"\".\n");
-  return (potential_energy + kinetic_energy);
+  prev_timestep = cvm::step_relative();
+
+  return COLVARS_OK;
 }
 
 
@@ -1966,6 +1964,10 @@ std::ostream & colvar::write_restart(std::ostream &os) {
 
   os << "}\n\n";
 
+  if (runave_os) {
+    cvm::main()->proxy->flush_output_stream(runave_os);
+  }
+
   return os;
 }
 
@@ -2075,55 +2077,61 @@ std::ostream & colvar::write_traj(std::ostream &os)
   return os;
 }
 
+
 int colvar::write_output_files()
 {
-  if (cvm::b_analysis) {
+  int error_code = COLVARS_OK;
 
+  if (is_enabled(f_cv_corrfunc)) {
     if (acf.size()) {
-      cvm::log("Writing acf to file \""+acf_outfile+"\".\n");
-
+      if (acf_outfile.size() == 0) {
+        acf_outfile = std::string(cvm::output_prefix()+"."+this->name+
+                                  ".corrfunc.dat");
+      }
+      cvm::log("Writing correlation function to file \""+acf_outfile+"\".\n");
       cvm::backup_file(acf_outfile.c_str());
       std::ostream *acf_os = cvm::proxy->output_stream(acf_outfile);
       if (!acf_os) return cvm::get_error();
-      write_acf(*acf_os);
+      error_code |= write_acf(*acf_os);
       cvm::proxy->close_output_stream(acf_outfile);
     }
-
-    if (runave_os) {
-      cvm::proxy->close_output_stream(runave_outfile);
-      runave_os = NULL;
-    }
   }
-  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
+
+  return error_code;
 }
 
 
 
 // ******************** ANALYSIS FUNCTIONS ********************
 
-void colvar::analyze()
+int colvar::analyze()
 {
+  int error_code = COLVARS_OK;
+
   if (is_enabled(f_cv_runave)) {
-    calc_runave();
+    error_code |= calc_runave();
   }
 
   if (is_enabled(f_cv_corrfunc)) {
-    calc_acf();
+    error_code |= calc_acf();
   }
+
+  return error_code;
 }
 
 
 inline void history_add_value(size_t const           &history_length,
-                               std::list<colvarvalue> &history,
-                               colvarvalue const      &new_value)
+                              std::list<colvarvalue> &history,
+                              colvarvalue const      &new_value)
 {
   history.push_front(new_value);
   if (history.size() > history_length)
     history.pop_back();
 }
 
+
 inline void history_incr(std::list< std::list<colvarvalue> >           &history,
-                          std::list< std::list<colvarvalue> >::iterator &history_p)
+                         std::list< std::list<colvarvalue> >::iterator &history_p)
 {
   if ((++history_p) == history.end())
     history_p = history.begin();
@@ -2133,18 +2141,21 @@ inline void history_incr(std::list< std::list<colvarvalue> >           &history,
 int colvar::calc_acf()
 {
   // using here an acf_stride-long list of vectors for either
-  // coordinates(acf_x_history) or velocities (acf_v_history); each vector can
+  // coordinates (acf_x_history) or velocities (acf_v_history); each vector can
   // contain up to acf_length values, which are contiguous in memory
   // representation but separated by acf_stride in the time series;
   // the pointer to each vector is changed at every step
 
+  colvar const *cfcv = cvm::colvar_by_name(acf_colvar_name);
+  if (cfcv == NULL) {
+    return cvm::error("Error: collective variable \""+acf_colvar_name+
+                      "\" is not defined at this time.\n", INPUT_ERROR);
+  }
+
   if (acf_x_history.empty() && acf_v_history.empty()) {
 
     // first-step operations
 
-    colvar *cfcv = (acf_colvar_name.size() ?
-                    cvm::colvar_by_name(acf_colvar_name) :
-                    this);
     if (colvarvalue::check_types(cfcv->value(), value())) {
       cvm::error("Error: correlation function between \""+cfcv->name+
                  "\" and \""+this->name+"\" cannot be calculated, "
@@ -2153,7 +2164,8 @@ int colvar::calc_acf()
     }
     acf_nframes = 0;
 
-    cvm::log("Colvar \""+this->name+"\": initializing ACF calculation.\n");
+    cvm::log("Colvar \""+this->name+"\": initializing correlation function "
+             "calculation.\n");
 
     if (acf.size() < acf_length+1)
       acf.resize(acf_length+1, 0.0);
@@ -2182,41 +2194,31 @@ int colvar::calc_acf()
       break;
     }
 
-  } else {
-
-    colvar *cfcv = (acf_colvar_name.size() ?
-                    cvm::colvar_by_name(acf_colvar_name) :
-                    this);
+  } else if (cvm::step_relative() > prev_timestep) {
 
     switch (acf_type) {
 
     case acf_vel:
 
-      if (is_enabled(f_cv_fdiff_velocity)) {
-        // calc() should do this already, but this only happens in a
-        // simulation; better do it again in case a trajectory is
-        // being read
-        v_reported = v_fdiff = fdiff_velocity(x_old, cfcv->value());
-      }
-
       calc_vel_acf((*acf_v_history_p), cfcv->velocity());
-      // store this value in the history
-      history_add_value(acf_length+acf_offset, *acf_v_history_p, cfcv->velocity());
-      // if stride is larger than one, cycle among different histories
+      history_add_value(acf_length+acf_offset, *acf_v_history_p,
+                        cfcv->velocity());
       history_incr(acf_v_history, acf_v_history_p);
       break;
 
     case acf_coor:
 
       calc_coor_acf((*acf_x_history_p), cfcv->value());
-      history_add_value(acf_length+acf_offset, *acf_x_history_p, cfcv->value());
+      history_add_value(acf_length+acf_offset, *acf_x_history_p,
+                        cfcv->value());
       history_incr(acf_x_history, acf_x_history_p);
       break;
 
     case acf_p2coor:
 
       calc_p2coor_acf((*acf_x_history_p), cfcv->value());
-      history_add_value(acf_length+acf_offset, *acf_x_history_p, cfcv->value());
+      history_add_value(acf_length+acf_offset, *acf_x_history_p,
+                        cfcv->value());
       history_incr(acf_x_history, acf_x_history_p);
       break;
 
@@ -2225,18 +2227,14 @@ int colvar::calc_acf()
     }
   }
 
-  if (is_enabled(f_cv_fdiff_velocity)) {
-    // set it for the next step
-    x_old = x;
-  }
-  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
+  return COLVARS_OK;
 }
 
 
-int colvar::calc_vel_acf(std::list<colvarvalue> &v_list,
-                           colvarvalue const      &v)
+void colvar::calc_vel_acf(std::list<colvarvalue> &v_list,
+                          colvarvalue const      &v)
 {
-  // loop over stored velocities and add to the ACF, but only the
+  // loop over stored velocities and add to the ACF, but only if the
   // length is sufficient to hold an entire row of ACF values
   if (v_list.size() >= acf_length+acf_offset) {
     std::list<colvarvalue>::iterator  vs_i = v_list.begin();
@@ -2255,7 +2253,6 @@ int colvar::calc_vel_acf(std::list<colvarvalue> &v_list,
 
     acf_nframes++;
   }
-  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
 }
 
 
@@ -2280,7 +2277,7 @@ void colvar::calc_coor_acf(std::list<colvarvalue> &x_list,
 
 
 void colvar::calc_p2coor_acf(std::list<colvarvalue> &x_list,
-                              colvarvalue const      &x)
+                             colvarvalue const      &x)
 {
   // same as above but with second order Legendre polynomial instead
   // of just the scalar product
@@ -2301,20 +2298,46 @@ void colvar::calc_p2coor_acf(std::list<colvarvalue> &x_list,
 }
 
 
-void colvar::write_acf(std::ostream &os)
+int colvar::write_acf(std::ostream &os)
 {
-  if (!acf_nframes)
-    cvm::log("Warning: ACF was not calculated (insufficient frames).\n");
+  if (!acf_nframes) {
+    return COLVARS_OK;
+  }
+
   os.setf(std::ios::scientific, std::ios::floatfield);
-  os << "# Autocorrelation function for collective variable \""
-     << this->name << "\"\n";
-  // one frame is used for normalization, the statistical sample is
-  // hence decreased
-  os << "# nframes = " << (acf_normalize ?
-                           acf_nframes - 1 :
-                           acf_nframes) << "\n";
+  os << "# ";
+  switch (acf_type) {
+  case acf_vel:
+    os << "Velocity";
+    break;
+  case acf_coor:
+    os << "Coordinate";
+    break;
+  case acf_p2coor:
+    os << "Coordinate (2nd Legendre poly)";
+    break;
+  }
+
+  if (acf_colvar_name == name) {
+    os << " autocorrelation function for variable \""
+       << this->name << "\"\n";
+  } else {
+    os << " correlation function between variables \"" //
+       << this->name << "\" and \"" << acf_colvar_name << "\"\n";
+  }
+
+  os << "# Number of samples = ";
+  if (acf_normalize) {
+    os << (acf_nframes-1) << " (one DoF is used for normalization)\n";
+  } else {
+    os << acf_nframes << "\n";
+  }
+
+  os << "# " << cvm::wrap_string("step", cvm::it_width-2) << " "
+     << cvm::wrap_string("corrfunc(step)", cvm::cv_width) << "\n";
 
   cvm::real const acf_norm = acf.front() / cvm::real(acf_nframes);
+
   std::vector<cvm::real>::iterator acf_i;
   size_t it = acf_offset;
   for (acf_i = acf.begin(); acf_i != acf.end(); ++acf_i) {
@@ -2325,11 +2348,15 @@ void colvar::write_acf(std::ostream &os)
             (*acf_i)/(acf_norm * cvm::real(acf_nframes)) :
             (*acf_i)/(cvm::real(acf_nframes)) ) << "\n";
   }
+
+  return os.good() ? COLVARS_OK : FILE_ERROR;
 }
 
 
-void colvar::calc_runave()
+int colvar::calc_runave()
 {
+  int error_code = COLVARS_OK;
+
   if (x_history.empty()) {
 
     runave.type(value().type());
@@ -2348,10 +2375,29 @@ void colvar::calc_runave()
 
   } else {
 
-    if ( (cvm::step_relative() % runave_stride) == 0) {
+    if ( (cvm::step_relative() % runave_stride) == 0 &&
+         (cvm::step_relative() > prev_timestep) ) {
 
       if ((*x_history_p).size() >= runave_length-1) {
 
+        if (runave_os == NULL) {
+          if (runave_outfile.size() == 0) {
+            runave_outfile = std::string(cvm::output_prefix()+"."+
+                                         this->name+".runave.traj");
+          }
+
+          size_t const this_cv_width = x.output_width(cvm::cv_width);
+          cvm::proxy->backup_file(runave_outfile);
+          runave_os = cvm::proxy->output_stream(runave_outfile);
+          runave_os->setf(std::ios::scientific, std::ios::floatfield);
+          *runave_os << "# " << cvm::wrap_string("step", cvm::it_width-2)
+                     << "   "
+                     << cvm::wrap_string("running average", this_cv_width)
+                     << " "
+                     << cvm::wrap_string("running stddev", this_cv_width)
+                     << "\n";
+        }
+
         runave = x;
         std::list<colvarvalue>::iterator xs_i;
         for (xs_i = (*x_history_p).begin();
@@ -2370,7 +2416,7 @@ void colvar::calc_runave()
         runave_variance *= 1.0 / cvm::real(runave_length-1);
 
         *runave_os << std::setw(cvm::it_width) << cvm::step_relative()
-                   << "  "
+                   << "   "
                    << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
                    << runave << " "
                    << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
@@ -2381,6 +2427,7 @@ void colvar::calc_runave()
     }
   }
 
+  return error_code;
 }
 
 // Static members
diff --git a/lib/colvars/colvar.h b/lib/colvars/colvar.h
index 989d55124f..a67749d577 100644
--- a/lib/colvars/colvar.h
+++ b/lib/colvars/colvar.h
@@ -291,6 +291,9 @@ public:
   /// \brief Calculate the colvar's value and related quantities
   int calc();
 
+  /// Carry out operations needed before next step is run
+  int end_of_step();
+
   /// \brief Calculate a subset of the colvar components (CVCs) currently active
   /// (default: all active CVCs)
   /// Note: both arguments refer to the sect of *active* CVCs, not all CVCs
@@ -410,8 +413,9 @@ public:
 
   /// Read the analysis tasks
   int parse_analysis(std::string const &conf);
+
   /// Perform analysis tasks
-  void analyze();
+  int analyze();
 
 
   /// Read the value from a collective variable trajectory file
@@ -489,23 +493,23 @@ protected:
   acf_type_e             acf_type;
 
   /// \brief Velocity ACF, scalar product between v(0) and v(t)
-  int calc_vel_acf(std::list<colvarvalue> &v_history,
-                     colvarvalue const      &v);
+  void calc_vel_acf(std::list<colvarvalue> &v_history,
+                    colvarvalue const      &v);
 
   /// \brief Coordinate ACF, scalar product between x(0) and x(t)
   /// (does not work with scalar numbers)
   void calc_coor_acf(std::list<colvarvalue> &x_history,
-                      colvarvalue const      &x);
+                     colvarvalue const      &x);
 
   /// \brief Coordinate ACF, second order Legendre polynomial between
   /// x(0) and x(t) (does not work with scalar numbers)
   void calc_p2coor_acf(std::list<colvarvalue> &x_history,
-                        colvarvalue const      &x);
+                       colvarvalue const      &x);
 
   /// Calculate the auto-correlation function (ACF)
   int calc_acf();
   /// Save the ACF to a file
-  void write_acf(std::ostream &os);
+  int write_acf(std::ostream &os);
 
   /// Length of running average series
   size_t         runave_length;
@@ -521,7 +525,7 @@ protected:
   cvm::real      runave_variance;
 
   /// Calculate the running average and its standard deviation
-  void calc_runave();
+  int calc_runave();
 
   /// If extended Lagrangian active: colvar energies (kinetic and harmonic potential)
   cvm::real kinetic_energy;
diff --git a/lib/colvars/colvarbias.cpp b/lib/colvars/colvarbias.cpp
index 29620fbee8..9363fcdcb6 100644
--- a/lib/colvars/colvarbias.cpp
+++ b/lib/colvars/colvarbias.cpp
@@ -236,6 +236,12 @@ void colvarbias::communicate_forces()
 }
 
 
+int colvarbias::end_of_step()
+{
+  return COLVARS_OK;
+}
+
+
 int colvarbias::change_configuration(std::string const &conf)
 {
   cvm::error("Error: change_configuration() not implemented.\n",
diff --git a/lib/colvars/colvarbias.h b/lib/colvars/colvarbias.h
index 083b9d7303..391826e79e 100644
--- a/lib/colvars/colvarbias.h
+++ b/lib/colvars/colvarbias.h
@@ -66,6 +66,9 @@ public:
   /// Send forces to the collective variables
   virtual void communicate_forces();
 
+  /// Carry out operations needed before next step is run
+  virtual int end_of_step();
+
   /// Load new configuration - force constant and/or centers only
   virtual int change_configuration(std::string const &conf);
 
diff --git a/lib/colvars/colvardeps.cpp b/lib/colvars/colvardeps.cpp
index 80bd6670d3..d20ee6e55c 100644
--- a/lib/colvars/colvardeps.cpp
+++ b/lib/colvars/colvardeps.cpp
@@ -567,6 +567,9 @@ void colvardeps::init_cv_requires() {
     // Most features are available, so we set them so
     // and list exceptions below
    }
+
+  feature_states[f_cv_fdiff_velocity].available =
+    cvm::main()->proxy->simulation_running();
 }
 
 
diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp
index 87b08b1ef8..d88a97a441 100644
--- a/lib/colvars/colvarmodule.cpp
+++ b/lib/colvars/colvarmodule.cpp
@@ -62,8 +62,6 @@ colvarmodule::colvarmodule(colvarproxy *proxy_in)
   use_scripted_forces = false;
   scripting_after_biases = false;
 
-  b_analysis = false;
-
   colvarmodule::debug_gradients_step_size = 1.0e-07;
 
   colvarmodule::rotation::monitor_crossings = false;
@@ -274,7 +272,12 @@ int colvarmodule::parse_global_params(std::string const &conf)
     }
   }
 
-  parse->get_keyval(conf, "analysis", b_analysis, b_analysis);
+  bool b_analysis = true;
+  if (parse->get_keyval(conf, "analysis", b_analysis, true,
+                        colvarparse::parse_silent)) {
+    cvm::log("Warning: keyword \"analysis\" is deprecated: it is now set "
+             "to true; individual analyses are performed only if requested.");
+  }
 
   parse->get_keyval(conf, "debugGradientsStepSize", debug_gradients_step_size,
                     debug_gradients_step_size,
@@ -715,9 +718,7 @@ int colvarmodule::calc()
   error_code |= calc_biases();
   error_code |= update_colvar_forces();
 
-  if (cvm::b_analysis) {
-    error_code |= analyze();
-  }
+  error_code |= analyze();
 
   // write trajectory files, if needed
   if (cv_traj_freq && cv_traj_name.size()) {
@@ -736,6 +737,8 @@ int colvarmodule::calc()
     write_output_files();
   }
 
+  error_code |= end_of_step();
+
   return error_code;
 }
 
@@ -1056,6 +1059,33 @@ int colvarmodule::analyze()
 }
 
 
+int colvarmodule::end_of_step()
+{
+  if (cvm::debug()) {
+    cvm::log("colvarmodule::end_of_step(), step = "+cvm::to_str(it)+".\n");
+  }
+
+  for (std::vector<colvar *>::iterator cvi = variables_active()->begin();
+       cvi != variables_active()->end();
+       cvi++) {
+    cvm::increase_depth();
+    (*cvi)->end_of_step();
+    cvm::decrease_depth();
+  }
+
+  // perform bias-specific analysis
+  for (std::vector<colvarbias *>::iterator bi = biases.begin();
+       bi != biases.end();
+       bi++) {
+    cvm::increase_depth();
+    (*bi)->end_of_step();
+    cvm::decrease_depth();
+  }
+
+  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
+}
+
+
 int colvarmodule::setup()
 {
   if (this->size() == 0) return cvm::get_error();
@@ -1895,7 +1925,6 @@ long      colvarmodule::it = 0;
 long      colvarmodule::it_restart = 0;
 size_t    colvarmodule::restart_out_freq = 0;
 size_t    colvarmodule::cv_traj_freq = 0;
-bool      colvarmodule::b_analysis = false;
 bool      colvarmodule::use_scripted_forces = false;
 bool      colvarmodule::scripting_after_biases = true;
 
diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h
index 3d93798e0a..99b797627e 100644
--- a/lib/colvars/colvarmodule.h
+++ b/lib/colvars/colvarmodule.h
@@ -422,6 +422,9 @@ public:
   /// Perform analysis
   int analyze();
 
+  /// Carry out operations needed before next step is run
+  int end_of_step();
+
   /// \brief Read a collective variable trajectory (post-processing
   /// only, not called at runtime)
   int read_traj(char const *traj_filename,
@@ -546,9 +549,6 @@ public:
   /// Frequency for collective variables trajectory output
   static size_t cv_traj_freq;
 
-  /// \brief True if only analysis is performed and not a run
-  static bool   b_analysis;
-
   /// Frequency for saving output restarts
   static size_t restart_out_freq;
   /// Output restart file name
diff --git a/lib/colvars/colvars_version.h b/lib/colvars/colvars_version.h
index a37d264a20..bd84d077d7 100644
--- a/lib/colvars/colvars_version.h
+++ b/lib/colvars/colvars_version.h
@@ -1,5 +1,5 @@
 #ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2018-10-16"
+#define COLVARS_VERSION "2018-11-16"
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars