diff --git a/doc/src/PDF/colvars-refman-lammps.pdf b/doc/src/PDF/colvars-refman-lammps.pdf
index 0dd82d4794..071bdb25b8 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/Makefile.fermi b/lib/colvars/Makefile.fermi
index e17798bbc3..906675ae12 100644
--- a/lib/colvars/Makefile.fermi
+++ b/lib/colvars/Makefile.fermi
@@ -93,8 +93,8 @@ colvarcomp_rotations.o: colvarcomp_rotations.cpp colvarmodule.h \
 colvar.o: colvar.cpp colvarmodule.h colvartypes.h colvarproxy.h \
  colvarvalue.h colvarparse.h colvardeps.h colvar.h colvarcomp.h \
  colvaratoms.h colvarscript.h colvarbias.h
-colvardeps.o: colvardeps.cpp colvarmodule.h colvartypes.h colvarproxy.h \
- colvarvalue.h colvardeps.h
+colvardeps.o: colvardeps.cpp colvardeps.h colvarmodule.h colvartypes.h \
+ colvarproxy.h colvarvalue.h colvarparse.h
 colvargrid.o: colvargrid.cpp colvarmodule.h colvartypes.h colvarproxy.h \
  colvarvalue.h colvarparse.h colvardeps.h colvar.h colvarcomp.h \
  colvaratoms.h colvargrid.h
@@ -104,12 +104,12 @@ colvarmodule.o: colvarmodule.cpp colvarmodule.h colvartypes.h \
  colvarbias_restraint.h colvarbias_histogram.h colvarbias_meta.h \
  colvarscript.h
 colvarparse.o: colvarparse.cpp colvarmodule.h colvartypes.h colvarproxy.h \
- colvarvalue.h colvarparse.h colvardeps.h
+ colvarvalue.h colvarparse.h
 colvarscript.o: colvarscript.cpp colvarscript.h colvarmodule.h \
  colvartypes.h colvarproxy.h colvarvalue.h colvarbias.h colvar.h \
  colvarparse.h colvardeps.h
 colvartypes.o: colvartypes.cpp colvarmodule.h colvartypes.h colvarproxy.h \
- colvarvalue.h colvarparse.h colvardeps.h
+ colvarvalue.h colvarparse.h
 colvarvalue.o: colvarvalue.cpp colvarmodule.h colvartypes.h colvarproxy.h \
  colvarvalue.h
 
diff --git a/lib/colvars/Makefile.g++ b/lib/colvars/Makefile.g++
index f1fa1734b9..c80fa1065e 100644
--- a/lib/colvars/Makefile.g++
+++ b/lib/colvars/Makefile.g++
@@ -92,8 +92,8 @@ colvarcomp_rotations.o: colvarcomp_rotations.cpp colvarmodule.h \
 colvar.o: colvar.cpp colvarmodule.h colvartypes.h colvarproxy.h \
  colvarvalue.h colvarparse.h colvardeps.h colvar.h colvarcomp.h \
  colvaratoms.h colvarscript.h colvarbias.h
-colvardeps.o: colvardeps.cpp colvarmodule.h colvartypes.h colvarproxy.h \
- colvarvalue.h colvardeps.h
+colvardeps.o: colvardeps.cpp colvardeps.h colvarmodule.h colvartypes.h \
+ colvarproxy.h colvarvalue.h colvarparse.h
 colvargrid.o: colvargrid.cpp colvarmodule.h colvartypes.h colvarproxy.h \
  colvarvalue.h colvarparse.h colvardeps.h colvar.h colvarcomp.h \
  colvaratoms.h colvargrid.h
@@ -103,12 +103,12 @@ colvarmodule.o: colvarmodule.cpp colvarmodule.h colvartypes.h \
  colvarbias_restraint.h colvarbias_histogram.h colvarbias_meta.h \
  colvarscript.h
 colvarparse.o: colvarparse.cpp colvarmodule.h colvartypes.h colvarproxy.h \
