diff --git a/doc/PDF/colvars-refman-lammps.pdf b/doc/PDF/colvars-refman-lammps.pdf
index d325bca62d..83b51844f2 100644
Binary files a/doc/PDF/colvars-refman-lammps.pdf and b/doc/PDF/colvars-refman-lammps.pdf differ
diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp
index 1ee35864bb..ebf03b814a 100644
--- a/lib/colvars/colvar.cpp
+++ b/lib/colvars/colvar.cpp
@@ -881,6 +881,7 @@ int colvar::calc_cvc_gradients(int first_cvc, size_t num_cvcs)
           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]);
         }
       }
diff --git a/lib/colvars/colvaratoms.cpp b/lib/colvars/colvaratoms.cpp
index 2ab20ca4ec..ad7d8fcaf0 100644
--- a/lib/colvars/colvaratoms.cpp
+++ b/lib/colvars/colvaratoms.cpp
@@ -681,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,9 +778,9 @@ 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 = cog;
+  cog_orig = this->center_of_geometry();
   if (ref_pos_group) {
-    ref_pos_group->cog_orig = ref_pos_group->cog;
+    ref_pos_group->cog_orig = ref_pos_group->center_of_geometry();
   }
 
   if (b_center) {
@@ -950,8 +950,6 @@ 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");
 
@@ -962,12 +960,10 @@ void cvm::atom_group::calc_fit_gradients()
     cvm::rvector atom_grad;
 
     for (size_t i = 0; i < this->size(); i++) {
-      atom_grad += (-1.0)/(cvm::real(group_for_fit->size())) * atoms[i].grad;
+      atom_grad += atoms[i].grad;
     }
-
-    // need to use the gradients in original frame
-    // we only rotate the sum for efficiency
     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;
@@ -981,17 +977,16 @@ void cvm::atom_group::calc_fit_gradients()
 
     for (size_t i = 0; i < this->size(); i++) {
 
-      // restore original position for this atom
+      // compute centered, unrotated position
       cvm::atom_pos const pos_orig =
-        rot_inv.rotate((b_center ? (atoms[i].pos - ref_pos_cog) : (atoms[i].pos))) +
-        (ref_pos_group ? ref_pos_group->cog_orig : cog_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++) {
-        // 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];
         }
@@ -1148,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/colvarcomp.cpp b/lib/colvars/colvarcomp.cpp
index 9e4d75ecaa..ef64989ad1 100644
--- a/lib/colvars/colvarcomp.cpp
+++ b/lib/colvars/colvarcomp.cpp
@@ -39,14 +39,17 @@ colvar::cvc::cvc(std::string const &conf)
   get_keyval(conf, "period", period, 0.0);
   get_keyval(conf, "wrapAround", wrap_center, 0.0);
 
-  get_keyval(conf, "debugGradients", b_debug_gradient, false, parse_silent);
+  // All cvcs implement this
+  provide(f_cvc_debug_gradient);
+  {
+    bool b_debug_gradient;
+    get_keyval(conf, "debugGradients", b_debug_gradient, false, parse_silent);
+    if (b_debug_gradient) enable(f_cvc_debug_gradient);
+  }
 
   // Attempt scalable calculations when in parallel? (By default yes, if available)
   get_keyval(conf, "scalable", b_try_scalable, true);
 
-  // All cvcs implement this
-  provide(f_cvc_debug_gradient);
-
   if (cvm::debug())
     cvm::log("Done initializing cvc base object.\n");
 }
@@ -107,7 +110,6 @@ int colvar::cvc::setup()
   if (b_try_scalable && is_available(f_cvc_scalable)) {
     enable(f_cvc_scalable);
   }
-  if (b_debug_gradient) enable(f_cvc_debug_gradient);
 
   return COLVARS_OK;
 }
@@ -176,30 +178,21 @@ void colvar::cvc::debug_gradients(cvm::atom_group *group)
   // cvm::log("gradients     = "+cvm::to_str (gradients)+"\n");
 
   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 original frame: we should print them in the rotated frame
-      if (group->b_rotate) {
-        for (j = 0; j < group_for_fit->fit_gradients.size(); j++) {
-          group_for_fit->fit_gradients[j] = rot_0.rotate(group_for_fit->fit_gradients[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) + "] = " + cvm::to_str(group_for_fit->fit_gradients[j]));
-      }
-
-      if (group->b_rotate) {
-        for (j = 0; j < group_for_fit->fit_gradients.size(); j++) {
-          group_for_fit->fit_gradients[j] = rot_inv.rotate(group_for_fit->fit_gradients[j]);
-        }
+                 "[" + 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])));
       }
     }
   }
@@ -208,9 +201,10 @@ void colvar::cvc::debug_gradients(cvm::atom_group *group)
   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
@@ -242,26 +236,26 @@ void colvar::cvc::debug_gradients(cvm::atom_group *group)
 
     for (size_t ia = 0; ia < ref_group->size(); ia++) {
 
-      // tests are best conducted in the unrotated (simulation) frame
+      // 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)| = "+
@@ -273,6 +267,10 @@ void colvar::cvc::debug_gradients(cvm::atom_group *group)
     }
   }
 
+  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.h b/lib/colvars/colvarcomp.h
index 027cd3499f..923b2c6bd1 100644
--- a/lib/colvars/colvarcomp.h
+++ b/lib/colvars/colvarcomp.h
@@ -224,9 +224,6 @@ public:
   /// \brief Whether or not this CVC will be computed in parallel whenever possible
   bool b_try_scalable;
 
-private:
-  bool b_debug_gradient;
-
 protected:
 
   /// \brief Cached value
diff --git a/lib/colvars/colvarcomp_distances.cpp b/lib/colvars/colvarcomp_distances.cpp
index a2a5102b8b..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"))
@@ -273,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)
@@ -373,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)
@@ -783,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;
   }
@@ -863,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());
   }
 }
@@ -1059,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());
diff --git a/lib/colvars/colvarcomp_rotations.cpp b/lib/colvars/colvarcomp_rotations.cpp
index 3bfab768ab..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));
diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp
index 0eef090793..b84f226ec5 100644
--- a/lib/colvars/colvarmodule.cpp
+++ b/lib/colvars/colvarmodule.cpp
@@ -1069,7 +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"
+     << "  version " << std::string(COLVARS_VERSION) << "\n"
      << "}\n\n";
 
   cvm::increase_depth();
@@ -1420,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 9bc5b4c886..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-20"
+#define COLVARS_VERSION "2016-04-27"
 #endif
 
 #ifndef COLVARS_DEBUG
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");
       }
     }
   }