diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp
index ebc3bf44f1..ebf03b814a 100644
--- a/lib/colvars/colvar.cpp
+++ b/lib/colvars/colvar.cpp
@@ -649,7 +649,7 @@ int colvar::parse_analysis(std::string const &conf)
     } else {
       cvm::log("Unknown type of correlation function, \""+
                         acf_type_str+"\".\n");
-      cvm::set_error_bits(INPUT_ERROR);
+      cvm::set_error_bit(INPUT_ERROR);
     }
 
     get_keyval(conf, "corrFuncOffset", acf_offset, 0);
@@ -720,11 +720,11 @@ int colvar::calc()
   // Note: if anything is added here, it should be added also in the SMP block of calc_colvars()
   int error_code = COLVARS_OK;
   if (is_enabled(f_cv_active)) {
-    error_code |= update_cvc_flags();
-    error_code |= calc_cvcs();
-    error_code |= collect_cvc_data();
+    cvm::combine_errors(error_code, update_cvc_flags());
+    cvm::combine_errors(error_code, calc_cvcs());
+    cvm::combine_errors(error_code, collect_cvc_data());
   }
-  return COLVARS_OK;
+  return error_code;
 }
 
 
@@ -735,15 +735,15 @@ int colvar::calc_cvcs(int first_cvc, size_t num_cvcs)
     cvm::log("Calculating colvar \""+this->name+"\", components "+
              cvm::to_str(first_cvc)+" through "+cvm::to_str(first_cvc+num_cvcs)+".\n");
 
-  error_code |= check_cvc_range(first_cvc, num_cvcs);
+  cvm::combine_errors(error_code, check_cvc_range(first_cvc, num_cvcs));
   if (error_code != COLVARS_OK) {
     return error_code;
   }
 
-  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_Jacobians(first_cvc, num_cvcs);
+  cvm::combine_errors(error_code, calc_cvc_values(first_cvc, num_cvcs));
+  cvm::combine_errors(error_code, calc_cvc_gradients(first_cvc, num_cvcs));
+  cvm::combine_errors(error_code, calc_cvc_sys_forces(first_cvc, num_cvcs));
+  cvm::combine_errors(error_code, calc_cvc_Jacobians(first_cvc, num_cvcs));
 
   if (cvm::debug())
     cvm::log("Done calculating colvar \""+this->name+"\".\n");
@@ -759,12 +759,12 @@ int colvar::collect_cvc_data()
 
   int error_code = COLVARS_OK;
 
-  error_code |= collect_cvc_values();
-  error_code |= collect_cvc_gradients();
-  error_code |= collect_cvc_sys_forces();
-  error_code |= collect_cvc_Jacobians();
+  cvm::combine_errors(error_code, collect_cvc_values());
+  cvm::combine_errors(error_code, collect_cvc_gradients());
+  cvm::combine_errors(error_code, collect_cvc_sys_forces());
+  cvm::combine_errors(error_code, collect_cvc_Jacobians());
 
-  error_code |= calc_colvar_properties();
+  cvm::combine_errors(error_code, calc_colvar_properties());
 
   if (cvm::debug())
     cvm::log("Done calculating colvar \""+this->name+"\"'s properties.\n");
@@ -879,6 +879,11 @@ int colvar::calc_cvc_gradients(int first_cvc, size_t num_cvcs)
       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();
+
+        if (cvcs[i]->is_enabled(f_cvc_debug_gradient)) {
+          cvm::log("Debugging gradients for " + cvcs[i]->description);
+          cvcs[i]->debug_gradients(cvcs[i]->atom_groups[ig]);
+        }
       }
     }
     cvm::decrease_depth();
@@ -1297,7 +1302,7 @@ bool colvar::periodic_boundaries(colvarvalue const &lb, colvarvalue const &ub) c
   if ( (!is_enabled(f_cv_lower_boundary)) || (!is_enabled(f_cv_upper_boundary)) ) {
     cvm::log("Error: checking periodicity for collective variable \""+this->name+"\" "
                     "requires lower and upper boundaries to be defined.\n");
-    cvm::set_error_bits(INPUT_ERROR);
+    cvm::set_error_bit(INPUT_ERROR);
   }
 
   if (period > 0.0) {
diff --git a/lib/colvars/colvaratoms.cpp b/lib/colvars/colvaratoms.cpp
index f13f6993b6..ad7d8fcaf0 100644
--- a/lib/colvars/colvaratoms.cpp
+++ b/lib/colvars/colvaratoms.cpp
@@ -297,7 +297,7 @@ int cvm::atom_group::parse(std::string const &conf)
     std::string numbers_conf = "";
     size_t pos = 0;
     while (key_lookup(group_conf, "atomNumbers", numbers_conf, pos)) {
-      parse_error |= add_atom_numbers(numbers_conf);
+      cvm::combine_errors(parse_error, add_atom_numbers(numbers_conf));
       numbers_conf = "";
     }
   }
@@ -306,7 +306,7 @@ int cvm::atom_group::parse(std::string const &conf)
     std::string index_group_name;
     if (get_keyval(group_conf, "indexGroup", index_group_name)) {
       // use an index group from the index file read globally
-      parse_error |= add_index_group(index_group_name);
+      cvm::combine_errors(parse_error, add_index_group(index_group_name));
     }
   }
 
@@ -315,7 +315,7 @@ int cvm::atom_group::parse(std::string const &conf)
     size_t pos = 0;
     while (key_lookup(group_conf, "atomNumbersRange",
                       range_conf, pos)) {
-      parse_error |= add_atom_numbers_range(range_conf);
+      cvm::combine_errors(parse_error, add_atom_numbers_range(range_conf));
       range_conf = "";
     }
   }
@@ -342,8 +342,8 @@ int cvm::atom_group::parse(std::string const &conf)
         cvm::error("Error: more instances of \"atomNameResidueRange\" than "
                    "values of \"psfSegID\".\n", INPUT_ERROR);
       } else {
-        parse_error |= add_atom_name_residue_range(psf_segids.size() ? *psii : std::string(""),
-                                                   range_conf);
+        cvm::combine_errors(parse_error, add_atom_name_residue_range(psf_segids.size() ?
+        *psii : std::string(""), range_conf));
         if (psf_segids.size()) psii++;
       }
       range_conf = "";
@@ -407,7 +407,7 @@ int cvm::atom_group::parse(std::string const &conf)
     index = (cvm::proxy)->init_atom_group(atoms_ids);
   }
 
-  parse_error |= parse_fitting_options(group_conf);
+  cvm::combine_errors(parse_error, parse_fitting_options(group_conf));
 
   // TODO move this to colvarparse object
   check_keywords(group_conf, key.c_str());
@@ -612,7 +612,6 @@ int cvm::atom_group::parse_fitting_options(std::string const &group_conf)
 
       // regardless of the configuration, fit gradients must be calculated by refPositionsGroup
       ref_pos_group->b_fit_gradients = this->b_fit_gradients;
-      this->b_fit_gradients = false;
     }
 
     atom_group *group_for_fit = ref_pos_group ? ref_pos_group : this;
@@ -682,7 +681,7 @@ int cvm::atom_group::parse_fitting_options(std::string const &group_conf)
                "If that happens, use refPositionsGroup (or a different definition for it if already defined) "
                "to align the coordinates.\n");
       // initialize rot member data