- colvarvalue.h colvarparse.h colvardeps.h
+ colvarvalue.h colvarparse.h
 colvarscript.o: colvarscript.cpp colvarscript.h colvarmodule.h \
  colvartypes.h colvarproxy.h colvarvalue.h colvarbias.h colvar.h \
  colvarparse.h colvardeps.h
 colvartypes.o: colvartypes.cpp colvarmodule.h colvartypes.h colvarproxy.h \
- colvarvalue.h colvarparse.h colvardeps.h
+ colvarvalue.h colvarparse.h
 colvarvalue.o: colvarvalue.cpp colvarmodule.h colvartypes.h colvarproxy.h \
  colvarvalue.h
 
diff --git a/lib/colvars/Makefile.mingw32-cross b/lib/colvars/Makefile.mingw32-cross
index cea3a53c25..eba83c555f 100644
--- a/lib/colvars/Makefile.mingw32-cross
+++ b/lib/colvars/Makefile.mingw32-cross
@@ -100,8 +100,8 @@ $(DIR)colvarcomp_rotations.o: colvarcomp_rotations.cpp colvarmodule.h \
 $(DIR)colvar.o: colvar.cpp colvarmodule.h colvartypes.h colvarproxy.h \
  colvarvalue.h colvarparse.h colvardeps.h colvar.h colvarcomp.h \
  colvaratoms.h colvarscript.h colvarbias.h
-$(DIR)colvardeps.o: colvardeps.cpp colvarmodule.h colvartypes.h colvarproxy.h \
- colvarvalue.h colvardeps.h
+$(DIR)colvardeps.o: colvardeps.cpp colvardeps.h colvarmodule.h colvartypes.h \
+ colvarproxy.h colvarvalue.h colvarparse.h
 $(DIR)colvargrid.o: colvargrid.cpp colvarmodule.h colvartypes.h colvarproxy.h \
  colvarvalue.h colvarparse.h colvardeps.h colvar.h colvarcomp.h \
  colvaratoms.h colvargrid.h
@@ -111,12 +111,12 @@ $(DIR)colvarmodule.o: colvarmodule.cpp colvarmodule.h colvartypes.h \
  colvarbias_restraint.h colvarbias_histogram.h colvarbias_meta.h \
  colvarscript.h
 $(DIR)colvarparse.o: colvarparse.cpp colvarmodule.h colvartypes.h colvarproxy.h \
- colvarvalue.h colvarparse.h colvardeps.h
+ colvarvalue.h colvarparse.h
 $(DIR)colvarscript.o: colvarscript.cpp colvarscript.h colvarmodule.h \
  colvartypes.h colvarproxy.h colvarvalue.h colvarbias.h colvar.h \
  colvarparse.h colvardeps.h
 $(DIR)colvartypes.o: colvartypes.cpp colvarmodule.h colvartypes.h colvarproxy.h \
- colvarvalue.h colvarparse.h colvardeps.h
+ colvarvalue.h colvarparse.h
 $(DIR)colvarvalue.o: colvarvalue.cpp colvarmodule.h colvartypes.h colvarproxy.h \
  colvarvalue.h
 
diff --git a/lib/colvars/Makefile.mingw64-cross b/lib/colvars/Makefile.mingw64-cross
index 7dd85b48f8..1d83b6a0a8 100644
--- a/lib/colvars/Makefile.mingw64-cross
+++ b/lib/colvars/Makefile.mingw64-cross
@@ -100,8 +100,8 @@ $(DIR)colvarcomp_rotations.o: colvarcomp_rotations.cpp colvarmodule.h \
 $(DIR)colvar.o: colvar.cpp colvarmodule.h colvartypes.h colvarproxy.h \
  colvarvalue.h colvarparse.h colvardeps.h colvar.h colvarcomp.h \
  colvaratoms.h colvarscript.h colvarbias.h
-$(DIR)colvardeps.o: colvardeps.cpp colvarmodule.h colvartypes.h colvarproxy.h \
- colvarvalue.h colvardeps.h
+$(DIR)colvardeps.o: colvardeps.cpp colvardeps.h colvarmodule.h colvartypes.h \
+ colvarproxy.h colvarvalue.h colvarparse.h
 $(DIR)colvargrid.o: colvargrid.cpp colvarmodule.h colvartypes.h colvarproxy.h \
  colvarvalue.h colvarparse.h colvardeps.h colvar.h colvarcomp.h \
  colvaratoms.h colvargrid.h
