diff --git a/doc/src/PDF/colvars-refman-lammps.pdf b/doc/src/PDF/colvars-refman-lammps.pdf
index 071bdb25b8..f1b4cae8bb 100644
Binary files a/doc/src/PDF/colvars-refman-lammps.pdf and b/doc/src/PDF/colvars-refman-lammps.pdf differ
diff --git a/doc/src/fix_colvars.txt b/doc/src/fix_colvars.txt
index 7e8f35922d..1991496d19 100644
--- a/doc/src/fix_colvars.txt
+++ b/doc/src/fix_colvars.txt
@@ -39,8 +39,8 @@ metadynamics, Steered Molecular Dynamics (SMD) and Umbrella Sampling
 (US) via a flexible harmonic restraint bias. The colvars library is
 hosted at "http://colvars.github.io/"_http://colvars.github.io/
 
-This documentation describes only the fix colvars command itself and 
-LAMMPS specific parts of the code.  The full documentation of the 
+This documentation describes only the fix colvars command itself and
+LAMMPS specific parts of the code.  The full documentation of the
 colvars library is available as "this supplementary PDF document"_PDF/colvars-refman-lammps.pdf
 
 A detailed discussion of the implementation of the portable collective
@@ -122,7 +122,7 @@ not a limitation of functionality.
 
 [Default:]
 
-The default options are input = NULL, output = out, seed = 1966, unwrap yes, 
+The default options are input = NULL, output = out, seed = 1966, unwrap yes,
 and tstat = NULL.
 
 :line
diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp
index 84880808af..dbf5fa5219 100644
--- a/lib/colvars/colvar.cpp
+++ b/lib/colvars/colvar.cpp
@@ -150,7 +150,7 @@ colvar::colvar(std::string const &conf)
     feature_states[f_cv_linear]->enabled = lin;
   }
 
-  // Colvar is homogeneous iff:
+  // Colvar is homogeneous if:
   // - it is linear (hence not scripted)
   // - all cvcs have coefficient 1 or -1
   // i.e. sum or difference of cvcs
@@ -375,28 +375,16 @@ colvar::colvar(std::string const &conf)
 
   {
     bool temp;
-    if (get_keyval(conf, "outputSystemForce", temp, false)) {
-      cvm::error("Colvar option outputSystemForce is deprecated.\n"
-        "Please use outputTotalForce, or outputSystemForce within an ABF bias.");
+    if (get_keyval(conf, "outputSystemForce", temp, false, colvarparse::parse_silent)) {
+      cvm::error("Option outputSystemForce is deprecated: only outputTotalForce is supported instead.\n"
+                 "The two are NOT identical: see http://colvars.github.io/totalforce.html.\n", INPUT_ERROR);
       return;
     }
   }
 
-  {
-    bool b_output_total_force;
-    get_keyval(conf, "outputTotalForce", b_output_total_force, false);
-    if (b_output_total_force) {
-      enable(f_cv_output_total_force);
-    }
-  }
-
-  {
-    bool b_output_applied_force;
-    get_keyval(conf, "outputAppliedForce", b_output_applied_force, false);
-    if (b_output_applied_force) {
-      enable(f_cv_output_applied_force);
-    }
-  }
+  get_keyval_feature(this, conf, "outputTotalForce", f_cv_output_total_force, false);
+  get_keyval_feature(this, conf, "outputAppliedForce", f_cv_output_applied_force, false);
+  get_keyval_feature(this, conf, "subtractAppliedForce", f_cv_subtract_applied_force, false);
 
   // Start in active state by default
   enable(f_cv_active);
@@ -409,6 +397,8 @@ colvar::colvar(std::string const &conf)
   fj.type(value());
   ft.type(value());
   ft_reported.type(value());
+  f_old.type(value());
+  f_old.reset();
 
   if (cvm::b_analysis)
     parse_analysis(conf);
@@ -519,6 +509,8 @@ int colvar::init_components(std::string const &conf)
     "number", "coordNum");
   error_code |= init_components_type<selfcoordnum>(conf, "self-coordination "
     "number", "selfCoordNum");
+  error_code |= init_components_type<groupcoordnum>(conf, "group-coordination "
+    "number", "groupCoord");
   error_code |= init_components_type<angle>(conf, "angle", "angle");
   error_code |= init_components_type<dipole_angle>(conf, "dipole angle", "dipoleAngle");
   error_code |= init_components_type<dihedral>(conf, "dihedral", "dihedral");
@@ -1104,6 +1096,14 @@ int colvar::calc_colvar_properties()
 
   } else {
 
+    if (is_enabled(f_cv_subtract_applied_force)) {
+      // correct the total force only if it has been measured
+      // TODO add a specific test instead of relying on sq norm
+      if (ft.norm2() > 0.0) {
+        ft -= f_old;
+      }
+    }
+
     x_reported = x;
     ft_reported = ft;
   }
@@ -1210,6 +1210,10 @@ cvm::real colvar::update_forces_energy()
     x_old = x;
   }
 
+  if (is_enabled(f_cv_subtract_applied_force)) {
+    f_old = f;
+  }
+
   if (cvm::debug())
     cvm::log("Done updating colvar \""+this->name+"\".\n");
   return (potential_energy + kinetic_energy);
@@ -1640,15 +1644,9 @@ std::ostream & colvar::write_traj(std::ostream &os)
   }
 
   if (is_enabled(f_cv_output_applied_force)) {
-    if (is_enabled(f_cv_extended_Lagrangian)) {
-      os << " "
-         << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
-         << fr;
-    } else {
-      os << " "
-         << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
-         << f;
-    }
+    os << " "
+       << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
+       << applied_force();
   }
 
   return os;
diff --git a/lib/colvars/colvar.h b/lib/colvars/colvar.h
index 0058b56a26..d7b1e6076f 100644
--- a/lib/colvars/colvar.h
+++ b/lib/colvars/colvar.h
@@ -175,6 +175,9 @@ public:
   /// (if defined) contribute to it
   colvarvalue f;
 
+  /// Applied force at the previous step (to be subtracted from total force if needed)
+  colvarvalue f_old;
+
   /// \brief Total force, as derived from the atomic trajectory;
   /// should equal the system force plus \link f \endlink
   colvarvalue ft;
@@ -272,10 +275,13 @@ public:
   /// \brief Calculate the quantities associated to the colvar (but not to the CVCs)
   int calc_colvar_properties();
 
-  /// Get the current biasing force
-  inline colvarvalue bias_force() const
+  /// Get the current applied force
+  inline colvarvalue const applied_force() const
   {
-    return fb;
+    if (is_enabled(f_cv_extended_Lagrangian)) {
+      return fr;
+    }
+    return f;
   }
 
   /// Set the total biasing force to zero
@@ -482,6 +488,7 @@ public:
   class dihedral;
   class coordnum;
   class selfcoordnum;
+  class groupcoordnum;
   class h_bond;
   class rmsd;
   class orientation_angle;
