diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp
index 74e52f507b..40f3c236ca 100644
--- a/lib/colvars/colvar.cpp
+++ b/lib/colvars/colvar.cpp
@@ -292,6 +292,11 @@ colvar::colvar (std::string const &conf)
     cvm::fatal_error ("Error: trying to expand boundaries that already "
                       "cover a whole period of a periodic colvar.\n");
   }
+  if (expand_boundaries && hard_lower_boundary && hard_upper_boundary) {
+    cvm::fatal_error ("Error: inconsistent configuration "
+                      "(trying to expand boundaries with both "
+                      "hardLowerBoundary and hardUpperBoundary enabled).\n");
+  }
 
   {
     bool b_extended_lagrangian;
@@ -676,12 +681,39 @@ void colvar::calc()
   if (cvm::debug())
     cvm::log ("Calculating colvar \""+this->name+"\".\n");
 
+  // prepare atom groups for calculation
+  for (size_t i = 0; i < cvcs.size(); i++) {
+    for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
+      cvm::atom_group &atoms = *(cvcs[i]->atom_groups[ig]);
+      atoms.reset_atoms_data();
+      atoms.read_positions();
+      if (atoms.b_center || atoms.b_rotate) {
+        atoms.calc_apply_roto_translation();
+      }
+      // each atom group will take care of its own ref_pos_group, if defined
+    }
+  }
+  if (tasks[task_output_velocity]) {
+    for (size_t i = 0; i < cvcs.size(); i++) {
+      for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
+        cvcs[i]->atom_groups[ig]->read_velocities();
+      }   
+    }
+  }
+  if (tasks[task_system_force]) {
+    for (size_t i = 0; i < cvcs.size(); i++) {
+      for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
+        cvcs[i]->atom_groups[ig]->read_system_forces();
+      }   
+    }
+  }
+
   // calculate the value of the colvar
 
   x.reset();
   if (x.type() == colvarvalue::type_scalar) {
+    // polynomial combination allowed
 
-    // scalar variable, polynomial combination allowed
     for (size_t i = 0; i < cvcs.size(); i++) {
       cvm::increase_depth();
       (cvcs[i])->calc_value();
@@ -697,6 +729,7 @@ void colvar::calc()
           (cvcs[i])->value().real_value );
     } 
   } else {
+    // only linear combination allowed
 
     for (size_t i = 0; i < cvcs.size(); i++) {
       cvm::increase_depth();
@@ -716,12 +749,29 @@ void colvar::calc()
               cvm::to_str (x, cvm::cv_width, cvm::cv_prec)+".\n");
 
   if (tasks[task_gradients]) {
-    // calculate the gradients
+
+    if (cvm::debug())
+      cvm::log ("Calculating gradients of colvar \""+this->name+"\".\n");
+
     for (size_t i = 0; i < cvcs.size(); i++) {
+      // calculate the gradients of each component
       cvm::increase_depth();
+
       (cvcs[i])->calc_gradients();
+
+      // if requested, propagate (via chain rule) the gradients above
+      // to the atoms used to define the roto-translation
+      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();
+      }   
+
       cvm::decrease_depth();
     }
+
+    if (cvm::debug())
+      cvm::log ("Done calculating gradients of colvar \""+this->name+"\".\n");
+
     if (tasks[task_collect_gradients]) {
       // Collect the atomic gradients inside colvar object
       for (int a = 0; a < atomic_gradients.size(); a++) {
@@ -760,6 +810,10 @@ void colvar::calc()
   }
 
   if (tasks[task_system_force]) {
+
+    if (cvm::debug())
+      cvm::log ("Calculating system force of colvar \""+this->name+"\".\n");
+
     ft.reset();
 
     if(!tasks[task_extended_lagrangian] && (cvm::step_relative() > 0)) {
@@ -778,6 +832,9 @@ void colvar::calc()
       // correction internally: biases such as colvarbias_abf will handle it
       ft += fj;
     }
+
+    if (cvm::debug())
+      cvm::log ("Done calculating system force of colvar \""+this->name+"\".\n");
   }
 
   if (tasks[task_fdiff_velocity]) {
diff --git a/lib/colvars/colvaratoms.cpp b/lib/colvars/colvaratoms.cpp
index 7a3ffa1903..1588fcee8b 100644
--- a/lib/colvars/colvaratoms.cpp
+++ b/lib/colvars/colvaratoms.cpp
@@ -3,29 +3,76 @@
 #include "colvaratoms.h"
 
 
-// atom member functions depend tightly on the MD interface, and are
+// member functions of the "atom" class depend tightly on the MD interface, and are
 // thus defined in colvarproxy_xxx.cpp
 
+// in this file only atom_group functions are defined
 
-// atom_group member functions
 
+// Note: "conf" is the configuration of the cvc who is using this atom group;
+// "key" is the name of the atom group (e.g. "atoms", "group1", "group2", ...)
 cvm::atom_group::atom_group (std::string const &conf,
-                             char const        *key,
-                             atom_group        *ref_pos_group_in)
+                             char const        *key)
   : b_center (false), b_rotate (false), b_user_defined_fit (false),
-    ref_pos_group (NULL), // this is always set within parse(),
-                          // regardless of ref_pos_group_in
+    ref_pos_group (NULL),
+    b_fit_gradients (false),
     noforce (false)
 {
   cvm::log ("Defining atom group \""+
             std::string (key)+"\".\n");
-  parse (conf, key, ref_pos_group_in);
+  // real work is done by parse
+  parse (conf, key);
+}
+
+
+cvm::atom_group::atom_group (std::vector<cvm::atom> const &atoms)
+  : b_dummy (false), b_center (false), b_rotate (false), 
+    ref_pos_group (NULL), b_fit_gradients (false), 
+    noforce (false)
+{
+  this->reserve (atoms.size());
+  for (size_t i = 0; i < atoms.size(); i++) {
+    this->push_back (atoms[i]);
+  }
+  total_mass = 0.0;
+  for (cvm::atom_iter ai = this->begin();
+       ai != this->end(); ai++) {
+    total_mass += ai->mass;
+  }
+}
+
+
+cvm::atom_group::atom_group()
+  : b_dummy (false), b_center (false), b_rotate (false), 
+    ref_pos_group (NULL), b_fit_gradients (false), 
+    noforce (false)
+{
+  total_mass = 0.0;
+}
+
+
+cvm::atom_group::~atom_group()
+{
+  if (ref_pos_group) {
+    delete ref_pos_group;
+    ref_pos_group = NULL;
+  }
+}
+
+
+void cvm::atom_group::add_atom (cvm::atom const &a)
+{
+  if (b_dummy) {
+    cvm::fatal_error ("Error: cannot add atoms to a dummy group.\n");
+  } else {
+    this->push_back (a);
+    total_mass += a.mass;
+  }
 }
 
 
 void cvm::atom_group::parse (std::string const &conf,
-                             char const        *key,
-                             atom_group        *ref_pos_group_in)
+                             char const        *key)
 {
   std::string group_conf;
 
@@ -226,6 +273,8 @@ void cvm::atom_group::parse (std::string const &conf,
   bool b_defined_center = get_keyval (group_conf, "centerReference", b_center, false, mode);
   bool b_defined_rotate = get_keyval (group_conf, "rotateReference", b_rotate, false, mode);
 
+  b_fit_gradients = get_keyval (group_conf, "fitGradients", b_fit_gradients, false, mode);
+
   // this cannot be shortened to one statement because lazy evaluation may prevent one
   // function from being called!
   b_user_defined_fit = b_defined_center || b_defined_rotate;
@@ -236,9 +285,9 @@ void cvm::atom_group::parse (std::string const &conf,
 
   if (b_center || b_rotate) {
 
-    // instead of this group, define another group (refPositionsGroup) to be the one
-    // used to fit the coordinates
     if (key_lookup (group_conf, "refPositionsGroup")) {
+      // instead of this group, define another group (refPositionsGroup)
+      // to be the one used to fit the coordinates
       if (ref_pos_group) {
         cvm::fatal_error ("Error: the atom group \""+
                           std::string (key)+"\" has already a reference group "
@@ -252,15 +301,7 @@ void cvm::atom_group::parse (std::string const &conf,
 
     atom_group *group_for_fit = ref_pos_group ? ref_pos_group : this;
 
-    if (get_keyval (group_conf, "refPositions", ref_pos, ref_pos, mode)) {
-      if (ref_pos.size() != group_for_fit->size()) {
-        cvm::fatal_error ("Error: the number of reference positions provided ("+
-                          cvm::to_str (ref_pos.size())+
-                          ") does not match the number of atoms of group \""+
-                          std::string (key)+
-                          "\" ("+cvm::to_str (group_for_fit->size())+").\n");
-      }
-    }
+    get_keyval (group_conf, "refPositions", ref_pos, ref_pos, mode);
 
     std::string ref_pos_file;
     if (get_keyval (group_conf, "refPositionsFile", ref_pos_file, std::string (""), mode)) {
@@ -310,16 +351,22 @@ void cvm::atom_group::parse (std::string const &conf,
 #endif
     }
 
+    if (b_fit_gradients) {
+      group_for_fit->fit_gradients.assign (group_for_fit->size(), cvm::atom_pos (0.0, 0.0, 0.0));
+      rot.request_group1_gradients (group_for_fit->size());
+    }
+
     if (b_rotate && !noforce) {
       cvm::log ("Warning: atom group \""+std::string (key)+
-                "\" will be fitted automatically onto a fixed orientation: "
-                "in few cases, torques applied on this group "
-                "may make the simulation unstable.  "
-                "If this happens, set \"disableForces\" to yes "
-                "for this group.\n");
+                "\" will be aligned to a fixed orientation given by the reference positions provided.  "
+                "If the internal structure of the group changes too much (i.e. its RMSD is comparable "
+                "to its radius of gyration), the optimal rotation and its gradients may become discontinuous.  "
+                "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());
     }
+
   }
 
   if (cvm::debug())
@@ -336,50 +383,6 @@ void cvm::atom_group::parse (std::string const &conf,
 }
 
 
-cvm::atom_group::atom_group (std::vector<cvm::atom> const &atoms)
-  : b_dummy (false), b_center (false), b_rotate (false), 
-    ref_pos_group (NULL), noforce (false)
-{
-  this->reserve (atoms.size());
-  for (size_t i = 0; i < atoms.size(); i++) {
-    this->push_back (atoms[i]);
-  }
-  total_mass = 0.0;
-  for (cvm::atom_iter ai = this->begin();
-       ai != this->end(); ai++) {
-    total_mass += ai->mass;
-  }
-}
-
-
-cvm::atom_group::atom_group()
-  : b_dummy (false), b_center (false), b_rotate (false), 
-    ref_pos_group (NULL), noforce (false)
-{
-  total_mass = 0.0;
-}
-
-
-cvm::atom_group::~atom_group()
-{
-  if (ref_pos_group) {
-    delete ref_pos_group;
-    ref_pos_group = NULL;
-  }
-}
-
-
-void cvm::atom_group::add_atom (cvm::atom const &a)
-{
-  if (b_dummy) {
-    cvm::fatal_error ("Error: cannot add atoms to a dummy group.\n");
-  } else {
-    this->push_back (a);
-    total_mass += a.mass;
-  }
-}
-
-
 void cvm::atom_group::create_sorted_ids (void)
 {
   // Only do the work if the vector is not yet populated
@@ -425,11 +428,6 @@ void cvm::atom_group::read_positions()
 {
   if (b_dummy) return;
 
-#if (! defined (COLVARS_STANDALONE))
-  if (!this->size())
-    cvm::fatal_error ("Error: no atoms defined in the requested group.\n");
-#endif
-
   for (cvm::atom_iter ai = this->begin();
        ai != this->end(); ai++) {
     ai->read_position();
@@ -437,13 +435,14 @@ void cvm::atom_group::read_positions()
 
   if (ref_pos_group)
     ref_pos_group->read_positions();
+}
 
+void cvm::atom_group::calc_apply_roto_translation()
+{
   atom_group *fit_group = ref_pos_group ? ref_pos_group : this;
 
   if (b_center) {
-    // store aside the current center of geometry (all positions will be
-    // set to the closest images to the first one) and then center on
-    // the origin
+    // center on the origin first
     cvm::atom_pos const cog = fit_group->center_of_geometry();
     for (cvm::atom_iter ai = this->begin();
          ai != this->end(); ai++) {
@@ -453,9 +452,7 @@ void cvm::atom_group::read_positions()
 
   if (b_rotate) {
     // rotate the group (around the center of geometry if b_center is
-    // true, around the origin otherwise); store the rotation, in
-    // order to bring back the forces to the original frame before
-    // applying them
+    // true, around the origin otherwise)
     rot.calc_optimal_rotation (fit_group->positions(), ref_pos);
 
     for (cvm::atom_iter ai = this->begin();
@@ -465,7 +462,7 @@ void cvm::atom_group::read_positions()
   }
 
   if (b_center) {
-    // use the center of geometry of ref_pos
+    // align with the center of geometry of ref_pos
     for (cvm::atom_iter ai = this->begin();
          ai != this->end(); ai++) {
       ai->pos += ref_pos_cog;
@@ -473,7 +470,6 @@ void cvm::atom_group::read_positions()
   }
 }
 
-
 void cvm::atom_group::apply_translation (cvm::rvector const &t) 
 {
   if (b_dummy) return;
@@ -484,7 +480,6 @@ void cvm::atom_group::apply_translation (cvm::rvector const &t)
   }
 }
 
-
 void cvm::atom_group::apply_rotation (cvm::rotation const &rot) 
 {
   if (b_dummy) return;
@@ -538,24 +533,6 @@ void cvm::atom_group::read_system_forces()
   }
 }
 
-/* This is deprecated.
-cvm::atom_pos cvm::atom_group::center_of_geometry (cvm::atom_pos const &ref_pos)
-{
-  if (b_dummy) {
-    cvm::select_closest_image (dummy_atom_pos, ref_pos);
-    return dummy_atom_pos;
-  }
-
-  cvm::atom_pos cog (0.0, 0.0, 0.0);
-  for (cvm::atom_iter ai = this->begin();
-       ai != this->end(); ai++) {
-    cvm::select_closest_image (ai->pos, ref_pos);
-    cog += ai->pos;
-  }
-  cog /= this->size();
-  return cog;
-} */
-
 cvm::atom_pos cvm::atom_group::center_of_geometry() const
 {
   if (b_dummy)
@@ -570,24 +547,6 @@ cvm::atom_pos cvm::atom_group::center_of_geometry() const
   return cog;
 }
 
-/* This is deprecated.
-cvm::atom_pos cvm::atom_group::center_of_mass (cvm::atom_pos const &ref_pos)
-{
-  if (b_dummy) {
-    cvm::select_closest_image (dummy_atom_pos, ref_pos);
-    return dummy_atom_pos;
-  }
-
-  cvm::atom_pos com (0.0, 0.0, 0.0);
-  for (cvm::atom_iter ai = this->begin();
-       ai != this->end(); ai++) {
-    cvm::select_closest_image (ai->pos, ref_pos);
-    com += ai->mass * ai->pos;
-  }
-  com /= this->total_mass;
-  return com;
-}*/
-
 cvm::atom_pos cvm::atom_group::center_of_mass() const
 {
   if (b_dummy)
@@ -614,6 +573,57 @@ void cvm::atom_group::set_weighted_gradient (cvm::rvector const &grad)
 }
 
 
+void cvm::atom_group::calc_fit_gradients()
+{
+  if (b_dummy) return;
+
+  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
+    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 ((*this)[i].grad) :
+        (*this)[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;
+      }
+    }
+  }
+
+  if (b_rotate) {
+    
+    // add the rotation matrix contribution to the gradients
+    cvm::rotation const rot_inv = rot.inverse();
+    cvm::atom_pos const cog = this->center_of_geometry();
+    
+    for (size_t i = 0; i < this->size(); i++) {
+
+      cvm::atom_pos const pos_orig = rot_inv.rotate ((b_center ? ((*this)[i].pos - cog) : ((*this)[i].pos)));
+
+      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, (*this)[i].grad);
+        // multiply by \cdot {\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");
+}
+
+
 std::vector<cvm::atom_pos> cvm::atom_group::positions() const
 {
   if (b_dummy)
@@ -699,7 +709,7 @@ void cvm::atom_group::apply_colvar_force (cvm::real const &force)
 
   if (b_rotate) {
 
-    // get the forces back from the rotated frame
+    // rotate forces back to the original frame
     cvm::rotation const rot_inv = rot.inverse();
     for (cvm::atom_iter ai = this->begin();
          ai != this->end(); ai++) {
@@ -708,13 +718,30 @@ void cvm::atom_group::apply_colvar_force (cvm::real const &force)
 
   } else {
 
-    // no need to manipulate gradients, they are still in the original
-    // frame
     for (cvm::atom_iter ai = this->begin();
          ai != this->end(); ai++) {
       ai->apply_force (force * ai->grad);
     }
   }
+
+  if (b_fit_gradients) {
+
+    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 * fit_gradients[j]));
+      }
+    } else {
+      for (size_t j = 0; j < group_for_fit->size(); j++) {
+        (*group_for_fit)[j].apply_force (force * fit_gradients[j]);
+      }
+    }
+  }
+
 }
 
 
@@ -778,4 +805,3 @@ void cvm::atom_group::apply_forces (std::vector<cvm::rvector> const &forces)
   }
 }
 
-
diff --git a/lib/colvars/colvaratoms.h b/lib/colvars/colvaratoms.h
index 6ffda6db56..d2be1dfe6e 100644
--- a/lib/colvars/colvaratoms.h
+++ b/lib/colvars/colvaratoms.h
@@ -140,39 +140,42 @@ public:
   void create_sorted_ids (void);
   
 
-  /// \brief Before calculating colvars, move the group to overlap the
-  /// center of mass of reference coordinates
+  /// \brief When updating atomic coordinates, translate them to align with the
+  /// center of mass of the reference coordinates
   bool b_center;
 
-  /// \brief Right after updating atom coordinates (and after
-  /// centering coordinates, if b_center is true), rotate the group to
-  /// overlap the reference coordinates.  You should not manipulate
-  /// atoms individually if you turn on this flag.
+  /// \brief When updating atom coordinates (and after
+  /// centering them if b_center is set), rotate the group to
+  /// align with the reference coordinates.
   ///
   /// Note: gradients will be calculated in the rotated frame: when
   /// forces will be applied, they will rotated back to the original
   /// frame
   bool b_rotate;
-  /// Rotation between the group and its reference coordinates
+  /// The rotation calculated automatically if b_rotate is defined
   cvm::rotation rot;
 
-  /// \brief Indicates that the user has specified centerReference or
-  /// rotateReference (and the related reference coordinates):
-  /// cvc's (eg rmsd, eigenvector) should not override these settings
+  /// \brief Indicates that the user has explicitly set centerReference or
+  /// rotateReference, and the corresponding reference:
+  /// cvc's (eg rmsd, eigenvector) will not override the user's choice
   bool b_user_defined_fit;
 
-  /// \brief use these reference coordinates, for b_center, b_rotate, or to compute
-  /// certain colvar components (orientation, rmsd, etc)
+  /// \brief Whether or not the derivatives of the roto-translation
+  /// should be included when calculating the colvar's gradients (default: no)
+  bool b_fit_gradients;
+
+  /// \brief use reference coordinates for b_center or b_rotate
   std::vector<cvm::atom_pos> ref_pos;
+
   /// \brief Center of geometry of the reference coordinates; regardless
   /// of whether b_center is true, ref_pos is centered to zero at
   /// initialization, and ref_pos_cog serves to center the positions
   cvm::atom_pos              ref_pos_cog;
-  /// \brief In case b_center or b_rotate is true, fit this group to
-  /// the reference positions (default: the parent group itself)
+
+  /// \brief If b_center or b_rotate is true, use this group to
+  /// define the transformation (default: this group itself)
   atom_group                *ref_pos_group;
-
-
+  
   /// Total mass of the atom group
   cvm::real total_mass;
 
@@ -186,14 +189,12 @@ public:
   /// which is a member function so that a group can be initialized
   /// also after construction
   atom_group (std::string const &conf,
-              char const        *key,
-              atom_group        *ref_pos_group = NULL);
+              char const        *key);
 
   /// \brief Initialize the group by looking up its configuration
   /// string in conf and parsing it
   void parse (std::string const &conf,
-              char const        *key,
-              atom_group        *ref_pos_group = NULL);
+              char const        *key);
 
   /// \brief Initialize the group after a temporary vector of atoms
   atom_group (std::vector<cvm::atom> const &atoms);
@@ -208,10 +209,12 @@ public:
   ~atom_group();
 
   /// \brief Get the current positions; if b_center or b_rotate are
-  /// true, center and/or rotate the coordinates right after reading
-  /// them
+  /// true, calc_apply_roto_translation() will be called too
   void read_positions();
 
+  /// \brief (Re)calculate the optimal roto-translation
+  void calc_apply_roto_translation();
+
   /// \brief Save center of geometry fo ref positions,
   /// then subtract it 
   void center_ref_pos();
@@ -238,6 +241,8 @@ public:
   {
     for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++)
       ai->reset_data();
+    if (ref_pos_group)
+      ref_pos_group->reset_atoms_data();
   }
 
   /// \brief Return a copy of the current atom positions
@@ -246,21 +251,15 @@ public:
   /// \brief Return a copy of the current atom positions, shifted by a constant vector
   std::vector<cvm::atom_pos> positions_shifted (cvm::rvector const &shift) const;
 
-  /// \brief Return the center of geometry of the positions \param ref_pos
-  /// Use the closest periodic images to this position
-  cvm::atom_pos center_of_geometry (cvm::atom_pos const &ref_pos);
   /// \brief Return the center of geometry of the positions, assuming
   /// that coordinates are already pbc-wrapped
   cvm::atom_pos center_of_geometry() const;
 
-  /// \brief Return the center of mass of the positions \param ref_pos
-  /// Use the closest periodic images to this position
-  cvm::atom_pos center_of_mass (cvm::atom_pos const &ref_pos);
   /// \brief Return the center of mass of the positions, assuming that
   /// coordinates are already pbc-wrapped
   cvm::atom_pos center_of_mass() const;
 
-  /// \brief Store atom positions from the previous step
+  /// \brief Atom positions at the previous step
   std::vector<cvm::atom_pos> old_pos;
 
 
@@ -279,6 +278,12 @@ public:
   /// \link center_of_mass() \endlink)
   void set_weighted_gradient (cvm::rvector const &grad);
 
+  /// \brief Calculate the derivatives of the fitting transformation
+  void calc_fit_gradients();
+
+  /// \brief Derivatives of the fitting transformation
+  std::vector<cvm::atom_pos> fit_gradients;
+
   /// \brief Used by a (scalar) colvar to apply its force on its \link
   /// atom_group \endlink members
   ///
@@ -290,7 +295,7 @@ public:
   /// colvar components, and use apply_force() or apply_forces()).  If
   /// the group is being rotated to a reference frame (e.g. to express
   /// the colvar independently from the solute rotation), the
-  /// gradients are temporarily to the original frame.
+  /// gradients are temporarily rotated to the original frame.
   void apply_colvar_force (cvm::real const &force);
 
   /// \brief Apply a force "to the center of mass", i.e. the force is
diff --git a/lib/colvars/colvarbias_meta.cpp b/lib/colvars/colvarbias_meta.cpp
index a3a483d606..2d03b5c8d8 100644
--- a/lib/colvars/colvarbias_meta.cpp
+++ b/lib/colvars/colvarbias_meta.cpp
@@ -49,6 +49,15 @@ colvarbias_meta::colvarbias_meta (std::string const &conf, char const *key)
 
   get_keyval (conf, "hillWidth", hill_width, std::sqrt (2.0 * PI) / 2.0);
 
+  {
+    bool b_replicas = false;
+    get_keyval (conf, "multipleReplicas", b_replicas, false);
+    if (b_replicas)
+      comm = multiple_replicas;
+    else 
+      comm = single_replica;
+  }
+
   get_keyval (conf, "useGrids", use_grids, true);
 
   if (use_grids) {
@@ -84,15 +93,6 @@ colvarbias_meta::colvarbias_meta (std::string const &conf, char const *key)
     dump_replica_fes = false;
   }
 
-  {
-    bool b_replicas = false;
-    get_keyval (conf, "multipleReplicas", b_replicas, false);
-    if (b_replicas)
-      comm = multiple_replicas;
-    else 
-      comm = single_replica;
-  }
-
   if (comm != single_replica) {
 
     if (expand_grids)
@@ -373,34 +373,36 @@ cvm::real colvarbias_meta::update()
         int       &new_size = new_sizes[i];
         bool changed_lb = false, changed_ub = false;
 
-        if (curr_bin[i] < min_buffer) {
-          int const extra_points = (min_buffer - curr_bin[i]);
-          new_lb -= extra_points * colvars[i]->width;
-          new_size += extra_points;
-          // changed offset in this direction => the pointer needs to
-          // be changed, too
-          curr_bin[i] += extra_points;
+        if (!colvars[i]->hard_lower_boundary)
+          if (curr_bin[i] < min_buffer) {
+            int const extra_points = (min_buffer - curr_bin[i]);
+            new_lb -= extra_points * colvars[i]->width;
+            new_size += extra_points;
+            // changed offset in this direction => the pointer needs to
+            // be changed, too
+            curr_bin[i] += extra_points;
 
-          changed_lb = true;
-          cvm::log ("Metadynamics bias \""+this->name+"\""+
-                    ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
-                    ": new lower boundary for colvar \""+
-                    colvars[i]->name+"\", at "+
-                    cvm::to_str (new_lower_boundaries[i])+".\n");
-        }
+            changed_lb = true;
+            cvm::log ("Metadynamics bias \""+this->name+"\""+
+                      ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
+                      ": new lower boundary for colvar \""+
+                      colvars[i]->name+"\", at "+
+                      cvm::to_str (new_lower_boundaries[i])+".\n");
+          }
 
-        if (curr_bin[i] > new_size - min_buffer - 1) {
-          int const extra_points = (curr_bin[i] - (new_size - 1) + min_buffer);
-          new_ub += extra_points * colvars[i]->width;
-          new_size += extra_points;
+        if (!colvars[i]->hard_upper_boundary)
+          if (curr_bin[i] > new_size - min_buffer - 1) {
+            int const extra_points = (curr_bin[i] - (new_size - 1) + min_buffer);
+            new_ub += extra_points * colvars[i]->width;
+            new_size += extra_points;
 
-          changed_ub = true;
-          cvm::log ("Metadynamics bias \""+this->name+"\""+
-                    ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
-                    ": new upper boundary for colvar \""+
-                    colvars[i]->name+"\", at "+
-                    cvm::to_str (new_upper_boundaries[i])+".\n");
-        }
+            changed_ub = true;
+            cvm::log ("Metadynamics bias \""+this->name+"\""+
+                      ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
+                      ": new upper boundary for colvar \""+
+                      colvars[i]->name+"\", at "+
+                      cvm::to_str (new_upper_boundaries[i])+".\n");
+          }
 
         if (changed_lb || changed_ub)
           changed_grids = true;
@@ -1521,7 +1523,8 @@ std::ostream & colvarbias_meta::write_restart (std::ostream& os)
 void colvarbias_meta::write_pmf()
 {
   // allocate a new grid to store the pmf
-  colvar_grid_scalar *pmf = new colvar_grid_scalar (colvars);
+  colvar_grid_scalar *pmf = new colvar_grid_scalar (*hills_energy);
+  pmf->create();
 
   std::string fes_file_name_prefix (cvm::output_prefix);
   
@@ -1544,17 +1547,18 @@ void colvarbias_meta::write_pmf()
       cvm::real const well_temper_scale = (bias_temperature + cvm::temperature()) / bias_temperature;
       pmf->multiply_constant (well_temper_scale);
     }
-    std::string const fes_file_name (fes_file_name_prefix +
-                                     ((comm != single_replica) ? ".partial" : "") +
-                                     (dump_fes_save ?
-                                      "."+cvm::to_str (cvm::step_absolute()) : "") +
-                                     ".pmf");
-    cvm::backup_file (fes_file_name.c_str());
-    std::ofstream fes_os (fes_file_name.c_str());
-    pmf->write_multicol (fes_os);
-    fes_os.close();
+    {
+      std::string const fes_file_name (fes_file_name_prefix +
+                                       ((comm != single_replica) ? ".partial" : "") +
+                                       (dump_fes_save ?
+                                        "."+cvm::to_str (cvm::step_absolute()) : "") +
+                                       ".pmf");
+      cvm::backup_file (fes_file_name.c_str());
+      std::ofstream fes_os (fes_file_name.c_str());
+      pmf->write_multicol (fes_os);
+      fes_os.close();
+    }
   }
-
   if (comm != single_replica) {
     // output the combined PMF from all replicas
     pmf->reset();
diff --git a/lib/colvars/colvarcomp.cpp b/lib/colvars/colvarcomp.cpp
index f1997a6dd9..b3e89527c3 100644
--- a/lib/colvars/colvarcomp.cpp
+++ b/lib/colvars/colvarcomp.cpp
@@ -71,3 +71,76 @@ void colvar::cvc::calc_Jacobian_derivative()
   cvm::fatal_error ("Error: calculation of inverse gradients is not implemented "
                     "for colvar components of type \""+function_type+"\".\n");
 }
+
+
+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")
+
+  // collect into a vector for convenience
+  std::vector<cvm::rvector> gradients (group.size());
+  for (size_t i = 0; i < group.size(); i++) {
+    gradients[i] = group[i].grad;
+  }
+
+  cvm::rotation const rot_0 = group.rot;
+  cvm::rotation const rot_inv = group.rot.inverse();
+
+  cvm::real const x_0 = x.real_value;
+
+  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_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]);
+      }
+    }
+  }
+
+  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;
+
+    for (size_t id = 0; id < 3; id++) {
+      // (re)read original positions
+      group.read_positions();
+      // change one coordinate 
+      group[ia].pos[id] += cvm::debug_gradients_step_size;
+      // (re)do the fit (if defined)
+      if (group.b_center || group.b_rotate) {
+        group.calc_apply_roto_translation();
+      }
+      calc_value();
+      cvm::real const x_1 = x.real_value;
+      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.b_fit_gradients && (group.ref_pos_group == NULL)) ?
+        (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");
+      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");
+
+    }
+  }
+}
diff --git a/lib/colvars/colvarcomp.h b/lib/colvars/colvarcomp.h
index 33bd73f645..a8a06814f4 100644
--- a/lib/colvars/colvarcomp.h
+++ b/lib/colvars/colvarcomp.h
@@ -1,5 +1,5 @@
-#ifndef COLVARDEF_H
-#define COLVARDEF_H
+#ifndef COLVARCOMP_H
+#define COLVARCOMP_H
 
 #include <fstream>
 #include <cmath>