@@ -111,12 +111,12 @@ $(DIR)colvarmodule.o: colvarmodule.cpp colvarmodule.h colvartypes.h \
  colvarbias_restraint.h colvarbias_histogram.h colvarbias_meta.h \
  colvarscript.h
 $(DIR)colvarparse.o: colvarparse.cpp colvarmodule.h colvartypes.h colvarproxy.h \
- colvarvalue.h colvarparse.h colvardeps.h
+ colvarvalue.h colvarparse.h
 $(DIR)colvarscript.o: colvarscript.cpp colvarscript.h colvarmodule.h \
  colvartypes.h colvarproxy.h colvarvalue.h colvarbias.h colvar.h \
  colvarparse.h colvardeps.h
 $(DIR)colvartypes.o: colvartypes.cpp colvarmodule.h colvartypes.h colvarproxy.h \
- colvarvalue.h colvarparse.h colvardeps.h
+ colvarvalue.h colvarparse.h
 $(DIR)colvarvalue.o: colvarvalue.cpp colvarmodule.h colvartypes.h colvarproxy.h \
  colvarvalue.h
 
diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp
index 1dd46524ac..84880808af 100644
--- a/lib/colvars/colvar.cpp
+++ b/lib/colvars/colvar.cpp
@@ -771,7 +771,7 @@ int colvar::calc_cvcs(int first_cvc, size_t num_cvcs)
 
   error_code |= calc_cvc_values(first_cvc, num_cvcs);
   error_code |= calc_cvc_gradients(first_cvc, num_cvcs);
-  error_code |= calc_cvc_sys_forces(first_cvc, num_cvcs);
+  error_code |= calc_cvc_total_force(first_cvc, num_cvcs);
   error_code |= calc_cvc_Jacobians(first_cvc, num_cvcs);
 
   if (cvm::debug())
@@ -790,7 +790,7 @@ int colvar::collect_cvc_data()
 
   error_code |= collect_cvc_values();
   error_code |= collect_cvc_gradients();
-  error_code |= collect_cvc_sys_forces();
+  error_code |= collect_cvc_total_forces();
   error_code |= collect_cvc_Jacobians();
   error_code |= calc_colvar_properties();
 
@@ -889,21 +889,22 @@ int colvar::calc_cvc_gradients(int first_cvc, size_t num_cvcs)
   size_t const cvc_max_count = num_cvcs ? num_cvcs : num_active_cvcs();
   size_t i, cvc_count;
 
-  if (is_enabled(f_cv_gradient)) {
+  if (cvm::debug())
+    cvm::log("Calculating gradients of colvar \""+this->name+"\".\n");
 
-    if (cvm::debug())
-      cvm::log("Calculating gradients of colvar \""+this->name+"\".\n");
+  // calculate the gradients of each component
+  cvm::increase_depth();
+  for (i = first_cvc, cvc_count = 0;
+      (i < cvcs.size()) && (cvc_count < cvc_max_count);
+      i++) {
+    if (!cvcs[i]->is_enabled()) continue;
+    cvc_count++;
 
-    // calculate the gradients of each component
-    cvm::increase_depth();
-    for (i = first_cvc, cvc_count = 0;
-        (i < cvcs.size()) && (cvc_count < cvc_max_count);
-        i++) {
-      if (!cvcs[i]->is_enabled()) continue;
-      cvc_count++;
+    if ((cvcs[i])->is_enabled(f_cvc_gradient)) {
       (cvcs[i])->calc_gradients();
       // if requested, propagate (via chain rule) the gradients above
       // to the atoms used to define the roto-translation
+      // This could be integrated in the CVC base class
       for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
         if (cvcs[i]->atom_groups[ig]->b_fit_gradients)
           cvcs[i]->atom_groups[ig]->calc_fit_gradients();
@@ -914,6 +915,7 @@ int colvar::calc_cvc_gradients(int first_cvc, size_t num_cvcs)
         }
       }
     }
+
     cvm::decrease_depth();
 
     if (cvm::debug())
@@ -976,7 +978,7 @@ int colvar::collect_cvc_gradients()
 }
 
 