-      rot.request_group1_gradients(this->size());
+      rot.request_group1_gradients(group_for_fit->size());
     }
   }
 
@@ -778,13 +777,19 @@ int cvm::atom_group::calc_required_properties()
 
 void cvm::atom_group::calc_apply_roto_translation()
 {
+  // store the laborarory-frame COGs for when they are needed later
+  cog_orig = this->center_of_geometry();
+  if (ref_pos_group) {
+    ref_pos_group->cog_orig = ref_pos_group->center_of_geometry();
+  }
+
   if (b_center) {
     // center on the origin first
-    cvm::atom_pos const cog = ref_pos_group ?
+    cvm::atom_pos const rpg_cog = ref_pos_group ?
       ref_pos_group->center_of_geometry() : this->center_of_geometry();
-    apply_translation(-1.0 * cog);
+    apply_translation(-1.0 * rpg_cog);
     if (ref_pos_group) {
-      ref_pos_group->apply_translation(-1.0 * cog);
+      ref_pos_group->apply_translation(-1.0 * rpg_cog);
     }
   }
 
@@ -945,26 +950,23 @@ void cvm::atom_group::calc_fit_gradients()
 {
   if (b_dummy) return;
 
-  if ((!b_center) && (!b_rotate)) return; // no fit
-
   if (cvm::debug())
     cvm::log("Calculating fit gradients.\n");
 
   atom_group *group_for_fit = ref_pos_group ? ref_pos_group : this;
-  group_for_fit->fit_gradients.assign(group_for_fit->size(), cvm::rvector(0.0, 0.0, 0.0));
 
   if (b_center) {
     // add the center of geometry contribution to the gradients
+    cvm::rvector atom_grad;
+
     for (size_t i = 0; i < this->size(); i++) {
-      // need to bring the gradients in original frame first
-      cvm::rvector const atom_grad = b_rotate ?
-        (rot.inverse()).rotate(atoms[i].grad) :
-        atoms[i].grad;
-      for (size_t j = 0; j < group_for_fit->size(); j++) {
-        group_for_fit->fit_gradients[j] +=
-          (-1.0)/(cvm::real(group_for_fit->size())) *
-          atom_grad;
-      }
+      atom_grad += atoms[i].grad;
+    }
+    if (b_rotate) atom_grad = (rot.inverse()).rotate(atom_grad);
+    atom_grad *= (-1.0)/(cvm::real(group_for_fit->size()));
+
+    for (size_t j = 0; j < group_for_fit->size(); j++) {
+      group_for_fit->fit_gradients[j] = atom_grad;
     }
   }
 
@@ -972,23 +974,26 @@ void cvm::atom_group::calc_fit_gradients()
 
     // add the rotation matrix contribution to the gradients
     cvm::rotation const rot_inv = rot.inverse();
-    cvm::atom_pos const cog = center_of_geometry();
 
     for (size_t i = 0; i < this->size(); i++) {
 
-      cvm::atom_pos const pos_orig = rot_inv.rotate((b_center ? (atoms[i].pos - cog) : (atoms[i].pos)));
+      // compute centered, unrotated position
+      cvm::atom_pos const pos_orig =
+        rot_inv.rotate((b_center ? (atoms[i].pos - ref_pos_cog) : (atoms[i].pos)));
+
+      // calculate \partial(R(q) \vec{x}_i)/\partial q) \cdot \partial\xi/\partial\vec{x}_i
+      cvm::quaternion const dxdq =
+        rot.q.position_derivative_inner(pos_orig, atoms[i].grad);
 
       for (size_t j = 0; j < group_for_fit->size(); j++) {
-        // calculate \partial(R(q) \vec{x}_i)/\partial q) \cdot \partial\xi/\partial\vec{x}_i
-        cvm::quaternion const dxdq =
-          rot.q.position_derivative_inner(pos_orig, atoms[i].grad);
-        // multiply by \cdot {\partial q}/\partial\vec{x}_j and add it to the fit gradients
+        // multiply by {\partial q}/\partial\vec{x}_j and add it to the fit gradients
         for (size_t iq = 0; iq < 4; iq++) {
           group_for_fit->fit_gradients[j] += dxdq[iq] * rot.dQ0_1[j][iq];
         }
       }
     }
   }
+
   if (cvm::debug())
     cvm::log("Done calculating fit gradients.\n");
 }
@@ -1138,17 +1143,9 @@ void cvm::atom_group::apply_colvar_force(cvm::real const &force)
 
     atom_group *group_for_fit = ref_pos_group ? ref_pos_group : this;
 
-    // add the contribution from the roto-translational fit to the gradients
-    if (b_rotate) {
-      // rotate forces back to the original frame
-      cvm::rotation const rot_inv = rot.inverse();
-      for (size_t j = 0; j < group_for_fit->size(); j++) {
-        (*group_for_fit)[j].apply_force(rot_inv.rotate(force * group_for_fit->fit_gradients[j]));
-      }
-    } else {
-      for (size_t j = 0; j < group_for_fit->size(); j++) {
-        (*group_for_fit)[j].apply_force(force * group_for_fit->fit_gradients[j]);
-      }
+    // Fit gradients are already calculated in "laboratory" frame
+    for (size_t j = 0; j < group_for_fit->size(); j++) {
+      (*group_for_fit)[j].apply_force(force * group_for_fit->fit_gradients[j]);
     }
   }
 
diff --git a/lib/colvars/colvaratoms.h b/lib/colvars/colvaratoms.h
index 827601fc4e..518ae53479 100644
--- a/lib/colvars/colvaratoms.h
+++ b/lib/colvars/colvaratoms.h
@@ -286,7 +286,7 @@ public:
   bool b_user_defined_fit;
 
   /// \brief Whether or not the derivatives of the roto-translation
-  /// should be included when calculating the colvar's gradients (default: no)
+  /// should be included when calculating the colvar's gradients (default: yes)
   bool b_fit_gradients;
 
   /// \brief use reference coordinates for b_center or b_rotate
@@ -359,10 +359,17 @@ public:
   /// \brief Calculate the center of geometry of the atomic positions, assuming
   /// that they are already pbc-wrapped
   int calc_center_of_geometry();
+
 private:
+
   /// \brief Center of geometry
   cvm::atom_pos cog;
+
+  /// \brief Center of geometry before any fitting
+  cvm::atom_pos cog_orig;
+
 public:
+
   /// \brief Return the center of geometry of the atomic positions
   inline cvm::atom_pos center_of_geometry() const
   {
diff --git a/lib/colvars/colvarbias_histogram.cpp b/lib/colvars/colvarbias_histogram.cpp
index bb55669614..ea91a5a55b 100644
--- a/lib/colvars/colvarbias_histogram.cpp
+++ b/lib/colvars/colvarbias_histogram.cpp
@@ -95,7 +95,7 @@ int colvarbias_histogram::update()
 {
   int error_code = COLVARS_OK;
   // update base class
-  error_code |= colvarbias::update();
+  cvm::combine_errors(error_code, colvarbias::update());
 
   if (cvm::debug()) {
     cvm::log("Updating histogram bias " + this->name);
@@ -148,7 +148,8 @@ int colvarbias_histogram::update()
     write_output_files();
   }
 
-  return error_code | cvm::get_error();
+  cvm::combine_errors(error_code, cvm::get_error());
+  return error_code;
 }
 
 
diff --git a/lib/colvars/colvarcomp.cpp b/lib/colvars/colvarcomp.cpp
index 2f9ea78398..ef64989ad1 100644
--- a/lib/colvars/colvarcomp.cpp
+++ b/lib/colvars/colvarcomp.cpp
@@ -26,7 +26,12 @@ colvar::cvc::cvc(std::string const &conf)
 
   init_cvc_requires();
 
-  get_keyval(conf, "name", this->name, std::string(""), parse_silent);
+  if (get_keyval(conf, "name", this->name, std::string(""), parse_silent)) {
+    // Temporary description until child object is initialized
+    description = "cvc " + name;
+  } else {
+    description = "uninitialized cvc";
+  }
 
   get_keyval(conf, "componentCoeff", sup_coeff, 1.0);
   get_keyval(conf, "componentExp", sup_np, 1);
@@ -34,6 +39,8 @@ colvar::cvc::cvc(std::string const &conf)
   get_keyval(conf, "period", period, 0.0);
   get_keyval(conf, "wrapAround", wrap_center, 0.0);
 
+  // All cvcs implement this
+  provide(f_cvc_debug_gradient);
   {
     bool b_debug_gradient;
     get_keyval(conf, "debugGradients", b_debug_gradient, false, parse_silent);
@@ -157,6 +164,8 @@ void colvar::cvc::debug_gradients(cvm::atom_group *group)
 {
   // this function should work for any scalar variable:
   // the only difference will be the name of the atom group (here, "group")
+  // NOTE: this assumes that groups for this cvc are non-overlapping,
+  // since atom coordinates are modified only within the current group
 
   if (group->b_dummy) return;
 
@@ -168,31 +177,34 @@ void colvar::cvc::debug_gradients(cvm::atom_group *group)
 
   // cvm::log("gradients     = "+cvm::to_str (gradients)+"\n");
 
-  // it only makes sense to debug the fit gradients
-  // when the fitting group is the same as this group
-  if (group->b_rotate || group->b_center)
-    if (group->b_fit_gradients && (group->ref_pos_group == NULL)) {
-      group->calc_fit_gradients();
-      if (group->b_rotate) {
-        // fit_gradients are in the original frame, we should print them in the rotated frame
-        for (size_t j = 0; j < group->fit_gradients.size(); j++) {
-          group->fit_gradients[j] = rot_0.rotate(group->fit_gradients[j]);
-        }
-      }
-      cvm::log("fit_gradients = "+cvm::to_str(group->fit_gradients)+"\n");
-      if (group->b_rotate) {
-        for (size_t j = 0; j < group->fit_gradients.size(); j++) {
-          group->fit_gradients[j] = rot_inv.rotate(group->fit_gradients[j]);
-        }
+  cvm::atom_group *group_for_fit = group->ref_pos_group ? group->ref_pos_group : group;
+  cvm::atom_pos fit_gradient_sum, gradient_sum;
+
+  // print the values of the fit gradients
+  if (group->b_rotate || group->b_center) {
+    if (group->b_fit_gradients) {
+      size_t j;
+
+      // fit_gradients are in the simulation frame: we should print them in the rotated frame
+      cvm::log("Fit gradients:\n");
+      for (j = 0; j < group_for_fit->fit_gradients.size(); j++) {
+        cvm::log((group->ref_pos_group ? std::string("refPosGroup") : group->key) +
+                 "[" + cvm::to_str(j) + "] = " +
+                 (group->b_rotate ?
+                  cvm::to_str(rot_0.rotate(group_for_fit->fit_gradients[j])) :
+                  cvm::to_str(group_for_fit->fit_gradients[j])));
       }
     }
+  }
 
+  // debug the gradients
   for (size_t ia = 0; ia < group->size(); ia++) {
 
     // tests are best conducted in the unrotated (simulation) frame
-    cvm::rvector const atom_grad = group->b_rotate ?
-      rot_inv.rotate((*group)[ia].grad) :
-      (*group)[ia].grad;
+    cvm::rvector const atom_grad = (group->b_rotate ?
+                                    rot_inv.rotate((*group)[ia].grad) :
+                                    (*group)[ia].grad);
+    gradient_sum += atom_grad;
 
     for (size_t id = 0; id < 3; id++) {
       // (re)read original positions
@@ -205,55 +217,59 @@ void colvar::cvc::debug_gradients(cvm::atom_group *group)
       if ((x.type() == colvarvalue::type_vector) && (x.size() == 1)) x_1 = x[0];
       cvm::log("Atom "+cvm::to_str(ia)+", component "+cvm::to_str(id)+":\n");
       cvm::log("dx(actual) = "+cvm::to_str(x_1 - x_0,
-                             21, 14)+"\n");
-      //cvm::real const dx_pred = (group->fit_gradients.size() && (group->ref_pos_group == NULL)) ?
+                            21, 14)+"\n");
       cvm::real const dx_pred = (group->fit_gradients.size()) ?
         (cvm::debug_gradients_step_size * (atom_grad[id] + group->fit_gradients[ia][id])) :
         (cvm::debug_gradients_step_size * atom_grad[id]);
       cvm::log("dx(interp) = "+cvm::to_str(dx_pred,
-                             21, 14)+"\n");
+                            21, 14)+"\n");
       cvm::log("|dx(actual) - dx(interp)|/|dx(actual)| = "+
                 cvm::to_str(std::fabs(x_1 - x_0 - dx_pred) /
-                             std::fabs(x_1 - x_0), 12, 5)+"\n");
+                            std::fabs(x_1 - x_0), 12, 5)+"\n");
     }
   }
 