diff --git a/lib/colvars/colvarbias_histogram.cpp b/lib/colvars/colvarbias_histogram.cpp
index 15ad5936aa..eef4b506e3 100644
--- a/lib/colvars/colvarbias_histogram.cpp
+++ b/lib/colvars/colvarbias_histogram.cpp
@@ -67,7 +67,7 @@ int colvarbias_histogram::init(std::string const &conf)
 
   if (colvar_array_size > 0) {
     weights.assign(colvar_array_size, 1.0);
-    get_keyval(conf, "weights", weights, weights, colvarparse::parse_silent);
+    get_keyval(conf, "weights", weights, weights);
   }
 
   for (i = 0; i < colvars.size(); i++) {
@@ -79,7 +79,7 @@ int colvarbias_histogram::init(std::string const &conf)
 
   {
     std::string grid_conf;
-    if (key_lookup(conf, "grid", grid_conf)) {
+    if (key_lookup(conf, "histogramGrid", grid_conf)) {
       grid->parse_params(grid_conf);
     }
   }
diff --git a/lib/colvars/colvarbias_meta.cpp b/lib/colvars/colvarbias_meta.cpp
index eb98a99259..596972253b 100644
--- a/lib/colvars/colvarbias_meta.cpp
+++ b/lib/colvars/colvarbias_meta.cpp
@@ -92,7 +92,11 @@ int colvarbias_meta::init(std::string const &conf)
     get_keyval(conf, "keepHills", keep_hills, false);
     if (! get_keyval(conf, "writeFreeEnergyFile", dump_fes, true))
       get_keyval(conf, "dumpFreeEnergyFile", dump_fes, true, colvarparse::parse_silent);
-    get_keyval(conf, "saveFreeEnergyFile", dump_fes_save, false);
+    if (get_keyval(conf, "saveFreeEnergyFile", dump_fes_save, false, colvarparse::parse_silent)) {
+      cvm::log("Option \"saveFreeEnergyFile\" is deprecated, "
+               "please use \"keepFreeEnergyFiles\" instead.");
+    }
+    get_keyval(conf, "keepFreeEnergyFiles", dump_fes_save, dump_fes_save);
 
     hills_energy           = new colvar_grid_scalar(colvars);
     hills_energy_gradients = new colvar_grid_gradient(colvars);
diff --git a/lib/colvars/colvarbias_restraint.cpp b/lib/colvars/colvarbias_restraint.cpp
index c6ef954e34..bd37a430a1 100644
--- a/lib/colvars/colvarbias_restraint.cpp
+++ b/lib/colvars/colvarbias_restraint.cpp
@@ -612,3 +612,250 @@ cvm::real colvarbias_restraint_linear::restraint_convert_k(cvm::real k,
 
 
 
+colvarbias_restraint_histogram::colvarbias_restraint_histogram(char const *key)
+  : colvarbias(key)
+{
+  lower_boundary = 0.0;
+  upper_boundary = 0.0;
+  width = 0.0;
+  gaussian_width = 0.0;
+}
+
+
+int colvarbias_restraint_histogram::init(std::string const &conf)
+{
+  colvarbias::init(conf);
+
+  get_keyval(conf, "lowerBoundary", lower_boundary, lower_boundary);
+  get_keyval(conf, "upperBoundary", upper_boundary, upper_boundary);
+  get_keyval(conf, "width", width, width);
+
+  if (width <= 0.0) {
+    cvm::error("Error: \"width\" must be positive.\n", INPUT_ERROR);
+  }
+
+  get_keyval(conf, "gaussianWidth", gaussian_width, 2.0 * width, colvarparse::parse_silent);
+  get_keyval(conf, "gaussianSigma", gaussian_width, 2.0 * width);
+
+  if (lower_boundary >= upper_boundary) {
+    cvm::error("Error: the upper boundary, "+
+               cvm::to_str(upper_boundary)+
+               ", is not higher than the lower boundary, "+
+               cvm::to_str(lower_boundary)+".\n",
+               INPUT_ERROR);
+  }
+
+  cvm::real const nbins = (upper_boundary - lower_boundary) / width;
+  int const nbins_round = (int)(nbins);
+
+  if (std::fabs(nbins - cvm::real(nbins_round)) > 1.0E-10) {
+    cvm::log("Warning: grid interval ("+
+             cvm::to_str(lower_boundary, cvm::cv_width, cvm::cv_prec)+" - "+
+             cvm::to_str(upper_boundary, cvm::cv_width, cvm::cv_prec)+
+             ") is not commensurate to its bin width ("+
+             cvm::to_str(width, cvm::cv_width, cvm::cv_prec)+").\n");
+  }
+
+  p.resize(nbins_round);
+  ref_p.resize(nbins_round);
+  p_diff.resize(nbins_round);
+
+  bool const inline_ref_p =
+    get_keyval(conf, "refHistogram", ref_p.data_array(), ref_p.data_array());
+  std::string ref_p_file;
+  get_keyval(conf, "refHistogramFile", ref_p_file, std::string(""));
+  if (ref_p_file.size()) {
+    if (inline_ref_p) {
+      cvm::error("Error: cannot specify both refHistogram and refHistogramFile at the same time.\n",
+                 INPUT_ERROR);
+    } else {
+      std::ifstream is(ref_p_file.c_str());
+      std::string data_s = "";
+      std::string line;
+      while (getline_nocomments(is, line)) {
+        data_s.append(line+"\n");
+      }
+      if (data_s.size() == 0) {
+        cvm::error("Error: file \""+ref_p_file+"\" empty or unreadable.\n", FILE_ERROR);
+      }
+      is.close();
+      cvm::vector1d<cvm::real> data;
+      if (data.from_simple_string(data_s) != 0) {
+        cvm::error("Error: could not read histogram from file \""+ref_p_file+"\".\n");
+      }
+      if (data.size() == 2*ref_p.size()) {
+        // file contains both x and p(x)
+        size_t i;
+        for (i = 0; i < ref_p.size(); i++) {
+          ref_p[i] = data[2*i+1];
+        }
+      } else if (data.size() == ref_p.size()) {
+        ref_p = data;
+      } else {
+        cvm::error("Error: file \""+ref_p_file+"\" contains a histogram of different length.\n",
+                   INPUT_ERROR);
+      }
+    }
+  }
+  cvm::real const ref_integral = ref_p.sum() * width;
+  if (std::fabs(ref_integral - 1.0) > 1.0e-03) {
+    cvm::log("Reference distribution not normalized, normalizing to unity.\n");
+    ref_p /= ref_integral;
+  }
+
+  get_keyval(conf, "writeHistogram", b_write_histogram, false);
+  get_keyval(conf, "forceConstant", force_k, 1.0);
+
+  return COLVARS_OK;
+}
+
+
+colvarbias_restraint_histogram::~colvarbias_restraint_histogram()
+{
+  p.resize(0);
+  ref_p.resize(0);
+  p_diff.resize(0);
+}
+
+
+int colvarbias_restraint_histogram::update()
+{
+  if (cvm::debug())
+    cvm::log("Updating the histogram restraint bias \""+this->name+"\".\n");
+
+  size_t vector_size = 0;
+  size_t icv;
+  for (icv = 0; icv < colvars.size(); icv++) {
+    vector_size += colvars[icv]->value().size();
+  }
+
+  cvm::real const norm = 1.0/(std::sqrt(2.0*PI)*gaussian_width*vector_size);
+
+  // calculate the histogram
+  p.reset();
+  for (icv = 0; icv < colvars.size(); icv++) {
+    colvarvalue const &cv = colvars[icv]->value();
+    if (cv.type() == colvarvalue::type_scalar) {
+      cvm::real const cv_value = cv.real_value;
+      size_t igrid;
+      for (igrid = 0; igrid < p.size(); igrid++) {
+        cvm::real const x_grid = (lower_boundary + (igrid+0.5)*width);
+        p[igrid] += norm * std::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) /
+                                    (2.0 * gaussian_width * gaussian_width));
+      }
+    } else if (cv.type() == colvarvalue::type_vector) {
+      size_t idim;
+      for (idim = 0; idim < cv.vector1d_value.size(); idim++) {
+        cvm::real const cv_value = cv.vector1d_value[idim];
+        size_t igrid;
+        for (igrid = 0; igrid < p.size(); igrid++) {
+          cvm::real const x_grid = (lower_boundary + (igrid+0.5)*width);
+          p[igrid] += norm * std::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) /
+                                      (2.0 * gaussian_width * gaussian_width));
+        }
+      }
+    } else {
+      // TODO
+    }
+  }
+
+  cvm::real const force_k_cv = force_k * vector_size;
+
+  // calculate the difference between current and reference
+  p_diff = p - ref_p;
+  bias_energy = 0.5 * force_k_cv * p_diff * p_diff;
+
+  // calculate the forces
+  for (icv = 0; icv < colvars.size(); icv++) {
+    colvarvalue const &cv = colvars[icv]->value();
+    colvarvalue &cv_force = colvar_forces[icv];
+    cv_force.type(cv);
+    cv_force.reset();
+
+    if (cv.type() == colvarvalue::type_scalar) {
+      cvm::real const cv_value = cv.real_value;
+      cvm::real &force = cv_force.real_value;
+      size_t igrid;
+      for (igrid = 0; igrid < p.size(); igrid++) {
+        cvm::real const x_grid = (lower_boundary + (igrid+0.5)*width);
+        force += force_k_cv * p_diff[igrid] *
+          norm * std::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) /
+                          (2.0 * gaussian_width * gaussian_width)) *
+          (-1.0 * (x_grid - cv_value) / (gaussian_width * gaussian_width));
+      }
+    } else if (cv.type() == colvarvalue::type_vector) {
+      size_t idim;
+      for (idim = 0; idim < cv.vector1d_value.size(); idim++) {
+        cvm::real const cv_value = cv.vector1d_value[idim];
+        cvm::real &force = cv_force.vector1d_value[idim];
+        size_t igrid;
+        for (igrid = 0; igrid < p.size(); igrid++) {
+          cvm::real const x_grid = (lower_boundary + (igrid+0.5)*width);
+          force += force_k_cv * p_diff[igrid] *
+            norm * std::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) /
+                            (2.0 * gaussian_width * gaussian_width)) *
+            (-1.0 * (x_grid - cv_value) / (gaussian_width * gaussian_width));
+        }
+      }
+    } else {
+      // TODO
+    }
+  }
+
+  return COLVARS_OK;
+}
+
+
+std::ostream & colvarbias_restraint_histogram::write_restart(std::ostream &os)
+{
+  if (b_write_histogram) {
+    std::string file_name(cvm::output_prefix+"."+this->name+".hist.dat");
+    std::ofstream os(file_name.c_str());
+    os << "# " << cvm::wrap_string(colvars[0]->name, cvm::cv_width)
+       << "  " << "p(" << cvm::wrap_string(colvars[0]->name, cvm::cv_width-3)
+       << ")\n";
+    size_t igrid;
+    for (igrid = 0; igrid < p.size(); igrid++) {
+      cvm::real const x_grid = (lower_boundary + (igrid+1)*width);
+      os << "  "
+         << std::setprecision(cvm::cv_prec)
+         << std::setw(cvm::cv_width)
+         << x_grid
+         << "  "
+         << std::setprecision(cvm::cv_prec)
+         << std::setw(cvm::cv_width)
+         << p[igrid] << "\n";
+    }
+    os.close();
+  }
+  return os;
+}
+
+
+std::istream & colvarbias_restraint_histogram::read_restart(std::istream &is)
+{
+  return is;
+}
+
+
+std::ostream & colvarbias_restraint_histogram::write_traj_label(std::ostream &os)
+{
+  os << " ";
+  if (b_output_energy) {
+    os << " E_"
+       << cvm::wrap_string(this->name, cvm::en_width-2);
+  }
+  return os;
+}
+
+
+std::ostream & colvarbias_restraint_histogram::write_traj(std::ostream &os)
+{
+  os << " ";
+  if (b_output_energy) {
+    os << " "
+       << std::setprecision(cvm::en_prec) << std::setw(cvm::en_width)
+       << bias_energy;
+  }
+  return os;
+}
diff --git a/lib/colvars/colvarbias_restraint.h b/lib/colvars/colvarbias_restraint.h
index a0f6da79ee..54097bd524 100644
--- a/lib/colvars/colvarbias_restraint.h
+++ b/lib/colvars/colvarbias_restraint.h
@@ -168,4 +168,52 @@ protected:
 };
 
 