@@ -106,10 +106,6 @@ public:
   /// Destructor
   virtual ~cvc();
 
-  /// \brief If true, calc_gradients() will calculate
-  /// finite-difference gradients alongside the analytical ones and
-  /// report their differences
-  bool b_debug_gradients;
 
   /// \brief If this flag is false (default), inverse gradients
   /// (derivatives of atom coordinates with respect to x) are
@@ -130,6 +126,12 @@ public:
   /// order to apply forces
   virtual void calc_gradients() = 0;
 
+  /// \brief If true, calc_gradients() will call debug_gradients() for every group needed
+  bool b_debug_gradients;
+
+  /// \brief Calculate finite-difference gradients alongside the analytical ones, for each Cartesian component
+  virtual void debug_gradients (cvm::atom_group &group);
+
   /// \brief Calculate the total force from the system using the
   /// inverse atomic gradients
   virtual void calc_force_invgrads();
@@ -527,11 +529,8 @@ public:
 /// \brief Colvar component: moment of inertia of an atom group
 /// (colvarvalue::type_scalar type, range [0:*))
 class colvar::inertia
-  : public colvar::cvc
+  : public colvar::gyration
 {
-protected:
-  /// Atoms involved
-  cvm::atom_group atoms;
 public:
   /// Constructor
   inertia (std::string const &conf);
@@ -591,7 +590,10 @@ protected:
   /// Reference coordinates
   std::vector<cvm::atom_pos>  ref_pos;
 
-  /// Eigenvector (of a normal or essential mode)
+  /// Geometric center of the reference coordinates
+  cvm::rvector                ref_pos_center;
+
+  /// Eigenvector (of a normal or essential mode): will always have zero center
   std::vector<cvm::rvector>   eigenvec;
 
   /// Inverse square norm of the eigenvector
@@ -1138,7 +1140,6 @@ public:
 };
 
 