-/*
- * The code below is WIP
- */
-//   if (group->ref_pos_group != NULL) {
-//     cvm::atom_group &ref = *group->ref_pos_group;
-//     group->calc_fit_gradients();
-//
-//     for (size_t ia = 0; ia < ref.size(); ia++) {
-//
-//       for (size_t id = 0; id < 3; id++) {
-//         // (re)read original positions
-//         group->read_positions();
-//         ref.read_positions();
-//         // change one coordinate
-//         ref[ia].pos[id] += cvm::debug_gradients_step_size;
-//         group->update();
-//         calc_value();
-//         cvm::real const x_1 = x.real_value;
-//         cvm::log("refPosGroup atom "+cvm::to_str(ia)+", component "+cvm::to_str (id)+":\n");
-//         cvm::log("dx(actual) = "+cvm::to_str (x_1 - x_0,
-//                                21, 14)+"\n");
-//         //cvm::real const dx_pred = (group->fit_gradients.size() && (group->ref_pos_group == NULL)) ?
-//         // cvm::real const dx_pred = (group->fit_gradients.size()) ?
-//         //  (cvm::debug_gradients_step_size * (atom_grad[id] + group->fit_gradients[ia][id])) :
-//         //  (cvm::debug_gradients_step_size * atom_grad[id]);
-//         cvm::real const dx_pred = cvm::debug_gradients_step_size * ref.fit_gradients[ia][id];
-//         cvm::log("dx(interp) = "+cvm::to_str (dx_pred,
-//                                21, 14)+"\n");
-//         cvm::log ("|dx(actual) - dx(interp)|/|dx(actual)| = "+
-//                   cvm::to_str(std::fabs (x_1 - x_0 - dx_pred) /
-//                                std::fabs (x_1 - x_0),
-//                                12, 5)+
-//                   ".\n");
-//       }
-//     }
-//   }
+  if ((group->b_fit_gradients) && (group->ref_pos_group != NULL)) {
+    cvm::atom_group *ref_group = group->ref_pos_group;
+    group->read_positions();
+    group->calc_required_properties();
+
+    for (size_t ia = 0; ia < ref_group->size(); ia++) {
+
+      // fit gradients are in the unrotated (simulation) frame
+      cvm::rvector const atom_grad = ref_group->fit_gradients[ia];
+      fit_gradient_sum += atom_grad;
+
+      for (size_t id = 0; id < 3; id++) {
+        // (re)read original positions
+        group->read_positions();
+        ref_group->read_positions();
+        // change one coordinate
+        (*ref_group)[ia].pos[id] += cvm::debug_gradients_step_size;
+        group->calc_required_properties();
+        calc_value();
+
+        cvm::real const x_1 = x.real_value;
+        cvm::log("refPosGroup atom "+cvm::to_str(ia)+", component "+cvm::to_str (id)+":\n");
+        cvm::log("dx(actual) = "+cvm::to_str (x_1 - x_0,
+                               21, 14)+"\n");
+
+        cvm::real const dx_pred = cvm::debug_gradients_step_size * atom_grad[id];
+
+        cvm::log("dx(interp) = "+cvm::to_str (dx_pred,
+                               21, 14)+"\n");
+        cvm::log ("|dx(actual) - dx(interp)|/|dx(actual)| = "+
+                  cvm::to_str(std::fabs (x_1 - x_0 - dx_pred) /
+                               std::fabs (x_1 - x_0),
+                               12, 5)+
+                  ".\n");
+      }
+    }
+  }
+
+  cvm::log("Gradient sum: " +  cvm::to_str(gradient_sum) +
+        "  Fit gradient sum: " + cvm::to_str(fit_gradient_sum) +
+        "  Total " + cvm::to_str(gradient_sum + fit_gradient_sum));
 
   return;
 }
diff --git a/lib/colvars/colvarcomp_angles.cpp b/lib/colvars/colvarcomp_angles.cpp
index 519288a01f..1f825a17c4 100644
--- a/lib/colvars/colvarcomp_angles.cpp
+++ b/lib/colvars/colvarcomp_angles.cpp
@@ -84,12 +84,6 @@ void colvar::angle::calc_gradients()
   group1->set_weighted_gradient(dxdr1);
   group2->set_weighted_gradient((dxdr1 + dxdr3) * (-1.0));
   group3->set_weighted_gradient(dxdr3);
-
-  if (is_enabled(f_cvc_debug_gradient)) {
-    debug_gradients(group1);
-    debug_gradients(group2);
-    debug_gradients(group3);
-  }
 }
 
 void colvar::angle::calc_force_invgrads()
diff --git a/lib/colvars/colvarcomp_distances.cpp b/lib/colvars/colvarcomp_distances.cpp
index 8d15c5d67a..a38fa7b91f 100644
--- a/lib/colvars/colvarcomp_distances.cpp
+++ b/lib/colvars/colvarcomp_distances.cpp
@@ -167,7 +167,7 @@ colvar::distance_z::distance_z(std::string const &conf)
   // this group is optional
   ref2 = parse_group(conf, "ref2", true);
 
-  if (ref2->size()) {
+  if (ref2 && ref2->size()) {
     cvm::log("Using axis joining the centers of mass of groups \"ref\" and \"ref2\"");
     fixed_axis = false;
     if (key_lookup(conf, "axis"))
@@ -249,15 +249,6 @@ void colvar::distance_z::calc_gradients()
         cvm::position_distance(main->center_of_mass(), ref1->center_of_mass()) + x.real_value * axis ));
     }
   }
-
-  if (is_enabled(f_cvc_debug_gradient)) {
-    cvm::log("Debugging gradients for group main:\n");
-    debug_gradients(main);
-    cvm::log("Debugging gradients for group ref1:\n");
-    debug_gradients(ref1);
-    cvm::log("Debugging gradients for group ref2:\n");
-    debug_gradients(ref2);
-  }
 }
 
 void colvar::distance_z::calc_force_invgrads()
@@ -282,7 +273,7 @@ void colvar::distance_z::apply_force(colvarvalue const &force)
   if (!ref1->noforce)
     ref1->apply_colvar_force(force.real_value);
 
-  if (ref2->size() && !ref2->noforce)
+  if (ref2 && ref2->size() && !ref2->noforce)
     ref2->apply_colvar_force(force.real_value);
 
   if (!main->noforce)
@@ -382,7 +373,7 @@ void colvar::distance_xy::apply_force(colvarvalue const &force)
   if (!ref1->noforce)
     ref1->apply_colvar_force(force.real_value);
 
-  if (ref2->size() && !ref2->noforce)
+  if (ref2 && ref2->size() && !ref2->noforce)
     ref2->apply_colvar_force(force.real_value);
 
   if (!main->noforce)
@@ -593,10 +584,6 @@ void colvar::distance_pairs::calc_value()
 void colvar::distance_pairs::calc_gradients()
 {
   // will be calculated on the fly in apply_force()
-  if (is_enabled(f_cvc_debug_gradient)) {
-    cvm::log("Debugging gradients:\n");
-    debug_gradients(group1);
-  }
 }
 
 void colvar::distance_pairs::apply_force(colvarvalue const &force)
@@ -667,11 +654,6 @@ void colvar::gyration::calc_gradients()
   for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) {
     ai->grad = drdx * ai->pos;
   }
-
-  if (is_enabled(f_cvc_debug_gradient)) {
-    cvm::log("Debugging gradients:\n");
-    debug_gradients(atoms);
-  }
 }
 
 
@@ -731,11 +713,6 @@ void colvar::inertia::calc_gradients()
   for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) {
     ai->grad = 2.0 * ai->pos;
   }
-
-  if (is_enabled(f_cvc_debug_gradient)) {
-    cvm::log("Debugging gradients:\n");
-    debug_gradients(atoms);
-  }
 }
 
 
@@ -786,11 +763,6 @@ void colvar::inertia_z::calc_gradients()
   for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) {
     ai->grad = 2.0 * (ai->pos * axis) * axis;
   }
-
-  if (is_enabled(f_cvc_debug_gradient)) {
-    cvm::log("Debugging gradients:\n");
-    debug_gradients(atoms);
-  }
 }
 
 
@@ -811,7 +783,7 @@ colvar::rmsd::rmsd(std::string const &conf)
 
   atoms = parse_group(conf, "atoms");
 
-  if (atoms->size() == 0) {
+  if (!atoms || atoms->size() == 0) {
     cvm::error("Error: \"atoms\" must contain at least 1 atom to compute RMSD.");
     return;
   }
@@ -891,18 +863,11 @@ colvar::rmsd::rmsd(std::string const &conf)
     cvm::log("This is a standard minimum RMSD, derivatives of the optimal rotation "
               "will not be computed as they cancel out in the gradients.");
     atoms->b_fit_gradients = false;
-  }
-
-  if (atoms->b_rotate) {
-    // TODO: finer-grained control of this would require exposing a
-    // "request_Jacobian_derivative()" method to the colvar, and the same
-    // from the colvar to biases
-    // TODO: this should not be enabled here anyway, as it is not specific of the
-    // component - instead it should be decided in a generic way by the atom group
 
     // request the calculation of the derivatives of the rotation defined by the atom group
     atoms->rot.request_group1_gradients(atoms->size());
     // request derivatives of optimal rotation wrt reference coordinates for Jacobian:
+    // this is only required for ABF, but we do both groups here for better caching
     atoms->rot.request_group2_gradients(atoms->size());
   }
 }