+/// Restrain the 1D histogram of a set of variables (or of a multidimensional one)
+// TODO this could be reimplemented more cleanly as a derived class of both restraint and histogram
+class colvarbias_restraint_histogram : public colvarbias {
+
+public:
+
+  colvarbias_restraint_histogram(char const *key);
+  int init(std::string const &conf);
+  ~colvarbias_restraint_histogram();
+
+  virtual int update();
+
+  virtual std::istream & read_restart(std::istream &is);
+  virtual std::ostream & write_restart(std::ostream &os);
+  virtual std::ostream & write_traj_label(std::ostream &os);
+  virtual std::ostream & write_traj(std::ostream &os);
+
+protected:
+
+  /// Probability density
+  cvm::vector1d<cvm::real> p;
+
+  /// Reference probability density
+  cvm::vector1d<cvm::real> ref_p;
+
+  /// Difference between probability density and reference
+  cvm::vector1d<cvm::real> p_diff;
+
+  /// Lower boundary of the grid
+  cvm::real lower_boundary;
+
+  /// Upper boundary of the grid
+  cvm::real upper_boundary;
+
+  /// Resolution of the grid
+  cvm::real width;
+
+  /// Width of the Gaussians
+  cvm::real gaussian_width;
+
+  /// Restraint force constant
+  cvm::real force_k;
+
+  /// Write the histogram to a file
+  bool b_write_histogram;
+};
+
+
 #endif
diff --git a/lib/colvars/colvarcomp.cpp b/lib/colvars/colvarcomp.cpp
index 534a927fda..73b1dbd81c 100644
--- a/lib/colvars/colvarcomp.cpp
+++ b/lib/colvars/colvarcomp.cpp
@@ -54,6 +54,21 @@ colvar::cvc::cvc(std::string const &conf)
 }
 
 
+int colvar::cvc::init_total_force_params(std::string const &conf)
+{
+  if (get_keyval_feature(this, conf, "oneSiteSystemForce",
+                         f_cvc_one_site_total_force, is_enabled(f_cvc_one_site_total_force))) {
+    cvm::log("Warning: keyword \"oneSiteSystemForce\" is deprecated: "
+             "please use \"oneSiteTotalForce\" instead.\n");
+  }
+  if (get_keyval_feature(this, conf, "oneSiteTotalForce",
+                         f_cvc_one_site_total_force, is_enabled(f_cvc_one_site_total_force))) {
+    cvm::log("Computing total force on group 1 only");
+  }
+  return COLVARS_OK;
+}
+
+
 cvm::atom_group *colvar::cvc::parse_group(std::string const &conf,
                                           char const *group_key,
                                           bool optional)
@@ -77,8 +92,6 @@ cvm::atom_group *colvar::cvc::parse_group(std::string const &conf,
 
       if (is_enabled(f_cvc_scalable)) {
         cvm::log("Will enable scalable calculation for group \""+group->key+"\".\n");
-      } else {
-        cvm::log("Scalable calculation is not available for group \""+group->key+"\" with the current configuration.\n");
       }
     }
 
diff --git a/lib/colvars/colvarcomp.h b/lib/colvars/colvarcomp.h
index f632754ab6..0649895e45 100644
--- a/lib/colvars/colvarcomp.h
+++ b/lib/colvars/colvarcomp.h
@@ -108,6 +108,9 @@ public:
                    char const *group_key,
                    bool optional = false);
 
+  /// \brief Parse options pertaining to total force calculation
+  virtual int init_total_force_params(std::string const &conf);
+
   /// \brief After construction, set data related to dependency handling
   int setup();
 
@@ -306,9 +309,6 @@ protected:
   cvm::rvector     dist_v;
   /// Use absolute positions, ignoring PBCs when present
   bool b_no_PBC;
-  /// Compute total force on first site only to avoid unwanted
-  /// coupling to other colvars (see e.g. Ciccotti et al., 2005)
-  bool b_1site_force;
 public:
   distance(std::string const &conf);
   distance();
@@ -388,9 +388,6 @@ protected:
   cvm::atom_group  *ref2;
   /// Use absolute positions, ignoring PBCs when present
   bool b_no_PBC;