-
 // metrics functions for cvc implementations with a periodicity
 
 inline cvm::real colvar::dihedral::dist2 (colvarvalue const &x1,
diff --git a/lib/colvars/colvarcomp_distances.cpp b/lib/colvars/colvarcomp_distances.cpp
index bc77c5a4b3..9519620b48 100644
--- a/lib/colvars/colvarcomp_distances.cpp
+++ b/lib/colvars/colvarcomp_distances.cpp
@@ -8,9 +8,6 @@
 
 
 
-/// \file cvc_distance.cpp \brief Collective variables
-/// determining various type of distances between two groups
-
 // "twogroup" flag defaults to true; set to false by selfCoordNum
 // (only distance-derived component based on only one group)
 
@@ -48,12 +45,6 @@ colvar::distance::distance()
 
 void colvar::distance::calc_value()
 {
-  group1.reset_atoms_data();
-  group2.reset_atoms_data();
-
-  group1.read_positions();
-  group2.read_positions();
-
   if (b_no_PBC) {
     dist_v = group2.center_of_mass() - group1.center_of_mass();
   } else {
@@ -113,12 +104,6 @@ colvar::distance_vec::distance_vec()
 
 void colvar::distance_vec::calc_value()
 {
-  group1.reset_atoms_data();
-  group2.reset_atoms_data();
-
-  group1.read_positions();
-  group2.read_positions();
-
   if (b_no_PBC) {
     x.rvector_value = group2.center_of_mass() - group1.center_of_mass();
   } else {
@@ -205,12 +190,6 @@ colvar::distance_z::distance_z()
 
 void colvar::distance_z::calc_value()
 {
-  main.reset_atoms_data();
-  ref1.reset_atoms_data();
-
-  main.read_positions();
-  ref1.read_positions();
-
   if (fixed_axis) {
     if (b_no_PBC) {
       dist_v = main.center_of_mass() - ref1.center_of_mass();
@@ -219,8 +198,6 @@ void colvar::distance_z::calc_value()
                                        main.center_of_mass());
     }
   } else {
-    ref2.reset_atoms_data();
-    ref2.read_positions();
 
     if (b_no_PBC) {
       dist_v = main.center_of_mass() -
@@ -310,12 +287,6 @@ colvar::distance_xy::distance_xy()
 
 void colvar::distance_xy::calc_value()
 {
-  ref1.reset_atoms_data();
-  main.reset_atoms_data();
-
-  ref1.read_positions();
-  main.read_positions();
-
   if (b_no_PBC) {
     dist_v = main.center_of_mass() - ref1.center_of_mass();
   } else {
@@ -323,9 +294,6 @@ void colvar::distance_xy::calc_value()
                                      main.center_of_mass());
   }
   if (!fixed_axis) {
-    ref2.reset_atoms_data();
-    ref2.read_positions();
-
     if (b_no_PBC) {
       v12 = ref2.center_of_mass() - ref1.center_of_mass();
     } else {
@@ -415,12 +383,6 @@ colvar::distance_dir::distance_dir()
 
 void colvar::distance_dir::calc_value()
 {
-  group1.reset_atoms_data();
-  group2.reset_atoms_data();
-
-  group1.read_positions();
-  group2.read_positions();
-
   if (b_no_PBC) {
     dist_v = group2.center_of_mass() - group1.center_of_mass();
   } else {
@@ -490,12 +452,6 @@ colvar::distance_inv::distance_inv()
 
 void colvar::distance_inv::calc_value()
 {
-  group1.reset_atoms_data();
-  group2.reset_atoms_data();
-
-  group1.read_positions();
-  group2.read_positions();
-
   x.real_value = 0.0;
   if (b_no_PBC) {
     for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) {
@@ -561,6 +517,17 @@ colvar::gyration::gyration (std::string const &conf)
   b_Jacobian_derivative = true;
   parse_group (conf, "atoms", atoms);
   atom_groups.push_back (&atoms);
+
+  if (atoms.b_user_defined_fit) {
+    cvm::log ("WARNING: explicit fitting parameters were provided for atom group \"atoms\".");
+  } else {
+    // default: fit everything
+    cvm::log ("No explicit fitting parameters: enabling centerReference by default to align with the origin\n");
+    cvm::log (" (rotateReference will not affect results for gyration and inertia, but it will for inertiaZ).\n");
+    atoms.b_center = true;
+    atoms.ref_pos.assign (1, cvm::rvector (0.0, 0.0, 0.0));
+  }
+
   x.type (colvarvalue::type_scalar);
 }
 
@@ -576,10 +543,6 @@ colvar::gyration::gyration()
 
 void colvar::gyration::calc_value()
 {
-  atoms.reset_atoms_data();
-  atoms.read_positions();
-  atoms.apply_translation ((-1.0) * atoms.center_of_geometry());
-
   x.real_value = 0.0;
   for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
     x.real_value += (ai->pos).norm2();
@@ -594,6 +557,11 @@ void colvar::gyration::calc_gradients()
   for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
     ai->grad = drdx * ai->pos;
   }
+
+  if (b_debug_gradients) {
+    cvm::log ("Debugging gradients:\n");
+    debug_gradients (atoms);
+  }
 }
 
 
@@ -625,11 +593,11 @@ void colvar::gyration::apply_force (colvarvalue const &force)
 
 
 colvar::inertia::inertia (std::string const &conf)
-  : cvc (conf)
+  : gyration (conf)
 {
   function_type = "inertia";
-  parse_group (conf, "atoms", atoms);
-  atom_groups.push_back (&atoms);
+  b_inverse_gradients = false;
+  b_Jacobian_derivative = false;
   x.type (colvarvalue::type_scalar);
 }
 
@@ -643,10 +611,6 @@ colvar::inertia::inertia()
 
 void colvar::inertia::calc_value()
 {
-  atoms.reset_atoms_data();
-  atoms.read_positions();
-  atoms.apply_translation ((-1.0) * atoms.center_of_mass());
-
   x.real_value = 0.0;
   for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
     x.real_value += ai->mass * (ai->pos).norm2();
@@ -659,6 +623,11 @@ void colvar::inertia::calc_gradients()
   for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
     ai->grad = 2.0 * ai->mass * ai->pos;
   }
+
+  if (b_debug_gradients) {
+    cvm::log ("Debugging gradients:\n");
+    debug_gradients (atoms);
+  }
 }
 
 
@@ -694,10 +663,6 @@ colvar::inertia_z::inertia_z()
 
 void colvar::inertia_z::calc_value()
 {
-  atoms.reset_atoms_data();
-  atoms.read_positions();
-  atoms.apply_translation ((-1.0) * atoms.center_of_mass());
-
   x.real_value = 0.0;
   for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
     cvm::real const iprod = ai->pos * axis;
@@ -711,6 +676,11 @@ void colvar::inertia_z::calc_gradients()
   for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
     ai->grad = 2.0 * ai->mass * (ai->pos * axis) * axis;
   }
+
+  if (b_debug_gradients) {
+    cvm::log ("Debugging gradients:\n");
+    debug_gradients (atoms);
+  }
 }
 
 
@@ -736,32 +706,15 @@ colvar::rmsd::rmsd (std::string const &conf)
   if (atoms.b_dummy) 
     cvm::fatal_error ("Error: \"atoms\" cannot be a dummy atom.");
 
-  if (atoms.b_user_defined_fit) {
-    cvm::log ("WARNING: explicit fitting parameters were provided for atom group \"atoms\".");
-    cvm::log ("Not computing standard minimum RMSD (if that is the desired result, leave atom group parameters at default values).");
-  } else {
-    // Default: fit everything
-    // NOTE: this won't work for class V cvc's
-    atoms.b_center = true;
-    atoms.b_rotate = true;
+  if (atoms.ref_pos_group != NULL) {
+    cvm::log ("The option \"refPositionsGroup\" (alternative group for fitting) was enabled: "
+              "Jacobian derivatives of the RMSD will not be calculated.\n");
+    b_Jacobian_derivative = false;
   }
 
-  if (atoms.b_center || atoms.b_rotate) {
-    // 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 2nd group for Jacobian: this is only
-    // required for ABF, but we do both groups here for better caching
-    atoms.rot.request_group2_gradients (atoms.size());
-
-    if (atoms.ref_pos_group != NULL) {
-      cvm::log ("The option \"refPositionsGroup\" (alternative group for fitting) was enabled: "
-                "Jacobian derivatives of the RMSD will not be calculated.\n");
-      b_Jacobian_derivative = false;
-    }
-  }
-
-  // the following is a simplified version of the corresponding atom group options: 
-  // to define the reference coordinates to compute this variable
+  // the following is a simplified version of the corresponding atom group options;
+  // we need this because the reference coordinates defined inside the atom group
+  // may be used only for fitting, and even more so if ref_pos_group is used
   if (get_keyval (conf, "refPositions", ref_pos, ref_pos)) {
     cvm::log ("Using reference positions from configuration file to calculate the variable.\n");
     if (ref_pos.size() != atoms.size()) {
@@ -800,41 +753,37 @@ colvar::rmsd::rmsd (std::string const &conf)
     }
   }
 
-  if (!atoms.b_user_defined_fit) {
-    // Default case: reference positions for calculating the rmsd are also those used
+  if (atoms.b_user_defined_fit) {
+    cvm::log ("WARNING: explicit fitting parameters were provided for atom group \"atoms\".");
+  } else {
+    // Default: fit everything
+    cvm::log ("No explicit fitting parameters: enabling both centerReference and rotateReference to compute standard minimum RMSD.");
+    atoms.b_center = true;
+    atoms.b_rotate = true;
+    // default case: reference positions for calculating the rmsd are also those used
     // for fitting
     atoms.ref_pos = ref_pos;
     atoms.center_ref_pos();
+
+    // 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());
   }
 }
 
   
 void colvar::rmsd::calc_value()
 {
-  atoms.reset_atoms_data();
-  atoms.read_positions();
-  // rotational fit is done internally
-
-  // atoms_cog = atoms.center_of_geometry();
-  // rot.calc_optimal_rotation (ref_pos, atoms.positions_shifted (-1.0 * atoms_cog));
-
-  // cvm::real group_pos_sum2 = 0.0;
-  // for (size_t i = 0; i < atoms.size(); i++) {
-  //   group_pos_sum2 += (atoms[i].pos - atoms_cog).norm2();
-  // }
-
-  // // value of the RMSD (Coutsias et al)
-  // cvm::real const MSD = 1.0/(cvm::real (atoms.size())) *
-  //   ( group_pos_sum2 + ref_pos_sum2 - 2.0 * rot.lambda );
-
-  // x.real_value = (MSD > 0.0) ? std::sqrt (MSD) : 0.0;
+  // rotational-translational fit is handled by the atom group
 
   x.real_value = 0.0;
   for (size_t ia = 0; ia < atoms.size(); ia++) {
     x.real_value += (atoms[ia].pos - ref_pos[ia]).norm2();
   }
   x.real_value /= cvm::real (atoms.size()); // MSD
-  x.real_value = (x.real_value > 0.0) ? std::sqrt (x.real_value) : 0.0;
+  x.real_value = std::sqrt (x.real_value);
 }
 
 
@@ -847,6 +796,11 @@ 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 (b_debug_gradients) {
+    cvm::log ("Debugging gradients:\n");
+    debug_gradients (atoms);
+  }
 }
 
 
@@ -908,11 +862,14 @@ void colvar::rmsd::calc_Jacobian_derivative()
       grad_rot_mat[1][2] =  2.0 * (g23 - g01); 
       grad_rot_mat[2][2] = -2.0 * (g11 + g22); 
 
-      cvm::atom_pos &x = atoms[ia].pos; 
+      cvm::atom_pos &y = ref_pos[ia];
 
-      for (size_t i = 0; i < 3; i++) {
-        for (size_t j = 0; j < 3; j++) {
-          divergence += grad_rot_mat[i][j][i] * x[j];
+      for (size_t alpha = 0; alpha < 3; alpha++) {
+        for (size_t beta = 0; beta < 3; beta++) {
+          divergence += grad_rot_mat[beta][alpha][alpha] * y[beta];
+        // Note: equation was derived for inverse rotation (see colvars paper)
+        // so here the matrix is transposed
+        // (eq would give   divergence += grad_rot_mat[alpha][beta][alpha] * y[beta];)
         }
       }
     }
@@ -935,23 +892,24 @@ colvar::eigenvector::eigenvector (std::string const &conf)
   parse_group (conf, "atoms", atoms);
   atom_groups.push_back (&atoms);
 
-  if (atoms.b_rotate) {
-    cvm::fatal_error ("Error: rotateReference should be disabled:"
-        "eigenvector component will set it internally.");
-  }
-
-  if (get_keyval (conf, "refPositions", ref_pos, ref_pos)) {
-    cvm::log ("Using reference positions from input file.\n");
-    if (ref_pos.size() != atoms.size()) {
-      cvm::fatal_error ("Error: reference positions do not "
-                        "match the number of requested atoms.\n");
-    }
-  }
 
   {
+    bool const b_inline = get_keyval (conf, "refPositions", ref_pos, ref_pos);
+
+    if (b_inline) {
+      cvm::log ("Using reference positions from input file.\n");
+      if (ref_pos.size() != atoms.size()) {
+        cvm::fatal_error ("Error: reference positions do not "
+                          "match the number of requested atoms.\n");
+      }
+    }
+
     std::string file_name;
     if (get_keyval (conf, "refPositionsFile", file_name)) {
 
+      if (b_inline)
+        cvm::fatal_error ("Error: refPositions and refPositionsFile cannot be specified at the same time.\n");
+
       std::string file_col;
       double file_col_value;
       if (get_keyval (conf, "refPositionsCol", file_col, std::string (""))) {
@@ -970,35 +928,59 @@ colvar::eigenvector::eigenvector (std::string const &conf)
     }
   }
 
-  // Set mobile frame of reference for atom group
-  atoms.b_center = true;
-  atoms.b_rotate = true;
-  atoms.ref_pos = ref_pos;
-  atoms.center_ref_pos();
+  // save for later the geometric center of the provided positions (may not be the origin)
+  cvm::rvector ref_pos_center (0.0, 0.0, 0.0);
+  for (size_t i = 0; i < atoms.size(); i++) {
+    ref_pos_center += ref_pos[i];
+  }
+  ref_pos_center *= 1.0 / atoms.size();
 
-  // now load the eigenvector
-  if (get_keyval (conf, "vector", eigenvec, eigenvec)) {
-    cvm::log ("Using vector components from input file.\n");
-    if (eigenvec.size() != atoms.size()) {
-      cvm::fatal_error ("Error: vector components do not "
-                        "match the number of requested atoms.\n");
-    }
+  if (atoms.b_user_defined_fit) {
+    cvm::log ("WARNING: explicit fitting parameters were provided for atom group \"atoms\".\n");
+  } else {
+    // default: fit everything
+    cvm::log ("Enabling centerReference and rotateReference (minimize RMSD before calculating the projection).\n");
+    atoms.b_center = true;
+    atoms.b_rotate = true;
+    atoms.ref_pos = ref_pos;
+    atoms.center_ref_pos();
+
+    // 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());
   }
 
   {
+    bool const b_inline = get_keyval (conf, "vector", eigenvec, eigenvec);
+    // now load the eigenvector
+    if (b_inline) {
+      cvm::log ("Using vector components from input file.\n");
+      if (eigenvec.size() != atoms.size()) {
+        cvm::fatal_error ("Error: vector components do not "
+                          "match the number of requested atoms.\n");
+      }
+    }
+
     std::string file_name;
     if (get_keyval (conf, "vectorFile", file_name)) {
 
-      std::string file_col;
-      if (!get_keyval (conf, "vectorCol", file_col, std::string (""))) {
-        cvm::fatal_error ("Error: parameter vectorCol is required if vectorFile is set.\n");
-      }
+      if (b_inline)
+        cvm::fatal_error ("Error: vector and vectorFile cannot be specified at the same time.\n");
 
+      std::string file_col;
       double file_col_value;
-      bool found = get_keyval (conf, "vectorColValue", file_col_value, 0.0);
-      if (found && !file_col_value)
-        cvm::fatal_error ("Error: eigenvectorColValue, "
-                          "if provided, must be non-zero.\n");
+      if (get_keyval (conf, "vectorCol", file_col, std::string (""))) {
+        // use PDB flags if column is provided
+        bool found = get_keyval (conf, "vectorColValue", file_col_value, 0.0);
+        if (found && !file_col_value)
+          cvm::fatal_error ("Error: vectorColValue, "
+                            "if provided, must be non-zero.\n");
+      } else {
+        // if not, use atom indices
+        atoms.create_sorted_ids();
+      }
 
       eigenvec.resize (atoms.size());
       cvm::load_coords (file_name.c_str(), eigenvec, atoms.sorted_ids, file_col, file_col_value);
@@ -1010,34 +992,64 @@ colvar::eigenvector::eigenvector (std::string const &conf)
                       "and eigenvector must be defined.\n");
   }
 
-  cvm::rvector center (0.0, 0.0, 0.0);
-  eigenvec_invnorm2 = 0.0;
-
+  cvm::rvector eig_center (0.0, 0.0, 0.0);
   for (size_t i = 0; i < atoms.size(); i++) {
-    center += eigenvec[i];
+    eig_center += eigenvec[i];
   }
-  center *= 1.0 / atoms.size();
+  eig_center *= 1.0 / atoms.size();
+  cvm::log ("Geometric center of the provided vector: "+cvm::to_str (eig_center)+"\n");
 
-  cvm::log ("Subtracting center of eigenvector components: " + cvm::to_str (center) + "\n");
+  bool b_vector_difference;
+  get_keyval (conf, "vectorDifference", b_vector_difference, false);
+  if (b_vector_difference) {
 
+    cvm::log ("Subtracting the reference positions from the provided vector.\n");
+
+    if (atoms.b_center) {
+      cvm::log ("centerReference is on: recentering the provided vector to the reference positions.\n");
+      // both sets should be centered on the origin for fitting
+      for (size_t i = 0; i < atoms.size(); i++) {
+        eigenvec[i] -= eig_center;
+        ref_pos[i]  -= ref_pos_center;
+      }
+    }
+
+    if (atoms.b_rotate) {
+      cvm::log ("rotateReference is on: aligning the provided vector to the reference positions.\n");
+      atoms.rot.calc_optimal_rotation (eigenvec, ref_pos);
+      for (size_t i = 0; i < atoms.size(); i++) {
+        eigenvec[i] = atoms.rot.rotate (eigenvec[i]);
+        eigenvec[i] -= ref_pos[i];
+      }
+    }
+
+    if (atoms.b_center) {
+      // bring back the ref positions to where they were
+      for (size_t i = 0; i < atoms.size(); i++) {
+        ref_pos[i] -= ref_pos_center;
+      }
+    }
+
+  } else {
+    cvm::log ("Centering the provided vector to zero.\n");
+    for (size_t i = 0; i < atoms.size(); i++) {
+      eigenvec[i] -= eig_center;
+    }
+  }
+
+  cvm::log ("The first three components (v1x, v1y, v1z) of the resulting vector are: "+cvm::to_str (eigenvec[0])+".\n");
+
+  // for inverse gradients
+  eigenvec_invnorm2 = 0.0;
   for (size_t i = 0; i < atoms.size(); i++) {
-    eigenvec[i] = eigenvec[i] - center;
     eigenvec_invnorm2 += eigenvec[i].norm2();
   }
   eigenvec_invnorm2 = 1.0 / eigenvec_invnorm2;
-
-  // request derivatives of optimal rotation wrt 2nd group
-  // for Jacobian
-  atoms.rot.request_group1_gradients(atoms.size());
-  atoms.rot.request_group2_gradients(atoms.size());
 }
 
   
 void colvar::eigenvector::calc_value()
 {
-  atoms.reset_atoms_data();
-  atoms.read_positions(); // this will also update atoms.rot
-
   x.real_value = 0.0;
   for (size_t i = 0; i < atoms.size(); i++) {
     x.real_value += (atoms[i].pos - ref_pos[i]) * eigenvec[i];
@@ -1060,7 +1072,14 @@ void colvar::eigenvector::calc_gradients()
     atoms[ia].grad = eigenvec[ia];
   }
 
-/*
+  std::vector<cvm::rvector> gradients (atoms.size());
+  for (size_t i = 0; i < atoms.size(); i++) {
+    gradients[i] = atoms[i].grad;
+  }
+  std::vector<cvm::rvector> gradients_jh_correct (atoms.size());
+
+  /*
+
   // Version B: complex, expensive, fragile. Right.
  
   // gradient of the rotation matrix
@@ -1080,7 +1099,6 @@ void colvar::eigenvector::calc_gradients()
   for (size_t ia = 0; ia < atoms.size(); ia++) {
 
     // Gradient of optimal quaternion wrt current Cartesian position
-    // WARNING: we want derivatives wrt the FIRST group here (unlike RMSD)
     cvm::vector1d< cvm::rvector, 4 >      &dq = atoms.rot.dQ0_1[ia];
 
     g11 = 2.0 * quat0[1]*dq[1];
@@ -1118,9 +1136,76 @@ void colvar::eigenvector::calc_gradients()
     }
 
     // Rotate correction term back to reference frame
-    atoms[ia].grad = eigenvec[ia] + quat0.rotate (rot_grad_term);
+    //    atoms[ia].grad = eigenvec[ia] + quat0.rotate (rot_grad_term);
+    //    atoms[ia].grad = eigenvec[ia] + quat0.rotate (rot_grad_term);
+
+    gradients_jh_correct[ia] = eigenvec[ia] + quat0.rotate (rot_grad_term);
+  }
+  */
+
+  // TODO replace this with a call to debug_gradients()
+  if (b_debug_gradients) {
+
+    // the following code should work for any scalar variable;
+    // the only difference is the name of the atom group object (here, "atoms")
+
+    cvm::rotation const rot_0 = atoms.rot;
+    cvm::rotation const rot_inv = atoms.rot.inverse();
+
+    cvm::real const x_0 = x.real_value;
+
+    cvm::log ("gradients     = "+cvm::to_str (gradients)+"\n");
+    /*    cvm::log ("gradients(JH-correct) = "+cvm::to_str (gradients_jh_correct)+"\n"); */
+
+    if (atoms.b_fit_gradients) {
+      atoms.calc_fit_gradients();
+      if (atoms.b_rotate) {
+        // fit_gradients are in the original frame, output them in the rotated frame
+        for (size_t j = 0; j < atoms.fit_gradients.size(); j++) {
+          atoms.fit_gradients[j] = rot_0.rotate (atoms.fit_gradients[j]);
+        }
+      }
+      cvm::log ("fit_gradients = "+cvm::to_str (atoms.fit_gradients)+"\n");
+      if (atoms.b_rotate) {
+        for (size_t j = 0; j < atoms.fit_gradients.size(); j++) {
+          atoms.fit_gradients[j] = rot_inv.rotate (atoms.fit_gradients[j]);
+        }
+      }
+    }
+
+    for (size_t ia = 0; ia < atoms.size(); ia++) {
+
+      // tests are best conducted in the unrotated (simulation) frame
+      cvm::rvector const atom_grad = atoms.b_rotate ?
+        rot_inv.rotate (atoms[ia].grad) :
+        atoms[ia].grad;
+
+      for (size_t id = 0; id < 3; id++) {
+        // (re)read original positions
+        atoms.read_positions();
+        // change one coordinate 
+        atoms[ia].pos[id] += cvm::debug_gradients_step_size;
+        // (re)do the fit (if defined)
+        if (atoms.b_center || atoms.b_rotate) {
+          atoms.calc_apply_roto_translation();
+        }
+        calc_value();
+        cvm::real const x_1 = x.real_value;
+        cvm::log ("Atom "+cvm::to_str (ia)+", component "+cvm::to_str (id)+":\n");
+        cvm::log ("dx(real) = "+cvm::to_str (x_1 - x_0)+"\n");
+        cvm::real const dx_pred = atoms.b_fit_gradients ?
+          (cvm::debug_gradients_step_size * (atom_grad[id] + atoms.fit_gradients[ia][id])) :
+          (cvm::debug_gradients_step_size * atom_grad[id]);
+        cvm::log ("dx(pred) = "+cvm::to_str (dx_pred)+"\n");
+        cvm::log ("|dx(real) - dx(pred)|/dx(real) = "+
+                  cvm::to_str (std::fabs (x_1 - x_0 - dx_pred) /
+                               std::fabs (x_1 - x_0),
+                               cvm::cv_width, cvm::cv_prec)+
+                  ".\n");
+
+      }
+    }
   }
-*/
 }
 
 
@@ -1196,3 +1281,4 @@ void colvar::eigenvector::calc_Jacobian_derivative()
 
 
 
+
diff --git a/lib/colvars/colvargrid.h b/lib/colvars/colvargrid.h
index c6e4d780a6..80a7b70b16 100644
--- a/lib/colvars/colvargrid.h
+++ b/lib/colvars/colvargrid.h
@@ -175,6 +175,7 @@ public:
   /// create() must be called after that;
   colvar_grid (colvar_grid<T> const &g) : has_data (false),
                                           nd (g.nd),
+                                          nx (g.nx),
                                           mult (g.mult),
                                           cv (g.cv),
                                           lower_boundaries (g.lower_boundaries),
diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp
index 00744740cf..67ef6cadd9 100644
--- a/lib/colvars/colvarmodule.cpp
+++ b/lib/colvars/colvarmodule.cpp
@@ -46,11 +46,11 @@ colvarmodule::colvarmodule (char const  *config_filename,
 
   parse->get_keyval (conf, "analysis", b_analysis, false);
 
-  parse->get_keyval (conf, "debugGradientsStepSize", debug_gradients_step_size, 1.0e-03,
+  parse->get_keyval (conf, "debugGradientsStepSize", debug_gradients_step_size, 1.0e-07,
                      colvarparse::parse_silent);
 
   parse->get_keyval (conf, "eigenvalueCrossingThreshold",
-                     colvarmodule::rotation::crossing_threshold, 1.0e-04,
+                     colvarmodule::rotation::crossing_threshold, 1.0e-02,
                      colvarparse::parse_silent);
 
   parse->get_keyval (conf, "colvarsTrajFrequency", cv_traj_freq, 100);
diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h
index 2375fafa39..e228a775d1 100644
--- a/lib/colvars/colvarmodule.h
+++ b/lib/colvars/colvarmodule.h
@@ -2,7 +2,7 @@
 #define COLVARMODULE_H
 
 #ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2012-11-20"
+#define COLVARS_VERSION "2013-01-23"
 #endif
 
 #ifndef COLVARS_DEBUG