-int colvar::calc_cvc_sys_forces(int first_cvc, size_t num_cvcs)
+int colvar::calc_cvc_total_force(int first_cvc, size_t num_cvcs)
 {
   size_t const cvc_max_count = num_cvcs ? num_cvcs : num_active_cvcs();
   size_t i, cvc_count;
@@ -1010,7 +1012,7 @@ int colvar::calc_cvc_sys_forces(int first_cvc, size_t num_cvcs)
 }
 
 
-int colvar::collect_cvc_sys_forces()
+int colvar::collect_cvc_total_forces()
 {
   if (is_enabled(f_cv_total_force_calc)) {
     ft.reset();
@@ -1097,10 +1099,8 @@ int colvar::calc_colvar_properties()
     // report the restraint center as "value"
     x_reported = xr;
     v_reported = vr;
-    // the "total force" with the extended Lagrangian is just the
-    // harmonic term acting on the extended coordinate
-    // Note: this is the force for current timestep
-    ft_reported = (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x);
+    // the "total force" with the extended Lagrangian is
+    // calculated in update_forces_energy() below
 
   } else {
 
@@ -1175,14 +1175,18 @@ cvm::real colvar::update_forces_energy()
 
     // 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)
+    // fr: bias force on extended variable (without harmonic spring), for output in trajectory
+    // f_ext: total force on extended variable (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);
 
+    // The total force acting on the extended variable is f_ext
+    // This will be used in the next timestep
+    ft_reported = f_ext;
+
     // 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
diff --git a/lib/colvars/colvar.h b/lib/colvars/colvar.h
index 665b3d6427..0058b56a26 100644
--- a/lib/colvars/colvar.h
+++ b/lib/colvars/colvar.h
@@ -254,7 +254,7 @@ public:
   /// \brief Same as \link colvar::calc_cvc_values \endlink but for gradients
   int calc_cvc_gradients(int first, size_t num_cvcs);
   /// \brief Same as \link colvar::calc_cvc_values \endlink but for total forces
-  int calc_cvc_sys_forces(int first, size_t num_cvcs);
+  int calc_cvc_total_force(int first, size_t num_cvcs);
   /// \brief Same as \link colvar::calc_cvc_values \endlink but for Jacobian derivatives/forces
   int calc_cvc_Jacobians(int first, size_t num_cvcs);
 
@@ -266,7 +266,7 @@ public:
   /// \brief Same as \link colvar::collect_cvc_values \endlink but for gradients
   int collect_cvc_gradients();
   /// \brief Same as \link colvar::collect_cvc_values \endlink but for total forces
-  int collect_cvc_sys_forces();
+  int collect_cvc_total_forces();
   /// \brief Same as \link colvar::collect_cvc_values \endlink but for Jacobian derivatives/forces
   int collect_cvc_Jacobians();
   /// \brief Calculate the quantities associated to the colvar (but not to the CVCs)
diff --git a/lib/colvars/colvarbias_abf.cpp b/lib/colvars/colvarbias_abf.cpp
index f2659951d4..870079a619 100644
--- a/lib/colvars/colvarbias_abf.cpp
+++ b/lib/colvars/colvarbias_abf.cpp
@@ -31,10 +31,8 @@ int colvarbias_abf::init(std::string const &conf)
 
   // ************* parsing general ABF options ***********************
 
-  get_keyval(conf, "applyBias",  apply_bias, true);
-  if (apply_bias) {
-    enable(f_cvb_apply_force);
-  } else {
+  get_keyval_feature((colvarparse *)this, conf, "applyBias",  f_cvb_apply_force, true);
+  if (!is_enabled(f_cvb_apply_force)){
     cvm::log("WARNING: ABF biases will *not* be applied!\n");
   }
 
@@ -84,9 +82,6 @@ int colvarbias_abf::init(std::string const &conf)
   // Request calculation of total force (which also checks for availability)
     if(enable(f_cvb_get_total_force)) return cvm::get_error();
   }
-  if (apply_bias) {
-    if(enable(f_cvb_apply_force)) return cvm::get_error();
-  }
 
   bool b_extended = false;
   for (size_t i = 0; i < colvars.size(); i++) {
@@ -249,6 +244,11 @@ int colvarbias_abf::update()
         // and subtract previous ABF force
         system_force[i] = colvars[i]->total_force().real_value
                         - colvar_forces[i].real_value;
+//         if (cvm::debug())
+//           cvm::log("ABF System force calc: cv " + cvm::to_str(i) +
+//                   " fs " + cvm::to_str(system_force[i]) +
+//                   " = ft " + cvm::to_str(colvars[i]->total_force().real_value) +
+//                   " - fa " + cvm::to_str(colvar_forces[i].real_value));
       }
       gradients->acc_force(force_bin, system_force);
     }
@@ -277,7 +277,7 @@ int colvarbias_abf::update()
   }
 
   // Compute and apply the new bias, if applicable