-  /// Compute total force on one site only to avoid unwanted
-  /// coupling to other colvars (see e.g. Ciccotti et al., 2005)
-  bool b_1site_force;
   /// Vector on which the distance vector is projected
   cvm::rvector axis;
   /// Norm of the axis
@@ -854,6 +851,62 @@ public:
                                   colvarvalue const &x2) const;
 };
 
+
+/// \brief Colvar component: coordination number between two groups
+/// (colvarvalue::type_scalar type, range [0:N1*N2])
+class colvar::groupcoordnum
+  : public colvar::distance
+{
+protected:
+  /// \brief "Cutoff" for isotropic calculation (default)
+  cvm::real     r0;
+  /// \brief "Cutoff vector" for anisotropic calculation
+  cvm::rvector  r0_vec;
+  /// \brief Wheter dist/r0 or \vec{dist}*\vec{1/r0_vec} should ne be
+  /// used
+  bool b_anisotropic;
+  /// Integer exponent of the function numerator
+  int en;
+  /// Integer exponent of the function denominator
+  int ed;
+public:
+  /// Constructor
+  groupcoordnum(std::string const &conf);
+  groupcoordnum();
+  virtual inline ~groupcoordnum() {}
+  virtual void calc_value();
+  virtual void calc_gradients();
+  virtual void apply_force(colvarvalue const &force);
+  template<bool b_gradients>
+  /// \brief Calculate a coordination number through the function
+  /// (1-x**n)/(1-x**m), x = |A1-A2|/r0 \param r0 "cutoff" for the
+  /// coordination number \param exp_num \i n exponent \param exp_den
+  /// \i m exponent \param A1 atom \param A2 atom
+  static cvm::real switching_function(cvm::real const &r0,
+                                      int const &exp_num, int const &exp_den,
+                                      cvm::atom &A1, cvm::atom &A2);
+
+  /*
+  template<bool b_gradients>
+  /// \brief Calculate a coordination number through the function
+  /// (1-x**n)/(1-x**m), x = |(A1-A2)*(r0_vec)^-|1 \param r0_vec
+  /// vector of different cutoffs in the three directions \param
+  /// exp_num \i n exponent \param exp_den \i m exponent \param A1
+  /// atom \param A2 atom
+  static cvm::real switching_function(cvm::rvector const &r0_vec,
+                                      int const &exp_num, int const &exp_den,
+                                      cvm::atom &A1, cvm::atom &A2);
+
+  virtual cvm::real dist2(colvarvalue const &x1,
+                          colvarvalue const &x2) const;
+  virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
+                                  colvarvalue const &x2) const;
+  virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
+                                  colvarvalue const &x2) const;
+  */
+};
+
+
 /// \brief Colvar component: hydrogen bond, defined as the product of
 /// a colvar::coordnum and 1/2*(1-cos((180-ang)/ang_tol))
 /// (colvarvalue::type_scalar type, range [0:1])
diff --git a/lib/colvars/colvarcomp_angles.cpp b/lib/colvars/colvarcomp_angles.cpp
index b1acb201ae..6b05fb16eb 100644
--- a/lib/colvars/colvarcomp_angles.cpp
+++ b/lib/colvars/colvarcomp_angles.cpp
@@ -18,9 +18,9 @@ colvar::angle::angle(std::string const &conf)
   group1 = parse_group(conf, "group1");
   group2 = parse_group(conf, "group2");
   group3 = parse_group(conf, "group3");
-  if (get_keyval(conf, "oneSiteSystemForce", b_1site_force, false)) {
-    cvm::log("Computing total force on group 1 only");
-  }
+
+  init_total_force_params(conf);
+
   x.type(colvarvalue::type_scalar);
 }
 