@@ -930,11 +895,6 @@ void colvar::rmsd::calc_gradients()
   for (size_t ia = 0; ia < atoms->size(); ia++) {
     (*atoms)[ia].grad = (drmsddx2 * 2.0 * ((*atoms)[ia].pos - ref_pos[ia]));
   }
-
-  if (is_enabled(f_cvc_debug_gradient)) {
-    cvm::log("Debugging gradients:\n");
-    debug_gradients(atoms);
-  }
 }
 
 
@@ -1031,7 +991,7 @@ colvar::eigenvector::eigenvector(std::string const &conf)
       cvm::log("Using reference positions from input file.\n");
       if (ref_pos.size() != atoms->size()) {
         cvm::error("Error: reference positions do not "
-                          "match the number of requested atoms->\n");
+                   "match the number of requested atoms.\n");
         return;
       }
     }
@@ -1064,9 +1024,14 @@ colvar::eigenvector::eigenvector(std::string const &conf)
     }
   }
 
+  if (ref_pos.size() == 0) {
+    cvm::error("Error: reference positions were not provided.\n", INPUT_ERROR);
+    return;
+  }
+
   if (ref_pos.size() != atoms->size()) {
-    cvm::error("Error: reference positions were not provided, or do not "
-                      "match the number of requested atoms->\n");
+    cvm::error("Error: reference positions do not "
+               "match the number of requested atoms.\n", INPUT_ERROR);
     return;
   }
 
@@ -1087,6 +1052,8 @@ colvar::eigenvector::eigenvector(std::string const &conf)
     atoms->b_rotate = true;
     atoms->ref_pos = ref_pos;
     atoms->center_ref_pos();
+    atoms->b_fit_gradients = false; // cancel out if group is fitted on itself
+                                    // and cvc is translationally invariant
 
     // request the calculation of the derivatives of the rotation defined by the atom group
     atoms->rot.request_group1_gradients(atoms->size());
@@ -1216,11 +1183,6 @@ void colvar::eigenvector::calc_gradients()
   for (size_t ia = 0; ia < atoms->size(); ia++) {
     (*atoms)[ia].grad = eigenvec[ia];
   }
-
-  if (is_enabled(f_cvc_debug_gradient)) {
-    cvm::log("Debugging gradients:\n");
-    debug_gradients(atoms);
-  }
 }
 
 
@@ -1313,6 +1275,7 @@ colvar::cartesian::cartesian(std::string const &conf)
 
   if (axes.size() == 0) {
     cvm::error("Error: a \"cartesian\" component was defined with all three axes disabled.\n");
+    return;
   }
 
   x.type(colvarvalue::type_vector);
diff --git a/lib/colvars/colvarcomp_rotations.cpp b/lib/colvars/colvarcomp_rotations.cpp
index 505f95e072..b318b0f7be 100644
--- a/lib/colvars/colvarcomp_rotations.cpp
+++ b/lib/colvars/colvarcomp_rotations.cpp
@@ -57,14 +57,14 @@ colvar::orientation::orientation(std::string const &conf)
   cvm::log("Centering the reference coordinates: it is "
             "assumed that each atom is the closest "
             "periodic image to the center of geometry.\n");
-  cvm::rvector cog(0.0, 0.0, 0.0);
+  cvm::rvector ref_cog(0.0, 0.0, 0.0);
   size_t i;
   for (i = 0; i < ref_pos.size(); i++) {
-    cog += ref_pos[i];
+    ref_cog += ref_pos[i];
   }
-  cog /= cvm::real(ref_pos.size());
+  ref_cog /= cvm::real(ref_pos.size());
   for (i = 0; i < ref_pos.size(); i++) {
-    ref_pos[i] -= cog;
+    ref_pos[i] -= ref_cog;
   }
 
   get_keyval(conf, "closestToQuaternion", ref_quat, cvm::quaternion(1.0, 0.0, 0.0, 0.0));
@@ -87,6 +87,7 @@ colvar::orientation::orientation()
 
 void colvar::orientation::calc_value()
 {
+  rot.b_debug_gradients = is_enabled(f_cvc_debug_gradient);
   atoms_cog = atoms->center_of_geometry();
 
   rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog));
@@ -163,10 +164,6 @@ void colvar::orientation_angle::calc_gradients()
   for (size_t ia = 0; ia < atoms->size(); ia++) {
     (*atoms)[ia].grad = (dxdq0 * (rot.dQ0_2[ia])[0]);
   }
-  if (is_enabled(f_cvc_debug_gradient)) {
-    cvm::log("Debugging orientationAngle component gradients:\n");
-    debug_gradients(atoms);
-  }
 }
 
 
@@ -210,10 +207,6 @@ void colvar::orientation_proj::calc_gradients()
   for (size_t ia = 0; ia < atoms->size(); ia++) {
     (*atoms)[ia].grad = (dxdq0 * (rot.dQ0_2[ia])[0]);
   }
-  if (is_enabled(f_cvc_debug_gradient)) {
-    cvm::log("Debugging orientationProj component gradients:\n");
-    debug_gradients(atoms);
-  }
 }
 
 
@@ -272,11 +265,6 @@ void colvar::tilt::calc_gradients()
       (*atoms)[ia].grad += (dxdq[iq] * (rot.dQ0_2[ia])[iq]);
     }
   }
-
-  if (is_enabled(f_cvc_debug_gradient)) {
-    cvm::log("Debugging tilt component gradients:\n");
-    debug_gradients(atoms);
-  }
 }
 
 
diff --git a/lib/colvars/colvardeps.cpp b/lib/colvars/colvardeps.cpp
index f7d829968b..5821c8b24a 100644
--- a/lib/colvars/colvardeps.cpp
+++ b/lib/colvars/colvardeps.cpp
@@ -355,6 +355,9 @@ void cvm::deps::init_cvc_requires() {
     f_description(f_cvc_inv_gradient, "inverse gradient");
     f_req_self(f_cvc_inv_gradient, f_cvc_gradient);
 
+    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_req_self(f_cvc_Jacobian, f_cvc_inv_gradient);
 
diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp
index 9991513aca..b84f226ec5 100644
--- a/lib/colvars/colvarmodule.cpp
+++ b/lib/colvars/colvarmodule.cpp
@@ -316,7 +316,8 @@ int colvarmodule::parse_biases(std::string const &conf)
 int colvarmodule::catch_input_errors(int result)
 {
   if (result != COLVARS_OK || get_error()) {
-    set_error_bits(result | INPUT_ERROR);
+    set_error_bit(result);
+    set_error_bit(INPUT_ERROR);
     parse->reset();
     return get_error();
   }
@@ -468,27 +469,27 @@ int colvarmodule::calc()
              cvm::to_str(cvm::step_absolute())+"\n");
   }
 
-  error_code |= calc_colvars();
+  combine_errors(error_code, calc_colvars());
   // set biasing forces to zero before biases are calculated and summed over
   for (std::vector<colvar *>::iterator cvi = colvars.begin();
        cvi != colvars.end(); cvi++) {
     (*cvi)->reset_bias_force();
   }