-  if ( apply_bias && samples->index_ok(bin) ) {
+  if (is_enabled(f_cvb_apply_force) && samples->index_ok(bin)) {
 
     size_t  count = samples->value(bin);
     cvm::real	fact = 1.0;
diff --git a/lib/colvars/colvarbias_abf.h b/lib/colvars/colvarbias_abf.h
index 8db0b75632..93079220ba 100644
--- a/lib/colvars/colvarbias_abf.h
+++ b/lib/colvars/colvarbias_abf.h
@@ -32,7 +32,6 @@ private:
   /// Base filename(s) for reading previous gradient data (replaces data from restart file)
   std::vector<std::string> input_prefix;
 
-  bool		apply_bias;
   bool		update_bias;
   bool		hide_Jacobian;
   size_t	full_samples;
diff --git a/lib/colvars/colvarcomp.cpp b/lib/colvars/colvarcomp.cpp
index 632e7fabda..534a927fda 100644
--- a/lib/colvars/colvarcomp.cpp
+++ b/lib/colvars/colvarcomp.cpp
@@ -43,7 +43,8 @@ colvar::cvc::cvc(std::string const &conf)
 
   // All cvcs implement this
   provide(f_cvc_debug_gradient);
-  get_keyval(conf, "debugGradients", set_feature(f_cvc_debug_gradient), false, parse_silent);
+  get_keyval_feature((colvarparse *)this, conf, "debugGradients",
+                     f_cvc_debug_gradient, false, parse_silent);
 
   // Attempt scalable calculations when in parallel? (By default yes, if available)
   get_keyval(conf, "scalable", b_try_scalable, true);
diff --git a/lib/colvars/colvardeps.cpp b/lib/colvars/colvardeps.cpp
index 7d8d05a2be..7dce7e31f6 100644
--- a/lib/colvars/colvardeps.cpp
+++ b/lib/colvars/colvardeps.cpp
@@ -1,4 +1,3 @@
-#include "colvarmodule.h"
 #include "colvardeps.h"
 
 
@@ -30,6 +29,18 @@ void colvardeps::provide(int feature_id) {
 }
 
 
+bool colvardeps::get_keyval_feature(colvarparse *cvp,
+                        std::string const &conf, char const *key,
+                        int feature_id, bool const &def_value,
+                        colvarparse::Parse_Mode const parse_mode)
+{
+  bool value;
+  bool const found = cvp->get_keyval(conf, key, value, def_value, parse_mode);
+  if (value) enable(feature_id);
+  return found;
+}
+
+
 int colvardeps::enable(int feature_id,
                       bool dry_run /* default: false */,
                       // dry_run: fail silently, do not enable if available
@@ -48,7 +59,7 @@ int colvardeps::enable(int feature_id,
   if (cvm::debug()) {
     cvm::log("DEPS: " + description +
       (dry_run ? " testing " : " requiring ") +
-      "\"" + f->description);
+      "\"" + f->description +"\"");
   }
 
   if (fs->enabled) {
@@ -140,11 +151,7 @@ int colvardeps::enable(int feature_id,
   // 4) solve deps in children
   for (i=0; i<f->requires_children.size(); i++) {
     int g = f->requires_children[i];
-    if (cvm::debug())
-      cvm::log("requires children " + features()[g]->description);
-//     cvm::log("children " + cvm::to_str(g));
     for (j=0; j<children.size(); j++) {
-//       cvm::log("child " +  children[j]->description);
       cvm::increase_depth();
       res = children[j]->enable(g, dry_run, false);
       cvm::decrease_depth();
@@ -215,7 +222,7 @@ void colvardeps::init_cvb_requires() {
   // Initialize feature_states for each instance
   feature_states.reserve(f_cvb_ntot);
   for (i = 0; i < f_cvb_ntot; i++) {
-    feature_states.push_back(new feature_state(this, feature_states.size(), true, false));
+    feature_states.push_back(new feature_state(true, false));
     // Most features are available, so we set them so
     // and list exceptions below
   }
@@ -319,7 +326,7 @@ void colvardeps::init_cv_requires() {
   // Initialize feature_states for each instance
   feature_states.reserve(f_cv_ntot);
   for (i = 0; i < f_cv_ntot; i++) {
-    feature_states.push_back(new feature_state(this, feature_states.size(), true, false));
+    feature_states.push_back(new feature_state(true, false));
     // Most features are available, so we set them so
     // and list exceptions below
    }
@@ -364,7 +371,7 @@ void colvardeps::init_cvc_requires() {
     f_description(f_cvc_debug_gradient, "debug gradient");
     f_req_self(f_cvc_debug_gradient, f_cvc_gradient);
 
-    f_description(f_cvc_Jacobian, "Jacobian");
+    f_description(f_cvc_Jacobian, "Jacobian derivative");
     f_req_self(f_cvc_Jacobian, f_cvc_inv_gradient);
 
     f_description(f_cvc_com_based, "depends on group centers of mass");
@@ -385,7 +392,7 @@ void colvardeps::init_cvc_requires() {
   // default as unavailable, not enabled
   feature_states.reserve(f_cvc_ntot);
   for (i = 0; i < colvardeps::f_cvc_ntot; i++) {
-    feature_states.push_back(new feature_state(this, feature_states.size(), false, false));
+    feature_states.push_back(new feature_state(false, false));
   }
 
   // Features that are implemented by all cvcs by default
@@ -429,7 +436,7 @@ void colvardeps::init_ag_requires() {
   // default as unavailable, not enabled
   feature_states.reserve(f_ag_ntot);
   for (i = 0; i < colvardeps::f_ag_ntot; i++) {
-    feature_states.push_back(new feature_state(this, feature_states.size(), false, false));
+    feature_states.push_back(new feature_state(false, false));
   }
 
   // Features that are implemented (or not) by all atom groups
diff --git a/lib/colvars/colvardeps.h b/lib/colvars/colvardeps.h
index a2f6f2d5f9..089461c532 100644
--- a/lib/colvars/colvardeps.h
+++ b/lib/colvars/colvardeps.h
@@ -1,10 +1,11 @@
 // -*- c++ -*-
 
-#include "colvarmodule.h"
-
 #ifndef COLVARDEPS_H
 #define COLVARDEPS_H
 
+#include "colvarmodule.h"
+#include "colvarparse.h"
+
 /// Parent class for a member object of a bias, cv or cvc etc. containing dependencies
 /// (features) and handling dependency resolution
 
@@ -27,16 +28,9 @@ public:
   std::string description; // reference to object name (cv, cvc etc.)
 
   /// This contains the current state of each feature for each object
-  class feature_state {
-  private:
-    colvardeps *const deps_object;
-    int const id;
-    operator int() { return 0; } // never cast as int
-  public:
-    inline colvardeps *object() const { return deps_object; }
-    inline int feature_id() const { return id; }
-    feature_state(colvardeps *o, int i, bool a, bool e)
-      : deps_object(o), id(i), available(a), enabled(e) {}
+  struct feature_state {
+    feature_state(bool a, bool e)
+    : available(a), enabled(e) {}
 
     /// Available means: supported, subject to dependencies as listed,
     /// MAY BE ENABLED AS A RESULT OF DEPENDENCY SOLVING
@@ -55,12 +49,6 @@ public:
   /// List of the state of all features
   std::vector<feature_state *> feature_states;
 
-  /// Allow setting a feature state while parsing its kewyord
-  inline feature_state * set_feature(int id)
-  {
-    return feature_states[id];
-  }
-
   /// Describes a feature and its dependecies
   /// used in a static array within each subclass
   class feature {
@@ -146,6 +134,12 @@ public:
 
   void provide(int feature_id); // set the feature's flag to available in local object
 
+  /// Parse a keyword and enable a feature accordingly
+  bool get_keyval_feature(colvarparse *cvp,
+                          std::string const &conf, char const *key,
+                          int feature_id, bool const &def_value,
+                          colvarparse::Parse_Mode const parse_mode = colvarparse::parse_normal);
+
   int enable(int f, bool dry_run = false, bool toplevel = true);  // enable a feature and recursively solve its dependencies
   // dry_run is set to true to recursively test if a feature is available, without enabling it
 //     int disable(int f);
diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h
index cb4272f0bd..8c61dcac2d 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-10"
+#define COLVARS_VERSION "2016-08-19"
 #endif
 
 #ifndef COLVARS_DEBUG
diff --git a/lib/colvars/colvarparse.cpp b/lib/colvars/colvarparse.cpp
index 47b1970220..8699391107 100644
--- a/lib/colvars/colvarparse.cpp
+++ b/lib/colvars/colvarparse.cpp
@@ -354,21 +354,6 @@ bool colvarparse::get_keyval(std::string const &conf,
 }
 
 
-bool colvarparse::get_keyval(std::string const &conf,
-                             char const *key,
-                             colvardeps::feature_state *value,
-                             bool const &def_value,
-                             Parse_Mode const parse_mode)
-{
-  bool feature_flag = def_value;
-  bool const b_found = get_keyval(conf, key, feature_flag, def_value, parse_mode);
-  if (feature_flag) {
-    value->object()->enable(value->feature_id());
-  }
-  return b_found;
-}
-
-
 // multiple-value keyword parsers
 
 bool colvarparse::get_keyval(std::string const &conf,
diff --git a/lib/colvars/colvarparse.h b/lib/colvars/colvarparse.h
index 7685b4eeb4..892f1632e1 100644
--- a/lib/colvars/colvarparse.h
+++ b/lib/colvars/colvarparse.h
@@ -8,7 +8,6 @@
 
 #include "colvarmodule.h"
 #include "colvarvalue.h"
-#include "colvardeps.h"
 
 
 /// \file colvarparse.h Parsing functions for collective variables
@@ -181,12 +180,6 @@ public:
                   bool &value,
                   bool const &def_value = false,
                   Parse_Mode const parse_mode = parse_normal);
-  bool get_keyval(std::string const &conf,
-                  char const *key,
-                  colvardeps::feature_state *value,
-                  bool const &def_value = false,
-                  Parse_Mode const parse_mode = parse_normal);
-
   bool get_keyval(std::string const &conf,
                   char const *key,
                   std::vector<int> &values,
diff --git a/lib/colvars/colvarscript.cpp b/lib/colvars/colvarscript.cpp
index 069a995135..dc3619fe5e 100644
--- a/lib/colvars/colvarscript.cpp
+++ b/lib/colvars/colvarscript.cpp
@@ -336,6 +336,11 @@ int colvarscript::proc_bias(int argc, char const *argv[]) {
     return COLVARS_OK;
   }
 
+  if (subcmd == "state") {
+    b->print_state();
+    return COLVARS_OK;
+  }
+
   // Subcommands for MW ABF
   if (subcmd == "bin") {
     int r = b->current_bin();