@@ -33,7 +33,6 @@ colvar::angle::angle(cvm::atom const &a1,
   provide(f_cvc_inv_gradient);
   provide(f_cvc_Jacobian);
   provide(f_cvc_com_based);
-  b_1site_force = false;
 
   group1 = new cvm::atom_group(std::vector<cvm::atom>(1, a1));
   group2 = new cvm::atom_group(std::vector<cvm::atom>(1, a2));
@@ -94,7 +93,7 @@ void colvar::angle::calc_force_invgrads()
   // centered on group2, which means group2 is kept fixed
   // when propagating changes in the angle)
 
-  if (b_1site_force) {
+  if (is_enabled(f_cvc_one_site_total_force)) {
     group1->read_total_forces();
     cvm::real norm_fact = 1.0 / dxdr1.norm2();
     ft.real_value = norm_fact * dxdr1 * group1->total_force();
@@ -140,9 +139,8 @@ colvar::dipole_angle::dipole_angle(std::string const &conf)
   group2 = parse_group(conf, "group2");
   group3 = parse_group(conf, "group3");
 
-  if (get_keyval(conf, "oneSiteSystemForce", b_1site_force, false)) {
-    cvm::log("Computing total force on group 1 only");
-  }
+  init_total_force_params(conf);
+
   x.type(colvarvalue::type_scalar);
 }
 
@@ -152,7 +150,6 @@ colvar::dipole_angle::dipole_angle(cvm::atom const &a1,
                       cvm::atom const &a3)
 {
   function_type = "dipole_angle";
-  b_1site_force = false;
 
   group1 = new cvm::atom_group(std::vector<cvm::atom>(1, a1));
   group2 = new cvm::atom_group(std::vector<cvm::atom>(1, a2));
@@ -250,14 +247,13 @@ colvar::dihedral::dihedral(std::string const &conf)
   provide(f_cvc_Jacobian);
   provide(f_cvc_com_based);
 
-  if (get_keyval(conf, "oneSiteSystemForce", b_1site_force, false)) {
-    cvm::log("Computing total force on group 1 only");
-  }
   group1 = parse_group(conf, "group1");
   group2 = parse_group(conf, "group2");
   group3 = parse_group(conf, "group3");
   group4 = parse_group(conf, "group4");
 
+  init_total_force_params(conf);
+
   x.type(colvarvalue::type_scalar);
 }
 
@@ -422,7 +418,7 @@ void colvar::dihedral::calc_force_invgrads()
   cvm::real const fact4 = d34 * std::sqrt(1.0 - dot4 * dot4);
 
   group1->read_total_forces();
-  if ( b_1site_force ) {
+  if (is_enabled(f_cvc_one_site_total_force)) {
     // This is only measuring the force on group 1
     ft.real_value = PI/180.0 * fact1 * (cross1 * group1->total_force());
   } else {
diff --git a/lib/colvars/colvarcomp_coordnums.cpp b/lib/colvars/colvarcomp_coordnums.cpp
index 390cdc1e59..d5e87e6b1e 100644
--- a/lib/colvars/colvarcomp_coordnums.cpp
+++ b/lib/colvars/colvarcomp_coordnums.cpp
@@ -338,3 +338,151 @@ void colvar::selfcoordnum::apply_force(colvarvalue const &force)
   }
 }
 
+
+// groupcoordnum member functions
+colvar::groupcoordnum::groupcoordnum(std::string const &conf)
+  : distance(conf), b_anisotropic(false)
+{
+  function_type = "groupcoordnum";
+  x.type(colvarvalue::type_scalar);
+
+  // group1 and group2 are already initialized by distance()
+  if (group1->b_dummy || group2->b_dummy)
+    cvm::fatal_error("Error: neither group can be a dummy atom\n");
+
+  bool const b_scale = get_keyval(conf, "cutoff", r0,
+                                  cvm::real(4.0 * cvm::unit_angstrom()));
+
+  if (get_keyval(conf, "cutoff3", r0_vec,
+                 cvm::rvector(4.0, 4.0, 4.0), parse_silent)) {
+
+    if (b_scale)
+      cvm::fatal_error("Error: cannot specify \"scale\" and "
+                       "\"scale3\" at the same time.\n");
+    b_anisotropic = true;
+    // remove meaningless negative signs
+    if (r0_vec.x < 0.0) r0_vec.x *= -1.0;
+    if (r0_vec.y < 0.0) r0_vec.y *= -1.0;
+    if (r0_vec.z < 0.0) r0_vec.z *= -1.0;
+  }
+
+  get_keyval(conf, "expNumer", en, int(6) );
+  get_keyval(conf, "expDenom", ed, int(12));
+
+  if ( (en%2) || (ed%2) ) {
+    cvm::fatal_error("Error: odd exponents provided, can only use even ones.\n");
+  }
+
+}
+
+
+colvar::groupcoordnum::groupcoordnum()
+  : b_anisotropic(false)
+{
+  function_type = "groupcoordnum";
+  x.type(colvarvalue::type_scalar);
+}
+
+
+template<bool calculate_gradients>
+cvm::real colvar::groupcoordnum::switching_function(cvm::real const &r0,
+                                                    int const &en,
+                                                    int const &ed,
+                                                    cvm::atom &A1,
+                                                    cvm::atom &A2)
+{
+  cvm::rvector const diff = cvm::position_distance(A1.pos, A2.pos);
+  cvm::real const l2 = diff.norm2()/(r0*r0);
+
+  // Assume en and ed are even integers, and avoid sqrt in the following
+  int const en2 = en/2;
+  int const ed2 = ed/2;
+
+  cvm::real const xn = std::pow(l2, en2);
+  cvm::real const xd = std::pow(l2, ed2);
+  cvm::real const func = (1.0-xn)/(1.0-xd);
+
+  if (calculate_gradients) {
+    cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) - func*ed2*(xd/l2))*(-1.0);
+    cvm::rvector const dl2dx = (2.0/(r0*r0))*diff;
+    A1.grad += (-1.0)*dFdl2*dl2dx;
+    A2.grad +=        dFdl2*dl2dx;
+  }
+
+  return func;
+}
+
+
+#if 0  // AMG: I don't think there's any reason to support anisotropic,
+       //      and I don't have those flags below in calc_value, but
+       //      if I need them, I'll also need to uncomment this method
+template<bool calculate_gradients>
+cvm::real colvar::groupcoordnum::switching_function(cvm::rvector const &r0_vec,
+                                                    int const &en,
+                                                    int const &ed,
+                                                    cvm::atom &A1,
+                                                    cvm::atom &A2)
+{
+  cvm::rvector const diff = cvm::position_distance(A1.pos, A2.pos);
+  cvm::rvector const scal_diff(diff.x/r0_vec.x, diff.y/r0_vec.y, diff.z/r0_vec.z);
+  cvm::real const l2 = scal_diff.norm2();
+
+  // Assume en and ed are even integers, and avoid sqrt in the following
+  int const en2 = en/2;
+  int const ed2 = ed/2;
+
+  cvm::real const xn = std::pow(l2, en2);
+  cvm::real const xd = std::pow(l2, ed2);
+  cvm::real const func = (1.0-xn)/(1.0-xd);
+
+  if (calculate_gradients) {
+    cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) - func*ed2*(xd/l2))*(-1.0);
+    cvm::rvector const dl2dx((2.0/(r0_vec.x*r0_vec.x))*diff.x,
+                             (2.0/(r0_vec.y*r0_vec.y))*diff.y,
+                             (2.0/(r0_vec.z*r0_vec.z))*diff.z);
+    A1.grad += (-1.0)*dFdl2*dl2dx;
+    A2.grad +=        dFdl2*dl2dx;
+  }
+  return func;
+}
+#endif
+
+
+
+void colvar::groupcoordnum::calc_value()
+{
+
+  // create fake atoms to hold the com coordinates
+  cvm::atom group1_com_atom;
+  cvm::atom group2_com_atom;
+  group1_com_atom.pos = group1->center_of_mass();
+  group2_com_atom.pos = group2->center_of_mass();
+
+  x.real_value = coordnum::switching_function<false>(r0, en, ed,
+                                                     group1_com_atom, group2_com_atom);
+
+}
+
+
+void colvar::groupcoordnum::calc_gradients()
+{
+  cvm::atom group1_com_atom;
+  cvm::atom group2_com_atom;
+  group1_com_atom.pos = group1->center_of_mass();
+  group2_com_atom.pos = group2->center_of_mass();
+
+  coordnum::switching_function<true>(r0, en, ed, group1_com_atom, group2_com_atom);
+  group1->set_weighted_gradient(group1_com_atom.grad);
+  group2->set_weighted_gradient(group2_com_atom.grad);
+
+}
+
+
+void colvar::groupcoordnum::apply_force(colvarvalue const &force)
+{
+  if (!group1->noforce)
+    group1->apply_colvar_force(force.real_value);
+
+  if (!group2->noforce)
+    group2->apply_colvar_force(force.real_value);
+}
diff --git a/lib/colvars/colvarcomp_distances.cpp b/lib/colvars/colvarcomp_distances.cpp
index f7bf2b87f0..15693778b2 100644
--- a/lib/colvars/colvarcomp_distances.cpp
+++ b/lib/colvars/colvarcomp_distances.cpp
@@ -18,14 +18,14 @@ colvar::distance::distance(std::string const &conf)
   provide(f_cvc_Jacobian);
   provide(f_cvc_com_based);
 
+  group1 = parse_group(conf, "group1");
+  group2 = parse_group(conf, "group2");
+
   if (get_keyval(conf, "forceNoPBC", b_no_PBC, false)) {
     cvm::log("Computing distance using absolute positions (not minimal-image)");
   }
-  if (get_keyval(conf, "oneSiteSystemForce", b_1site_force, false)) {
-    cvm::log("Computing total force on group 1 only");
-  }
-  group1 = parse_group(conf, "group1");
-  group2 = parse_group(conf, "group2");
+
+  init_total_force_params(conf);
 
   x.type(colvarvalue::type_scalar);
 }
@@ -38,7 +38,6 @@ colvar::distance::distance()
   provide(f_cvc_inv_gradient);
   provide(f_cvc_Jacobian);
   provide(f_cvc_com_based);
-  b_1site_force = false;
   b_no_PBC = false;
   x.type(colvarvalue::type_scalar);
 }
