diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp
index eec49d1fad..344a9a1301 100644
--- a/lib/colvars/colvar.cpp
+++ b/lib/colvars/colvar.cpp
@@ -1187,40 +1187,6 @@ cvm::real colvar::update_forces_energy()
       f -= fj;
   }
 
-  if (is_enabled(f_cv_extended_Lagrangian)) {
-
-    cvm::real dt = cvm::dt();
-    cvm::real f_ext;
-
-    // the total force is applied to the fictitious mass, while the
-    // atoms only feel the harmonic force
-    // fr: bias force on extended coordinate (without harmonic spring), for output in trajectory
-    // f_ext: total force on extended coordinate (including harmonic spring)
-    // f: - initially, external biasing force
-    //    - after this code block, colvar force to be applied to atomic coordinates, ie. spring force
-    //      (note: wall potential is added to f after this block)
-    fr    = f;
-    f_ext = f + (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x);
-    f     =     (-0.5 * ext_force_k) * this->dist2_rgrad(xr, x);
-
-    // leapfrog: starting from x_i, f_i, v_(i-1/2)
-    vr  += (0.5 * dt) * f_ext / ext_mass;
-    // Because of leapfrog, kinetic energy at time i is approximate
-    kinetic_energy = 0.5 * ext_mass * vr * vr;
-    potential_energy = 0.5 * ext_force_k * this->dist2(xr, x);
-    // leap to v_(i+1/2)
-    if (is_enabled(f_cv_Langevin)) {
-      vr -= dt * ext_gamma * vr.real_value;
-      vr += dt * ext_sigma * cvm::rand_gaussian() / ext_mass;
-    }
-    vr  += (0.5 * dt) * f_ext / ext_mass;
-    xr  += dt * vr;
-    xr.apply_constraints();
-    if (this->b_periodic) this->wrap(xr);
-  }
-
-
-  // Adding wall potential to "true" colvar force, whether or not an extended coordinate is in use
   if (is_enabled(f_cv_lower_wall) || is_enabled(f_cv_upper_wall)) {
 
     // Wall force
@@ -1257,6 +1223,37 @@ cvm::real colvar::update_forces_energy()
     }
   }
 
+  if (is_enabled(f_cv_extended_Lagrangian)) {
+
+    cvm::real dt = cvm::dt();
+    cvm::real f_ext;
+
+    // the total force is applied to the fictitious mass, while the
+    // atoms only feel the harmonic force
+    // fr: bias force on extended coordinate (without harmonic spring), for output in trajectory
+    // f_ext: total force on extended coordinate (including harmonic spring)
+    // f: - initially, external biasing force
+    //    - after this code block, colvar force to be applied to atomic coordinates, ie. spring force
+    fr    = f;
+    f_ext = f + (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x);
+    f     =     (-0.5 * ext_force_k) * this->dist2_rgrad(xr, x);
+
+    // leapfrog: starting from x_i, f_i, v_(i-1/2)
+    vr  += (0.5 * dt) * f_ext / ext_mass;
+    // Because of leapfrog, kinetic energy at time i is approximate
+    kinetic_energy = 0.5 * ext_mass * vr * vr;
+    potential_energy = 0.5 * ext_force_k * this->dist2(xr, x);
+    // leap to v_(i+1/2)
+    if (is_enabled(f_cv_Langevin)) {
+      vr -= dt * ext_gamma * vr.real_value;
+      vr += dt * ext_sigma * cvm::rand_gaussian() / ext_mass;
+    }
+    vr  += (0.5 * dt) * f_ext / ext_mass;
+    xr  += dt * vr;
+    xr.apply_constraints();
+    if (this->b_periodic) this->wrap(xr);
+  }
+
   f_accumulated += f;
 
   if (is_enabled(f_cv_fdiff_velocity)) {
diff --git a/lib/colvars/colvarbias_abf.cpp b/lib/colvars/colvarbias_abf.cpp
index 648ef9e721..475e4ecef0 100644
--- a/lib/colvars/colvarbias_abf.cpp
+++ b/lib/colvars/colvarbias_abf.cpp
@@ -12,6 +12,7 @@ colvarbias_abf::colvarbias_abf(char const *key)
     samples(NULL),
     z_gradients(NULL),
     z_samples(NULL),
+    czar_gradients(NULL),
     last_gradients(NULL),
     last_samples(NULL)
 {
@@ -144,6 +145,7 @@ int colvarbias_abf::init(std::string const &conf)
     z_gradients->request_actual_value();
     z_gradients->samples = z_samples;
     z_samples->has_parent_data = true;
+    czar_gradients = new colvar_grid_gradient(colvars);
   }
 
   // For shared ABF, we store a second set of grids.
@@ -178,6 +180,21 @@ colvarbias_abf::~colvarbias_abf()
     gradients = NULL;
   }
 