-  error_code |= calc_biases();
-  error_code |= update_colvar_forces();
+  combine_errors(error_code, calc_biases());
+  combine_errors(error_code, update_colvar_forces());
 
   if (cvm::b_analysis) {
-    error_code |= analyze();
+    combine_errors(error_code, analyze());
   }
 
   // write trajectory files, if needed
   if (cv_traj_freq && cv_traj_name.size()) {
-    error_code |= write_traj_files();
+    combine_errors(error_code, write_traj_files());
   }
 
   // write restart files, if needed
   if (restart_out_freq && restart_out_name.size()) {
-    error_code |= write_restart_files();
+    combine_errors(error_code, write_restart_files());
   }
 
   return error_code;
@@ -520,7 +521,7 @@ int colvarmodule::calc_colvars()
     cvm::increase_depth();
     for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
 
-      error_code |= (*cvi)->update_cvc_flags();
+      combine_errors(error_code, (*cvi)->update_cvc_flags());
 
       size_t num_items = (*cvi)->num_active_cvcs();
       colvars_smp.reserve(colvars_smp.size() + num_items);
@@ -535,11 +536,11 @@ int colvarmodule::calc_colvars()
     cvm::decrease_depth();
 
     // calculate colvar components in parallel
-    error_code |= proxy->smp_colvars_loop();
+    combine_errors(error_code, proxy->smp_colvars_loop());
 
     cvm::increase_depth();
     for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
-      error_code |= (*cvi)->collect_cvc_data();
+      combine_errors(error_code, (*cvi)->collect_cvc_data());
     }
     cvm::decrease_depth();
 
@@ -548,7 +549,7 @@ int colvarmodule::calc_colvars()
     // calculate colvars one at a time
     cvm::increase_depth();
     for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
-      error_code |= (*cvi)->calc();
+      combine_errors(error_code, (*cvi)->calc());
       if (cvm::get_error()) {
         return COLVARS_ERROR;
       }
@@ -556,7 +557,8 @@ int colvarmodule::calc_colvars()
     cvm::decrease_depth();
   }
 
-  return error_code | cvm::get_error();
+  combine_errors(error_code, cvm::get_error());
+  return error_code;
 }
 
 
@@ -575,21 +577,21 @@ int colvarmodule::calc_biases()
 
     if (use_scripted_forces && !scripting_after_biases) {
       // calculate biases and scripted forces in parallel
-      error_code |= proxy->smp_biases_script_loop();
+      combine_errors(error_code, proxy->smp_biases_script_loop());
     } else {
       // calculate biases in parallel
-      error_code |= proxy->smp_biases_loop();
+      combine_errors(error_code, proxy->smp_biases_loop());
     }
 
   } else {
 
     if (use_scripted_forces && !scripting_after_biases) {
-      error_code |= calc_scripted_forces();
+      combine_errors(error_code, calc_scripted_forces());
     }
 
     cvm::increase_depth();
     for (bi = biases.begin(); bi != biases.end(); bi++) {
-      error_code |= (*bi)->update();
+      combine_errors(error_code, (*bi)->update());
       if (cvm::get_error()) {
         return COLVARS_ERROR;
       }
@@ -627,7 +629,7 @@ int colvarmodule::update_colvar_forces()
   cvm::decrease_depth();
 
   if (use_scripted_forces && scripting_after_biases) {
-    error_code |= calc_scripted_forces();
+    combine_errors(error_code, calc_scripted_forces());
   }
 
   cvm::real total_colvar_energy = 0.0;
@@ -875,17 +877,17 @@ int colvarmodule::setup_output()
      std::string(""));
 
   if (cv_traj_freq && cv_traj_name.size()) {
-    error_code |= open_traj_file(cv_traj_name);
+    combine_errors(error_code, open_traj_file(cv_traj_name));
   }
 
   for (std::vector<colvarbias *>::iterator bi = biases.begin();
        bi != biases.end();
        bi++) {
-    error_code |= (*bi)->setup_output();
+    combine_errors(error_code, (*bi)->setup_output());
   }
 
   if (error_code != COLVARS_OK || cvm::get_error()) {
-    set_error_bits(FILE_ERROR);
+    set_error_bit(FILE_ERROR);
   }
 
   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
@@ -904,6 +906,13 @@ std::istream & colvarmodule::read_restart(std::istream &is)
                           colvarparse::parse_silent);
         it = it_restart;
       }
+      std::string restart_version;
+      parse->get_keyval(restart_conf, "version",
+                        restart_version, std::string(""),
+                        colvarparse::parse_silent);
+      if (restart_version.size() && (restart_version != std::string(COLVARS_VERSION))) {
+        cvm::log("This state file was generated with version "+restart_version+"\n");
+      }
     }
     is.clear();
     parse->clear_keyword_registry();
@@ -1060,6 +1069,7 @@ std::ostream & colvarmodule::write_restart(std::ostream &os)
      << "  step " << std::setw(it_width)
      << it << "\n"
      << "  dt " << dt() << "\n"
+     << "  version " << std::string(COLVARS_VERSION) << "\n"
      << "}\n\n";
 
   cvm::increase_depth();
@@ -1213,14 +1223,30 @@ size_t & cvm::depth()
 }
 
 
-void colvarmodule::set_error_bits(int code)
+void colvarmodule::set_error_bit(int code)
 {
+  if (code > 0) {
+    cvm::fatal_error("Error: set_error_bit() received positive error code.\n");
+    return;
+  }
   proxy->smp_lock();
-  errorCode |= code;
-  errorCode |= COLVARS_ERROR;
+  errorCode = -1 * ((-1 * errorCode) | (-1 * code) | (-1 * COLVARS_ERROR));
   proxy->smp_unlock();
 }
 