@@ -67,7 +66,7 @@ void colvar::distance::calc_gradients()
 void colvar::distance::calc_force_invgrads()
 {
   group1->read_total_forces();
-  if ( b_1site_force ) {
+  if (is_enabled(f_cvc_one_site_total_force)) {
     ft.real_value = -1.0 * (group1->total_force() * dist_v.unit());
   } else {
     group2->read_total_forces();
@@ -97,6 +96,7 @@ colvar::distance_vec::distance_vec(std::string const &conf)
   : distance(conf)
 {
   function_type = "distance_vec";
+  provide(f_cvc_com_based);
   x.type(colvarvalue::type_3vector);
 }
 
@@ -105,6 +105,7 @@ colvar::distance_vec::distance_vec()
   : distance()
 {
   function_type = "distance_vec";
+  provide(f_cvc_com_based);
   x.type(colvarvalue::type_3vector);
 }
 
@@ -185,9 +186,9 @@ colvar::distance_z::distance_z(std::string const &conf)
   if (get_keyval(conf, "forceNoPBC", b_no_PBC, false)) {
     cvm::log("Computing distance using absolute positions (not minimal-image)");
   }
-  if (get_keyval(conf, "oneSiteSystemForce", b_1site_force, false)) {
-    cvm::log("Computing total force on group \"main\" only");
-  }
+
+  init_total_force_params(conf);
+
 }
 
 colvar::distance_z::distance_z()
@@ -251,7 +252,7 @@ void colvar::distance_z::calc_force_invgrads()
 {
   main->read_total_forces();
 
-  if (fixed_axis && !b_1site_force) {
+  if (fixed_axis && !is_enabled(f_cvc_one_site_total_force)) {
     ref1->read_total_forces();
     ft.real_value = 0.5 * ((main->total_force() - ref1->total_force()) * axis);
   } else {
@@ -351,7 +352,7 @@ void colvar::distance_xy::calc_force_invgrads()
 {
   main->read_total_forces();
 
-  if (fixed_axis && !b_1site_force) {
+  if (fixed_axis && !is_enabled(f_cvc_one_site_total_force)) {
     ref1->read_total_forces();
     ft.real_value = 0.5 / x.real_value * ((main->total_force() - ref1->total_force()) * dist_v_ortho);
   } else {
@@ -382,6 +383,7 @@ colvar::distance_dir::distance_dir(std::string const &conf)
   : distance(conf)
 {
   function_type = "distance_dir";
+  provide(f_cvc_com_based);
   x.type(colvarvalue::type_unit3vector);
 }
 
@@ -390,6 +392,7 @@ colvar::distance_dir::distance_dir()
   : distance()
 {
   function_type = "distance_dir";
+  provide(f_cvc_com_based);
   x.type(colvarvalue::type_unit3vector);
 }
 
@@ -461,7 +464,6 @@ colvar::distance_inv::distance_inv()
 {
   function_type = "distance_inv";
   exponent = 6;
-  b_1site_force = false;
   x.type(colvarvalue::type_scalar);
 }
 
diff --git a/lib/colvars/colvardeps.cpp b/lib/colvars/colvardeps.cpp
index 7dce7e31f6..7736ea14bd 100644
--- a/lib/colvars/colvardeps.cpp
+++ b/lib/colvars/colvardeps.cpp
@@ -293,6 +293,9 @@ void colvardeps::init_cv_requires() {
     f_description(f_cv_output_total_force, "output total force");
     f_req_self(f_cv_output_total_force, f_cv_total_force);
 
+    f_description(f_cv_subtract_applied_force, "subtract applied force from total force");
+    f_req_self(f_cv_subtract_applied_force, f_cv_total_force);
+
     f_description(f_cv_lower_boundary, "lower boundary");
     f_req_self(f_cv_lower_boundary, f_cv_scalar);
 
@@ -376,6 +379,11 @@ void colvardeps::init_cvc_requires() {
 
     f_description(f_cvc_com_based, "depends on group centers of mass");
 
+    // Compute total force on first site only to avoid unwanted
+    // coupling to other colvars (see e.g. Ciccotti et al., 2005)
+    f_description(f_cvc_one_site_total_force, "compute total collective force only from one group center");
+    f_req_self(f_cvc_one_site_total_force, f_cvc_com_based);
+
     f_description(f_cvc_scalable, "scalable calculation");
     f_req_self(f_cvc_scalable, f_cvc_scalable_com);
 
diff --git a/lib/colvars/colvardeps.h b/lib/colvars/colvardeps.h
index 089461c532..31dfd0e302 100644
--- a/lib/colvars/colvardeps.h
+++ b/lib/colvars/colvardeps.h
@@ -176,6 +176,8 @@ public:
     f_cv_total_force,
     /// \brief Calculate total force from atomic forces
     f_cv_total_force_calc,
+    /// \brief Subtract the applied force from the total force
+    f_cv_subtract_applied_force,
     /// \brief Estimate Jacobian derivative
     f_cv_Jacobian,
     /// \brief Do not report the Jacobian force as part of the total force
@@ -236,6 +238,7 @@ public:
     /// \brief If enabled, calc_gradients() will call debug_gradients() for every group needed
     f_cvc_debug_gradient,
     f_cvc_Jacobian,
+    f_cvc_one_site_total_force,
     f_cvc_com_based,
     f_cvc_scalable,
     f_cvc_scalable_com,
diff --git a/lib/colvars/colvargrid.h b/lib/colvars/colvargrid.h
index 473348d3e7..91d8be9573 100644
--- a/lib/colvars/colvargrid.h
+++ b/lib/colvars/colvargrid.h
@@ -382,8 +382,8 @@ public:
   inline int current_bin_scalar(int const i, int const iv) const
   {
     return value_to_bin_scalar(actual_value[i] ?
-        cv[i]->actual_value().vector1d_value[iv] :
-        cv[i]->value().vector1d_value[iv], i);
+                               cv[i]->actual_value().vector1d_value[iv] :
+                               cv[i]->value().vector1d_value[iv], i);
   }
 
   /// \brief Use the lower boundary and the width to report which bin
@@ -395,8 +395,8 @@ public:
 
   /// \brief Same as the standard version, but uses another grid definition
   inline int value_to_bin_scalar(colvarvalue const &value,
-                                  colvarvalue const &new_offset,
-                                  cvm::real   const &new_width) const
+                                 colvarvalue const &new_offset,
+                                 cvm::real   const &new_width) const
   {
     return (int) std::floor( (value.real_value - new_offset.real_value) / new_width );
   }
@@ -410,22 +410,22 @@ public:
 
   /// \brief Same as the standard version, but uses different parameters
   inline colvarvalue bin_to_value_scalar(int const &i_bin,
-                                          colvarvalue const &new_offset,
-                                          cvm::real const &new_width) const
+                                         colvarvalue const &new_offset,
+                                         cvm::real const &new_width) const
   {
     return new_offset.real_value + new_width * (0.5 + i_bin);
   }
 
   /// Set the value at the point with index ix
   inline void set_value(std::vector<int> const &ix,
-                         T const &t,
-                         size_t const &imult = 0)
+                        T const &t,
+                        size_t const &imult = 0)
   {
     data[this->address(ix)+imult] = t;
     has_data = true;
   }
 
- /// \brief Get the change from this to other_grid
+  /// \brief Get the change from this to other_grid
   /// and store the result in this.
   /// this_grid := other_grid - this_grid
   /// Grids must have the same dimensions.
@@ -434,13 +434,13 @@ public:
 
     if (other_grid.multiplicity() != this->multiplicity()) {
       cvm::error("Error: trying to subtract two grids with "
-                        "different multiplicity.\n");
+                 "different multiplicity.\n");
       return;
     }
 
     if (other_grid.data.size() != this->data.size()) {
       cvm::error("Error: trying to subtract two grids with "
-                        "different size.\n");
+                 "different size.\n");
       return;
     }
 
@@ -457,13 +457,13 @@ public:
   {
     if (other_grid.multiplicity() != this->multiplicity()) {
       cvm::error("Error: trying to copy two grids with "
-                        "different multiplicity.\n");
+                 "different multiplicity.\n");
       return;
     }
 
     if (other_grid.data.size() != this->data.size()) {
       cvm::error("Error: trying to copy two grids with "
-                        "different size.\n");
+                 "different size.\n");
       return;
     }
 
@@ -493,7 +493,7 @@ public:
   /// \brief Get the binned value indexed by ix, or the first of them
   /// if the multiplicity is larger than 1
   inline T const & value(std::vector<int> const &ix,
-                          size_t const &imult = 0) const
+                         size_t const &imult = 0) const
   {
     return data[this->address(ix) + imult];
   }
@@ -541,7 +541,7 @@ public:
   /// boundaries; a negative number is returned if the given point is
   /// off-grid
   inline cvm::real bin_distance_from_boundaries(std::vector<colvarvalue> const &values,
-                                                 bool skip_hard_boundaries = false)
+                                                bool skip_hard_boundaries = false)
   {
     cvm::real minimum = 1.0E+16;
     for (size_t i = 0; i < nd; i++) {
@@ -574,7 +574,7 @@ public:
   {
     if (other_grid.multiplicity() != this->multiplicity()) {
       cvm::error("Error: trying to merge two grids with values of "
-                        "different multiplicity.\n");
+                 "different multiplicity.\n");
       return;
     }
 
@@ -593,8 +593,8 @@ public:
       for (size_t i = 0; i < nd; i++) {
         oix[i] =
           value_to_bin_scalar(bin_to_value_scalar(ix[i], gb[i], gw[i]),
-                               ogb[i],
-                               ogw[i]);
+                              ogb[i],
+                              ogw[i]);
       }
 
       if (! other_grid.index_ok(oix)) {
@@ -614,11 +614,11 @@ public:
   /// \brief Add data from another grid of the same type, AND
   /// identical definition (boundaries, widths)
   void add_grid(colvar_grid<T> const &other_grid,
-                 cvm::real scale_factor = 1.0)
+                cvm::real scale_factor = 1.0)
   {
     if (other_grid.multiplicity() != this->multiplicity()) {
       cvm::error("Error: trying to sum togetehr two grids with values of "
-                        "different multiplicity.\n");
+                 "different multiplicity.\n");
       return;
     }
     if (scale_factor != 1.0)
@@ -636,7 +636,7 @@ public:
   /// \brief Return the value suitable for output purposes (so that it
   /// may be rescaled or manipulated without changing it permanently)
   virtual inline T value_output(std::vector<int> const &ix,
-                                 size_t const &imult = 0)
+                                size_t const &imult = 0)
   {
     return value(ix, imult);
   }
@@ -645,9 +645,9 @@ public:
   /// into the internal representation (the two may be different,
   /// e.g. when using colvar_grid_count)
   virtual inline void value_input(std::vector<int> const &ix,
-                                   T const &t,
-                                   size_t const &imult = 0,
-                                   bool add = false)
+                                  T const &t,
+                                  size_t const &imult = 0,
+                                  bool add = false)
   {
     if ( add )
       data[address(ix) + imult] += t;
@@ -737,7 +737,8 @@ public:
   }
 
   /// Read a grid definition from a config string