+  if (z_samples) {
+    delete z_samples;
+    z_samples = NULL;
+  }
+
+  if (z_gradients) {
+    delete z_gradients;
+    z_gradients = NULL;
+  }
+
+  if (czar_gradients) {
+    delete czar_gradients;
+    czar_gradients = NULL;
+  }
+
   // shared ABF
   // We used to only do this if "shared" was defined,
   // but now we can call shared externally
@@ -237,7 +254,11 @@ int colvarbias_abf::update()
         z_bin[i] = z_samples->current_bin_scalar(i);
       }
       if ( z_samples->index_ok(z_bin) ) {
-        // Set increment flag to 0 to only increment
+        for (size_t i=0; i<colvars.size(); i++) {
+          // If we are outside the range of xi, the force has not been obtained above
+          // the function is just an accessor, so cheap to call again anyway
+          force[i] = colvars[i]->system_force();
+        }
         z_gradients->acc_force(z_bin, force);
       }
     }
@@ -404,14 +425,10 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
 {
   std::string  samples_out_name = prefix + ".count";
   std::string  gradients_out_name = prefix + ".grad";
-  std::string  z_gradients_out_name = prefix + ".zgrad";
-  std::string  z_samples_out_name = prefix + ".zcount";
   std::ios::openmode mode = (append ? std::ios::app : std::ios::out);
 
   cvm::ofstream samples_os;
   cvm::ofstream gradients_os;
-  cvm::ofstream z_samples_os;
-  cvm::ofstream z_gradients_os;
 
   if (!append) cvm::backup_file(samples_out_name.c_str());
   samples_os.open(samples_out_name.c_str(), mode);
@@ -442,6 +459,14 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
   }
 
   if (z_gradients) {
+    // Write eABF-related quantities
+    std::string  z_samples_out_name = prefix + ".zcount";
+    std::string  z_gradients_out_name = prefix + ".zgrad";
+    std::string  czar_gradients_out_name = prefix + ".czar";
+    cvm::ofstream z_samples_os;
+    cvm::ofstream z_gradients_os;
+    cvm::ofstream czar_gradients_os;
+
     if (!append) cvm::backup_file(z_samples_out_name.c_str());
     z_samples_os.open(z_samples_out_name.c_str(), mode);
     if (!z_samples_os.is_open()) {
@@ -458,6 +483,24 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
     z_gradients->write_multicol(z_gradients_os);
     z_gradients_os.close();
 
+    // Calculate CZAR estimator of gradients
+    for (std::vector<int> ix = czar_gradients->new_index();
+          czar_gradients->index_ok(ix); czar_gradients->incr(ix)) {
+      for (size_t n = 0; n < czar_gradients->multiplicity(); n++) {
+        czar_gradients->set_value(ix, z_gradients->value_output(ix, n)
+            - cvm::temperature() * cvm::boltzmann() * z_samples->log_gradient_finite_diff(ix, n),
+            n);
+      }
+    }
+
+    if (!append) cvm::backup_file(czar_gradients_out_name.c_str());
+    czar_gradients_os.open(czar_gradients_out_name.c_str(), mode);
+    if (!czar_gradients_os.is_open()) {
+      cvm::error("Error opening CZAR gradient file " + czar_gradients_out_name + " for writing");
+    }
+    czar_gradients->write_multicol(czar_gradients_os);
+    czar_gradients_os.close();
+
     if (colvars.size() == 1) {
       std::string  z_pmf_out_name = prefix + ".zpmf";
       if (!append) cvm::backup_file(z_pmf_out_name.c_str());
@@ -468,6 +511,16 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
       z_gradients->write_1D_integral(z_pmf_os);
       z_pmf_os << std::endl;
       z_pmf_os.close();
+
+      std::string  czar_pmf_out_name = prefix + ".czarpmf";
+      if (!append) cvm::backup_file(czar_pmf_out_name.c_str());
+      cvm::ofstream czar_pmf_os;
+      // Do numerical integration and output a PMF
+      czar_pmf_os.open(czar_pmf_out_name.c_str(), mode);
+      if (!czar_pmf_os.is_open())  cvm::error("Error opening CZAR pmf file " + czar_pmf_out_name + " for writing");
+      czar_gradients->write_1D_integral(czar_pmf_os);
+      czar_pmf_os << std::endl;
+      czar_pmf_os.close();
     }
   }
 
@@ -545,24 +598,27 @@ std::ostream & colvarbias_abf::write_restart(std::ostream& os)
 {
 
   std::ios::fmtflags flags(os.flags());
-  os.setf(std::ios::fmtflags(0), std::ios::floatfield); // default floating-point format
 
   os << "abf {\n"
      << "  configuration {\n"
      << "    name " << this->name << "\n";
   os << "  }\n";
 
-  os << "samples\n";
+  os.setf(std::ios::fmtflags(0), std::ios::floatfield); // default floating-point format
+  os << "\nsamples\n";
   samples->write_raw(os, 8);
+  os.flags(flags);
 
   os << "\ngradient\n";
-  gradients->write_raw(os);
+  gradients->write_raw(os, 8);
 
   if (z_gradients) {
-    os << "z_samples\n";
+    os.setf(std::ios::fmtflags(0), std::ios::floatfield); // default floating-point format
+    os << "\nz_samples\n";
     z_samples->write_raw(os, 8);
+    os.flags(flags);
     os << "\nz_gradient\n";
-    z_gradients->write_raw(os);
+    z_gradients->write_raw(os, 8);
   }
 
   os << "}\n\n";
@@ -638,7 +694,7 @@ std::istream & colvarbias_abf::read_restart(std::istream& is)
   }
 
   if (z_gradients) {
-    if ( !(is >> key)   || !(key == "z_samples")) {
+    if ( !(is >> key) || !(key == "z_samples")) {
       cvm::log("Error: in reading restart configuration for ABF bias \""+
                 this->name+"\" at position "+
                 cvm::to_str(is.tellg())+" in stream.\n");
diff --git a/lib/colvars/colvarbias_abf.h b/lib/colvars/colvarbias_abf.h
index 3527eeb63a..18b9a463b5 100644
--- a/lib/colvars/colvarbias_abf.h
+++ b/lib/colvars/colvarbias_abf.h
@@ -60,6 +60,8 @@ private:
   colvar_grid_gradient  *z_gradients;
   /// n-dim grid of number of samples on "real" coordinate for eABF z-based estimator
   colvar_grid_count     *z_samples;
+  /// n-dim grid contining CZAR estimator of "real" free energy gradients
+  colvar_grid_gradient  *czar_gradients;
 
   // shared ABF
   bool     shared_on;
diff --git a/lib/colvars/colvarbias_histogram.cpp b/lib/colvars/colvarbias_histogram.cpp
index 0d52c18361..542dde2017 100644
--- a/lib/colvars/colvarbias_histogram.cpp
+++ b/lib/colvars/colvarbias_histogram.cpp
@@ -75,13 +75,12 @@ int colvarbias_histogram::init(std::string const &conf)
   }
 
   grid = new colvar_grid_scalar();
+  grid->init_from_colvars(colvars);
 
   {
     std::string grid_conf;
     if (key_lookup(conf, "grid", grid_conf)) {
       grid->parse_params(grid_conf);
-    } else {
-      grid->init_from_colvars(colvars);
     }
   }
 
@@ -256,6 +255,9 @@ std::istream & colvarbias_histogram::read_restart(std::istream& is)
 
 std::ostream & colvarbias_histogram::write_restart(std::ostream& os)
 {
+  std::ios::fmtflags flags(os.flags());
+  os.setf(std::ios::fmtflags(0), std::ios::floatfield);
+
   os << "histogram {\n"
      << "  configuration {\n"
      << "    name " << this->name << "\n";
@@ -266,5 +268,6 @@ std::ostream & colvarbias_histogram::write_restart(std::ostream& os)
 
   os << "}\n\n";
 
+  os.flags(flags);
   return os;
 }
diff --git a/lib/colvars/colvargrid.cpp b/lib/colvars/colvargrid.cpp
index 3b07a0b3d5..ca2d935e1c 100644
--- a/lib/colvars/colvargrid.cpp
+++ b/lib/colvars/colvargrid.cpp
@@ -35,13 +35,13 @@ colvar_grid_scalar::colvar_grid_scalar(colvar_grid_scalar const &g)
 }
 
 colvar_grid_scalar::colvar_grid_scalar(std::vector<int> const &nx_i)
-  : colvar_grid<cvm::real>(nx_i, 0.0, 1), samples(NULL)
+  : colvar_grid<cvm::real>(nx_i, 0.0, 1), samples(NULL), grad(NULL)
 {
   grad = new cvm::real[nd];
 }
 
 colvar_grid_scalar::colvar_grid_scalar(std::vector<colvar *> &colvars, bool margin)
-  : colvar_grid<cvm::real>(colvars, 0.0, 1, margin), samples(NULL)
+  : colvar_grid<cvm::real>(colvars, 0.0, 1, margin), samples(NULL), grad(NULL)
 {
   grad = new cvm::real[nd];
 }
diff --git a/lib/colvars/colvargrid.h b/lib/colvars/colvargrid.h
index e00405d278..473348d3e7 100644
--- a/lib/colvars/colvargrid.h
+++ b/lib/colvars/colvargrid.h
@@ -1160,6 +1160,50 @@ public:
     }
     has_data = true;
   }
