diff --git a/lib/colvars/colvarcomp.h b/lib/colvars/colvarcomp.h
index 078c9b4cd..964a73396 100644
--- a/lib/colvars/colvarcomp.h
+++ b/lib/colvars/colvarcomp.h
@@ -1,1445 +1,1445 @@
 #ifndef COLVARCOMP_H
 #define COLVARCOMP_H
 
 // Declaration of colvar::cvc base class and derived ones.
 //
 // Future cvc's could be declared on additional header files.
 // After the declaration of a new derived class, its metric
 // functions must be reimplemented as well.
 // If the new cvc has no symmetry or periodicity,
 // this can be done straightforwardly by using the macro:
 // simple_scalar_dist_functions (derived_class)
 
 
 #include <fstream>
 #include <cmath>
 
 
 #include "colvarmodule.h"
 #include "colvar.h"
 #include "colvaratoms.h"
 
 
 /// \brief Colvar component (base class); most implementations of
 /// \link cvc \endlink utilize one or more \link
 /// colvarmodule::atom \endlink or \link colvarmodule::atom_group
 /// \endlink objects to access atoms.
 ///
 /// A \link cvc \endlink object (or an object of a
 /// cvc-derived class) specifies how to calculate a collective
 /// variable, its gradients and other related physical quantities
 /// which do not depend only on the numeric value (the \link colvar
 /// \endlink class already serves this purpose).
 ///
 /// No restriction is set to what kind of calculation a \link
 /// cvc \endlink object performs (usually calculate an
 /// analytical function of atomic coordinates).  The only constraint
 /// is that the value calculated is implemented as a \link colvarvalue
 /// \endlink object.  This serves to provide a unique way to calculate
 /// scalar and non-scalar collective variables, and specify if and how
 /// they can be combined together by the parent \link colvar \endlink
 /// object.
 ///
 /// <b> If you wish to implement a new collective variable component, you
 /// should write your own class by inheriting directly from \link
 /// cvc \endlink, or one of its derived classes (for instance,
 /// \link distance \endlink is frequently used, because it provides
 /// useful data and function members for any colvar based on two
 /// atom groups).  The steps are: \par
 /// 1. add the name of this class under colvar.h \par
 /// 2. add a call to the parser in colvar.C, within the function colvar::colvar() \par
 /// 3. declare the class in colvarcomp.h \par
 /// 4. implement the class in one of the files colvarcomp_*.C
 ///
 /// </b>
 /// The cvm::atom and cvm::atom_group classes are available to
 /// transparently communicate with the simulation program.  However,
 /// they are not strictly needed, as long as all the degrees of
 /// freedom associated to the cvc are properly evolved by a simple
 /// call to e.g. apply_force().
 
 class colvar::cvc
   : public colvarparse
 {
 public:
 
   /// \brief The name of the object (helps to identify this
   /// cvc instance when debugging)
   std::string name;
 
   /// \brief Description of the type of collective variable
   ///
   /// Normally this string is set by the parent \link colvar \endlink
   /// object within its constructor, when all \link cvc \endlink
   /// objects are initialized; therefore the main "config string"
   /// constructor does not need to define it.  If a \link cvc
   /// \endlink is initialized and/or a different constructor is used,
   /// this variable definition should be set within the constructor.
   std::string function_type;
 
   /// \brief Type of \link colvarvalue \endlink that this cvc
   /// provides
   colvarvalue::Type type() const;
 
   /// \brief Coefficient in the polynomial combination (default: 1.0)
   cvm::real sup_coeff;
   /// \brief Exponent in the polynomial combination (default: 1)
   int       sup_np;
 
   /// \brief Period of this cvc value, (default: 0.0, non periodic)
   cvm::real period;
 
   /// \brief If the component is periodic, wrap around this value (default: 0.0)
   cvm::real wrap_center;
 
   bool b_periodic;
 
   /// \brief Constructor
   ///
   /// At least one constructor which reads a string should be defined
   /// for every class inheriting from cvc \param conf Contents
   /// of the configuration file pertaining to this \link cvc
   /// \endlink
   cvc (std::string const &conf);
 
   /// \brief Within the constructor, make a group parse its own
   /// options from the provided configuration string
   void parse_group (std::string const &conf,
                     char const *group_key,
                     cvm::atom_group &group,
                     bool optional = false);
 
   /// \brief Default constructor (used when \link cvc \endlink
   /// objects are declared within other ones)
   cvc();
 
   /// Destructor
   virtual ~cvc();
 
 
   /// \brief If this flag is false (default), inverse gradients
   /// (derivatives of atom coordinates with respect to x) are
   /// unavailable; it should be set to true by the constructor of each
   /// derived object capable of calculating them
   bool b_inverse_gradients;
 
   /// \brief If this flag is false (default), the Jacobian derivative
   /// (divergence of the inverse gradients) is unavailable; it should
   /// be set to true by the constructor of each derived object capable
   /// of calculating it
   bool b_Jacobian_derivative;
 
   /// \brief Calculate the variable
   virtual void calc_value() = 0;
 
   /// \brief Calculate the atomic gradients, to be reused later in
   /// 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();
 
   /// \brief Calculate the divergence of the inverse atomic gradients
   virtual void calc_Jacobian_derivative();
 
 
   /// \brief Return the previously calculated value
   virtual colvarvalue value() const;
 
   /// \brief Return the previously calculated system force
   virtual colvarvalue system_force() const;
 
   /// \brief Return the previously calculated divergence of the
   /// inverse atomic gradients
   virtual colvarvalue Jacobian_derivative() const;
 
   /// \brief Apply the collective variable force, by communicating the
   /// atomic forces to the simulation program (\b Note: the \link ft
   /// \endlink member is not altered by this function)
   ///
   /// Note: multiple calls to this function within the same simulation
   /// step will add the forces altogether \param cvforce The
   /// collective variable force, usually coming from the biases and
   /// eventually manipulated by the parent \link colvar \endlink
   /// object
   virtual void apply_force (colvarvalue const &cvforce) = 0;
 
   /// \brief Square distance between x1 and x2 (can be redefined to
   /// transparently implement constraints, symmetries and
   /// periodicities)
   ///
   /// colvar::cvc::dist2() and the related functions are
   /// declared as "const" functions, but not "static", because
   /// additional parameters defining the metrics (e.g. the
   /// periodicity) may be specific to each colvar::cvc object.
   ///
   /// If symmetries or periodicities are present, the
   /// colvar::cvc::dist2() should be redefined to return the
   /// "closest distance" value and colvar::cvc::dist2_lgrad(),
   /// colvar::cvc::dist2_rgrad() to return its gradients.
   ///
   /// If constraints are present (and not already implemented by any
   /// of the \link colvarvalue \endlink types), the
   /// colvar::cvc::dist2_lgrad() and
   /// colvar::cvc::dist2_rgrad() functions should be redefined
   /// to provide a gradient which is compatible with the constraint,
   /// i.e. already deprived of its component normal to the constraint
   /// hypersurface.
   ///
   /// Finally, another useful application, if you are performing very
   /// many operations with these functions, could be to override the
   /// \link colvarvalue \endlink member functions and access directly
   /// its member data.  For instance: to define dist2(x1,x2) as
   /// (x2.real_value-x1.real_value)*(x2.real_value-x1.real_value) in
   /// case of a scalar \link colvarvalue \endlink type.
   virtual cvm::real dist2 (colvarvalue const &x1,
                            colvarvalue const &x2) const;
 
   /// \brief Gradient (with respect to x1) of the square distance (can
   /// be redefined to transparently implement constraints, symmetries
   /// and periodicities)
   virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
 
   /// \brief Gradient (with respect to x2) of the square distance (can
   /// be redefined to transparently implement constraints, symmetries
   /// and periodicities)
   virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
 
   /// \brief Return a positive number if x2>x1, zero if x2==x1,
   /// negative otherwise (can be redefined to transparently implement
   /// constraints, symmetries and periodicities) \b Note: \b it \b
   /// only \b works \b with \b scalar \b variables, otherwise raises
   /// an error
   virtual cvm::real compare (colvarvalue const &x1,
                              colvarvalue const &x2) const;
 
   /// \brief Wrapp value (for periodic/symmetric cvcs)
   virtual void wrap (colvarvalue &x) const;
 
   /// \brief Pointers to all atom groups, to let colvars collect info
   /// e.g. atomic gradients
   std::vector<cvm::atom_group *> atom_groups;
 
 protected:
 
   /// \brief Cached value
   colvarvalue x;
 
   /// \brief Value at the previous step
   colvarvalue x_old;
 
   /// \brief Calculated system force (\b Note: this is calculated from
   /// the total atomic forces read from the program, subtracting fromt
   /// the "internal" forces of the system the "external" forces from
   /// the colvar biases)
   colvarvalue ft;
 
   /// \brief Calculated Jacobian derivative (divergence of the inverse
   /// gradients): serves to calculate the phase space correction
   colvarvalue jd;
 };
 
 
 
 
 inline colvarvalue::Type colvar::cvc::type() const
 {
   return x.type();
 }
 
 inline colvarvalue colvar::cvc::value() const
 {
   return x;
 }
 
 inline colvarvalue colvar::cvc::system_force() const
 {
   return ft;
 }
 
 inline colvarvalue colvar::cvc::Jacobian_derivative() const
 {
   return jd;
 }
 
 
 inline cvm::real colvar::cvc::dist2 (colvarvalue const &x1,
                                      colvarvalue const &x2) const
 {
   return x1.dist2 (x2);
 }
 
 inline colvarvalue colvar::cvc::dist2_lgrad (colvarvalue const &x1,
                                              colvarvalue const &x2) const
 {
   return x1.dist2_grad (x2);
 }
 
 inline colvarvalue colvar::cvc::dist2_rgrad (colvarvalue const &x1,
                                              colvarvalue const &x2) const
 {
   return x2.dist2_grad (x1);
 }
 
 inline cvm::real colvar::cvc::compare (colvarvalue const &x1,
                                        colvarvalue const &x2) const
 {
   if (this->type() == colvarvalue::type_scalar) {
     return cvm::real (x1 - x2);
   } else {
     cvm::fatal_error ("Error: you requested an operation which requires "
                       "comparison between two non-scalar values.\n");
     return 0.0;
   }
 }
 
 inline void colvar::cvc::wrap (colvarvalue &x) const
 {
   return;
 }
 
 
 /// \brief Colvar component: distance between the centers of mass of
 /// two groups (colvarvalue::type_scalar type, range [0:*))
 ///
 /// This class also serves as the template for many collective
 /// variables with two atom groups: in this case, the
 /// distance::distance() constructor should be called on the same
 /// configuration string, to make the same member data and functions
 /// available to the derived object
 class colvar::distance
   : public colvar::cvc
 {
 protected:
   /// First atom group
   cvm::atom_group  group1;
   /// Second atom group
   cvm::atom_group  group2;
   /// Vector distance, cached to be recycled
   cvm::rvector     dist_v;
   /// Use absolute positions, ignoring PBCs when present
   bool b_no_PBC;
   /// Compute system force on first site only to avoid unwanted
   /// coupling to other colvars (see e.g. Ciccotti et al., 2005)
   bool b_1site_force;
 public:
   distance (std::string const &conf, bool twogroups = true);
   distance();
   virtual inline ~distance() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void calc_force_invgrads();
   virtual void calc_Jacobian_derivative();
   virtual void apply_force (colvarvalue const &force);
   virtual cvm::real dist2 (colvarvalue const &x1,
                            colvarvalue const &x2) const;
   virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual cvm::real compare (colvarvalue const &x1,
                              colvarvalue const &x2) const;
 };
 
 
 // \brief Colvar component: distance vector between centers of mass
 // of two groups (\link colvarvalue::type_vector \endlink type,
 // range (-*:*)x(-*:*)x(-*:*))
 class colvar::distance_vec
   : public colvar::distance
 {
 public:
   distance_vec (std::string const &conf);
   distance_vec();
   virtual inline ~distance_vec() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force (colvarvalue const &force);
   /// Redefined to handle the box periodicity
   virtual cvm::real dist2 (colvarvalue const &x1,
                            colvarvalue const &x2) const;
   /// Redefined to handle the box periodicity
   virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   /// Redefined to handle the box periodicity
   virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   /// Redefined to handle the box periodicity
   virtual cvm::real compare (colvarvalue const &x1,
                              colvarvalue const &x2) const;
 };
 
 
 /// \brief Colvar component: distance unit vector (direction) between
 /// centers of mass of two groups (colvarvalue::type_unitvector type,
 /// range [-1:1]x[-1:1]x[-1:1])
 class colvar::distance_dir
   : public colvar::distance
 {
 public:
   distance_dir (std::string const &conf);
   distance_dir();
   virtual inline ~distance_dir() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force (colvarvalue const &force);
   virtual cvm::real dist2 (colvarvalue const &x1,
                            colvarvalue const &x2) const;
   virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual cvm::real compare (colvarvalue const &x1,
                              colvarvalue const &x2) const;
 };
 
 
 /// \brief Colvar component: projection of the distance vector along
 /// an axis (colvarvalue::type_scalar type, range (-*:*))
 class colvar::distance_z
   : public colvar::cvc
 {
 protected:
   /// Main atom group
   cvm::atom_group  main;
   /// Reference atom group
   cvm::atom_group  ref1;
   /// Optional, second ref atom group
   cvm::atom_group  ref2;
   /// Use absolute positions, ignoring PBCs when present
   bool b_no_PBC;
   /// Compute system force on one site only to avoid unwanted
   /// coupling to other colvars (see e.g. Ciccotti et al., 2005)
   bool b_1site_force;
   /// Vector on which the distance vector is projected
   cvm::rvector axis;
   /// Norm of the axis
   cvm::real axis_norm;
   /// Vector distance, cached to be recycled
   cvm::rvector     dist_v;
   /// Flag: using a fixed axis vector?
   bool fixed_axis;
 public:
   distance_z (std::string const &conf);
   distance_z();
   virtual inline ~distance_z() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void calc_force_invgrads();
   virtual void calc_Jacobian_derivative();
   virtual void apply_force (colvarvalue const &force);
   virtual cvm::real dist2 (colvarvalue const &x1,
                            colvarvalue const &x2) const;
   virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual cvm::real compare (colvarvalue const &x1,
                              colvarvalue const &x2) const;
   /// \brief Redefined to make use of the user-provided period
   virtual void wrap (colvarvalue &x) const;
 };
 
 
 /// \brief Colvar component: projection of the distance vector on a
 /// plane (colvarvalue::type_scalar type, range [0:*))
 class colvar::distance_xy
   : public colvar::distance_z
 {
 protected:
   /// Components of the distance vector orthogonal to the axis
   cvm::rvector dist_v_ortho;
   /// Vector distances
   cvm::rvector v12, v13;
 public:
   distance_xy (std::string const &conf);
   distance_xy();
   virtual inline ~distance_xy() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void calc_force_invgrads();
   virtual void calc_Jacobian_derivative();
   virtual void apply_force (colvarvalue const &force);
   virtual cvm::real dist2 (colvarvalue const &x1,
                            colvarvalue const &x2) const;
   virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual cvm::real compare (colvarvalue const &x1,
                              colvarvalue const &x2) const;
 };
 
 
 /// \brief Colvar component: average distance between two groups of atoms, weighted as the sixth power,
 /// as in NMR refinements (colvarvalue::type_scalar type, range (0:*))
 class colvar::distance_inv
   : public colvar::distance
 {
 protected:
   /// Components of the distance vector orthogonal to the axis
   int exponent;
 public:
   distance_inv (std::string const &conf);
   distance_inv();
   virtual inline ~distance_inv() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force (colvarvalue const &force);
   virtual cvm::real dist2 (colvarvalue const &x1,
                            colvarvalue const &x2) const;
   virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual cvm::real compare (colvarvalue const &x1,
                              colvarvalue const &x2) const;
 };
 
 
 
 /// \brief Colvar component: Radius of gyration of an atom group
 /// (colvarvalue::type_scalar type, range [0:*))
 class colvar::gyration
   : public colvar::cvc
 {
 protected:
   /// Atoms involved
   cvm::atom_group atoms;
 public:
   /// Constructor
   gyration (std::string const &conf);
   gyration();
   virtual inline ~gyration() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void calc_force_invgrads();
   virtual void calc_Jacobian_derivative();
   virtual void apply_force (colvarvalue const &force);
   virtual cvm::real dist2 (colvarvalue const &x1,
                            colvarvalue const &x2) const;
   virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual cvm::real compare (colvarvalue const &x1,
                              colvarvalue const &x2) const;
 };
 
 
 /// \brief Colvar component: moment of inertia of an atom group
 /// (colvarvalue::type_scalar type, range [0:*))
 class colvar::inertia
   : public colvar::gyration
 {
 public:
   /// Constructor
   inertia (std::string const &conf);
   inertia();
   virtual inline ~inertia() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force (colvarvalue const &force);
   virtual cvm::real dist2 (colvarvalue const &x1,
                            colvarvalue const &x2) const;
   virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual cvm::real compare (colvarvalue const &x1,
                              colvarvalue const &x2) const;
 };
 
 
 /// \brief Colvar component: moment of inertia of an atom group
 /// around a user-defined axis (colvarvalue::type_scalar type, range [0:*))
 class colvar::inertia_z
   : public colvar::inertia
 {
 protected:
   /// Vector on which the inertia tensor is projected
   cvm::rvector axis;
 public:
   /// Constructor
   inertia_z (std::string const &conf);
   inertia_z();
   virtual inline ~inertia_z() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force (colvarvalue const &force);
   virtual cvm::real dist2 (colvarvalue const &x1,
                            colvarvalue const &x2) const;
   virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual cvm::real compare (colvarvalue const &x1,
                              colvarvalue const &x2) const;
 };
 
 
 /// \brief Colvar component: projection of 3N coordinates onto an
 /// eigenvector (colvarvalue::type_scalar type, range (-*:*))
 class colvar::eigenvector
   : public colvar::cvc
 {
 protected:
 
   /// Atom group
   cvm::atom_group             atoms;
 
   /// Reference coordinates
   std::vector<cvm::atom_pos>  ref_pos;
 
   /// Geometric center of the reference coordinates
-  cvm::rvector                ref_pos_center;
+  cvm::atom_pos               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
   cvm::real                   eigenvec_invnorm2;
 
 public:
 
   /// Constructor
   eigenvector (std::string const &conf);
   virtual inline ~eigenvector() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void calc_force_invgrads();
   virtual void calc_Jacobian_derivative();
   virtual void apply_force (colvarvalue const &force);
   virtual cvm::real dist2 (colvarvalue const &x1,
                            colvarvalue const &x2) const;
   virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual cvm::real compare (colvarvalue const &x1,
                              colvarvalue const &x2) const;
 };
 
 
 
 /// \brief Colvar component: angle between the centers of mass of
 /// three groups (colvarvalue::type_scalar type, range [0:PI])
 class colvar::angle
   : public colvar::cvc
 {
 protected:
 
   /// Atom group
   cvm::atom_group group1;
   /// Atom group
   cvm::atom_group group2;
   /// Atom group
   cvm::atom_group group3;
 
   /// Inter site vectors
   cvm::rvector r21, r23;
   /// Inter site vector norms
   cvm::real r21l, r23l;
   /// Derivatives wrt group centers of mass
   cvm::rvector dxdr1, dxdr3;
 
   /// Compute system force on first site only to avoid unwanted
   /// coupling to other colvars (see e.g. Ciccotti et al., 2005)
   /// (or to allow dummy atoms)
   bool b_1site_force;
 public:
 
   /// Initialize by parsing the configuration
   angle (std::string const &conf);
   /// \brief Initialize the three groups after three atoms
   angle (cvm::atom const &a1, cvm::atom const &a2, cvm::atom const &a3);
   angle();
   virtual inline ~angle() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void calc_force_invgrads();
   virtual void calc_Jacobian_derivative();
   virtual void apply_force (colvarvalue const &force);
   virtual cvm::real dist2 (colvarvalue const &x1,
                            colvarvalue const &x2) const;
   virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual cvm::real compare (colvarvalue const &x1,
                              colvarvalue const &x2) const;
 };
 
 
 /// \brief Colvar component: dihedral between the centers of mass of
 /// four groups (colvarvalue::type_scalar type, range [-PI:PI])
 class colvar::dihedral
   : public colvar::cvc
 {
 protected:
 
   /// Atom group
   cvm::atom_group group1;
   /// Atom group
   cvm::atom_group group2;
   /// Atom group
   cvm::atom_group group3;
   /// Atom group
   cvm::atom_group group4;
   /// Inter site vectors
   cvm::rvector r12, r23, r34;
 
   /// \brief Compute system force on first site only to avoid unwanted
   /// coupling to other colvars (see e.g. Ciccotti et al., 2005)
   bool b_1site_force;
 
 public:
 
   /// Initialize by parsing the configuration
   dihedral (std::string  const &conf);
   /// \brief Initialize the four groups after four atoms
   dihedral (cvm::atom const &a1, cvm::atom const &a2, cvm::atom const &a3, cvm::atom const &a4);
   dihedral();
   virtual inline ~dihedral() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void calc_force_invgrads();
   virtual void calc_Jacobian_derivative();
   virtual void apply_force (colvarvalue const &force);
 
   /// Redefined to handle the 2*PI periodicity
   virtual cvm::real dist2 (colvarvalue const &x1,
                            colvarvalue const &x2) const;
   /// Redefined to handle the 2*PI periodicity
   virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   /// Redefined to handle the 2*PI periodicity
   virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   /// Redefined to handle the 2*PI periodicity
   virtual cvm::real compare (colvarvalue const &x1,
                              colvarvalue const &x2) const;
   /// Redefined to handle the 2*PI periodicity
   virtual void wrap (colvarvalue &x) const;
 };
 
 
 /// \brief Colvar component: coordination number between two groups
 /// (colvarvalue::type_scalar type, range [0:N1*N2])
 class colvar::coordnum
   : public colvar::distance
 {
 protected:
   /// \brief "Cutoff" for isotropic calculation (default)
   cvm::real     r0;
   /// \brief "Cutoff vector" for anisotropic calculation
   cvm::rvector  r0_vec;
   /// \brief Wheter dist/r0 or \vec{dist}*\vec{1/r0_vec} should ne be
   /// used
   bool b_anisotropic;
   /// Integer exponent of the function numerator
   int en;
   /// Integer exponent of the function denominator
   int ed;
   /// \brief If true, group2 will be treated as a single atom
   /// (default: loop over all pairs of atoms in group1 and group2)
   bool b_group2_center_only;
 public:
   /// Constructor
   coordnum (std::string const &conf);
   coordnum();
   virtual inline ~coordnum() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force (colvarvalue const &force);
   template<bool b_gradients>
   /// \brief Calculate a coordination number through the function
   /// (1-x**n)/(1-x**m), x = |A1-A2|/r0 \param r0 "cutoff" for the
   /// coordination number \param exp_num \i n exponent \param exp_den
   /// \i m exponent \param A1 atom \param A2 atom
   static cvm::real switching_function (cvm::real const &r0,
                                        int const &exp_num, int const &exp_den,
                                        cvm::atom &A1, cvm::atom &A2);
 
   template<bool b_gradients>
   /// \brief Calculate a coordination number through the function
   /// (1-x**n)/(1-x**m), x = |(A1-A2)*(r0_vec)^-|1 \param r0_vec
   /// vector of different cutoffs in the three directions \param
   /// exp_num \i n exponent \param exp_den \i m exponent \param A1
   /// atom \param A2 atom
   static cvm::real switching_function (cvm::rvector const &r0_vec,
                                        int const &exp_num, int const &exp_den,
                                        cvm::atom &A1, cvm::atom &A2);
 
   virtual cvm::real dist2 (colvarvalue const &x1,
                            colvarvalue const &x2) const;
   virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual cvm::real compare (colvarvalue const &x1,
                              colvarvalue const &x2) const;
 };
 
 /// \brief Colvar component: self-coordination number within a group
 /// (colvarvalue::type_scalar type, range [0:N*(N-1)/2])
 class colvar::selfcoordnum
   : public colvar::distance
 {
 protected:
   /// \brief "Cutoff" for isotropic calculation (default)
   cvm::real     r0;
   /// Integer exponent of the function numerator
   int en;
   /// Integer exponent of the function denominator
   int ed;
 public:
   /// Constructor
   selfcoordnum (std::string const &conf);
   selfcoordnum();
   virtual inline ~selfcoordnum() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force (colvarvalue const &force);
   template<bool b_gradients>
   /// \brief Calculate a coordination number through the function
   /// (1-x**n)/(1-x**m), x = |A1-A2|/r0 \param r0 "cutoff" for the
   /// coordination number \param exp_num \i n exponent \param exp_den
   /// \i m exponent \param A1 atom \param A2 atom
   static cvm::real switching_function (cvm::real const &r0,
                                        int const &exp_num, int const &exp_den,
                                        cvm::atom &A1, cvm::atom &A2);
 
   virtual cvm::real dist2 (colvarvalue const &x1,
                            colvarvalue const &x2) const;
   virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual cvm::real compare (colvarvalue const &x1,
                              colvarvalue const &x2) const;
 };
 
 /// \brief Colvar component: hydrogen bond, defined as the product of
 /// a colvar::coordnum and 1/2*(1-cos((180-ang)/ang_tol))
 /// (colvarvalue::type_scalar type, range [0:1])
 class colvar::h_bond
   : public colvar::cvc
 {
 protected:
   /// Atoms involved in the component
   cvm::atom     acceptor, donor;
   /// \brief "Cutoff" distance between acceptor and donor
   cvm::real     r0;
   /// Integer exponent of the function numerator
   int en;
   /// Integer exponent of the function denominator
   int ed;
 public:
   h_bond (std::string const &conf);
   /// Constructor for atoms already allocated
   h_bond (cvm::atom const &acceptor,
           cvm::atom const &donor,
           cvm::real r0, int en, int ed);
   h_bond();
   virtual ~h_bond();
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force (colvarvalue const &force);
 
   virtual cvm::real dist2 (colvarvalue const &x1,
                            colvarvalue const &x2) const;
   virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual cvm::real compare (colvarvalue const &x1,
                              colvarvalue const &x2) const;
 };
 
 
 /// \brief Colvar component: alpha helix content of a contiguous
 /// segment of 5 or more residues, implemented as a sum of phi/psi
 /// dihedral angles and hydrogen bonds (colvarvalue::type_scalar type,
 /// range [0:1])
 // class colvar::alpha_dihedrals
 //   : public colvar::cvc
 // {
 // protected:
 
 //   /// Alpha-helical reference phi value
 //   cvm::real phi_ref;
 
 //   /// Alpha-helical reference psi value
 //   cvm::real psi_ref;
 
 //   /// List of phi dihedral angles
 //   std::vector<dihedral *> phi;
 
 //   /// List of psi dihedral angles
 //   std::vector<dihedral *> psi;
 
 //   /// List of hydrogen bonds
 //   std::vector<h_bond *>   hb;
 
 // public:
 
 //   alpha_dihedrals (std::string const &conf);
 //   alpha_dihedrals();
 //   virtual inline ~alpha_dihedrals() {}
 //   virtual void calc_value();
 //   virtual void calc_gradients();
 //   virtual void apply_force (colvarvalue const &force);
 //   virtual cvm::real dist2 (colvarvalue const &x1,
 //                            colvarvalue const &x2) const;
 //   virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
 //                                    colvarvalue const &x2) const;
 //   virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
 //                                    colvarvalue const &x2) const;
 //   virtual cvm::real compare (colvarvalue const &x1,
 //                              colvarvalue const &x2) const;
 // };
 
 
 /// \brief Colvar component: alpha helix content of a contiguous
 /// segment of 5 or more residues, implemented as a sum of Ca-Ca-Ca
 /// angles and hydrogen bonds (colvarvalue::type_scalar type, range
 /// [0:1])
 class colvar::alpha_angles
   : public colvar::cvc
 {
 protected:
 
   /// Reference Calpha-Calpha angle (default: 88 degrees)
   cvm::real theta_ref;
 
   /// Tolerance on the Calpha-Calpha angle
   cvm::real theta_tol;
 
   /// List of Calpha-Calpha angles
   std::vector<angle *> theta;
 
   /// List of hydrogen bonds
   std::vector<h_bond *>   hb;
 
   /// Contribution of the hb terms
   cvm::real hb_coeff;
 
 public:
 
   alpha_angles (std::string const &conf);
   alpha_angles();
   virtual ~alpha_angles();
   void calc_value();
   void calc_gradients();
   void apply_force (colvarvalue const &force);
   virtual cvm::real dist2 (colvarvalue const &x1,
                            colvarvalue const &x2) const;
   virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual cvm::real compare (colvarvalue const &x1,
                              colvarvalue const &x2) const;
 };
 
 /// \brief Colvar component: dihedPC
 /// Projection of the config onto a dihedral principal component
 /// See e.g. Altis et al., J. Chem. Phys 126, 244111 (2007)
 /// Based on a set of 'dihedral' cvcs
 class colvar::dihedPC
   : public colvar::cvc
 {
 protected:
 
   std::vector<dihedral *> theta;
   std::vector<cvm::real> coeffs;
 
 public:
 
   dihedPC (std::string const &conf);
   dihedPC();
   virtual  ~dihedPC();
   void calc_value();
   void calc_gradients();
   void apply_force (colvarvalue const &force);
   virtual cvm::real dist2 (colvarvalue const &x1,
                            colvarvalue const &x2) const;
   virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual cvm::real compare (colvarvalue const &x1,
                              colvarvalue const &x2) const;
 };
 
 /// \brief Colvar component: orientation in space of an atom group,
 /// with respect to a set of reference coordinates
 /// (colvarvalue::type_quaternion type, range
 /// [-1:1]x[-1:1]x[-1:1]x[-1:1])
 class colvar::orientation
   : public colvar::cvc
 {
 protected:
 
   /// Atom group
   cvm::atom_group            atoms;
   /// Center of geometry of the group
   cvm::atom_pos              atoms_cog;
 
   /// Reference coordinates
   std::vector<cvm::atom_pos> ref_pos;
 
   /// Rotation object
   cvm::rotation              rot;
 
   /// \brief This is used to remove jumps in the sign of the
   /// quaternion, which may be annoying in the colvars trajectory
   cvm::quaternion            ref_quat;
 
 public:
 
   orientation (std::string const &conf);
   orientation();
   virtual inline ~orientation() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force (colvarvalue const &force);
   virtual cvm::real dist2 (colvarvalue const &x1,
                            colvarvalue const &x2) const;
   virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual cvm::real compare (colvarvalue const &x1,
                              colvarvalue const &x2) const;
 };
 
 
 /// \brief Colvar component: angle of rotation with respect to a set
 /// of reference coordinates (colvarvalue::type_scalar type, range
 /// [0:PI))
 class colvar::orientation_angle
   : public colvar::orientation
 {
 public:
 
   orientation_angle (std::string const &conf);
   orientation_angle();
   virtual inline ~orientation_angle() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force (colvarvalue const &force);
   virtual cvm::real dist2 (colvarvalue const &x1,
                            colvarvalue const &x2) const;
   virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual cvm::real compare (colvarvalue const &x1,
                              colvarvalue const &x2) const;
 };
 
 
 /// \brief Colvar component: projection of the orientation vector onto
 /// a predefined axis (colvarvalue::type_scalar type, range [-1:1])
 class colvar::tilt
   : public colvar::orientation
 {
 protected:
 
   cvm::rvector axis;
 
 public:
 
   tilt (std::string const &conf);
   tilt();
   virtual inline ~tilt() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force (colvarvalue const &force);
   virtual cvm::real dist2 (colvarvalue const &x1,
                            colvarvalue const &x2) const;
   virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual cvm::real compare (colvarvalue const &x1,
                              colvarvalue const &x2) const;
 };
 
 
 
 /// \brief Colvar component: angle of rotation around a predefined
 /// axis (colvarvalue::type_scalar type, range [-PI:PI])
 class colvar::spin_angle
   : public colvar::orientation
 {
 protected:
 
   cvm::rvector axis;
 
 public:
 
   spin_angle (std::string const &conf);
   spin_angle();
   virtual inline ~spin_angle() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force (colvarvalue const &force);
   /// Redefined to handle the 2*PI periodicity
   virtual cvm::real dist2 (colvarvalue const &x1,
                            colvarvalue const &x2) const;
   /// Redefined to handle the 2*PI periodicity
   virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   /// Redefined to handle the 2*PI periodicity
   virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   /// Redefined to handle the 2*PI periodicity
   virtual cvm::real compare (colvarvalue const &x1,
                              colvarvalue const &x2) const;
   /// Redefined to handle the 2*PI periodicity
   virtual void wrap (colvarvalue &x) const;
 };
 
 
 
 /// \brief Colvar component: root mean square deviation (RMSD) of a
 /// group with respect to a set of reference coordinates; uses \link
 /// colvar::orientation \endlink to calculate the rotation matrix
 /// (colvarvalue::type_scalar type, range [0:*))
 class colvar::rmsd
   : public colvar::cvc
 {
 protected:
 
   /// Atom group
   cvm::atom_group             atoms;
 
   /// Reference coordinates (for RMSD calculation only)
   std::vector<cvm::atom_pos>  ref_pos;
 
 public:
 
   /// Constructor
   rmsd (std::string const &conf);
   virtual inline ~rmsd() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void calc_force_invgrads();
   virtual void calc_Jacobian_derivative();
   virtual void apply_force (colvarvalue const &force);
   virtual cvm::real dist2 (colvarvalue const &x1,
                            colvarvalue const &x2) const;
   virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                    colvarvalue const &x2) const;
   virtual cvm::real compare (colvarvalue const &x1,
                              colvarvalue const &x2) const;
 };
 
 
 
 // metrics functions for cvc implementations
 
 // simple definitions of the distance functions; these are useful only
 // for optimization (the type check performed in the default
 // colvarcomp functions is skipped)
 
 // definitions assuming the scalar type
 
 #define simple_scalar_dist_functions(TYPE)                              \
                                                                         \
   inline cvm::real colvar::TYPE::dist2 (colvarvalue const &x1,          \
                                         colvarvalue const &x2) const    \
   {                                                                     \
     return (x1.real_value - x2.real_value)*(x1.real_value - x2.real_value); \
   }                                                                     \
                                                                         \
   inline colvarvalue colvar::TYPE::dist2_lgrad (colvarvalue const &x1,  \
                                                 colvarvalue const &x2) const \
   {                                                                     \
     return 2.0 * (x1.real_value - x2.real_value);                       \
   }                                                                     \
                                                                         \
   inline colvarvalue colvar::TYPE::dist2_rgrad (colvarvalue const &x1,  \
                                                 colvarvalue const &x2) const \
   {                                                                     \
     return this->dist2_lgrad (x2, x1);                                  \
   }                                                                     \
                                                                         \
   inline cvm::real colvar::TYPE::compare (colvarvalue const &x1,        \
                                           colvarvalue const &x2) const  \
   {                                                                     \
     return this->dist2_lgrad (x1, x2);                                  \
   }                                                                     \
                                                                         \
 
   simple_scalar_dist_functions (distance)
   // NOTE: distance_z has explicit functions, see below
   simple_scalar_dist_functions (distance_xy)
   simple_scalar_dist_functions (distance_inv)
   simple_scalar_dist_functions (angle)
   simple_scalar_dist_functions (coordnum)
   simple_scalar_dist_functions (selfcoordnum)
   simple_scalar_dist_functions (h_bond)
   simple_scalar_dist_functions (gyration)
   simple_scalar_dist_functions (inertia)
   simple_scalar_dist_functions (inertia_z)
   simple_scalar_dist_functions (rmsd)
   simple_scalar_dist_functions (orientation_angle)
   simple_scalar_dist_functions (tilt)
   simple_scalar_dist_functions (eigenvector)
   //  simple_scalar_dist_functions (alpha_dihedrals)
   simple_scalar_dist_functions (alpha_angles)
   simple_scalar_dist_functions (dihedPC)
 
 
 // metrics functions for cvc implementations with a periodicity
 
 inline cvm::real colvar::dihedral::dist2 (colvarvalue const &x1,
                                           colvarvalue const &x2) const
 {
   cvm::real diff = x1.real_value - x2.real_value;
   diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
   return diff * diff;
 }
 
 inline colvarvalue colvar::dihedral::dist2_lgrad (colvarvalue const &x1,
                                                   colvarvalue const &x2) const
 {
   cvm::real diff = x1.real_value - x2.real_value;
   diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
   return 2.0 * diff;
 }
 
 inline colvarvalue colvar::dihedral::dist2_rgrad (colvarvalue const &x1,
                                                   colvarvalue const &x2) const
 {
   cvm::real diff = x1.real_value - x2.real_value;
   diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
   return (-2.0) * diff;
 }
 
 inline cvm::real colvar::dihedral::compare (colvarvalue const &x1,
                                             colvarvalue const &x2) const
 {
   return dist2_lgrad (x1, x2);
 }
 
 inline void colvar::dihedral::wrap (colvarvalue &x) const
 {
   if ((x.real_value - wrap_center) >= 180.0) {
     x.real_value -= 360.0;
     return;
   }
 
   if ((x.real_value - wrap_center) < -180.0) {
     x.real_value += 360.0;
     return;
   }
 
   return;
 }
 
 inline cvm::real colvar::spin_angle::dist2 (colvarvalue const &x1,
                                           colvarvalue const &x2) const
 {
   cvm::real diff = x1.real_value - x2.real_value;
   diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
   return diff * diff;
 }
 
 inline colvarvalue colvar::spin_angle::dist2_lgrad (colvarvalue const &x1,
                                                   colvarvalue const &x2) const
 {
   cvm::real diff = x1.real_value - x2.real_value;
   diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
   return 2.0 * diff;
 }
 
 inline colvarvalue colvar::spin_angle::dist2_rgrad (colvarvalue const &x1,
                                                   colvarvalue const &x2) const
 {
   cvm::real diff = x1.real_value - x2.real_value;
   diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
   return (-2.0) * diff;
 }
 
 inline cvm::real colvar::spin_angle::compare (colvarvalue const &x1,
                                             colvarvalue const &x2) const
 {
   return dist2_lgrad (x1, x2);
 }
 
 inline void colvar::spin_angle::wrap (colvarvalue &x) const
 {
   if ((x.real_value - wrap_center) >= 180.0) {
     x.real_value -= 360.0;
     return;
   }
 
   if ((x.real_value - wrap_center) < -180.0) {
     x.real_value += 360.0;
     return;
   }
 
   return;
 }
 
 
 // Projected distance
 // Differences should always be wrapped around 0 (ignoring wrap_center)
 inline cvm::real colvar::distance_z::dist2 (colvarvalue const &x1,
                                             colvarvalue const &x2) const
 {
   cvm::real diff = x1.real_value - x2.real_value;
   if (period != 0.0) {
     cvm::real shift = std::floor (diff/period + 0.5);
     diff -= shift * period;
   }
   return diff * diff;
 }
 
 inline colvarvalue colvar::distance_z::dist2_lgrad (colvarvalue const &x1,
                                                     colvarvalue const &x2) const
 {
   cvm::real diff = x1.real_value - x2.real_value;
   if (period != 0.0) {
     cvm::real shift = std::floor (diff/period + 0.5);
     diff -= shift * period;
   }
   return 2.0 * diff;
 }
 
 inline colvarvalue colvar::distance_z::dist2_rgrad (colvarvalue const &x1,
                                                     colvarvalue const &x2) const
 {
   cvm::real diff = x1.real_value - x2.real_value;
   if (period != 0.0) {
     cvm::real shift = std::floor (diff/period + 0.5);
     diff -= shift * period;
   }
   return (-2.0) * diff;
 }
 
 inline cvm::real colvar::distance_z::compare (colvarvalue const &x1,
                                               colvarvalue const &x2) const
 {
   return dist2_lgrad (x1, x2);
 }
 
 inline void colvar::distance_z::wrap (colvarvalue &x) const
 {
   if (! this->b_periodic) {
     // don't wrap if the period has not been set
     return;
   }
 
   cvm::real shift = std::floor ((x.real_value - wrap_center) / period + 0.5);
   x.real_value -= shift * period;
   return;
 }
 
 
 // distance between three dimensional vectors
 //
 // TODO apply PBC to distance_vec
 // Note: differences should be centered around (0, 0, 0)!
 
 inline cvm::real colvar::distance_vec::dist2 (colvarvalue const &x1,
                                               colvarvalue const &x2) const
 {
   return cvm::position_dist2 (x1.rvector_value, x2.rvector_value);
 }
 
 inline colvarvalue colvar::distance_vec::dist2_lgrad (colvarvalue const &x1,
                                                       colvarvalue const &x2) const
 {
   return 2.0 * cvm::position_distance(x2.rvector_value, x1.rvector_value);
 }
 
 inline colvarvalue colvar::distance_vec::dist2_rgrad (colvarvalue const &x1,
                                                       colvarvalue const &x2) const
 {
   return 2.0 * cvm::position_distance(x2.rvector_value, x1.rvector_value);
 }
 
 inline cvm::real colvar::distance_vec::compare (colvarvalue const &x1,
                                                 colvarvalue const &x2) const
 {
   cvm::fatal_error ("Error: cannot compare() two distance vectors.\n");
   return 0.0;
 }
 
 inline cvm::real colvar::distance_dir::dist2 (colvarvalue const &x1,
                                               colvarvalue const &x2) const
 {
   return (x1.rvector_value - x2.rvector_value).norm2();
 }
 
 inline colvarvalue colvar::distance_dir::dist2_lgrad (colvarvalue const &x1,
                                                       colvarvalue const &x2) const
 {
   return colvarvalue ((x1.rvector_value - x2.rvector_value), colvarvalue::type_unitvector);
 }
 
 inline colvarvalue colvar::distance_dir::dist2_rgrad (colvarvalue const &x1,
                                                       colvarvalue const &x2) const
 {
   return colvarvalue ((x2.rvector_value - x1.rvector_value), colvarvalue::type_unitvector);
 }
 
 inline cvm::real colvar::distance_dir::compare (colvarvalue const &x1,
                                                 colvarvalue const &x2) const
 {
   cvm::fatal_error ("Error: cannot compare() two distance directions.\n");
   return 0.0;
 }
 
 // distance between quaternions
 
 inline cvm::real colvar::orientation::dist2 (colvarvalue const &x1,
                                              colvarvalue const &x2) const
 {
   return x1.quaternion_value.dist2 (x2);
 }
 
 inline colvarvalue colvar::orientation::dist2_lgrad (colvarvalue const &x1,
                                                      colvarvalue const &x2) const
 {
   return x1.quaternion_value.dist2_grad (x2);
 }
 
 inline colvarvalue colvar::orientation::dist2_rgrad (colvarvalue const &x1,
                                                      colvarvalue const &x2) const
 {
   return x2.quaternion_value.dist2_grad (x1);
 }
 
 inline cvm::real colvar::orientation::compare (colvarvalue const &x1,
                                                colvarvalue const &x2) const
 {
   cvm::fatal_error ("Error: cannot compare() two quaternions.\n");
   return 0.0;
 }
 
 
 #endif
 
 
 // Emacs
 // Local Variables:
 // mode: C++
 // End:
diff --git a/lib/colvars/colvarcomp_distances.cpp b/lib/colvars/colvarcomp_distances.cpp
index f59bab8a9..339ba95f7 100644
--- a/lib/colvars/colvarcomp_distances.cpp
+++ b/lib/colvars/colvarcomp_distances.cpp
@@ -1,1151 +1,1151 @@
 #include <cmath>
 
 #include "colvarmodule.h"
 #include "colvarvalue.h"
 #include "colvarparse.h"
 #include "colvar.h"
 #include "colvarcomp.h"
 
 
 
 // "twogroup" flag defaults to true; set to false by selfCoordNum
 // (only distance-derived component based on only one group)
 
 colvar::distance::distance (std::string const &conf, bool twogroups)
   : cvc (conf)
 {
   function_type = "distance";
   b_inverse_gradients = true;
   b_Jacobian_derivative = true;
   if (get_keyval (conf, "forceNoPBC", b_no_PBC, false)) {
     cvm::log ("Computing distance using absolute positions (not minimal-image)");
   }
   if (twogroups && get_keyval (conf, "oneSiteSystemForce", b_1site_force, false)) {
     cvm::log ("Computing system force on group 1 only");
   }
   parse_group (conf, "group1", group1);
   atom_groups.push_back (&group1);
   if (twogroups) {
     parse_group (conf, "group2", group2);
     atom_groups.push_back (&group2);
   }
   x.type (colvarvalue::type_scalar);
 }
 
 
 colvar::distance::distance()
   : cvc ()
 {
   function_type = "distance";
   b_inverse_gradients = true;
   b_Jacobian_derivative = true;
   b_1site_force = false;
   x.type (colvarvalue::type_scalar);
 }
 
 void colvar::distance::calc_value()
 {
   if (b_no_PBC) {
     dist_v = group2.center_of_mass() - group1.center_of_mass();
   } else {
     dist_v = cvm::position_distance (group1.center_of_mass(),
                                      group2.center_of_mass());
   }
   x.real_value = dist_v.norm();
 }
 
 void colvar::distance::calc_gradients()
 {
   cvm::rvector const u = dist_v.unit();
   group1.set_weighted_gradient (-1.0 * u);
   group2.set_weighted_gradient (       u);
 }
 
 void colvar::distance::calc_force_invgrads()
 {
   group1.read_system_forces();
   if ( b_1site_force ) {
     ft.real_value = -1.0 * (group1.system_force() * dist_v.unit());
   } else {
     group2.read_system_forces();
     ft.real_value = 0.5 * ((group2.system_force() - group1.system_force()) * dist_v.unit());
   }
 }
 
 void colvar::distance::calc_Jacobian_derivative()
 {
   jd.real_value = x.real_value ? (2.0 / x.real_value) : 0.0;
 }
 
 void colvar::distance::apply_force (colvarvalue const &force)
 {
   if (!group1.noforce)
     group1.apply_colvar_force (force.real_value);
 
   if (!group2.noforce)
     group2.apply_colvar_force (force.real_value);
 }
 
 
 
 colvar::distance_vec::distance_vec (std::string const &conf)
   : distance (conf)
 {
   function_type = "distance_vec";
   x.type (colvarvalue::type_vector);
 }
 
 colvar::distance_vec::distance_vec()
   : distance()
 {
   function_type = "distance_vec";
   x.type (colvarvalue::type_vector);
 }
 
 void colvar::distance_vec::calc_value()
 {
   if (b_no_PBC) {
     x.rvector_value = group2.center_of_mass() - group1.center_of_mass();
   } else {
     x.rvector_value = cvm::position_distance (group1.center_of_mass(),
                                               group2.center_of_mass());
   }
 }
 
 void colvar::distance_vec::calc_gradients()
 {
   // gradients are not stored: a 3x3 matrix for each atom would be
   // needed to store just the identity matrix
 }
 
 void colvar::distance_vec::apply_force (colvarvalue const &force)
 {
   if (!group1.noforce)
     group1.apply_force (-1.0 * force.rvector_value);
 
   if (!group2.noforce)
     group2.apply_force (       force.rvector_value);
 }
 
 
 
 colvar::distance_z::distance_z (std::string const &conf)
   : cvc (conf)
 {
   function_type = "distance_z";
   b_inverse_gradients = true;
   b_Jacobian_derivative = true;
   x.type (colvarvalue::type_scalar);
 
   // TODO detect PBC from MD engine (in simple cases)
   // and then update period in real time
   if (period != 0.0)
     b_periodic = true;
 
   if ((wrap_center != 0.0) && (period == 0.0)) {
     cvm::fatal_error ("Error: wrapAround was defined in a distanceZ component,"
                       " but its period has not been set.\n");
   }
 
   parse_group (conf, "main", main);
   parse_group (conf, "ref", ref1);
   atom_groups.push_back (&main);
   atom_groups.push_back (&ref1);
   // this group is optional
   parse_group (conf, "ref2", ref2, true);
 
   if (ref2.size()) {
     atom_groups.push_back (&ref2);
     cvm::log ("Using axis joining the centers of mass of groups \"ref\" and \"ref2\"");
     fixed_axis = false;
     if (key_lookup (conf, "axis"))
       cvm::log ("Warning: explicit axis definition will be ignored!");
   } else {
     if (get_keyval (conf, "axis", axis, cvm::rvector (0.0, 0.0, 1.0))) {
       if (axis.norm2() == 0.0)
         cvm::fatal_error ("Axis vector is zero!");
       if (axis.norm2() != 1.0) {
         axis = axis.unit();
         cvm::log ("The normalized axis is: "+cvm::to_str (axis)+".\n");
       }
     }
     fixed_axis = true;
   }
 
   if (get_keyval (conf, "forceNoPBC", b_no_PBC, false)) {
     cvm::log ("Computing distance using absolute positions (not minimal-image)");
   }
   if (get_keyval (conf, "oneSiteSystemForce", b_1site_force, false)) {
     cvm::log ("Computing system force on group \"main\" only");
   }
 }
 
 colvar::distance_z::distance_z()
 {
   function_type = "distance_z";
   b_inverse_gradients = true;
   b_Jacobian_derivative = true;
   x.type (colvarvalue::type_scalar);
 }
 
 void colvar::distance_z::calc_value()
 {
   if (fixed_axis) {
     if (b_no_PBC) {
       dist_v = main.center_of_mass() - ref1.center_of_mass();
     } else {
       dist_v = cvm::position_distance (ref1.center_of_mass(),
                                        main.center_of_mass());
     }
   } else {
 
     if (b_no_PBC) {
       dist_v = main.center_of_mass() -
                (0.5 * (ref1.center_of_mass() + ref2.center_of_mass()));
       axis = ref2.center_of_mass() - ref1.center_of_mass();
     } else {
       dist_v = cvm::position_distance (0.5 * (ref1.center_of_mass() +
                ref2.center_of_mass()), main.center_of_mass());
       axis = cvm::position_distance (ref1.center_of_mass(), ref2.center_of_mass());
     }
     axis_norm = axis.norm();
     axis = axis.unit();
   }
   x.real_value = axis * dist_v;
   this->wrap (x);
 }
 
 void colvar::distance_z::calc_gradients()
 {
   main.set_weighted_gradient ( axis );
 
   if (fixed_axis) {
     ref1.set_weighted_gradient (-1.0 * axis);
   } else {
     if (b_no_PBC) {
       ref1.set_weighted_gradient ( 1.0 / axis_norm * (main.center_of_mass() - ref2.center_of_mass() -
                                    x.real_value * axis ));
       ref2.set_weighted_gradient ( 1.0 / axis_norm * (ref1.center_of_mass() - main.center_of_mass() +
                                    x.real_value * axis ));
     } else {
       ref1.set_weighted_gradient ( 1.0 / axis_norm * (
         cvm::position_distance (ref2.center_of_mass(), main.center_of_mass()) - x.real_value * axis ));
       ref2.set_weighted_gradient ( 1.0 / axis_norm * (
         cvm::position_distance (main.center_of_mass(), ref1.center_of_mass()) + x.real_value * axis ));
     }
   }
 }
 
 void colvar::distance_z::calc_force_invgrads()
 {
   main.read_system_forces();
 
   if (fixed_axis && !b_1site_force) {
     ref1.read_system_forces();
     ft.real_value = 0.5 * ((main.system_force() - ref1.system_force()) * axis);
   } else {
     ft.real_value = main.system_force() * axis;
   }
 }
 
 void colvar::distance_z::calc_Jacobian_derivative()
 {
   jd.real_value = 0.0;
 }
 
 void colvar::distance_z::apply_force (colvarvalue const &force)
 {
   if (!ref1.noforce)
     ref1.apply_colvar_force (force.real_value);
 
   if (ref2.size() && !ref2.noforce)
     ref2.apply_colvar_force (force.real_value);
 
   if (!main.noforce)
     main.apply_colvar_force (force.real_value);
 }
 
 
 
 colvar::distance_xy::distance_xy (std::string const &conf)
   : distance_z (conf)
 {
   function_type = "distance_xy";
   b_inverse_gradients = true;
   b_Jacobian_derivative = true;
   x.type (colvarvalue::type_scalar);
 }
 
 colvar::distance_xy::distance_xy()
   : distance_z()
 {
   function_type = "distance_xy";
   b_inverse_gradients = true;
   b_Jacobian_derivative = true;
   x.type (colvarvalue::type_scalar);
 }
 
 void colvar::distance_xy::calc_value()
 {
   if (b_no_PBC) {
     dist_v = main.center_of_mass() - ref1.center_of_mass();
   } else {
     dist_v = cvm::position_distance (ref1.center_of_mass(),
                                      main.center_of_mass());
   }
   if (!fixed_axis) {
     if (b_no_PBC) {
       v12 = ref2.center_of_mass() - ref1.center_of_mass();
     } else {
       v12 = cvm::position_distance (ref1.center_of_mass(), ref2.center_of_mass());
     }
     axis_norm = v12.norm();
     axis = v12.unit();
   }
 
   dist_v_ortho = dist_v - (dist_v * axis) * axis;
   x.real_value = dist_v_ortho.norm();
 }
 
 void colvar::distance_xy::calc_gradients()
 {
   // Intermediate quantity (r_P3 / r_12 where P is the projection
   // of 3 (main) on the plane orthogonal to 12, containing 1 (ref1))
   cvm::real A;
   cvm::real x_inv;
 
   if (x.real_value == 0.0) return;
   x_inv = 1.0 / x.real_value;
 
   if (fixed_axis) {
     ref1.set_weighted_gradient (-1.0 * x_inv * dist_v_ortho);
     main.set_weighted_gradient (       x_inv * dist_v_ortho);
   } else {
     if (b_no_PBC) {
       v13 = main.center_of_mass() - ref1.center_of_mass();
     } else {
       v13 = cvm::position_distance (ref1.center_of_mass(), main.center_of_mass());
     }
     A = (dist_v * axis) / axis_norm;
 
     ref1.set_weighted_gradient ( (A - 1.0) * x_inv * dist_v_ortho);
     ref2.set_weighted_gradient ( -A        * x_inv * dist_v_ortho);
     main.set_weighted_gradient (      1.0  * x_inv * dist_v_ortho);
   }
 }
 
 void colvar::distance_xy::calc_force_invgrads()
 {
   main.read_system_forces();
 
   if (fixed_axis && !b_1site_force) {
     ref1.read_system_forces();
     ft.real_value = 0.5 / x.real_value * ((main.system_force() - ref1.system_force()) * dist_v_ortho);
   } else {
     ft.real_value = 1.0 / x.real_value * main.system_force() * dist_v_ortho;
   }
 }
 
 void colvar::distance_xy::calc_Jacobian_derivative()
 {
   jd.real_value = x.real_value ? (1.0 / x.real_value) : 0.0;
 }
 
 void colvar::distance_xy::apply_force (colvarvalue const &force)
 {
   if (!ref1.noforce)
     ref1.apply_colvar_force (force.real_value);
 
   if (ref2.size() && !ref2.noforce)
     ref2.apply_colvar_force (force.real_value);
 
   if (!main.noforce)
     main.apply_colvar_force (force.real_value);
 }
 
 
 
 colvar::distance_dir::distance_dir (std::string const &conf)
   : distance (conf)
 {
   function_type = "distance_dir";
   x.type (colvarvalue::type_unitvector);
 }
 
 
 colvar::distance_dir::distance_dir()
   : distance()
 {
   function_type = "distance_dir";
   x.type (colvarvalue::type_unitvector);
 }
 
 
 void colvar::distance_dir::calc_value()
 {
   if (b_no_PBC) {
     dist_v = group2.center_of_mass() - group1.center_of_mass();
   } else {
     dist_v = cvm::position_distance (group1.center_of_mass(),
                                      group2.center_of_mass());
   }
   x.rvector_value = dist_v.unit();
 }
 
 
 void colvar::distance_dir::calc_gradients()
 {
   // gradients are computed on the fly within apply_force()
   // Note: could be a problem if a future bias relies on gradient
   // calculations...
 }
 
 
 void colvar::distance_dir::apply_force (colvarvalue const &force)
 {
   // remove the radial force component
   cvm::real const iprod = force.rvector_value * x.rvector_value;
   cvm::rvector const force_tang = force.rvector_value - iprod * x.rvector_value;
 
   if (!group1.noforce)
     group1.apply_force (-1.0 * force_tang);
 
   if (!group2.noforce)
     group2.apply_force (       force_tang);
 }
 
 
 
 colvar::distance_inv::distance_inv (std::string const &conf)
   : distance (conf)
 {
   function_type = "distance_inv";
   get_keyval (conf, "exponent", exponent, 6);
   if (exponent%2) {
     cvm::fatal_error ("Error: odd exponent provided, can only use even ones.\n");
   }
   if (exponent <= 0) {
     cvm::fatal_error ("Error: negative or zero exponent provided.\n");
   }
 
   for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) {
     for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
       if (ai1->id == ai2->id)
         cvm::fatal_error ("Error: group1 and group1 have some atoms in common: this is not allowed for distanceInv.\n");
     }
   }
 
   b_inverse_gradients = false;
   b_Jacobian_derivative = false;
   x.type (colvarvalue::type_scalar);
 }
 
 colvar::distance_inv::distance_inv()
 {
   function_type = "distance_inv";
   exponent = 6;
   b_inverse_gradients = false;
   b_Jacobian_derivative = false;
   b_1site_force = false;
   x.type (colvarvalue::type_scalar);
 }
 
 void colvar::distance_inv::calc_value()
 {
   x.real_value = 0.0;
   if (b_no_PBC) {
     for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) {
       for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
         cvm::rvector const dv = ai2->pos - ai1->pos;
         cvm::real const d2 = dv.norm2();
         cvm::real dinv = 1.0;
         for (int ne = 0; ne < exponent/2; ne++)
           dinv *= 1.0/d2;
         x.real_value += dinv;
         cvm::rvector const dsumddv = -(cvm::real (exponent)) * dinv/d2 * dv;
         ai1->grad += -1.0 * dsumddv;
         ai2->grad +=        dsumddv;
       }
     }
   } else {
     for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) {
       for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
         cvm::rvector const dv = cvm::position_distance (ai1->pos, ai2->pos);
         cvm::real const d2 = dv.norm2();
         cvm::real dinv = 1.0;
         for (int ne = 0; ne < exponent/2; ne++)
           dinv *= 1.0/d2;
         x.real_value += dinv;
         cvm::rvector const dsumddv = -(cvm::real (exponent)) * dinv/d2 * dv;
         ai1->grad += -1.0 * dsumddv;
         ai2->grad +=        dsumddv;
       }
     }
   }
 
   x.real_value *= 1.0 / cvm::real (group1.size() * group2.size());
   x.real_value = std::pow (x.real_value, -1.0/(cvm::real (exponent)));
 }
 
 void colvar::distance_inv::calc_gradients()
 {
   cvm::real const dxdsum = (-1.0/(cvm::real (exponent))) * std::pow (x.real_value, exponent+1) / cvm::real (group1.size() * group2.size());
   for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) {
     ai1->grad *= dxdsum;
   }
   for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
     ai2->grad *= dxdsum;
   }
 }
 
 void colvar::distance_inv::apply_force (colvarvalue const &force)
 {
   if (!group1.noforce)
     group1.apply_colvar_force (force.real_value);
 
   if (!group2.noforce)
     group2.apply_colvar_force (force.real_value);
 }
 
 
 
 colvar::gyration::gyration (std::string const &conf)
   : cvc (conf)
 {
   function_type = "gyration";
   b_inverse_gradients = true;
   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 {
     atoms.b_center = true;
-    atoms.ref_pos.assign (1, cvm::rvector (0.0, 0.0, 0.0));
+    atoms.ref_pos.assign (1, cvm::atom_pos (0.0, 0.0, 0.0));
   }
 
   x.type (colvarvalue::type_scalar);
 }
 
 
 colvar::gyration::gyration()
 {
   function_type = "gyration";
   b_inverse_gradients = true;
   b_Jacobian_derivative = true;
   x.type (colvarvalue::type_scalar);
 }
 
 
 void colvar::gyration::calc_value()
 {
   x.real_value = 0.0;
   for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
     x.real_value += (ai->pos).norm2();
   }
   x.real_value = std::sqrt (x.real_value / cvm::real (atoms.size()));
 }
 
 
 void colvar::gyration::calc_gradients()
 {
   cvm::real const drdx = 1.0/(cvm::real (atoms.size()) * x.real_value);
   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);
   }
 }
 
 
 void colvar::gyration::calc_force_invgrads()
 {
   atoms.read_system_forces();
 
   cvm::real const dxdr = 1.0/x.real_value;
   ft.real_value = 0.0;
 
   for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
     ft.real_value += dxdr * ai->pos * ai->system_force;
   }
 }
 
 
 void colvar::gyration::calc_Jacobian_derivative()
 {
   jd = x.real_value ? (3.0 * cvm::real (atoms.size()) - 4.0) / x.real_value : 0.0;
 }
 
 
 void colvar::gyration::apply_force (colvarvalue const &force)
 {
   if (!atoms.noforce)
     atoms.apply_colvar_force (force.real_value);
 }
 
 
 
 colvar::inertia::inertia (std::string const &conf)
   : gyration (conf)
 {
   function_type = "inertia";
   b_inverse_gradients = false;
   b_Jacobian_derivative = false;
   x.type (colvarvalue::type_scalar);
 }
 
 
 colvar::inertia::inertia()
 {
   function_type = "inertia";
   x.type (colvarvalue::type_scalar);
 }
 
 
 void colvar::inertia::calc_value()
 {
   x.real_value = 0.0;
   for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
     x.real_value += (ai->pos).norm2();
   }
 }
 
 
 void colvar::inertia::calc_gradients()
 {
   for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
     ai->grad = 2.0 * ai->pos;
   }
 
   if (b_debug_gradients) {
     cvm::log ("Debugging gradients:\n");
     debug_gradients (atoms);
   }
 }
 
 
 void colvar::inertia::apply_force (colvarvalue const &force)
 {
   if (!atoms.noforce)
     atoms.apply_colvar_force (force.real_value);
 }
 
 
 colvar::inertia_z::inertia_z (std::string const &conf)
   : inertia (conf)
 {
   function_type = "inertia_z";
   if (get_keyval (conf, "axis", axis, cvm::rvector (0.0, 0.0, 1.0))) {
     if (axis.norm2() == 0.0)
       cvm::fatal_error ("Axis vector is zero!");
     if (axis.norm2() != 1.0) {
       axis = axis.unit();
       cvm::log ("The normalized axis is: "+cvm::to_str (axis)+".\n");
     }
   }
   x.type (colvarvalue::type_scalar);
 }
 
 
 colvar::inertia_z::inertia_z()
 {
   function_type = "inertia_z";
   x.type (colvarvalue::type_scalar);
 }
 
 
 void colvar::inertia_z::calc_value()
 {
   x.real_value = 0.0;
   for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
     cvm::real const iprod = ai->pos * axis;
     x.real_value += iprod * iprod;
   }
 }
 
 
 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 (b_debug_gradients) {
     cvm::log ("Debugging gradients:\n");
     debug_gradients (atoms);
   }
 }
 
 
 void colvar::inertia_z::apply_force (colvarvalue const &force)
 {
   if (!atoms.noforce)
     atoms.apply_colvar_force (force.real_value);
 }
 
 
 
 colvar::rmsd::rmsd (std::string const &conf)
   : cvc (conf)
 {
   b_inverse_gradients = true;
   b_Jacobian_derivative = true;
   function_type = "rmsd";
   x.type (colvarvalue::type_scalar);
 
   parse_group (conf, "atoms", atoms);
   atom_groups.push_back (&atoms);
 
   if (atoms.b_dummy)
     cvm::fatal_error ("Error: \"atoms\" cannot be a dummy atom.");
 
   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;
   // 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()) {
       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 \"atoms\" ("+
                         cvm::to_str (atoms.size())+").\n");
     }
   }
   {
     std::string ref_pos_file;
     if (get_keyval (conf, "refPositionsFile", ref_pos_file, std::string (""))) {
 
       if (ref_pos.size()) {
         cvm::fatal_error ("Error: cannot specify \"refPositionsFile\" and "
                           "\"refPositions\" at the same time.\n");
       }
 
       std::string ref_pos_col;
       double ref_pos_col_value;
 
       if (get_keyval (conf, "refPositionsCol", ref_pos_col, std::string (""))) {
         // if provided, use PDB column to select coordinates
         bool found = get_keyval (conf, "refPositionsColValue", ref_pos_col_value, 0.0);
         if (found && !ref_pos_col_value)
           cvm::fatal_error ("Error: refPositionsColValue, "
                             "if provided, must be non-zero.\n");
       } else {
         // if not, rely on existing atom indices for the group
         atoms.create_sorted_ids();
       }
 
       ref_pos.resize (atoms.size());
       cvm::load_coords (ref_pos_file.c_str(), ref_pos, atoms.sorted_ids,
                         ref_pos_col, ref_pos_col_value);
     }
   }
 
   if (atoms.b_user_defined_fit) {
     cvm::log ("WARNING: explicit fitting parameters were provided for atom group \"atoms\".");
   } else {
     // Default: fit everything
     cvm::log ("Enabling \"centerReference\" and \"rotateReference\", to minimize RMSD before calculating it as a variable: "
               "if this is not the desired behavior, disable them explicitly within the \"atoms\" block.\n");
     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()
 {
   // 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 = std::sqrt (x.real_value);
 }
 
 
 void colvar::rmsd::calc_gradients()
 {
   cvm::real const drmsddx2 = (x.real_value > 0.0) ?
     0.5 / (x.real_value * cvm::real (atoms.size())) :
     0.0;
 
   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);
   }
 }
 
 
 void colvar::rmsd::apply_force (colvarvalue const &force)
 {
   if (!atoms.noforce)
     atoms.apply_colvar_force (force.real_value);
 }
 
 
 void colvar::rmsd::calc_force_invgrads()
 {
   atoms.read_system_forces();
   ft.real_value = 0.0;
 
   // Note: gradient square norm is 1/N_atoms
 
   for (size_t ia = 0; ia < atoms.size(); ia++) {
     ft.real_value += atoms[ia].grad * atoms[ia].system_force;
   }
   ft.real_value *= atoms.size();
 }
 
 
 void colvar::rmsd::calc_Jacobian_derivative()
 {
   // divergence of the rotated coordinates (including only derivatives of the rotation matrix)
   cvm::real divergence = 0.0;
 
   if (atoms.b_rotate) {
 
     // gradient of the rotation matrix
     cvm::matrix2d <cvm::rvector, 3, 3> grad_rot_mat;
     // gradients of products of 2 quaternion components
     cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23;
     for (size_t ia = 0; ia < atoms.size(); ia++) {
 
       // Gradient of optimal quaternion wrt current Cartesian position
       cvm::vector1d< cvm::rvector, 4 >      &dq = atoms.rot.dQ0_1[ia];
 
       g11 = 2.0 * (atoms.rot.q)[1]*dq[1];
       g22 = 2.0 * (atoms.rot.q)[2]*dq[2];
       g33 = 2.0 * (atoms.rot.q)[3]*dq[3];
       g01 = (atoms.rot.q)[0]*dq[1] + (atoms.rot.q)[1]*dq[0];
       g02 = (atoms.rot.q)[0]*dq[2] + (atoms.rot.q)[2]*dq[0];
       g03 = (atoms.rot.q)[0]*dq[3] + (atoms.rot.q)[3]*dq[0];
       g12 = (atoms.rot.q)[1]*dq[2] + (atoms.rot.q)[2]*dq[1];
       g13 = (atoms.rot.q)[1]*dq[3] + (atoms.rot.q)[3]*dq[1];
       g23 = (atoms.rot.q)[2]*dq[3] + (atoms.rot.q)[3]*dq[2];
 
       // Gradient of the rotation matrix wrt current Cartesian position
       grad_rot_mat[0][0] = -2.0 * (g22 + g33);
       grad_rot_mat[1][0] =  2.0 * (g12 + g03);
       grad_rot_mat[2][0] =  2.0 * (g13 - g02);
       grad_rot_mat[0][1] =  2.0 * (g12 - g03);
       grad_rot_mat[1][1] = -2.0 * (g11 + g33);
       grad_rot_mat[2][1] =  2.0 * (g01 + g23);
       grad_rot_mat[0][2] =  2.0 * (g02 + g13);
       grad_rot_mat[1][2] =  2.0 * (g23 - g01);
       grad_rot_mat[2][2] = -2.0 * (g11 + g22);
 
       cvm::atom_pos &y = ref_pos[ia];
 
       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];)
         }
       }
     }
   }
 
   jd.real_value = x.real_value > 0.0 ? (3.0 * atoms.size() - 4.0 - divergence) / x.real_value : 0.0;
 }
 
 
 
 
 colvar::eigenvector::eigenvector (std::string const &conf)
   : cvc (conf)
 {
   b_inverse_gradients = true;
   b_Jacobian_derivative = true;
   function_type = "eigenvector";
   x.type (colvarvalue::type_scalar);
 
   parse_group (conf, "atoms", atoms);
   atom_groups.push_back (&atoms);
 
 
   {
     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 (""))) {
         // use PDB flags if column is provided
         bool found = get_keyval (conf, "refPositionsColValue", file_col_value, 0.0);
         if (found && !file_col_value)
           cvm::fatal_error ("Error: refPositionsColValue, "
                             "if provided, must be non-zero.\n");
       } else {
         // if not, use atom indices
         atoms.create_sorted_ids();
       }
 
       ref_pos.resize (atoms.size());
       cvm::load_coords (file_name.c_str(), ref_pos, atoms.sorted_ids, file_col, file_col_value);
     }
   }
 
   // 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();
 
   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\", to minimize RMSD before calculating the vector projection: "
               "if this is not the desired behavior, disable them explicitly within the \"atoms\" block.\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)) {
 
       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;
       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);
     }
   }
 
   if (!ref_pos.size() || !eigenvec.size()) {
     cvm::fatal_error ("Error: both reference coordinates"
                       "and eigenvector must be defined.\n");
   }
 