-  int parse_params(std::string const &conf)
+  int parse_params(std::string const &conf,
+                   colvarparse::Parse_Mode const parse_mode = colvarparse::parse_normal)
   {
     if (cvm::debug()) cvm::log("Reading grid configuration from string.\n");
 
@@ -746,30 +747,33 @@ public:
 
     {
       size_t nd_in = 0;
+      // this is only used in state files
       colvarparse::get_keyval(conf, "n_colvars", nd_in, nd, colvarparse::parse_silent);
       if (nd_in != nd) {
         cvm::error("Error: trying to read data for a grid "
-                    "that contains a different number of colvars ("+
-                    cvm::to_str(nd_in)+") than the grid defined "
-                    "in the configuration file("+cvm::to_str(nd)+
-                    ").\n");
+                   "that contains a different number of colvars ("+
+                   cvm::to_str(nd_in)+") than the grid defined "
+                   "in the configuration file("+cvm::to_str(nd)+
+                   ").\n");
         return COLVARS_ERROR;
       }
     }
 
+    // underscore keywords are used in state file
     colvarparse::get_keyval(conf, "lower_boundaries",
-                             lower_boundaries, lower_boundaries, colvarparse::parse_silent);
+                            lower_boundaries, lower_boundaries, colvarparse::parse_silent);
     colvarparse::get_keyval(conf, "upper_boundaries",
-                             upper_boundaries, upper_boundaries, colvarparse::parse_silent);
+                            upper_boundaries, upper_boundaries, colvarparse::parse_silent);
 
-    // support also camel case
+    // camel case keywords are used in config file
     colvarparse::get_keyval(conf, "lowerBoundaries",
-                             lower_boundaries, lower_boundaries, colvarparse::parse_silent);
+                            lower_boundaries, lower_boundaries, parse_mode);
     colvarparse::get_keyval(conf, "upperBoundaries",
-                             upper_boundaries, upper_boundaries, colvarparse::parse_silent);
+                            upper_boundaries, upper_boundaries, parse_mode);
 
-    colvarparse::get_keyval(conf, "widths", widths, widths, colvarparse::parse_silent);
+    colvarparse::get_keyval(conf, "widths", widths, widths, parse_mode);
 
+    // only used in state file
     colvarparse::get_keyval(conf, "sizes", nx, nx, colvarparse::parse_silent);
 
     if (nd < lower_boundaries.size()) nd = lower_boundaries.size();
@@ -808,13 +812,13 @@ public:
   {
     for (size_t i = 0; i < nd; i++) {
       if ( (std::sqrt(cv[i]->dist2(cv[i]->lower_boundary,
-                                     lower_boundaries[i])) > 1.0E-10) ||
+                                   lower_boundaries[i])) > 1.0E-10) ||
            (std::sqrt(cv[i]->dist2(cv[i]->upper_boundary,
-                                     upper_boundaries[i])) > 1.0E-10) ||
+                                   upper_boundaries[i])) > 1.0E-10) ||
            (std::sqrt(cv[i]->dist2(cv[i]->width,
-                                     widths[i])) > 1.0E-10) ) {
+                                   widths[i])) > 1.0E-10) ) {
         cvm::error("Error: restart information for a grid is "
-                    "inconsistent with that of its colvars.\n");
+                   "inconsistent with that of its colvars.\n");
         return;
       }
     }
@@ -830,19 +834,19 @@ public:
       // matter: boundaries should be EXACTLY the same (otherwise,
       // map_grid() should be used)
       if ( (std::fabs(other_grid.lower_boundaries[i] -
-                       lower_boundaries[i]) > 1.0E-10) ||
+                      lower_boundaries[i]) > 1.0E-10) ||
            (std::fabs(other_grid.upper_boundaries[i] -
-                       upper_boundaries[i]) > 1.0E-10) ||
+                      upper_boundaries[i]) > 1.0E-10) ||
            (std::fabs(other_grid.widths[i] -
-                       widths[i]) > 1.0E-10) ||
+                      widths[i]) > 1.0E-10) ||
            (data.size() != other_grid.data.size()) ) {
-      cvm::error("Error: inconsistency between "
-                  "two grids that are supposed to be equal, "
-                  "aside from the data stored.\n");
-      return;
+        cvm::error("Error: inconsistency between "
+                   "two grids that are supposed to be equal, "
+                   "aside from the data stored.\n");
+        return;
+      }
     }
   }
-}
 
 
   /// \brief Read grid entry in restart file
@@ -853,7 +857,7 @@ public:
     if ((is >> key) && (key == std::string("grid_parameters"))) {
       is.seekg(start_pos, std::ios::beg);
       is >> colvarparse::read_block("grid_parameters", conf);
-      parse_params(conf);
+      parse_params(conf, colvarparse::parse_silent);
     } else {
       cvm::log("Grid parameters are missing in the restart file, using those from the configuration.\n");
       is.seekg(start_pos, std::ios::beg);
@@ -871,11 +875,11 @@ public:
   }
 
 
-/// \brief Write the grid data without labels, as they are
-/// represented in memory
-/// \param buf_size Number of values per line
+  /// \brief Write the grid data without labels, as they are
+  /// represented in memory
+  /// \param buf_size Number of values per line
   std::ostream & write_raw(std::ostream &os,
-                            size_t const buf_size = 3)
+                           size_t const buf_size = 3)
   {
     std::streamsize const w = os.width();
     std::streamsize const p = os.precision();
@@ -935,10 +939,10 @@ public:
     os << std::setw(2) << "# " << nd << "\n";
     for (size_t i = 0; i < nd; i++) {
       os << "# "
-        << std::setw(10) << lower_boundaries[i]
-        << std::setw(10) << widths[i]
-        << std::setw(10) << nx[i] << "  "
-        << periodic[i] << "\n";
+         << std::setw(10) << lower_boundaries[i]
+         << std::setw(10) << widths[i]
+         << std::setw(10) << nx[i] << "  "
+         << periodic[i] << "\n";
     }
 
 
@@ -951,14 +955,14 @@ public:
 
       for (size_t i = 0; i < nd; i++) {
         os << " "
-          << std::setw(w) << std::setprecision(p)
-          << bin_to_value_scalar(ix[i], i);
+           << std::setw(w) << std::setprecision(p)
+           << bin_to_value_scalar(ix[i], i);
       }
       os << " ";
       for (size_t imult = 0; imult < mult; imult++) {
         os << " "
-          << std::setw(w) << std::setprecision(p)
-          << value_output(ix, imult);
+           << std::setw(w) << std::setprecision(p)
+           << value_output(ix, imult);
       }
       os << "\n";
     }