+
+  /// \brief Return the log-gradient from finite differences
+  /// on the *same* grid for dimension n
+  inline const cvm::real log_gradient_finite_diff( const std::vector<int> &ix0,
+    int n = 0)
+  {
+    cvm::real A0, A1;
+    std::vector<int> ix;
+
+    // factor for mesh width, 2.0 for central finite difference
+    // but only 1.0 on edges for non-PBC coordinates
+    cvm::real factor;
+
+    if (periodic[n]) {
+      factor = 2.;
+      ix = ix0;
+      ix[n]--; wrap(ix);
+      A0 = data[address(ix)];
+      ix = ix0;
+      ix[n]++; wrap(ix);
+      A1 = data[address(ix)];
+    } else {
+      factor = 0.;
+      ix = ix0;
+      if (ix[n] > 0) { // not left edge
+        ix[n]--;
+        factor += 1.;
+      }
+      A0 = data[address(ix)];
+      ix = ix0;
+      if (ix[n]+1 < nx[n]) { // not right edge
+        ix[n]++;
+        factor += 1.;
+      }
+      A1 = data[address(ix)];
+    }
+    if (A0 == 0 || A1 == 0) {
+      // can't handle empty bins
+      return 0.;
+    } else {
+      return (std::log((cvm::real)A1) - std::log((cvm::real)A0))
+        / (widths[n] * factor);
+    }
+  }
 };
 
 
diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h
index e6d21af85b..b8fc01be9e 100644
--- a/lib/colvars/colvarmodule.h
+++ b/lib/colvars/colvarmodule.h
@@ -4,7 +4,7 @@
 #define COLVARMODULE_H
 
 #ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2016-08-01"
+#define COLVARS_VERSION "2016-08-05"
 #endif
 
 #ifndef COLVARS_DEBUG