-  cvm::rvector eig_center (0.0, 0.0, 0.0);
+  cvm::atom_pos eig_center (0.0, 0.0, 0.0);
   for (size_t i = 0; i < atoms.size(); i++) {
     eig_center += eigenvec[i];
   }
   eig_center *= 1.0 / atoms.size();
   cvm::log ("Geometric center of the provided vector: "+cvm::to_str (eig_center)+"\n");
 
   bool b_difference_vector = false;
   get_keyval (conf, "differenceVector", b_difference_vector, false);
 
   if (b_difference_vector) {
 
     if (atoms.b_center) {
       // 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) {
       atoms.rot.calc_optimal_rotation (eigenvec, ref_pos);
       for (size_t i = 0; i < atoms.size(); i++) {
         eigenvec[i] = atoms.rot.rotate (eigenvec[i]);
       }
     }
     cvm::log ("\"differenceVector\" is on: subtracting the reference positions from the provided vector: v = v - x0.\n");
     for (size_t i = 0; i < atoms.size(); 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_invnorm2 += eigenvec[i].norm2();
   }
   eigenvec_invnorm2 = 1.0 / eigenvec_invnorm2;
 
   if (b_difference_vector) {
     cvm::log ("\"differenceVector\" is on: normalizing the vector.\n");
     for (size_t i = 0; i < atoms.size(); i++) {
       eigenvec[i] *= eigenvec_invnorm2;
     }
   } else {
     cvm::log ("The norm of the vector is |v| = "+cvm::to_str (eigenvec_invnorm2)+".\n");
   }
 }
 
 
 void colvar::eigenvector::calc_value()
 {
   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];
   }
 }
 
 
 void colvar::eigenvector::calc_gradients()
 {
   for (size_t ia = 0; ia < atoms.size(); ia++) {
     atoms[ia].grad = eigenvec[ia];
   }
 
   if (b_debug_gradients) {
     cvm::log ("Debugging gradients:\n");
     debug_gradients (atoms);
   }
 }
 
 
 void colvar::eigenvector::apply_force (colvarvalue const &force)
 {
   if (!atoms.noforce)
     atoms.apply_colvar_force (force.real_value);
 }
 
 
 void colvar::eigenvector::calc_force_invgrads()
 {
   atoms.read_system_forces();
   ft.real_value = 0.0;
 
   for (size_t ia = 0; ia < atoms.size(); ia++) {
     ft.real_value += eigenvec_invnorm2 * atoms[ia].grad *
       atoms[ia].system_force;
   }
 }
 
 
 void colvar::eigenvector::calc_Jacobian_derivative()
 {
   // gradient of the rotation matrix
   cvm::matrix2d <cvm::rvector, 3, 3> grad_rot_mat;
   cvm::quaternion &quat0 = atoms.rot.q;
 
   // gradients of products of 2 quaternion components
   cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23;
 
   cvm::atom_pos x_relative;
   cvm::real sum = 0.0;
 
   for (size_t ia = 0; ia < atoms.size(); ia++) {
 
     // Gradient of optimal quaternion wrt current Cartesian position
     // trick: d(R^-1)/dx = d(R^t)/dx = (dR/dx)^t
     // we can just transpose the derivatives of the direct matrix
     cvm::vector1d< cvm::rvector, 4 >      &dq_1 = atoms.rot.dQ0_1[ia];
 
     g11 = 2.0 * quat0[1]*dq_1[1];
     g22 = 2.0 * quat0[2]*dq_1[2];
     g33 = 2.0 * quat0[3]*dq_1[3];
     g01 = quat0[0]*dq_1[1] + quat0[1]*dq_1[0];
     g02 = quat0[0]*dq_1[2] + quat0[2]*dq_1[0];
     g03 = quat0[0]*dq_1[3] + quat0[3]*dq_1[0];
     g12 = quat0[1]*dq_1[2] + quat0[2]*dq_1[1];
     g13 = quat0[1]*dq_1[3] + quat0[3]*dq_1[1];
     g23 = quat0[2]*dq_1[3] + quat0[3]*dq_1[2];
 
     // Gradient of the inverse rotation matrix wrt current Cartesian position
     // (transpose of the gradient of the direct rotation)
     grad_rot_mat[0][0] = -2.0 * (g22 + g33);
     grad_rot_mat[0][1] =  2.0 * (g12 + g03);
     grad_rot_mat[0][2] =  2.0 * (g13 - g02);
     grad_rot_mat[1][0] =  2.0 * (g12 - g03);
     grad_rot_mat[1][1] = -2.0 * (g11 + g33);
     grad_rot_mat[1][2] =  2.0 * (g01 + g23);
     grad_rot_mat[2][0] =  2.0 * (g02 + g13);
     grad_rot_mat[2][1] =  2.0 * (g23 - g01);
     grad_rot_mat[2][2] = -2.0 * (g11 + g22);
 
     for (size_t i = 0; i < 3; i++) {
       for (size_t j = 0; j < 3; j++) {
         sum += grad_rot_mat[i][j][i] * eigenvec[ia][j];
       }
     }
   }
 
   jd.real_value = sum * std::sqrt (eigenvec_invnorm2);
 }