@@ -986,7 +990,7 @@ public:
 
     if ( !(is >> hash) || (hash != "#") ) {
       cvm::error("Error reading grid at position "+
-                  cvm::to_str(is.tellg())+" in stream(read \"" + hash + "\")\n");
+                 cvm::to_str(is.tellg())+" in stream(read \"" + hash + "\")\n");
       return is;
     }
 
@@ -1008,7 +1012,7 @@ public:
     for (size_t i = 0; i < nd; i++ ) {
       if ( !(is >> hash) || (hash != "#") ) {
         cvm::error("Error reading grid at position "+
-                    cvm::to_str(is.tellg())+" in stream(read \"" + hash + "\")\n");
+                   cvm::to_str(is.tellg())+" in stream(read \"" + hash + "\")\n");
         return is;
       }
 
@@ -1016,10 +1020,10 @@ public:
 
 
       if ( (std::fabs(lower - lower_boundaries[i].real_value) > 1.0e-10) ||
-          (std::fabs(width - widths[i] ) > 1.0e-10) ||
-          (nx_read[i] != nx[i]) ) {
+           (std::fabs(width - widths[i] ) > 1.0e-10) ||
+           (nx_read[i] != nx[i]) ) {
         cvm::log("Warning: reading from different grid definition (colvar "
-                  + cvm::to_str(i+1) + "); remapping data on new grid.\n");
+                 + cvm::to_str(i+1) + "); remapping data on new grid.\n");
         remap = true;
       }
     }
@@ -1063,7 +1067,6 @@ public:
 
   /// \brief Write the grid data without labels, as they are
   /// represented in memory
-  /// \param buf_size Number of values per line
   std::ostream & write_opendx(std::ostream &os)
   {
     // write the header
@@ -1122,11 +1125,11 @@ public:
 
   /// Constructor
   colvar_grid_count(std::vector<int> const &nx_i,
-                     size_t const           &def_count = 0);
+                    size_t const           &def_count = 0);
 
   /// Constructor from a vector of colvars
   colvar_grid_count(std::vector<colvar *>  &colvars,
-                     size_t const           &def_count = 0);
+                    size_t const           &def_count = 0);
 
   /// Increment the counter at given position
   inline void incr_count(std::vector<int> const &ix)
@@ -1136,7 +1139,7 @@ public:
 
   /// \brief Get the binned count indexed by ix from the newly read data
   inline size_t const & new_count(std::vector<int> const &ix,
-                                   size_t const &imult = 0)
+                                  size_t const &imult = 0)
   {
     return new_data[address(ix) + imult];
   }
@@ -1145,9 +1148,9 @@ public:
   /// into the internal representation (it may have been rescaled or
   /// manipulated)
   virtual inline void value_input(std::vector<int> const &ix,
-                                   size_t const &t,
-                                   size_t const &imult = 0,
-                                   bool add = false)
+                                  size_t const &t,
+                                  size_t const &imult = 0,
+                                  bool add = false)
   {
     if (add) {
       data[address(ix)] += t;
@@ -1164,7 +1167,7 @@ public:
   /// \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)
+                                                   int n = 0)
   {
     cvm::real A0, A1;
     std::vector<int> ix;
@@ -1377,7 +1380,7 @@ public:
   /// \brief Return the value of the function at ix divided by its
   /// number of samples (if the count grid is defined)
   virtual inline cvm::real value_output(std::vector<int> const &ix,
-                                         size_t const &imult = 0)
+                                        size_t const &imult = 0)
   {
     if (samples)
       return (samples->value(ix) > 0) ?
@@ -1391,9 +1394,9 @@ public:
   /// into the internal representation (it may have been rescaled or
   /// manipulated)
   virtual inline void value_input(std::vector<int> const &ix,
-                                   cvm::real const &new_value,
-                                   size_t const &imult = 0,
-                                   bool add = false)
+                                  cvm::real const &new_value,
+                                  size_t const &imult = 0,
+                                  bool add = false)
   {
     if (add) {
       if (samples)
diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp
index 49a192c3ca..aaaa6fac0d 100644
--- a/lib/colvars/colvarmodule.cpp
+++ b/lib/colvars/colvarmodule.cpp
@@ -293,6 +293,9 @@ int colvarmodule::parse_biases(std::string const &conf)
   /// initialize histograms
   parse_biases_type<colvarbias_histogram>(conf, "histogram", n_histo_biases);
 
+  /// initialize histogram restraints
+  parse_biases_type<colvarbias_restraint_histogram>(conf, "histogramRestraint", n_rest_biases);
+
   /// initialize linear restraints
   parse_biases_type<colvarbias_restraint_linear>(conf, "linear", n_rest_biases);
 
diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h
index 7f40e2548e..b33d3ef1a7 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-09-14"
+#define COLVARS_VERSION "2016-09-30"
 #endif
 
 #ifndef COLVARS_DEBUG
diff --git a/lib/colvars/colvarscript.cpp b/lib/colvars/colvarscript.cpp
index dc3619fe5e..86b52ac819 100644
--- a/lib/colvars/colvarscript.cpp
+++ b/lib/colvars/colvarscript.cpp
@@ -243,11 +243,17 @@ int colvarscript::proc_colvar(int argc, char const *argv[]) {
   }
 
   if (subcmd == "getappliedforce") {
-    result = (cv->bias_force()).to_simple_string();
+    result = (cv->applied_force()).to_simple_string();
     return COLVARS_OK;
   }
 
   if (subcmd == "getsystemforce") {
+    // TODO warning here
+    result = (cv->total_force()).to_simple_string();
+    return COLVARS_OK;
+  }
+
+  if (subcmd == "gettotalforce") {
     result = (cv->total_force()).to_simple_string();
     return COLVARS_OK;
   }
diff --git a/lib/colvars/colvartypes.h b/lib/colvars/colvartypes.h
index 65259c7220..2c3d49f686 100644
--- a/lib/colvars/colvartypes.h
+++ b/lib/colvars/colvartypes.h
@@ -57,6 +57,12 @@ public:
     }
   }
 
+  /// Return a reference to the data
+  inline std::vector<T> &data_array()
+  {
+    return data;
+  }
+
   inline ~vector1d()
   {
     data.clear();
@@ -203,6 +209,16 @@ public:
     return std::sqrt(this->norm2());
   }
 
+  inline cvm::real sum() const
+  {
+    cvm::real result = 0.0;
+    size_t i;
+    for (i = 0; i < this->size(); i++) {
+      result += (*this)[i];
+    }
+    return result;
+  }
+
   /// Slicing
   inline vector1d<T> const slice(size_t const i1, size_t const i2) const
   {
@@ -295,11 +311,23 @@ public:
   {
     std::stringstream stream(s);
     size_t i = 0;
-    while ((stream >> (*this)[i]) && (i < this->size())) {
-      i++;
-    }
-    if (i < this->size()) {
-      return COLVARS_ERROR;
+    if (this->size()) {
+      while ((stream >> (*this)[i]) && (i < this->size())) {
+        i++;
+      }
+      if (i < this->size()) {
+        return COLVARS_ERROR;
+      }
+    } else {
+      T input;
+      while (stream >> input) {
+        if ((i % 100) == 0) {
+          data.reserve(data.size()+100);
+        }
+        data.resize(data.size()+1);
+        data[i] = input;
+        i++;
+      }
     }
     return COLVARS_OK;
   }
@@ -434,6 +462,12 @@ public:
     this->clear();
   }
 
+  /// Return a reference to the data
+  inline std::vector<T> &data_array()
+  {
+    return data;
+  }
+
   inline row & operator [] (size_t const i)
   {
     return rows[i];