+bool colvarmodule::get_error_bit(int code)
+{
+  return bool((-1 * errorCode) & (-1 * code));
+}
+
+void colvarmodule::combine_errors(int &target, int const code)
+{
+  if (code > 0 || target > 0) {
+    cvm::fatal_error("Error: combine_errors() received positive error code.\n");
+    return;
+  }
+  target = -1 * ((-1 * target) | (-1 * code));
+}
 
 void colvarmodule::clear_error()
 {
@@ -1232,7 +1258,7 @@ void colvarmodule::clear_error()
 
 void cvm::error(std::string const &message, int code)
 {
-  set_error_bits(code);
+  set_error_bit(code);
   proxy->error(message);
 }
 
@@ -1241,7 +1267,7 @@ void cvm::fatal_error(std::string const &message)
 {
   // TODO once all non-fatal errors have been set to be handled by error(),
   // set DELETE_COLVARS here for VMD to handle it
-  set_error_bits(FATAL_ERROR);
+  set_error_bit(FATAL_ERROR);
   proxy->fatal_error(message);
 }
 
@@ -1394,7 +1420,7 @@ colvarproxy              *colvarmodule::proxy = NULL;
 
 
 // static runtime data
-cvm::real colvarmodule::debug_gradients_step_size = 1.0e-03;
+cvm::real colvarmodule::debug_gradients_step_size = 1.0e-07;
 int       colvarmodule::errorCode = 0;
 long      colvarmodule::it = 0;
 long      colvarmodule::it_restart = 0;
diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h
index eac071805c..6f60814133 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-04-14"
+#define COLVARS_VERSION "2016-04-27"
 #endif
 
 #ifndef COLVARS_DEBUG
@@ -20,19 +20,19 @@
 /// shared between all object instances) to be accessed from other
 /// objects.
 
-// Error codes
+// Error codes are the negative of powers-of-two
+// as a result, error codes should NOT be bitwise-ORed but only
+// accessed through set_error_bit() and get_error_bit()
 #define COLVARS_OK 0
 #define COLVARS_ERROR   -1
-#define GENERAL_ERROR   -1  // TODO this can be simply merged with COLVARS_ERROR
-#define COLVARS_NOT_IMPLEMENTED (-1<<1)
-#define INPUT_ERROR     (-1<<2) // out of bounds or inconsistent input
-#define BUG_ERROR       (-1<<3) // Inconsistent state indicating bug
-#define FILE_ERROR      (-1<<4)
-#define MEMORY_ERROR    (-1<<5)
-#define FATAL_ERROR     (-1<<6) // Should be set, or not, together with other bits
-#define DELETE_COLVARS  (-1<<7) // Instruct the caller to delete cvm
-#define COLVARS_NO_SUCH_FRAME (-1<<8) // Cannot load the requested frame
-
+#define COLVARS_NOT_IMPLEMENTED (-1*(1<<1))
+#define INPUT_ERROR     (-1*(1<<2)) // out of bounds or inconsistent input
+#define BUG_ERROR       (-1*(1<<3)) // Inconsistent state indicating bug
+#define FILE_ERROR      (-1*(1<<4))
+#define MEMORY_ERROR    (-1*(1<<5))
+#define FATAL_ERROR     (-1*(1<<6)) // Should be set, or not, together with other bits
+#define DELETE_COLVARS  (-1*(1<<7)) // Instruct the caller to delete cvm
+#define COLVARS_NO_SUCH_FRAME (-1*(1<<8)) // Cannot load the requested frame
 
 #include <iostream>
 #include <iomanip>
@@ -115,7 +115,11 @@ protected:
 
 public:
 
-  static void set_error_bits(int code);
+  static void set_error_bit(int code);
+
+  static bool get_error_bit(int code);
+
+  static void combine_errors(int &target, int const code);
 
   static inline int get_error()
   {
@@ -395,7 +399,7 @@ public:
   static void fatal_error(std::string const &message);
 
   /// Print a message to the main log and set global error code
-  static void error(std::string const &message, int code = GENERAL_ERROR);
+  static void error(std::string const &message, int code = COLVARS_ERROR);
 
   /// Print a message to the main log and exit normally
   static void exit(std::string const &message);
diff --git a/lib/colvars/colvarproxy.h b/lib/colvars/colvarproxy.h
index 585a40cdfc..b80b575729 100644
--- a/lib/colvars/colvarproxy.h
+++ b/lib/colvars/colvarproxy.h
@@ -69,10 +69,11 @@ public:
   virtual cvm::real rand_gaussian(void) = 0;
 
   /// \brief Get the current frame number
+  // Negative return values indicate error
   virtual int frame() { return COLVARS_NOT_IMPLEMENTED; }
 
   /// \brief Set the current frame number
-  // return 0 on success, -1 on failure
+  // return frame number on success, COLVARS_NO_SUCH_FRAME otherwise
   virtual int frame(int) { return COLVARS_NOT_IMPLEMENTED; }
 
   /// \brief Prefix to be used for input files (restarts, not
diff --git a/lib/colvars/colvarscript.cpp b/lib/colvars/colvarscript.cpp
index bd1c8878fb..cabb520248 100644
--- a/lib/colvars/colvarscript.cpp
+++ b/lib/colvars/colvarscript.cpp
@@ -53,17 +53,16 @@ int colvarscript::run(int argc, char const *argv[]) {
   }
 
   if (cmd == "delete") {
-    colvars->reset();
     // Note: the delete bit may be ignored by some backends
     // it is mostly useful in VMD
-    colvars->set_error_bits(DELETE_COLVARS);
+    colvars->set_error_bit(DELETE_COLVARS);
     return COLVARS_OK;
   }
 
   if (cmd == "update") {
-    error_code |= proxy->update_input();
-    error_code |= colvars->calc();
-    error_code |= proxy->update_output();
+    cvm::combine_errors(error_code, proxy->update_input());
+    cvm::combine_errors(error_code, colvars->calc());
+    cvm::combine_errors(error_code, proxy->update_output());
     if (error_code) {
       result += "Error updating the colvars module.\n";
     }
@@ -143,8 +142,8 @@ int colvarscript::run(int argc, char const *argv[]) {
     }
     proxy->output_prefix_str = argv[2];
     int error = 0;
-    error |= colvars->setup_output();
-    error |= colvars->write_output_files();
+    cvm::combine_errors(error, colvars->setup_output());
+    cvm::combine_errors(error, colvars->write_output_files());
     return error ? COLVARSCRIPT_ERROR : COLVARS_OK;
   }
 
@@ -406,7 +405,7 @@ Input and output:\n\
   list biases                 -- return a list of all biases\n\
   load <file name>            -- load a state file (requires configuration)\n\
   save <file name>            -- save a state file (requires configuration)\n\
-  update                      -- recalculate colvars and biases based\n\
+  update                      -- recalculate colvars and biases\n\
   printframe                  -- return a summary of the current frame\n\
   printframelabels            -- return labels to annotate printframe's output\n";
 
diff --git a/lib/colvars/colvarscript.h b/lib/colvars/colvarscript.h
index 25138387ca..ccd3ff82c5 100644
--- a/lib/colvars/colvarscript.h
+++ b/lib/colvars/colvarscript.h
@@ -9,7 +9,7 @@
 #include "colvarbias.h"
 #include "colvarproxy.h"
 
-// TODO merge these into colvarmodule.h
+// Only these error values are part of the scripting interface
 #define COLVARSCRIPT_ERROR -1
 #define COLVARSCRIPT_OK 0
 
diff --git a/lib/colvars/colvartypes.cpp b/lib/colvars/colvartypes.cpp
index ad297a6ac9..6f773bd5f7 100644
--- a/lib/colvars/colvartypes.cpp
+++ b/lib/colvars/colvartypes.cpp
@@ -312,10 +312,8 @@ void colvarmodule::rotation::calc_optimal_rotation(std::vector<cvm::atom_pos> co
   S_backup.resize(4,4);
   S_backup = S;
 
-  if (cvm::debug()) {
-    if (b_debug_gradients) {
-      cvm::log("S     = "+cvm::to_str(cvm::to_str(S_backup), cvm::cv_width, cvm::cv_prec)+"\n");
-    }
+  if (b_debug_gradients) {
+    cvm::log("S     = "+cvm::to_str(cvm::to_str(S_backup), cvm::cv_width, cvm::cv_prec)+"\n");
   }
 
   diagonalize_matrix(S, S_eigval, S_eigvec);
@@ -344,25 +342,23 @@ void colvarmodule::rotation::calc_optimal_rotation(std::vector<cvm::atom_pos> co
     q_old = q;
   }
 
-  if (cvm::debug()) {
-    if (b_debug_gradients) {
-      cvm::log("L0 = "+cvm::to_str(L0, cvm::cv_width, cvm::cv_prec)+
-               ", Q0 = "+cvm::to_str(Q0, cvm::cv_width, cvm::cv_prec)+
-               ", Q0*Q0 = "+cvm::to_str(Q0.inner(Q0), cvm::cv_width, cvm::cv_prec)+
-               "\n");
-      cvm::log("L1 = "+cvm::to_str(L1, cvm::cv_width, cvm::cv_prec)+
-               ", Q1 = "+cvm::to_str(Q1, cvm::cv_width, cvm::cv_prec)+
-               ", Q0*Q1 = "+cvm::to_str(Q0.inner(Q1), cvm::cv_width, cvm::cv_prec)+
-               "\n");
-      cvm::log("L2 = "+cvm::to_str(L2, cvm::cv_width, cvm::cv_prec)+
-               ", Q2 = "+cvm::to_str(Q2, cvm::cv_width, cvm::cv_prec)+
-               ", Q0*Q2 = "+cvm::to_str(Q0.inner(Q2), cvm::cv_width, cvm::cv_prec)+
-               "\n");
-      cvm::log("L3 = "+cvm::to_str(L3, cvm::cv_width, cvm::cv_prec)+
-               ", Q3 = "+cvm::to_str(Q3, cvm::cv_width, cvm::cv_prec)+
-               ", Q0*Q3 = "+cvm::to_str(Q0.inner(Q3), cvm::cv_width, cvm::cv_prec)+
-               "\n");
-    }
+  if (b_debug_gradients) {
+    cvm::log("L0 = "+cvm::to_str(L0, cvm::cv_width, cvm::cv_prec)+
+             ", Q0 = "+cvm::to_str(Q0, cvm::cv_width, cvm::cv_prec)+
+             ", Q0*Q0 = "+cvm::to_str(Q0.inner(Q0), cvm::cv_width, cvm::cv_prec)+
+             "\n");
+    cvm::log("L1 = "+cvm::to_str(L1, cvm::cv_width, cvm::cv_prec)+
+             ", Q1 = "+cvm::to_str(Q1, cvm::cv_width, cvm::cv_prec)+
+             ", Q0*Q1 = "+cvm::to_str(Q0.inner(Q1), cvm::cv_width, cvm::cv_prec)+
+             "\n");
+    cvm::log("L2 = "+cvm::to_str(L2, cvm::cv_width, cvm::cv_prec)+
+             ", Q2 = "+cvm::to_str(Q2, cvm::cv_width, cvm::cv_prec)+
+             ", Q0*Q2 = "+cvm::to_str(Q0.inner(Q2), cvm::cv_width, cvm::cv_prec)+
+             "\n");
+    cvm::log("L3 = "+cvm::to_str(L3, cvm::cv_width, cvm::cv_prec)+
+             ", Q3 = "+cvm::to_str(Q3, cvm::cv_width, cvm::cv_prec)+
+             ", Q0*Q3 = "+cvm::to_str(Q0.inner(Q3), cvm::cv_width, cvm::cv_prec)+
+             "\n");
   }
 
   // calculate derivatives of L0 and Q0 with respect to each atom in
@@ -472,46 +468,43 @@ void colvarmodule::rotation::calc_optimal_rotation(std::vector<cvm::atom_pos> co
       }
     }
 
-    if (cvm::debug()) {
+    if (b_debug_gradients) {
 
-      if (b_debug_gradients) {
+      cvm::matrix2d<cvm::real> S_new(4, 4);
+      cvm::vector1d<cvm::real> S_new_eigval(4);
+      cvm::matrix2d<cvm::real> S_new_eigvec(4, 4);
 
-        cvm::matrix2d<cvm::real> S_new(4, 4);
-        cvm::vector1d<cvm::real> S_new_eigval(4);
-        cvm::matrix2d<cvm::real> S_new_eigvec(4, 4);
+      // make an infitesimal move along each cartesian coordinate of
+      // this atom, and solve again the eigenvector problem
+      for (size_t comp = 0; comp < 3; comp++) {
 
-        // make an infitesimal move along each cartesian coordinate of
-        // this atom, and solve again the eigenvector problem
-        for (size_t comp = 0; comp < 3; comp++) {
-
-          S_new = S_backup;
-          // diagonalize the new overlap matrix
-          for (size_t i = 0; i < 4; i++) {
-            for (size_t j = 0; j < 4; j++) {
-              S_new[i][j] +=
-                colvarmodule::debug_gradients_step_size * ds_2[i][j][comp];
-            }
+        S_new = S_backup;
+        // diagonalize the new overlap matrix
+        for (size_t i = 0; i < 4; i++) {
+          for (size_t j = 0; j < 4; j++) {
+            S_new[i][j] +=
+              colvarmodule::debug_gradients_step_size * ds_2[i][j][comp];
           }
-
-          //           cvm::log("S_new = "+cvm::to_str(cvm::to_str (S_new), cvm::cv_width, cvm::cv_prec)+"\n");
-
-          diagonalize_matrix(S_new, S_new_eigval, S_new_eigvec);
-
-          cvm::real const &L0_new = S_new_eigval[0];
-          cvm::quaternion const Q0_new(S_new_eigvec[0]);
-
-          cvm::real const DL0 = (dl0_2[comp]) * colvarmodule::debug_gradients_step_size;
-          cvm::quaternion const DQ0(dq0_2[0][comp] * colvarmodule::debug_gradients_step_size,
-                                    dq0_2[1][comp] * colvarmodule::debug_gradients_step_size,
-                                    dq0_2[2][comp] * colvarmodule::debug_gradients_step_size,
-                                    dq0_2[3][comp] * colvarmodule::debug_gradients_step_size);
-
-          cvm::log(  "|(l_0+dl_0) - l_0^new|/l_0 = "+
-                     cvm::to_str(std::fabs(L0+DL0 - L0_new)/L0, cvm::cv_width, cvm::cv_prec)+
-                     ", |(q_0+dq_0) - q_0^new| = "+
-                     cvm::to_str((Q0+DQ0 - Q0_new).norm(), cvm::cv_width, cvm::cv_prec)+
-                     "\n");
         }
+
+        //           cvm::log("S_new = "+cvm::to_str(cvm::to_str (S_new), cvm::cv_width, cvm::cv_prec)+"\n");
+
+        diagonalize_matrix(S_new, S_new_eigval, S_new_eigvec);
+
+        cvm::real const &L0_new = S_new_eigval[0];
+        cvm::quaternion const Q0_new(S_new_eigvec[0]);
+
+        cvm::real const DL0 = (dl0_2[comp]) * colvarmodule::debug_gradients_step_size;
+        cvm::quaternion const DQ0(dq0_2[0][comp] * colvarmodule::debug_gradients_step_size,
+                                  dq0_2[1][comp] * colvarmodule::debug_gradients_step_size,
+                                  dq0_2[2][comp] * colvarmodule::debug_gradients_step_size,
+                                  dq0_2[3][comp] * colvarmodule::debug_gradients_step_size);
+
+        cvm::log(  "|(l_0+dl_0) - l_0^new|/l_0 = "+
+                   cvm::to_str(std::fabs(L0+DL0 - L0_new)/L0, cvm::cv_width, cvm::cv_prec)+
+                   ", |(q_0+dq_0) - q_0^new| = "+
+                   cvm::to_str((Q0+DQ0 - Q0_new).norm(), cvm::cv_width, cvm::cv_prec)+
+                   "\n");
       }
     }
   }