initial import of colvars code from svn archive

lib/colvars/colvar.h

@@ -0,0 +1,561 @@
+// -*- c++ -*-
+
+#ifndef COLVAR_H
+#define COLVAR_H
+
+#include <iostream>
+#include <iomanip>
+#include <cmath>
+
+#include "colvarmodule.h"
+#include "colvarvalue.h"
+#include "colvarparse.h"
+
+
+/// \brief A collective variable (main class); to be defined, it needs
+/// at least one object of a derived class of colvar::cvc; it
+/// calculates and returns a \link colvarvalue \endlink object
+///
+/// This class parses the configuration, defines the behaviour and
+/// stores the value (\link colvar::x \endlink) and all related data
+/// of a collective variable.  How the value is calculated is defined
+/// in \link colvar::cvc \endlink and its derived classes.  The
+/// \link colvar \endlink object contains pointers to multiple \link
+/// colvar::cvc \endlink derived objects, which can be combined
+/// together into one collective variable.  This makes possible to
+/// implement new collective variables at runtime based on the
+/// existing ones.  Currently, this possibility is limited to a
+/// polynomial, using the coefficients cvc::sup_coeff and the
+/// exponents cvc::sup_np.  In case of non-scalar variables,
+/// only exponents equal to 1 are accepted.
+///
+/// Please note that most of its members are \link colvarvalue
+/// \endlink objects, i.e. they can handle different data types
+/// together, and must all be set to the same type of colvar::x by
+/// using the colvarvalue::type() member function before using them
+/// together in assignments or other operations; this is usually done
+/// automatically in the constructor.  If you add a new member of
+/// \link colvarvalue \endlink type, you should also add its
+/// initialization line in the \link colvar \endlink constructor.
+
+class colvar : public colvarparse {
+
+public:
+
+  /// Name
+  std::string name;
+
+  /// Type of value
+  colvarvalue::Type type() const;
+
+  /// \brief Current value (previously obtained from calc() or read_traj())
+  colvarvalue const & value() const;
+
+  /// \brief Current velocity (previously obtained from calc() or read_traj())
+  colvarvalue const & velocity() const;
+
+  /// \brief Current system force (previously obtained from calc() or
+  /// read_traj()).  Note: this is calculated using the atomic forces
+  /// from the last simulation step.
+  ///
+  /// Total atom forces are read from the MD program, the total force
+  /// acting on the collective variable is calculated summing those
+  /// from all colvar components, the bias and walls forces are
+  /// subtracted.
+  colvarvalue const & system_force() const;
+
+  /// \brief
+
+  /// \brief Typical fluctuation amplitude for this collective
+  /// variable (e.g. local width of a free energy basin)
+  ///
+  /// In metadynamics calculations, \link colvarbias_meta \endlink,
+  /// this value is used to calculate the width of a gaussian.  In ABF
+  /// calculations, \link colvarbias_abf \endlink, it is used to
+  /// calculate the grid spacing in the direction of this collective
+  /// variable.
+  cvm::real width;
+
+  /// \brief True if this \link colvar \endlink is a linear
+  /// combination of \link cvc \endlink elements
+  bool b_linear;
+
+  /// \brief True if all \link cvc \endlink objects are capable
+  /// of calculating inverse gradients
+  bool b_inverse_gradients;
+
+  /// \brief True if all \link cvc \endlink objects are capable
+  /// of calculating Jacobian forces
+  bool b_Jacobian_force;
+
+  /// \brief Options controlling the behaviour of the colvar during
+  /// the simulation, which are set from outside the cvcs
+  enum task {
+    /// \brief Gradients are calculated and temporarily stored, so
+    /// that external forces can be applied
+    task_gradients,
+    /// \brief Collect atomic gradient data from all cvcs into vector
+    /// atomic_gradients
+    task_collect_gradients,
+    /// \brief Calculate the velocity with finite differences
+    task_fdiff_velocity,
+    /// \brief The system force is calculated, projecting the atomic
+    /// forces on the inverse gradients
+    task_system_force,
+    /// \brief The variable has a harmonic restraint around a moving
+    /// center with fictitious mass; bias forces will be applied to
+    /// the center
+    task_extended_lagrangian,
+    /// \brief The extended system coordinate undergoes Langevin
+    /// dynamics
+    task_langevin,
+    /// \brief Output the potential and kinetic energies
+    /// (for extended Lagrangian colvars only)
+    task_output_energy,
+    /// \brief Compute analytically the "force" arising from the
+    /// largest entropy component (for example, that from the angular
+    /// states orthogonal to a distance vector)
+    task_Jacobian_force,
+    /// \brief Report the Jacobian force as part of the system force
+    /// if disabled, apply a correction internally to cancel it
+    task_report_Jacobian_force,
+    /// \brief Output the value to the trajectory file (on by default)
+    task_output_value,
+    /// \brief Output the velocity to the trajectory file
+    task_output_velocity,
+    /// \brief Output the applied force to the trajectory file
+    task_output_applied_force,
+    /// \brief Output the system force to the trajectory file
+    task_output_system_force,
+    /// \brief A lower boundary is defined
+    task_lower_boundary,
+    /// \brief An upper boundary is defined
+    task_upper_boundary,
+    /// \brief Provide a discretization of the values of the colvar to
+    /// be used by the biases or in analysis (needs lower and upper
+    /// boundary)
+    task_grid,
+    /// \brief Apply a restraining potential (|x-xb|^2) to the colvar
+    /// when it goes below the lower wall
+    task_lower_wall,
+    /// \brief Apply a restraining potential (|x-xb|^2) to the colvar
+    /// when it goes above the upper wall
+    task_upper_wall,
+    /// \brief Compute running average
+    task_runave,
+    /// \brief Compute time correlation function
+    task_corrfunc,
+    /// \brief Number of possible tasks
+    task_ntot
+  };
+
+  /// Tasks performed by this colvar
+  bool tasks[task_ntot];
+
+protected:
+
+
+  /*
+    extended:
+    H = H_{0} + \sum_{i} 1/2*\lambda*(S_i(x(t))-s_i(t))^2 \\
+    + \sum_{i} 1/2*m_i*(ds_i(t)/dt)^2 \\
+    + \sum_{t'<t} W * exp (-1/2*\sum_{i} (s_i(t')-s_i(t))^2/(\delta{}s_i)^2) \\
+    + \sum_{w} (\sum_{i}c_{w,i}s_i(t) - D_w)^M/(\Sigma_w)^M
+
+    normal:
+    H = H_{0} + \sum_{t'<t} W * exp (-1/2*\sum_{i} (S_i(x(t'))-S_i(x(t)))^2/(\delta{}S_i)^2) \\
+    + \sum_{w} (\sum_{i}c_{w,i}S_i(t) - D_w)^M/(\Sigma_w)^M
+
+    output:
+    H = H_{0}   (only output S(x), no forces)
+
+    Here:
+    S(x(t)) = x
+    s(t)    = xr
+    DS = Ds = delta
+  */
+
+
+  /// Value of the colvar
+  colvarvalue x;
+
+  /// Cached reported value (x may be manipulated)
+  colvarvalue x_reported;
+
+  /// Finite-difference velocity
+  colvarvalue v_fdiff;
+
+  inline colvarvalue fdiff_velocity (colvarvalue const &xold,
+                                     colvarvalue const &xnew)
+  {
+    // using the gradient of the square distance to calculate the
+    // velocity (non-scalar variables automatically taken into
+    // account)
+    cvm::real dt = cvm::dt();
+    return ( ( (dt > 0.0) ? (1.0/dt) : 1.0 ) *
+             0.5 * dist2_lgrad (xnew, xold) );
+  }
+
+  /// Cached reported velocity
+  colvarvalue v_reported;
+
+  // Options for task_extended_lagrangian
+  /// Restraint center
+  colvarvalue xr;
+  /// Velocity of the restraint center
+  colvarvalue vr;
+  /// Mass of the restraint center
+  cvm::real ext_mass;
+  /// Restraint force constant
+  cvm::real ext_force_k;
+  /// Friction coefficient for Langevin extended dynamics
+  cvm::real ext_gamma;
+  /// Amplitude of Gaussian white noise for Langevin extended dynamics
+  cvm::real ext_sigma;
+  
+  /// \brief Harmonic restraint force
+  colvarvalue fr;
+
+  /// \brief Jacobian force, when task_Jacobian_force is enabled
+  colvarvalue fj;
+
+  /// Cached reported system force
+  colvarvalue ft_reported;
+
+public:
+
+
+  /// \brief Bias force; reset_bias_force() should be called before
+  /// the biases are updated
+  colvarvalue fb;
+
+  /// \brief Total \em applied force; fr (if task_extended_lagrangian
+  /// is defined), fb (if biases are applied) and the walls' forces
+  /// (if defined) contribute to it
+  colvarvalue f;
+
+  /// \brief Total force, as derived from the atomic trajectory;
+  /// should equal the total system force plus \link f \endlink
+  colvarvalue ft;
+
+
+  /// Period, if it is a constant
+  cvm::real period;
+
+  /// \brief Same as above, but also takes into account components
+  /// with a variable period, such as distanceZ
+  bool b_periodic;
+
+
+  /// \brief Expand the boundaries of multiples of width, to keep the
+  /// value always within range
+  bool expand_boundaries;
+
+  /// \brief Location of the lower boundary
+  colvarvalue lower_boundary;
+  /// \brief Location of the lower wall
+  colvarvalue lower_wall;
+  /// \brief Force constant for the lower boundary potential (|x-xb|^2)
+  cvm::real   lower_wall_k;
+
+  /// \brief Location of the upper boundary
+  colvarvalue upper_boundary;
+  /// \brief Location of the upper wall
+  colvarvalue upper_wall;
+  /// \brief Force constant for the upper boundary potential (|x-xb|^2)
+  cvm::real   upper_wall_k;
+
+  /// \brief Is the interval defined by the two boundaries periodic?
+  bool periodic_boundaries() const;
+
+  /// \brief Is the interval defined by the two boundaries periodic?
+  bool periodic_boundaries (colvarvalue const &lb, colvarvalue const &ub) const;
+
+
+  /// Constructor
+  colvar (std::string const &conf);
+
+  /// Enable the specified task
+  void enable (colvar::task const &t);
+
+  /// Disable the specified task
+  void disable (colvar::task const &t);
+
+  /// Destructor
+  ~colvar();
+
+
+  /// \brief Calculate the colvar value and all the other requested
+  /// quantities
+  void calc();
+
+  /// \brief Calculate just the value, and store it in the argument
+  void calc_value (colvarvalue &xn);
+
+  /// Set the total biasing force to zero
+  void reset_bias_force();
+
+  /// Add to the total force from biases
+  void add_bias_force (colvarvalue const &force);
+
+  /// \brief Collect all forces on this colvar, integrate internal
+  /// equations of motion of internal degrees of freedom; see also
+  /// colvar::communicate_forces()
+  /// return colvar energy if extended Lagrandian active
+  cvm::real update();
+
+  /// \brief Communicate forces (previously calculated in
+  /// colvar::update()) to the external degrees of freedom
+  void communicate_forces();
+
+
+  /// \brief Use the internal metrics (as from \link cvc
+  /// \endlink objects) to calculate square distances and gradients
+  ///
+  /// Handles correctly symmetries and periodic boundary conditions
+  cvm::real dist2 (colvarvalue const &x1,
+                   colvarvalue const &x2) const;
+
+  /// \brief Use the internal metrics (as from \link cvc
+  /// \endlink objects) to calculate square distances and gradients
+  ///
+  /// Handles correctly symmetries and periodic boundary conditions
+  colvarvalue dist2_lgrad (colvarvalue const &x1,
+                           colvarvalue const &x2) const;
+
+  /// \brief Use the internal metrics (as from \link cvc
+  /// \endlink objects) to calculate square distances and gradients
+  ///
+  /// Handles correctly symmetries and periodic boundary conditions
+  colvarvalue dist2_rgrad (colvarvalue const &x1,
+                           colvarvalue const &x2) const;
+
+  /// \brief Use the internal metrics (as from \link cvc
+  /// \endlink objects) to compare colvar values
+  ///
+  /// Handles correctly symmetries and periodic boundary conditions
+  cvm::real compare (colvarvalue const &x1,
+                     colvarvalue const &x2) const;
+
+  /// \brief Use the internal metrics (as from \link cvc
+  /// \endlink objects) to wrap a value into a standard interval
+  ///
+  /// Handles correctly symmetries and periodic boundary conditions
+  void wrap (colvarvalue &x) const;
+
+
+  /// Read the analysis tasks
+  void parse_analysis (std::string const &conf);
+  /// Perform analysis tasks
+  void analyse();
+
+
+  /// Read the value from a collective variable trajectory file
+  std::istream & read_traj (std::istream &is);
+  /// Output formatted values to the trajectory file
+  std::ostream & write_traj (std::ostream &os);
+  /// Write a label to the trajectory file (comment line)
+  std::ostream & write_traj_label (std::ostream &os);
+
+
+  /// Read the collective variable from a restart file
+  std::istream & read_restart (std::istream &is);
+  /// Write the collective variable to a restart file
+  std::ostream & write_restart (std::ostream &os);
+
+
+protected:
+
+  /// Previous value (to calculate velocities during analysis)
+  colvarvalue            x_old;
+
+  /// Time series of values and velocities used in correlation
+  /// functions 
+  std::list< std::list<colvarvalue> > acf_x_history, acf_v_history;
+  /// Time series of values and velocities used in correlation
+  /// functions (pointers)x
+  std::list< std::list<colvarvalue> >::iterator acf_x_history_p, acf_v_history_p;
+
+  /// Time series of values and velocities used in running averages
+  std::list< std::list<colvarvalue> > x_history;
+  /// Time series of values and velocities used in correlation
+  /// functions (pointers)x
+  std::list< std::list<colvarvalue> >::iterator x_history_p;
+
+  /// \brief Collective variable with which the correlation is
+  /// calculated (default: itself)
+  std::string            acf_colvar_name;
+  /// Length of autocorrelation function (ACF)
+  size_t                 acf_length;
+  /// After how many steps the ACF starts
+  size_t                 acf_offset;
+  /// How many timesteps separate two ACF values
+  size_t                 acf_stride;
+  /// Number of frames for each ACF point
+  size_t                 acf_nframes;
+  /// Normalize the ACF to a maximum value of 1?
+  bool                   acf_normalize;
+  /// ACF values
+  std::vector<cvm::real> acf;
+  /// Name of the file to write the ACF
+  std::string            acf_outfile;
+
+  /// Type of autocorrelation function (ACF)
+  enum acf_type_e {
+    /// Unset type
+    acf_notset,
+    /// Velocity ACF, scalar product between v(0) and v(t)
+    acf_vel,
+    /// Coordinate ACF, scalar product between x(0) and x(t)
+    acf_coor,
+    /// \brief Coordinate ACF, second order Legendre polynomial
+    /// between x(0) and x(t) (does not work with scalar numbers)
+    acf_p2coor
+  };
+
+  /// Type of autocorrelation function (ACF)
+  acf_type_e             acf_type;
+
+  /// \brief Velocity ACF, scalar product between v(0) and v(t)
+  void calc_vel_acf (std::list<colvarvalue> &v_history,
+                     colvarvalue const      &v);
+
+  /// \brief Coordinate ACF, scalar product between x(0) and x(t)
+  /// (does not work with scalar numbers)
+  void calc_coor_acf (std::list<colvarvalue> &x_history,
+                      colvarvalue const      &x);
+
+  /// \brief Coordinate ACF, second order Legendre polynomial between
+  /// x(0) and x(t) (does not work with scalar numbers)
+  void calc_p2coor_acf (std::list<colvarvalue> &x_history,
+                        colvarvalue const      &x);
+
+  /// Calculate the auto-correlation function (ACF)
+  void calc_acf();
+  /// Save the ACF to a file
+  void write_acf (std::ostream &os);
+
+  /// Length of running average series
+  size_t         runave_length;
+  /// Timesteps to skip between two values in the running average series
+  size_t         runave_stride;
+  /// Name of the file to write the running average
+  std::ofstream  runave_os;
+  /// Current value of the running average
+  colvarvalue    runave;
+  /// Current value of the square deviation from the running average
+  cvm::real      runave_variance;
+
+  /// Calculate the running average and its standard deviation
+  void calc_runave();
+
+  /// If extended Lagrangian active: colvar energies (kinetic and harmonic potential)
+  cvm::real kinetic_energy;
+  cvm::real potential_energy;
+public:
+
+
+  // collective variable component base class
+  class cvc;
+
+  // currently available collective variable components
+
+  // scalar colvar components
+  class distance;
+  class distance_z;
+  class distance_xy;
+  class min_distance;
+  class angle;
+  class dihedral;
+  class coordnum;
+  class selfcoordnum;
+  class h_bond;
+  class rmsd;
+  class logmsd;
+  class orientation_angle;
+  class tilt;
+  class spin_angle;
+  class gyration;
+  class eigenvector;
+  class alpha_dihedrals;
+  class alpha_angles;
+  class dihedPC;
+
+  // non-scalar components
+  class distance_vec;
+  class distance_dir;
+  class orientation;
+
+protected:
+
+  /// \brief Array of \link cvc \endlink objects
+  std::vector<cvc *> cvcs;
+
+  /// \brief Initialize the sorted list of atom IDs for atoms involved
+  /// in all cvcs (called when enabling task_collect_gradients)
+  void build_atom_list (void);
+
+public:
+  /// \brief Sorted array of (zero-based) IDs for all atoms involved
+  std::vector<int> atom_ids;
+
+  /// \brief Array of atomic gradients collected from all cvcs
+  /// with appropriate components, rotations etc.
+  /// For scalar variables only!
+  std::vector<cvm::rvector> atomic_gradients;
+
+  inline size_t n_components () const {
+    return cvcs.size();
+  }
+};
+
+
+inline colvar * cvm::colvar_p (std::string const &name)
+{
+  for (std::vector<colvar *>::iterator cvi = cvm::colvars.begin();
+       cvi != cvm::colvars.end();
+       cvi++) {
+    if ((*cvi)->name == name) {
+      return (*cvi);
+    }
+  }
+  return NULL;
+}
+
+
+inline colvarvalue::Type colvar::type() const
+{
+  return x.type();
+}
+
+
+inline colvarvalue const & colvar::value() const
+{
+  return x_reported;
+}
+
+
+inline colvarvalue const & colvar::velocity() const
+{
+  return v_reported;
+}
+
+
+inline colvarvalue const & colvar::system_force() const
+{
+  return ft_reported;
+}
+
+
+inline void colvar::add_bias_force (colvarvalue const &force)
+{
+  fb += force;
+}
+
+
+inline void colvar::reset_bias_force() {
+  fb.reset();
+}
+
+#endif
+

lib/colvars/colvaratoms.C

@@ -0,0 +1,761 @@
+#include "colvarmodule.h"
+#include "colvarparse.h"
+#include "colvaratoms.h"
+
+
+// atom member functions depend tightly on the MD interface, and are
+// thus defined in colvarproxy_xxx.cpp
+
+
+// atom_group member functions
+
+cvm::atom_group::atom_group (std::string const &conf,
+                             char const        *key,
+                             atom_group        *ref_pos_group_in)
+  : b_center (false), b_rotate (false),
+    ref_pos_group (NULL), // this is always set within parse(),
+                          // regardless of ref_pos_group_in
+    noforce (false)
+{
+  cvm::log ("Defining atom group \""+
+            std::string (key)+"\".\n");
+  parse (conf, key, ref_pos_group_in);
+}
+
+
+void cvm::atom_group::parse (std::string const &conf,
+                             char const        *key,
+                             atom_group        *ref_pos_group_in)
+{
+  std::string group_conf;
+
+  // save_delimiters is set to false for this call, because "conf" is
+  // not the config string of this group, but of its parent object
+  // (which has already taken care of the delimiters)
+  save_delimiters = false;
+  key_lookup (conf, key, group_conf, dummy_pos);
+  // restoring the normal value, because we do want keywords checked
+  // inside "group_conf"
+  save_delimiters = true;
+
+  if (group_conf.size() == 0) {
+    cvm::fatal_error ("Error: atom group \""+
+                      std::string (key)+"\" is set, but "
+                      "has no definition.\n");
+  }
+
+  cvm::increase_depth();
+
+  cvm::log ("Initializing atom group \""+std::string (key)+"\".\n");
+
+  // whether or not to include messages in the log
+  colvarparse::Parse_Mode mode = parse_silent;
+  {
+    bool b_verbose;
+    get_keyval (group_conf, "verboseOutput", b_verbose, false, parse_silent);
+    if (b_verbose) mode = parse_normal;
+  }
+
+  {
+    // get the atoms by numbers
+    std::vector<int> atom_indexes;
+    if (get_keyval (group_conf, "atomNumbers", atom_indexes, atom_indexes, mode)) {
+      if (atom_indexes.size()) {
+        this->reserve (this->size()+atom_indexes.size());
+        for (size_t i = 0; i < atom_indexes.size(); i++) {
+          this->push_back (cvm::atom (atom_indexes[i]));
+        }
+      } else
+        cvm::fatal_error ("Error: no numbers provided for \""
+                          "atomNumbers\".\n");
+    }
+  }
+
+  {
+    std::string range_conf = "";
+    size_t pos = 0;
+    while (key_lookup (group_conf, "atomNumbersRange",
+                       range_conf, pos)) {
+
+      if (range_conf.size()) {
+        std::istringstream is (range_conf);
+        int initial, final;
+        char dash;
+        if ( (is >> initial) && (initial > 0) &&
+             (is >> dash) && (dash == '-') &&
+             (is >> final) && (final > 0) ) {
+          for (int anum = initial; anum <= final; anum++) {
+            this->push_back (cvm::atom (anum));
+          }
+          range_conf = "";
+          continue;
+        }
+      }
+
+      cvm::fatal_error ("Error: no valid definition for \""
+                        "atomNumbersRange\", \""+
+                        range_conf+"\".\n");
+    }
+  }
+
+  std::vector<std::string> psf_segids;
+  get_keyval (group_conf, "psfSegID", psf_segids, std::vector<std::string> (), mode);
+  for (std::vector<std::string>::iterator psii = psf_segids.begin();
+       psii < psf_segids.end(); psii++) {
+
+    if ( (psii->size() == 0) || (psii->size() > 4) ) {
+      cvm::fatal_error ("Error: invalid segmend identifier provided, \""+
+                        (*psii)+"\".\n");
+    }
+  }
+
+  {
+    std::string range_conf = "";
+    size_t pos = 0;
+    size_t range_count = 0;
+    std::vector<std::string>::iterator psii = psf_segids.begin();
+    while (key_lookup (group_conf, "atomNameResidueRange",
+                       range_conf, pos)) {
+
+      if (psii > psf_segids.end()) {
+        cvm::fatal_error ("Error: more instances of \"atomNameResidueRange\" than "
+                          "values of \"psfSegID\".\n");
+      }
+
+      std::string const &psf_segid = psf_segids.size() ? *psii : std::string ("");
+
+      if (range_conf.size()) {
+        
+        std::istringstream is (range_conf);
+        std::string atom_name;
+        int initial, final;
+        char dash;
+        if ( (is >> atom_name) && (atom_name.size())  &&
+             (is >> initial) && (initial > 0) &&
+             (is >> dash) && (dash == '-') &&
+             (is >> final) && (final > 0) ) {
+          for (int resid = initial; resid <= final; resid++) {
+            this->push_back (cvm::atom (resid, atom_name, psf_segid));
+          }
+          range_conf = "";
+        } else {
+          cvm::fatal_error ("Error: cannot parse definition for \""
+                            "atomNameResidueRange\", \""+
+                            range_conf+"\".\n");
+        }
+
+      } else {
+        cvm::fatal_error ("Error: atomNameResidueRange with empty definition.\n");
+      }
+
+      if (psf_segid.size())
+        psii++;
+    }
+  }
+
+  {
+    // read the atoms from a file
+    std::string atoms_file_name;
+    if (get_keyval (group_conf, "atomsFile", atoms_file_name, std::string (""), mode)) {
+
+      std::string atoms_col;
+      if (!get_keyval (group_conf, "atomsCol", atoms_col, std::string (""), mode)) {
+        cvm::fatal_error ("Error: parameter atomsCol is required if atomsFile is set.\n");
+      }
+
+      double atoms_col_value;
+      bool const atoms_col_value_defined = get_keyval (group_conf, "atomsColValue", atoms_col_value, 0.0, mode);
+      if (atoms_col_value_defined && (!atoms_col_value))
+        cvm::fatal_error ("Error: atomsColValue, "
+                          "if provided, must be non-zero.\n");
+
+      cvm::load_atoms (atoms_file_name.c_str(), *this, atoms_col, atoms_col_value);
+    }
+  }
+
+  for (std::vector<cvm::atom>::iterator a1 = this->begin(); 
+       a1 != this->end(); a1++) {
+    std::vector<cvm::atom>::iterator a2 = a1;
+    ++a2;
+    for ( ; a2 != this->end(); a2++) {
+      if (a1->id == a2->id) {
+        if (cvm::debug())
+          cvm::log ("Discarding doubly counted atom with number "+
+                    cvm::to_str (a1->id+1)+".\n");
+        a2 = this->erase (a2);
+        if (a2 == this->end())
+          break;
+      }
+    }
+  }
+
+  if (get_keyval (group_conf, "dummyAtom", dummy_atom_pos, cvm::atom_pos(), mode)) {
+    b_dummy = true;
+    this->total_mass = 1.0;
+  } else 
+    b_dummy = false;
+
+  if (b_dummy && (this->size())) 
+    cvm::fatal_error ("Error: cannot set up group \""+
+                      std::string (key)+"\" as a dummy atom "
+                      "and provide it with atom definitions.\n");
+
+#if (! defined (COLVARS_STANDALONE))
+  if ( (!b_dummy) && (!cvm::b_analysis) && (!(this->size())) ) {
+    cvm::fatal_error ("Error: no atoms defined for atom group \""+
+                      std::string (key)+"\".\n");
+  }
+#endif
+
+  if (!b_dummy) {
+    this->total_mass = 0.0;
+    for (cvm::atom_iter ai = this->begin();
+         ai != this->end(); ai++) {
+      this->total_mass += ai->mass;
+    }
+  }
+
+  get_keyval (group_conf, "disableForces",   noforce,  false, mode);
+
+  get_keyval (group_conf, "centerReference", b_center, false, mode);
+  get_keyval (group_conf, "rotateReference", b_rotate, false, mode);
+
+  if (b_center || b_rotate) {
+
+    if (b_dummy)
+      cvm::fatal_error ("Error: cannot set \"centerReference\" or "
+                        "\"rotateReference\" with \"dummyAtom\".\n");
+
+    // use refPositionsGroup instead of this group as the one which is
+    // used to fit the coordinates
+    if (key_lookup (group_conf, "refPositionsGroup")) {
+      if (ref_pos_group) {
+        cvm::fatal_error ("Error: the atom group \""+
+                          std::string (key)+"\" has already a reference group "
+                          "for the rototranslational fit, which was communicated by the "
+                          "colvar component.  You should not use refPositionsGroup "
+                          "in this case.\n");
+      }
+      cvm::log ("Within atom group \""+std::string (key)+"\":\n");
+      ref_pos_group = new atom_group (group_conf, "refPositionsGroup");
+    }
+
+    atom_group *ag = ref_pos_group ? ref_pos_group : this;
+
+    if (get_keyval (group_conf, "refPositions", ref_pos, ref_pos, mode)) {
+      cvm::log ("Using reference positions from input file.\n");
+      if (ref_pos.size() != ag->size()) {
+        cvm::fatal_error ("Error: the number of reference positions provided ("+
+                          cvm::to_str (ref_pos.size())+
+                          ") does not match the number of atoms within \""+
+                          std::string (key)+
+                          "\" ("+cvm::to_str (ag->size())+").\n");
+      }
+    }
+
+    std::string ref_pos_file;
+    if (get_keyval (group_conf, "refPositionsFile", ref_pos_file, std::string (""), mode)) {
+
+      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 (group_conf, "refPositionsCol", ref_pos_col, std::string (""), mode)) {
+        // if provided, use PDB column to select coordinates
+        bool found = get_keyval (group_conf, "refPositionsColValue", ref_pos_col_value, 0.0, mode);
+        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
+        ag->create_sorted_ids();
+      }
+      cvm::load_coords (ref_pos_file.c_str(), ref_pos, ag->sorted_ids,
+                        ref_pos_col, ref_pos_col_value);
+    }
+
+    if (ref_pos.size()) {
+
+      if (b_rotate) {
+        if (ref_pos.size() != ag->size())
+          cvm::fatal_error ("Error: the number of reference positions provided ("+
+                            cvm::to_str (ref_pos.size())+
+                            ") does not match the number of atoms within \""+
+                            std::string (key)+
+                            "\" ("+cvm::to_str (ag->size())+").\n");
+      }
+
+      // save the center of mass of ref_pos and then subtract it from
+      // them; in this way it is possible to use the coordinates for
+      // the rotational fit, if needed
+      ref_pos_cog = cvm::atom_pos (0.0, 0.0, 0.0);
+      std::vector<cvm::atom_pos>::iterator pi = ref_pos.begin();
+      for ( ; pi != ref_pos.end(); pi++) {
+        ref_pos_cog += *pi;
+      }
+      ref_pos_cog /= (cvm::real) ref_pos.size();
+
+      for (std::vector<cvm::atom_pos>::iterator pi = ref_pos.begin();
+           pi != ref_pos.end(); pi++) {
+        (*pi) -= ref_pos_cog;
+      }
+    } else {
+#if (! defined (COLVARS_STANDALONE))
+      if (!cvm::b_analysis)
+        cvm::fatal_error ("Error: no reference positions provided.\n");
+#endif
+    }
+
+    if (b_rotate && !noforce) {
+      cvm::log ("Warning: atom group \""+std::string (key)+
+                "\" is set to be rotated to a reference orientation: "
+                "a torque different than zero on this group "
+                "could make the simulation unstable.  "
+                "If this happens, set \"disableForces\" to yes "
+                "for this group.\n");
+    }
+
+  }
+
+
+  if (cvm::debug())
+    cvm::log ("Done initializing atom group with name \""+
+              std::string (key)+"\".\n");
+
+  this->check_keywords (group_conf, key);
+
+  cvm::log ("Atom group \""+std::string (key)+"\" defined, "+
+            cvm::to_str (this->size())+" initialized: total mass = "+
+            cvm::to_str (this->total_mass)+".\n");
+
+  cvm::decrease_depth();
+}
+
+
+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
+  if (sorted_ids.size())
+    return;
+
+  std::list<int> temp_id_list;
+  for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
+    temp_id_list.push_back (ai->id);
+  }
+  temp_id_list.sort();
+  temp_id_list.unique();
+  if (temp_id_list.size() != this->size()) {
+    cvm::fatal_error ("Error: duplicate atom IDs in atom group? (found " +
+                      cvm::to_str(temp_id_list.size()) +
+                      " unique atom IDs instead of" +
+                      cvm::to_str(this->size()) + ").\n");
+  }
+  sorted_ids = std::vector<int> (temp_id_list.begin(), temp_id_list.end());
+  return;
+}
+
+
+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();
+  }
+
+  if (ref_pos_group)
+    ref_pos_group->read_positions();
+
+  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
+    cvm::atom_pos const cog = fit_group->center_of_geometry();
+    for (cvm::atom_iter ai = this->begin();
+         ai != this->end(); ai++) {
+      ai->pos -= cog;
+    }
+  }
+
+  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
+    rot.calc_optimal_rotation (fit_group->positions(), ref_pos);
+
+    for (cvm::atom_iter ai = this->begin();
+         ai != this->end(); ai++) {
+      ai->pos = rot.rotate (ai->pos);
+    }
+  }
+
+  if (b_center) {
+    // use the center of geometry of ref_pos
+    for (cvm::atom_iter ai = this->begin();
+         ai != this->end(); ai++) {
+      ai->pos += ref_pos_cog;
+    }
+  }
+}
+
+
+void cvm::atom_group::apply_translation (cvm::rvector const &t) 
+{
+  if (b_dummy) return;
+
+  for (cvm::atom_iter ai = this->begin();
+       ai != this->end(); ai++) {
+    ai->pos += t;
+  }
+}
+
+
+void cvm::atom_group::apply_rotation (cvm::rotation const &rot) 
+{
+  if (b_dummy) return;
+
+  for (cvm::atom_iter ai = this->begin();
+       ai != this->end(); ai++) {
+    ai->pos = rot.rotate (ai->pos);
+  }
+}
+
+
+void cvm::atom_group::read_velocities()
+{
+  if (b_dummy) return;
+
+  if (b_rotate) {
+
+    for (cvm::atom_iter ai = this->begin();
+         ai != this->end(); ai++) {
+      ai->read_velocity();
+      ai->vel = rot.rotate (ai->vel);
+    }
+
+  } else {
+
+    for (cvm::atom_iter ai = this->begin();
+         ai != this->end(); ai++) {
+      ai->read_velocity();
+    }
+  }
+}
+
+void cvm::atom_group::read_system_forces()
+{
+  if (b_dummy) return;
+
+  if (b_rotate) {
+
+    for (cvm::atom_iter ai = this->begin();
+         ai != this->end(); ai++) {
+      ai->read_system_force();
+      ai->system_force = rot.rotate (ai->system_force);
+    }
+
+  } else {
+
+    for (cvm::atom_iter ai = this->begin();
+         ai != this->end(); ai++) {
+      ai->read_system_force();
+    }
+  }
+}
+
+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)
+    return dummy_atom_pos;
+
+  cvm::atom_pos cog (0.0, 0.0, 0.0);
+  for (cvm::atom_const_iter ai = this->begin();
+       ai != this->end(); ai++) {
+    cog += ai->pos;
+  }
+  cog /= this->size();
+  return cog;
+}
+
+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)
+    return dummy_atom_pos;
+
+  cvm::atom_pos com (0.0, 0.0, 0.0);
+  for (cvm::atom_const_iter ai = this->begin();
+       ai != this->end(); ai++) {
+    com += ai->mass * ai->pos;
+  }
+  com /= this->total_mass;
+  return com;
+}
+
+
+void cvm::atom_group::set_weighted_gradient (cvm::rvector const &grad)
+{
+  if (b_dummy) return;
+
+  for (cvm::atom_iter ai = this->begin();
+       ai != this->end(); ai++) {
+    ai->grad = (ai->mass/this->total_mass) * grad;
+  }
+}
+
+
+std::vector<cvm::atom_pos> cvm::atom_group::positions() const
+{
+  if (b_dummy)
+    cvm::fatal_error ("Error: positions are not available "
+                      "from a dummy atom group.\n");
+
+  std::vector<cvm::atom_pos> x (this->size(), 0.0);
+  cvm::atom_const_iter ai = this->begin();
+  std::vector<cvm::atom_pos>::iterator xi = x.begin();
+  for ( ; ai != this->end(); xi++, ai++) {
+    *xi = ai->pos;
+  }
+  return x;
+}
+
+std::vector<cvm::atom_pos> cvm::atom_group::positions_shifted (cvm::rvector const &shift) const
+{
+  if (b_dummy)
+    cvm::fatal_error ("Error: positions are not available "
+                      "from a dummy atom group.\n");
+
+  std::vector<cvm::atom_pos> x (this->size(), 0.0);
+  cvm::atom_const_iter ai = this->begin();
+  std::vector<cvm::atom_pos>::iterator xi = x.begin();
+  for ( ; ai != this->end(); xi++, ai++) {
+    *xi = (ai->pos + shift);
+  }
+  return x;
+}
+
+std::vector<cvm::rvector> cvm::atom_group::velocities() const
+{
+  if (b_dummy)
+    cvm::fatal_error ("Error: velocities are not available "
+                      "from a dummy atom group.\n");
+
+  std::vector<cvm::rvector> v (this->size(), 0.0);
+  cvm::atom_const_iter ai = this->begin();
+  std::vector<cvm::atom_pos>::iterator vi = v.begin();
+  for ( ; ai != this->end(); vi++, ai++) {
+    *vi = ai->vel;
+  }
+  return v;
+}
+
+std::vector<cvm::rvector> cvm::atom_group::system_forces() const
+{
+  if (b_dummy)
+    cvm::fatal_error ("Error: system forces are not available "
+                      "from a dummy atom group.\n");
+
+  std::vector<cvm::rvector> f (this->size(), 0.0);
+  cvm::atom_const_iter ai = this->begin();
+  std::vector<cvm::atom_pos>::iterator fi = f.begin();
+  for ( ; ai != this->end(); fi++, ai++) {
+    *fi = ai->system_force;
+  }
+  return f;
+}
+
+cvm::rvector cvm::atom_group::system_force() const
+{
+  if (b_dummy)
+    cvm::fatal_error ("Error: system forces are not available "
+                      "from a dummy atom group.\n");
+
+  cvm::rvector f (0.0);
+  for (cvm::atom_const_iter ai = this->begin(); ai != this->end(); ai++) {
+    f += ai->system_force;
+  }
+  return f;
+}
+
+
+void cvm::atom_group::apply_colvar_force (cvm::real const &force)
+{
+  if (b_dummy)
+    return;
+
+  if (noforce)
+    cvm::fatal_error ("Error: sending a force to a group that has "
+                      "\"disableForces\" defined.\n");
+
+  if (b_rotate) {
+
+    // get the forces back from the rotated frame
+    cvm::rotation const rot_inv = rot.inverse();
+    for (cvm::atom_iter ai = this->begin();
+         ai != this->end(); ai++) {
+      ai->apply_force (rot_inv.rotate (force * ai->grad));
+    }
+
+  } 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);
+    }
+  }
+}
+
+
+void cvm::atom_group::apply_force (cvm::rvector const &force)
+{
+  if (b_dummy)
+    return;
+
+  if (noforce)
+    cvm::fatal_error ("Error: sending a force to a group that has "
+                      "\"disableForces\" defined.\n");
+
+  if (b_rotate) {
+
+    cvm::rotation const rot_inv = rot.inverse();
+    for (cvm::atom_iter ai = this->begin();
+         ai != this->end(); ai++) {
+      ai->apply_force (rot_inv.rotate ((ai->mass/this->total_mass) * force));
+    }
+
+  } else {
+
+    for (cvm::atom_iter ai = this->begin();
+         ai != this->end(); ai++) {
+      ai->apply_force ((ai->mass/this->total_mass) * force);
+    } 
+  }
+}
+
+
+void cvm::atom_group::apply_forces (std::vector<cvm::rvector> const &forces)
+{
+  if (b_dummy)
+    return;
+
+  if (noforce)
+    cvm::fatal_error ("Error: sending a force to a group that has "
+                      "\"disableForces\" defined.\n");
+
+  if (forces.size() != this->size()) {
+    cvm::fatal_error ("Error: trying to apply an array of forces to an atom "
+                      "group which does not have the same length.\n");
+  }
+
+  if (b_rotate) {
+
+    cvm::rotation const rot_inv = rot.inverse();
+    cvm::atom_iter ai = this->begin();
+    std::vector<cvm::rvector>::const_iterator fi = forces.begin();
+    for ( ; ai != this->end(); fi++, ai++) {
+      ai->apply_force (rot_inv.rotate (*fi));
+    }
+
+  } else {
+
+    cvm::atom_iter ai = this->begin();
+    std::vector<cvm::rvector>::const_iterator fi = forces.begin();
+    for ( ; ai != this->end(); fi++, ai++) {
+      ai->apply_force (*fi);
+    }
+  }
+}
+
+

lib/colvars/colvaratoms.h

@@ -0,0 +1,318 @@
+#ifndef COLVARATOMS_H
+#define COLVARATOMS_H
+
+#include "colvarmodule.h"
+#include "colvarparse.h"
+
+/// \brief Stores numeric id, mass and all mutable data for an atom,
+/// mostly used by a \link cvc \endlink
+/// 
+/// This class may be used (although not necessarily) to keep atomic
+/// data (id, mass, position and collective variable derivatives)
+/// altogether.  There may be multiple instances with identical
+/// numeric id, all acting independently: forces communicated through
+/// these instances will be summed together.
+///
+/// Read/write operations depend on the underlying code: hence, some
+/// member functions are defined in colvarproxy_xxx.h.
+class colvarmodule::atom {
+
+protected:
+
+  /// \brief Index in the list of atoms involved by the colvars (\b
+  /// NOT in the global topology!)
+  size_t        index;
+
+public:
+
+  /// Internal identifier (zero-based)
+  int              id;
+
+  /// Mass
+  cvm::real      mass;
+
+  /// \brief Current position (copied from the program, can be
+  /// manipulated)
+  cvm::atom_pos   pos;
+
+  /// \brief Current velocity (copied from the program, can be
+  /// manipulated)
+  cvm::rvector    vel;
+
+  /// \brief System force at the previous step (copied from the
+  /// program, can be manipulated)
+  cvm::rvector    system_force;
+
+  /// \brief Gradient of a scalar collective variable with respect
+  /// to this atom
+  /// 
+  /// This can only handle a scalar collective variable (i.e. when
+  /// the \link colvarvalue::real_value \endlink member is used
+  /// from the \link colvarvalue \endlink class), which is also the
+  /// most frequent case. For more complex types of \link
+  /// colvarvalue \endlink objects, atomic gradients should be
+  /// defined within the specific \link cvc \endlink
+  /// implementation
+  cvm::rvector   grad;
+
+  /// \brief Default constructor, setting id to a non-valid value
+  inline atom() {}
+
+  /// \brief Initialize an atom for collective variable calculation
+  /// and get its internal identifier \param atom_number Atom index in
+  /// the system topology (starting from 1)
+  atom (int const &atom_number);
+
+  /// \brief Initialize an atom for collective variable calculation
+  /// and get its internal identifier \param residue Residue number
+  /// \param atom_name Name of the atom in the residue \param
+  /// segment_id For PSF topologies, the segment identifier; for other
+  /// type of topologies, may not be required
+  atom (cvm::residue_id const &residue,
+        std::string const     &atom_name,
+        std::string const     &segment_id = std::string (""));
+
+  /// Copy constructor
+  atom (atom const &a);
+
+  /// Destructor
+  ~atom();
+
+  /// Set non-constant data (everything except id and mass) to zero
+  inline void reset_data() {
+    pos = atom_pos (0.0);
+    vel = grad = system_force = rvector (0.0);
+  }
+
+  /// Get the current position
+  void read_position();
+
+  /// Get the current velocity
+  void read_velocity();
+
+  /// Get the system force
+  void read_system_force();
+
+  /// \brief Apply a force to the atom
+  ///
+  /// The force will be used later by the MD integrator, the
+  /// collective variables module does not integrate equations of
+  /// motion.  Multiple calls to this function by either the same
+  /// \link atom \endlink object or different objects with identical
+  /// \link id \endlink, will all add to the existing MD force.
+  void apply_force (cvm::rvector const &new_force);
+};
+
+
+
+
+/// \brief Group of \link atom \endlink objects, mostly used by a
+/// \link cvc \endlink
+///
+/// This class inherits from \link colvarparse \endlink and from
+/// std::vector<colvarmodule::atom>, and hence all functions and
+/// operators (including the bracket operator, group[i]) can be used
+/// on an \link atom_group \endlink object.  It can be initialized as
+/// a vector, or by parsing a keyword in the configuration.
+class colvarmodule::atom_group
+  : public std::vector<cvm::atom>,
+    public colvarparse
+{
+public:
+  // Note: all members here are kept public, to make possible to any
+  // object accessing and manipulating them
+
+
+  /// \brief If this option is on, this group merely acts as a wrapper
+  /// for a fixed position; any calls to atoms within or to
+  /// functions that return disaggregated data will fail
+  bool b_dummy;
+  /// \brief dummy atom position
+  cvm::atom_pos dummy_atom_pos;
+
+  /// Sorted list of zero-based (internal) atom ids
+  /// (populated on-demand by create_sorted_ids)
+  std::vector<int> sorted_ids;
+
+  /// Allocates and populates the sorted list of atom ids
+  void create_sorted_ids (void);
+  
+
+  /// \brief Before calculating colvars, move the group to overlap the
+  /// center of mass of 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.
+  ///
+  /// 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
+  cvm::rotation rot;
+
+  /// \brief In case b_center or b_rotate is true, use these reference
+  /// coordinates
+  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)
+  atom_group                *ref_pos_group;
+
+
+  /// Total mass of the atom group
+  cvm::real total_mass;
+
+  /// \brief Don't apply any force on this group (use its coordinates
+  /// only to calculate a colvar)
+  bool        noforce;
+
+
+  /// \brief Initialize the group by looking up its configuration
+  /// string in conf and parsing it; this is actually done by parse(),
+  /// 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);
+
+  /// \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);
+
+  /// \brief Initialize the group after a temporary vector of atoms
+  atom_group (std::vector<cvm::atom> const &atoms);
+
+  /// \brief Add an atom to this group  
+  void add_atom (cvm::atom const &a);
+
+  /// \brief Default constructor
+  atom_group();
+
+  /// \brief Destructor
+  ~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
+  void read_positions();
+
+  /// \brief Move all positions
+  void apply_translation (cvm::rvector const &t);
+
+  /// \brief Rotate all positions
+  void apply_rotation (cvm::rotation const &q);
+
+
+  /// \brief Get the current velocities; this must be called always
+  /// *after* read_positions(); if b_rotate is defined, the same
+  /// rotation applied to the coordinates will be used
+  void read_velocities();
+
+  /// \brief Get the current system_forces; this must be called always
+  /// *after* read_positions(); if b_rotate is defined, the same
+  /// rotation applied to the coordinates will be used
+  void read_system_forces();
+
+  /// Call reset_data() for each atom
+  inline void reset_atoms_data()
+  {
+    for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++)
+      ai->reset_data();
+  }
+
+  /// \brief Return a copy of the current atom positions
+  std::vector<cvm::atom_pos> positions() const;
+
+  /// \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
+  std::vector<cvm::atom_pos> old_pos;
+
+
+  /// \brief Return a copy of the current atom velocities
+  std::vector<cvm::rvector> velocities() const;
+
+
+  /// \brief Return a copy of the system forces
+  std::vector<cvm::rvector> system_forces() const;
+  /// \brief Return a copy of the aggregated total force on the group
+  cvm::rvector system_force() const;
+
+
+  /// \brief Shorthand: save the specified gradient on each atom,
+  /// weighting with the atom mass (mostly used in combination with
+  /// \link center_of_mass() \endlink)
+  void set_weighted_gradient (cvm::rvector const &grad);
+
+  /// \brief Used by a (scalar) colvar to apply its force on its \link
+  /// atom_group \endlink members
+  ///
+  /// The (scalar) force is multiplied by the colvar gradient for each
+  /// atom; this should be used when a colvar with scalar \link
+  /// colvarvalue \endlink type is used (this is the most frequent
+  /// case: for colvars with a non-scalar type, the most convenient
+  /// solution is to sum together the Cartesian forces from all the
+  /// 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.
+  void apply_colvar_force (cvm::real const &force);
+
+  /// \brief Apply a force "to the center of mass", i.e. the force is
+  /// distributed on each atom according to its mass
+  ///
+  /// If the group is being rotated to a reference frame (e.g. to
+  /// express the colvar independently from the solute rotation), the
+  /// force is rotated back to the original frame.  Colvar gradients
+  /// are not used, either because they were not defined (e.g because
+  /// the colvar has not a scalar value) or the biases require to
+  /// micromanage the force.
+  void apply_force (cvm::rvector const &force);
+
+  /// \brief Apply an array of forces directly on the individual
+  /// atoms; the length of the specified vector must be the same of
+  /// this \link atom_group \endlink.
+  ///
+  /// If the group is being rotated to a reference frame (e.g. to
+  /// express the colvar independently from the solute rotation), the
+  /// forces are rotated back to the original frame.  Colvar gradients
+  /// are not used, either because they were not defined (e.g because
+  /// the colvar has not a scalar value) or the biases require to
+  /// micromanage the forces.
+  void apply_forces (std::vector<cvm::rvector> const &forces);
+};
+
+
+
+
+#endif
+
+
+// Emacs
+// Local Variables:
+// mode: C++
+// End:

lib/colvars/colvarbias.C

@@ -0,0 +1,428 @@
+#include "colvarmodule.h"
+#include "colvarvalue.h"
+#include "colvarbias.h"
+
+
+colvarbias::colvarbias (std::string const &conf, char const *key)
+  : colvarparse(), has_data (false)
+{
+  cvm::log ("Initializing a new \""+std::string (key)+"\" instance.\n");
+
+  size_t rank = 1;
+  std::string const key_str (key);
+
+  if (to_lower_cppstr (key_str) == std::string ("abf")) {
+    rank = cvm::n_abf_biases+1;
+  }
+  if (to_lower_cppstr (key_str) == std::string ("harmonic")) {
+    rank = cvm::n_harm_biases+1;
+  }
+  if (to_lower_cppstr (key_str) == std::string ("histogram")) {
+    rank = cvm::n_histo_biases+1;
+  }
+  if (to_lower_cppstr (key_str) == std::string ("metadynamics")) {
+    rank = cvm::n_meta_biases+1;
+  }
+
+  get_keyval (conf, "name", name, key_str+cvm::to_str (rank));
+
+  // lookup the associated colvars
+  std::vector<std::string> colvars_str;
+  if (get_keyval (conf, "colvars", colvars_str)) {
+    for (size_t i = 0; i < colvars_str.size(); i++) {
+      add_colvar (colvars_str[i]);
+    }
+  }
+
+  if (!colvars.size()) {
+    cvm::fatal_error ("Error: no collective variables specified.\n");
+  }
+}
+
+
+colvarbias::colvarbias()
+  : colvarparse(), has_data (false)
+{}
+
+
+void colvarbias::add_colvar (std::string const &cv_name)
+{
+  if (colvar *cvp = cvm::colvar_p (cv_name)) {
+    cvp->enable (colvar::task_gradients);
+    if (cvm::debug())
+      cvm::log ("Applying this bias to collective variable \""+
+                cvp->name+"\".\n"); 
+    colvars.push_back (cvp);
+    colvar_forces.push_back (colvarvalue (cvp->type()));
+  } else {
+    cvm::fatal_error ("Error: cannot find a colvar named \""+
+                      cv_name+"\".\n");
+  }
+}
+
+
+void colvarbias::communicate_forces()
+{
+  for (size_t i = 0; i < colvars.size(); i++) {
+    if (cvm::debug()) {
+      cvm::log ("Communicating a force to colvar \""+
+                colvars[i]->name+"\", of type \""+
+                colvarvalue::type_desc[colvars[i]->type()]+"\".\n");
+    }
+    colvars[i]->add_bias_force (colvar_forces[i]);
+  }
+}    
+
+
+
+colvarbias_harmonic::colvarbias_harmonic (std::string const &conf,
+                                          char const *key)
+  : colvarbias (conf, key), 
+    target_nsteps (0)
+{
+  get_keyval (conf, "forceConstant", force_k, 1.0);
+  for (size_t i = 0; i < colvars.size(); i++) {
+    if (colvars[i]->width != 1.0)
+      cvm::log ("The force constant for colvar \""+colvars[i]->name+
+                "\" will be rescaled to "+
+                cvm::to_str (force_k/(colvars[i]->width*colvars[i]->width))+
+                " according to the specified width.\n");
+  }
+
+  // get the initial restraint centers
+  colvar_centers.resize (colvars.size());
+  colvar_centers_raw.resize (colvars.size());
+  for (size_t i = 0; i < colvars.size(); i++) {
+    colvar_centers[i].type (colvars[i]->type());
+    colvar_centers_raw[i].type (colvars[i]->type());
+  }
+  if (get_keyval (conf, "centers", colvar_centers, colvar_centers)) {
+    for (size_t i = 0; i < colvars.size(); i++) {
+      colvar_centers[i].apply_constraints();
+      colvar_centers_raw[i] = colvar_centers[i];
+    }
+  } else {
+    colvar_centers.clear();
+    cvm::fatal_error ("Error: must define the initial centers of the restraints.\n");
+  }
+
+  if (colvar_centers.size() != colvars.size())
+    cvm::fatal_error ("Error: number of harmonic centers does not match "
+                      "that of collective variables.\n");
+
+  if (get_keyval (conf, "targetCenters", target_centers, colvar_centers)) {
+    b_chg_centers = true;
+    target_nstages = 0;
+    for (size_t i = 0; i < target_centers.size(); i++) {
+      target_centers[i].apply_constraints();
+    }
+  } else {
+    b_chg_centers = false;
+    target_centers.clear();
+  }
+
+  if (get_keyval (conf, "targetForceConstant", target_force_k, 0.0)) {
+    if (b_chg_centers)
+      cvm::fatal_error ("Error: cannot specify both targetCenters and targetForceConstant.\n");
+
+    starting_force_k = force_k;
+    b_chg_force_k = true;
+
+    get_keyval (conf, "lambdaSchedule", lambda_schedule, lambda_schedule);
+    if (lambda_schedule.size()) {
+      // There is one more lambda-point than stages
+      target_nstages = lambda_schedule.size() - 1;
+    } else {
+      target_nstages = 0;
+    }
+  } else {
+    b_chg_force_k = false;
+  }
+
+  if (b_chg_centers || b_chg_force_k) {
+    get_keyval (conf, "targetNumSteps", target_nsteps, 0);
+    if (!target_nsteps)
+      cvm::fatal_error ("Error: targetNumSteps must be non-zero.\n");
+
+    if (get_keyval (conf, "targetNumStages", target_nstages, target_nstages) &&
+        lambda_schedule.size()) {
+      cvm::fatal_error ("Error: targetNumStages and lambdaSchedule are incompatible.\n");
+    }
+
+    if (target_nstages) {
+      // This means that either numStages of lambdaSchedule has been provided
+      stage = -1;
+      restraint_FE = 0.0;
+    }
+
+    if (get_keyval (conf, "targetForceExponent", force_k_exp, 1.0)) {
+      if (! b_chg_force_k)
+        cvm::log ("Warning: not changing force constant: targetForceExponent will be ignored\n");
+      if (force_k_exp < 1.0)
+        cvm::log ("Warning: for all practical purposes, targetForceExponent should be 1.0 or greater.\n");
+    }
+  }
+
+  if (cvm::debug())
+    cvm::log ("Done initializing a new harmonic restraint bias.\n");
+}
+
+
+cvm::real colvarbias_harmonic::update()
+{
+  bias_energy = 0.0;
+
+  if (cvm::debug())
+    cvm::log ("Updating the harmonic bias \""+this->name+"\".\n");
+
+  // Setup first stage of staged variable force constant calculation
+  if (b_chg_force_k && target_nstages && stage == -1) {
+    stage = 0;
+    cvm::real lambda;
+    if (lambda_schedule.size()) {
+      lambda = lambda_schedule[0];
+    } else {
+      lambda = 0.0;
+    }
+    force_k = starting_force_k + (target_force_k - starting_force_k)
+              * std::pow (lambda, force_k_exp);
+    cvm::log ("Harmonic restraint, stage " + cvm::to_str(stage) +
+        " : lambda = " + cvm::to_str(lambda));
+    cvm::log ("Setting force constant to " + cvm::to_str (force_k));
+  }
+  
+  // Force and energy calculation
+  for (size_t i = 0; i < colvars.size(); i++) {
+    colvar_forces[i] =
+      (-0.5) * force_k /
+      (colvars[i]->width * colvars[i]->width) *
+      colvars[i]->dist2_lgrad (colvars[i]->value(),
+                               colvar_centers[i]);
+    bias_energy += 0.5 * force_k / (colvars[i]->width * colvars[i]->width) *
+              colvars[i]->dist2(colvars[i]->value(), colvar_centers[i]);
+    if (cvm::debug())
+      cvm::log ("dist_grad["+cvm::to_str (i)+
+                "] = "+cvm::to_str (colvars[i]->dist2_lgrad (colvars[i]->value(),
+                               colvar_centers[i]))+"\n");
+  }
+
+  if (cvm::debug())
+    cvm::log ("Current forces for the harmonic bias \""+
+              this->name+"\": "+cvm::to_str (colvar_forces)+".\n");
+
+  if (b_chg_centers) {
+
+    if (!centers_incr.size()) {
+      // if this is the first calculation, calculate the advancement
+      // at each simulation step (or stage, if applicable)
+      // (take current stage into account: it can be non-zero
+      //  if we are restarting a staged calculation)
+      centers_incr.resize (colvars.size());
+      for (size_t i = 0; i < colvars.size(); i++) {
+        centers_incr[i].type (colvars[i]->type());
+        centers_incr[i] = (target_centers[i] - colvar_centers[i]) /
+          cvm::real ( target_nstages ? (target_nstages - stage) :
+                                      (target_nsteps - cvm::step_absolute()));
+      }
+      if (cvm::debug())
+        cvm::log ("Center increment for the harmonic bias \""+
+                  this->name+"\": "+cvm::to_str (centers_incr)+".\n");
+    }
+
+    if (cvm::debug())
+      cvm::log ("Current centers for the harmonic bias \""+
+                this->name+"\": "+cvm::to_str (colvar_centers)+".\n");
+
+    if (target_nstages) {
+      if (cvm::step_absolute() > 0
+            && (cvm::step_absolute() % target_nsteps) == 0
+            && stage < target_nstages) {
+
+          // any per-stage calculation, e.g. free energy stuff
+          // should be done here
+          for (size_t i = 0; i < colvars.size(); i++) {
+            colvar_centers_raw[i] += centers_incr[i];
+            colvar_centers[i] = colvar_centers_raw[i];
+            colvars[i]->wrap(colvar_centers[i]);
+            colvar_centers[i].apply_constraints();
+          }
+          stage++;
+          cvm::log ("Moving restraint stage " + cvm::to_str(stage) +
+              " : setting centers to " + cvm::to_str (colvar_centers));
+      }
+    } else if (cvm::step_absolute() < target_nsteps) {
+      // move the restraint centers in the direction of the targets
+      // (slow growth)
+      for (size_t i = 0; i < colvars.size(); i++) {
+        colvar_centers_raw[i] += centers_incr[i];
+        colvar_centers[i] = colvar_centers_raw[i];
+        colvars[i]->wrap(colvar_centers[i]);
+        colvar_centers[i].apply_constraints();
+      }
+    }
+  }
+
+  if (b_chg_force_k) {
+    // Coupling parameter, between 0 and 1
+    cvm::real lambda;
+
+    if (target_nstages) {
+      // Square distance for restraint energy calculation
+      // normalized by square colvar width
+      cvm::real dist_sq = 0.0;
+      for (size_t i = 0; i < colvars.size(); i++) {
+        dist_sq += colvars[i]->dist2 (colvars[i]->value(), colvar_centers[i])
+          / (colvars[i]->width * colvars[i]->width);
+      }
+
+      // TI calculation: estimate free energy derivative
+      // need current lambda
+      if (lambda_schedule.size()) {
+        lambda = lambda_schedule[stage];
+      } else {
+        lambda = cvm::real(stage) / cvm::real(target_nstages);
+      }
+      restraint_FE += 0.5 * force_k_exp * std::pow(lambda, force_k_exp - 1.0)
+        * (target_force_k - starting_force_k) * dist_sq;
+
+      // Finish current stage...
+      if (cvm::step_absolute() % target_nsteps == 0 &&
+          cvm::step_absolute() > 0) {
+
+          cvm::log ("Lambda= " + cvm::to_str (lambda) + " dA/dLambda= "
+              + cvm::to_str (restraint_FE / cvm::real(target_nsteps)));
+      
+        //  ...and move on to the next one
+        if (stage < target_nstages) {
+
+          restraint_FE = 0.0;
+          stage++;
+          if (lambda_schedule.size()) {
+            lambda = lambda_schedule[stage];
+          } else {
+            lambda = cvm::real(stage) / cvm::real(target_nstages);
+          }
+          force_k = starting_force_k + (target_force_k - starting_force_k)
+                    * std::pow (lambda, force_k_exp);
+          cvm::log ("Harmonic restraint, stage " + cvm::to_str(stage) +
+              " : lambda = " + cvm::to_str(lambda));
+          cvm::log ("Setting force constant to " + cvm::to_str (force_k));
+        }
+      }
+    } else if (cvm::step_absolute() <= target_nsteps) {
+      // update force constant (slow growth)
+      lambda = cvm::real(cvm::step_absolute()) / cvm::real(target_nsteps);
+      force_k = starting_force_k + (target_force_k - starting_force_k)
+          * std::pow (lambda, force_k_exp);
+    }
+  }
+
+  if (cvm::debug())
+    cvm::log ("Done updating the harmonic bias \""+this->name+"\".\n");
+  return bias_energy;
+}
+
+
+std::istream & colvarbias_harmonic::read_restart (std::istream &is)
+{
+  size_t const start_pos = is.tellg();
+
+  cvm::log ("Restarting harmonic bias \""+
+            this->name+"\".\n");
+
+  std::string key, brace, conf;
+  if ( !(is >> key)   || !(key == "harmonic") ||
+       !(is >> brace) || !(brace == "{") ||
+       !(is >> colvarparse::read_block ("configuration", conf)) ) {
+
+    cvm::log ("Error: in reading restart configuration for harmonic bias \""+
+              this->name+"\" at position "+
+              cvm::to_str (is.tellg())+" in stream.\n");
+    is.clear();
+    is.seekg (start_pos, std::ios::beg);
+    is.setstate (std::ios::failbit);
+    return is;
+  }
+
+//   int id = -1; 
+  std::string name = "";
+//   if ( ( (colvarparse::get_keyval (conf, "id", id, -1, colvarparse::parse_silent)) &&
+//          (id != this->id) ) ||
+  if ( (colvarparse::get_keyval (conf, "name", name, std::string (""), colvarparse::parse_silent)) &&
+       (name != this->name) )
+    cvm::fatal_error ("Error: in the restart file, the "
+                      "\"harmonic\" block has a wrong name\n");
+//   if ( (id == -1) && (name == "") ) {
+  if (name.size() == 0) {
+    cvm::fatal_error ("Error: \"harmonic\" block in the restart file "
+                      "has no identifiers.\n");
+  }
+
+  if (b_chg_centers) {
+    cvm::log ("Reading the updated restraint centers from the restart.\n");
+    if (!get_keyval (conf, "centers", colvar_centers))
+      cvm::fatal_error ("Error: restraint centers are missing from the restart.\n");
+    if (!get_keyval (conf, "centers_raw", colvar_centers_raw))
+      cvm::fatal_error ("Error: \"raw\" restraint centers are missing from the restart.\n");
+  }
+
+  if (b_chg_force_k) {
+    cvm::log ("Reading the updated force constant from the restart.\n");
+    if (!get_keyval (conf, "forceConstant", force_k))
+      cvm::fatal_error ("Error: force constant is missing from the restart.\n");
+  }
+
+  if (target_nstages) {
+    cvm::log ("Reading current stage from the restart.\n");
+    if (!get_keyval (conf, "stage", stage))
+      cvm::fatal_error ("Error: current stage is missing from the restart.\n");
+  }
+
+  is >> brace;
+  if (brace != "}") {
+    cvm::fatal_error ("Error: corrupt restart information for harmonic bias \""+
+                      this->name+"\": no matching brace at position "+
+                      cvm::to_str (is.tellg())+" in the restart file.\n");
+    is.setstate (std::ios::failbit);
+  }
+  return is;
+}
+
+
+std::ostream & colvarbias_harmonic::write_restart (std::ostream &os)
+{
+  os << "harmonic {\n"
+     << "  configuration {\n"
+    //      << "    id " << this->id << "\n"
+     << "    name " << this->name << "\n";
+
+  if (b_chg_centers) {
+    os << "    centers ";
+    for (size_t i = 0; i < colvars.size(); i++) {
+      os << " " << colvar_centers[i];
+    }
+    os << "\n";
+    os << "    centers_raw ";
+    for (size_t i = 0; i < colvars.size(); i++) {
+      os << " " << colvar_centers_raw[i];
+    }
+    os << "\n";
+  }
+
+  if (b_chg_force_k) {
+    os << "    forceConstant "
+       << std::setprecision (cvm::en_prec)
+       << std::setw (cvm::en_width) << force_k << "\n";
+  }
+
+  if (target_nstages) {
+    os << "    stage " << std::setw (cvm::it_width)
+       << stage << "\n";
+  }
+
+  os << "  }\n"
+     << "}\n\n";
+
+  return os;
+}
+

lib/colvars/colvarbias.h

@@ -0,0 +1,152 @@
+#ifndef COLVARBIAS_H
+#define COLVARBIAS_H
+
+#include "colvar.h"
+#include "colvarparse.h"
+
+
+/// \brief Collective variable bias, base class
+class colvarbias : public colvarparse {
+public:
+
+  /// Numeric id of this bias
+  int            id;
+
+  /// Name of this bias
+  std::string    name;
+  
+  /// Add a new collective variable to this bias
+  void add_colvar (std::string const &cv_name);
+
+  /// Retrieve colvar values and calculate their biasing forces
+  /// Return bias energy
+  virtual cvm::real update() = 0;
+
+  /// Perform analysis tasks
+  virtual inline void analyse() {}
+
+  /// Send forces to the collective variables
+  void communicate_forces();
+
+  /// \brief Constructor
+  /// 
+  /// The constructor of the colvarbias base class is protected, so
+  /// that it can only be called from inherited classes
+  colvarbias (std::string const &conf, char const *key);
+
+  /// Default constructor
+  colvarbias();
+
+  /// Destructor
+  virtual inline ~colvarbias() {}
+
+  /// Read the bias configuration from a restart file
+  virtual std::istream & read_restart (std::istream &is) = 0;
+
+  /// Write the bias configuration to a restart file
+  virtual std::ostream & write_restart (std::ostream &os) = 0;
+
+protected:
+
+  /// \brief Pointers to collective variables to which the bias is
+  /// applied; current values and metric functions will be obtained
+  /// through each colvar object
+  std::vector<colvar *>    colvars;
+
+  /// \brief Current forces from this bias to the colvars
+  std::vector<colvarvalue> colvar_forces;
+
+  /// \brief Current energy of this bias (colvar_forces should be
+  /// obtained by deriving this)
+  cvm::real                bias_energy;
+
+  /// \brief Whether this bias has already accumulated information
+  /// (when relevant)
+  bool                     has_data;
+
+};
+
+
+/// \brief Harmonic restraint, optionally moving towards a target
+/// (implementation of \link colvarbias \endlink)
+class colvarbias_harmonic : public colvarbias {
+
+public:
+
+  /// Retrieve colvar values and calculate their biasing forces
+  virtual cvm::real update();
+
+  /// Read the bias configuration from a restart file
+  virtual std::istream & read_restart (std::istream &is);
+
+  /// Write the bias configuration to a restart file
+  virtual std::ostream & write_restart (std::ostream &os);
+
+  /// \brief Constructor
+  colvarbias_harmonic (std::string const &conf, char const *key);
+
+  /// Destructor
+  virtual inline ~colvarbias_harmonic() {}
+
+
+protected:
+
+  /// \brief Restraint centers
+  std::vector<colvarvalue> colvar_centers;
+
+  /// \brief Restraint centers without wrapping or constraints applied
+  std::vector<colvarvalue> colvar_centers_raw;
+
+  /// \brief Restraint force constant
+  cvm::real force_k;
+
+  /// \brief Moving target?
+  bool b_chg_centers;
+
+  /// \brief Changing force constant?
+  bool b_chg_force_k;
+
+  /// \brief Restraint force constant (target value)
+  cvm::real target_force_k;
+
+  /// \brief Restraint force constant (starting value)
+  cvm::real starting_force_k;
+
+  /// \brief Lambda-schedule for custom varying force constant
+  std::vector<cvm::real> lambda_schedule;
+
+  /// \brief New restraint centers
+  std::vector<colvarvalue> target_centers;
+
+  /// \brief Amplitude of the restraint centers' increment at each step
+  /// (or stage) towards the new values (calculated from target_nsteps)
+  std::vector<colvarvalue> centers_incr;
+
+  /// \brief Exponent for varying the force constant
+  cvm::real force_k_exp;
+
+  /// \brief Number of steps required to reach the target force constant
+  /// or restraint centers
+  size_t target_nsteps;
+
+  /// \brief Number of stages over which to perform the change
+  /// If zero, perform a continuous change
+  size_t target_nstages;
+
+  /// \brief Number of current stage of the perturbation
+  size_t stage;
+  
+  /// \brief Intermediate quantity to compute the restraint free energy
+  /// (in TI, would be the accumulating FE derivative)
+  cvm::real restraint_FE;
+};
+
+
+#endif
+
+
+
+// Emacs
+// Local Variables:
+// mode: C++
+// End:

lib/colvars/colvarbias_abf.C

@@ -0,0 +1,492 @@
+/********************************************************************************
+ * Implementation of the ABF and histogram biases                               *
+ ********************************************************************************/
+
+#include "colvarmodule.h"
+#include "colvar.h"
+#include "colvarbias_abf.h"
+
+/// ABF bias constructor; parses the config file
+
+colvarbias_abf::colvarbias_abf (std::string const &conf, char const *key)
+  : colvarbias (conf, key),
+    gradients (NULL),
+    samples (NULL)
+{
+  if (cvm::temperature() == 0.0)
+    cvm::log ("WARNING: ABF should not be run without a thermostat or at 0 Kelvin!\n");
+  
+  // ************* parsing general ABF options ***********************
+
+  get_keyval (conf, "applyBias",  apply_bias, true);
+  if (!apply_bias) cvm::log ("WARNING: ABF biases will *not* be applied!\n");
+
+  get_keyval (conf, "updateBias",  update_bias, true);
+  if (!update_bias) cvm::log ("WARNING: ABF biases will *not* be updated!\n");
+
+  get_keyval (conf, "hideJacobian", hide_Jacobian, false);
+  if (hide_Jacobian) {
+    cvm::log ("Jacobian (geometric) forces will be handled internally.\n");
+  } else {
+    cvm::log ("Jacobian (geometric) forces will be included in reported free energy gradients.\n");
+  }
+
+  get_keyval (conf, "fullSamples", full_samples, 200);
+  if ( full_samples <= 1 ) full_samples = 1; 
+  min_samples = full_samples / 2;
+  // full_samples - min_samples >= 1 is guaranteed
+
+  get_keyval (conf, "inputPrefix",  input_prefix, std::vector<std::string> ());
+  get_keyval (conf, "outputFreq", output_freq, cvm::restart_out_freq);
+  get_keyval (conf, "historyFreq", history_freq, 0);
+  b_history_files = (history_freq > 0);
+
+  // ************* checking the associated colvars *******************
+
+  if (colvars.size() == 0) {
+    cvm::fatal_error ("Error: no collective variables specified for the ABF bias.\n");
+  }
+
+  for (size_t i = 0; i < colvars.size(); i++) {
+
+    if (colvars[i]->type() != colvarvalue::type_scalar) {
+      cvm::fatal_error ("Error: ABF bias can only use scalar-type variables.\n");
+    }
+
+    colvars[i]->enable (colvar::task_gradients);
+
+    if (update_bias) {
+      // Request calculation of system force (which also checks for availability)
+      colvars[i]->enable (colvar::task_system_force);
+
+      if (!colvars[i]->tasks[colvar::task_extended_lagrangian]) {
+        // request computation of Jacobian force
+        colvars[i]->enable (colvar::task_Jacobian_force);
+
+        // request Jacobian force as part as system force
+        // except if the user explicitly requires the "silent" Jacobian
+        // correction AND the colvar has a single component
+        if (hide_Jacobian) {
+          if (colvars[i]->n_components() > 1) {
+            cvm::log ("WARNING: colvar \"" + colvars[i]->name
+            + "\" has multiple components; reporting its Jacobian forces\n");
+            colvars[i]->enable (colvar::task_report_Jacobian_force);
+          }
+        } else {
+          colvars[i]->enable (colvar::task_report_Jacobian_force);
+        }
+      }
+    }
+
+    // Here we could check for orthogonality of the Cartesian coordinates
+    // and make it just a warning if some parameter is set? 
+  }
+
+  bin.assign (colvars.size(), 0);
+  prev_bin.assign (colvars.size(), 0);
+  prev_force = new cvm::real [colvars.size()];
+
+  // Construct empty grids based on the colvars
+  samples   = new colvar_grid_count    (colvars);
+  gradients = new colvar_grid_gradient (colvars);
+  gradients->samples = samples;
+  samples->has_parent_data = true;
+
+  // If custom grids are provided, read them
+  if ( input_prefix.size() > 0 ) {
+    read_gradients_samples ();
+  }
+  
+  cvm::log ("Finished ABF setup.\n");
+}
+
+/// Destructor
+colvarbias_abf::~colvarbias_abf()
+{
+  if (samples) {
+    delete samples;
+    samples = NULL;
+  }
+
+  if (gradients) {
+    delete gradients;
+    gradients = NULL;
+  }
+
+  delete [] prev_force;
+}
+
+
+/// Update the FE gradient, compute and apply biasing force
+/// also output data to disk if needed
+
+cvm::real colvarbias_abf::update()
+{
+  if (cvm::debug()) cvm::log ("Updating ABF bias " + this->name);
+
+  if (cvm::step_relative() == 0) {
+	
+    // At first timestep, do only:
+    // initialization stuff (file operations relying on n_abf_biases
+    // compute current value of colvars
+
+    if ( cvm::n_abf_biases == 1 && cvm::n_meta_biases == 0 ) {
+      // This is the only ABF bias
+      output_prefix = cvm::output_prefix;
+    } else {
+      output_prefix = cvm::output_prefix + "." + this->name;
+    }
+
+    for (size_t i=0; i<colvars.size(); i++) {
+      bin[i] = samples->current_bin_scalar(i);
+    }
+
+  } else {
+
+    for (size_t i=0; i<colvars.size(); i++) {
+      bin[i] = samples->current_bin_scalar(i);
+    }
+
+    if ( update_bias && samples->index_ok (prev_bin) ) {
+      // Only if requested and within bounds of the grid...
+
+      for (size_t i=0; i<colvars.size(); i++) {	  // get forces (lagging by 1 timestep) from colvars
+        prev_force[i] = colvars[i]->system_force();
+      }
+      gradients->acc_force (prev_bin, prev_force);
+    }
+  }
+
+  // save bin for next timestep
+  prev_bin = bin;  
+
+  // Reset biasing forces from previous timestep
+  for (size_t i=0; i<colvars.size(); i++) {
+    colvar_forces[i].reset();
+  }
+
+  // Compute and apply the new bias, if applicable
+  if ( apply_bias && samples->index_ok (bin) ) {
+
+    size_t  count = samples->value (bin);
+    cvm::real	fact = 1.0;
+
+    // Factor that ensures smooth introduction of the force
+    if ( count < full_samples ) {
+      fact = ( count < min_samples) ? 0.0 :
+        (cvm::real (count - min_samples)) / (cvm::real (full_samples - min_samples));
+    }
+	
+    const cvm::real * grad  = &(gradients->value (bin));
+
+    if ( fact != 0.0 ) {
+
+      if ( (colvars.size() == 1) && colvars[0]->periodic_boundaries() ) {
+        // Enforce a zero-mean bias on periodic, 1D coordinates
+        colvar_forces[0].real_value += fact * (grad[0] / cvm::real (count) - gradients->average ());
+      } else {
+        for (size_t i=0; i<colvars.size(); i++) {
+          // subtracting the mean force (opposite of the FE gradient) means adding the gradient
+          colvar_forces[i].real_value += fact * grad[i] / cvm::real (count);
+          // without .real_value, the above would do (cheap) runtime type checking
+        }
+      }
+    }
+  }
+
+  if (output_freq && (cvm::step_absolute() % output_freq) == 0) {
+    if (cvm::debug()) cvm::log ("ABF bias trying to write gradients and samples to disk");
+    write_gradients_samples (output_prefix);
+  }
+  if (b_history_files && (cvm::step_absolute() % history_freq) == 0) {
+    // append to existing file only if cvm::step_absolute() > 0
+    // otherwise, backup and replace
+    write_gradients_samples (output_prefix + ".hist", (cvm::step_absolute() > 0));
+  }
+  return 0.0; // TODO compute bias energy whenever possible (i.e. 1D with updateBias off)
+}
+
+
+void colvarbias_abf::write_gradients_samples (const std::string &prefix, bool append)
+{
+  std::string  samples_out_name = prefix + ".count";
+  std::string  gradients_out_name = prefix + ".grad";
+  std::ios::openmode mode = (append ? std::ios::app : std::ios::out);
+
+  std::ofstream samples_os;
+  std::ofstream gradients_os;
+
+  if (!append) cvm::backup_file (samples_out_name.c_str());
+  samples_os.open (samples_out_name.c_str(), mode);
+  if (!samples_os.good()) cvm::fatal_error ("Error opening ABF samples file " + samples_out_name + " for writing");
+  samples->write_multicol (samples_os);
+  samples_os.close ();
+
+  if (!append) cvm::backup_file (gradients_out_name.c_str());
+  gradients_os.open (gradients_out_name.c_str(), mode);
+  if (!gradients_os.good())	cvm::fatal_error ("Error opening ABF gradient file " + gradients_out_name + " for writing");
+  gradients->write_multicol (gradients_os);
+  gradients_os.close ();
+
+  if (colvars.size () == 1) {
+    std::string  pmf_out_name = prefix + ".pmf";
+    if (!append) cvm::backup_file (pmf_out_name.c_str());
+    std::ofstream pmf_os;
+    // Do numerical integration and output a PMF
+    pmf_os.open (pmf_out_name.c_str(), mode);
+    if (!pmf_os.good())	cvm::fatal_error ("Error opening pmf file " + pmf_out_name + " for writing");
+    gradients->write_1D_integral (pmf_os);
+    pmf_os << std::endl;
+    pmf_os.close ();
+  }
+  return;
+}
+
+void colvarbias_abf::read_gradients_samples ()
+{
+  std::string samples_in_name, gradients_in_name;
+
+  for ( size_t i = 0; i < input_prefix.size(); i++ ) {
+    samples_in_name = input_prefix[i] + ".count";
+    gradients_in_name = input_prefix[i] + ".grad";
+    // For user-provided files, the per-bias naming scheme may not apply
+    
+    std::ifstream is;
+
+    cvm::log ("Reading sample count from " + samples_in_name + " and gradients from " + gradients_in_name);
+    is.open (samples_in_name.c_str());
+    if (!is.good()) cvm::fatal_error ("Error opening ABF samples file " + samples_in_name + " for reading");
+    samples->read_multicol (is, true);
+    is.close ();
+    is.clear();
+
+    is.open (gradients_in_name.c_str());
+    if (!is.good())	cvm::fatal_error ("Error opening ABF gradient file " + gradients_in_name + " for reading");
+    gradients->read_multicol (is, true);
+    is.close ();
+  }
+  return;
+}
+
+
+std::ostream & colvarbias_abf::write_restart (std::ostream& os)
+{
+
+  std::ios::fmtflags flags (os.flags ()); 
+  os.setf(std::ios::fmtflags (0), std::ios::floatfield); // default floating-point format
+
+  os << "abf {\n"
+     << "  configuration {\n"
+     << "    name " << this->name << "\n";
+  os << "  }\n";
+
+  os << "samples\n";
+  samples->write_raw (os, 8);
+
+  os << "\ngradient\n";
+  gradients->write_raw (os);
+
+  os << "}\n\n";
+
+  os.flags (flags);
+  return os;
+}
+
+
+std::istream & colvarbias_abf::read_restart (std::istream& is)
+{
+  if ( input_prefix.size() > 0 ) {
+    cvm::fatal_error ("ERROR: cannot provide both inputPrefix and restart information (colvarsInput)");
+  }
+
+  size_t const start_pos = is.tellg();
+
+  cvm::log ("Restarting ABF bias \""+
+            this->name+"\".\n");
+  std::string key, brace, conf;
+
+  if ( !(is >> key)   || !(key == "abf") ||
+       !(is >> brace) || !(brace == "{") ||
+       !(is >> colvarparse::read_block ("configuration", conf)) ) {
+    cvm::log ("Error: in reading restart configuration for ABF bias \""+
+              this->name+"\" at position "+
+              cvm::to_str (is.tellg())+" in stream.\n");
+    is.clear();
+    is.seekg (start_pos, std::ios::beg);
+    is.setstate (std::ios::failbit);
+    return is;
+  }
+
+  std::string name = "";
+  if ( (colvarparse::get_keyval (conf, "name", name, std::string (""), colvarparse::parse_silent)) &&
+         (name != this->name) )
+    cvm::fatal_error ("Error: in the restart file, the "
+                      "\"abf\" block has wrong name (" + name + ")\n");
+  if ( name == "" ) {
+    cvm::fatal_error ("Error: \"abf\" block in the restart file has no name.\n");
+  }
+  
+  if ( !(is >> key)   || !(key == "samples")) {
+    cvm::log ("Error: in reading restart configuration for ABF bias \""+
+              this->name+"\" at position "+
+              cvm::to_str (is.tellg())+" in stream.\n");
+    is.clear();
+    is.seekg (start_pos, std::ios::beg);
+    is.setstate (std::ios::failbit);
+    return is;
+  }
+  if (! samples->read_raw (is)) {
+    samples->read_raw_error();
+  }
+
+  if ( !(is >> key)   || !(key == "gradient")) {
+    cvm::log ("Error: in reading restart configuration for ABF bias \""+
+              this->name+"\" at position "+
+              cvm::to_str (is.tellg())+" in stream.\n");
+    is.clear();
+    is.seekg (start_pos, std::ios::beg);
+    is.setstate (std::ios::failbit);
+    return is;
+  }
+  if (! gradients->read_raw (is)) {
+    gradients->read_raw_error();
+  }
+
+  is >> brace;
+  if (brace != "}") {
+    cvm::fatal_error ("Error: corrupt restart information for ABF bias \""+
+                      this->name+"\": no matching brace at position "+
+                      cvm::to_str (is.tellg())+" in the restart file.\n");
+    is.setstate (std::ios::failbit);
+  }
+  return is;
+}
+
+
+
+
+/// Histogram "bias" constructor
+
+colvarbias_histogram::colvarbias_histogram (std::string const &conf, char const *key)
+  : colvarbias (conf, key),
+    grid (NULL)
+{
+  get_keyval (conf, "outputfreq", output_freq, cvm::restart_out_freq);
+
+  if ( output_freq == 0 ) {
+    cvm::fatal_error ("User required histogram with zero output frequency");
+  }
+
+  grid   = new colvar_grid_count    (colvars);
+  bin.assign (colvars.size(), 0);
+
+  out_name = cvm::output_prefix + "." + this->name + ".dat";
+  cvm::log ("Histogram will be written to file " + out_name); 
+
+  cvm::log ("Finished histogram setup.\n");
+}
+
+/// Destructor
+colvarbias_histogram::~colvarbias_histogram()
+{
+  if (grid_os.good())	grid_os.close();
+
+  if (grid) {
+    delete grid;
+    grid = NULL;
+  }
+}
+
+/// Update the grid
+cvm::real colvarbias_histogram::update()
+{
+  if (cvm::debug()) cvm::log ("Updating Grid bias " + this->name);
+
+  for (size_t i=0; i<colvars.size(); i++) {
+    bin[i] = grid->current_bin_scalar(i);
+  }
+
+  if ( grid->index_ok (bin) ) {	  // Only within bounds of the grid...
+    grid->incr_count (bin);
+  }
+
+  if (output_freq && (cvm::step_absolute() % output_freq) == 0) {
+    if (cvm::debug()) cvm::log ("Histogram bias trying to write grid to disk");
+
+    grid_os.open (out_name.c_str());
+    if (!grid_os.good()) cvm::fatal_error ("Error opening histogram file " + out_name + " for writing");
+    grid->write_multicol (grid_os);
+    grid_os.close ();
+  }
+  return 0.0; // no bias energy for histogram
+}
+
+
+std::istream & colvarbias_histogram::read_restart (std::istream& is)
+{
+  size_t const start_pos = is.tellg();
+
+  cvm::log ("Restarting collective variable histogram \""+
+            this->name+"\".\n");
+  std::string key, brace, conf;
+
+  if ( !(is >> key)   || !(key == "histogram") ||
+       !(is >> brace) || !(brace == "{") ||
+       !(is >> colvarparse::read_block ("configuration", conf)) ) {
+    cvm::log ("Error: in reading restart configuration for histogram \""+
+              this->name+"\" at position "+
+              cvm::to_str (is.tellg())+" in stream.\n");
+    is.clear();
+    is.seekg (start_pos, std::ios::beg);
+    is.setstate (std::ios::failbit);
+    return is;
+  }
+
+  int id = -1;
+  std::string name = "";
+  if ( (colvarparse::get_keyval (conf, "name", name, std::string (""), colvarparse::parse_silent)) &&
+         (name != this->name) )
+    cvm::fatal_error ("Error: in the restart file, the "
+                      "\"histogram\" block has a wrong name: different system?\n");
+  if ( (id == -1) && (name == "") ) {
+    cvm::fatal_error ("Error: \"histogram\" block in the restart file "
+                      "has no name.\n");
+  }
+  
+  if ( !(is >> key)   || !(key == "grid")) {
+    cvm::log ("Error: in reading restart configuration for histogram \""+
+              this->name+"\" at position "+
+              cvm::to_str (is.tellg())+" in stream.\n");
+    is.clear();
+    is.seekg (start_pos, std::ios::beg);
+    is.setstate (std::ios::failbit);
+    return is;
+  }
+  if (! grid->read_raw (is)) {
+    grid->read_raw_error();
+  }
+
+  is >> brace;
+  if (brace != "}") {
+    cvm::fatal_error ("Error: corrupt restart information for ABF bias \""+
+                      this->name+"\": no matching brace at position "+
+                      cvm::to_str (is.tellg())+" in the restart file.\n");
+    is.setstate (std::ios::failbit);
+  }
+  return is;
+}
+
+std::ostream & colvarbias_histogram::write_restart (std::ostream& os)
+{
+  os << "histogram {\n"
+     << "  configuration {\n"
+     << "    name " << this->name << "\n";
+  os << "  }\n";
+
+  os << "grid\n";
+  grid->write_raw (os, 8);
+
+  os << "}\n\n";
+
+  return os;
+}

lib/colvars/colvarbias_abf.h

@@ -0,0 +1,95 @@
+/************************************************************************
+ * Headers for the ABF and histogram biases                             *
+ ************************************************************************/
+
+#ifndef COLVARBIAS_ABF_H
+#define COLVARBIAS_ABF_H
+
+#include <vector>
+#include <list>
+#include <sstream>
+#include <iomanip>
+//#include <cmath>
+
+#include "colvarbias.h"
+#include "colvargrid.h"
+
+typedef cvm::real* gradient_t;
+
+
+/// ABF bias
+class colvarbias_abf : public colvarbias {
+
+public:
+
+  colvarbias_abf (std::string const &conf, char const *key);
+  ~colvarbias_abf ();
+
+  cvm::real update ();
+
+private:
+
+  /// Filename prefix for human-readable gradient/sample count output
+  std::string	output_prefix;
+
+  /// Base filename(s) for reading previous gradient data (replaces data from restart file)
+  std::vector<std::string> input_prefix;
+
+  bool		apply_bias;
+  bool		update_bias;
+  bool		hide_Jacobian;
+  size_t	full_samples;
+  size_t	min_samples;
+  /// frequency for updating output files (default: same as restartFreq?)
+  int		output_freq;
+  /// Write combined files with a history of all output data?
+  bool      b_history_files;
+  size_t    history_freq;
+
+  // Internal data and methods
+
+  std::vector<int>  bin, prev_bin;
+  gradient_t	    prev_force;
+
+  /// n-dim grid of free energy gradients
+  colvar_grid_gradient  *gradients;
+  /// n-dim grid of number of samples
+  colvar_grid_count     *samples;
+
+  /// Write human-readable FE gradients and sample count
+  void		  write_gradients_samples (const std::string &prefix, bool append = false);
+
+  /// Read human-readable FE gradients and sample count (if not using restart)
+  void		  read_gradients_samples ();
+
+  std::istream& read_restart  (std::istream&);
+  std::ostream& write_restart (std::ostream&);
+};
+
+
+/// Histogram "bias" (does as the name says)
+class colvarbias_histogram : public colvarbias {
+
+public:
+
+  colvarbias_histogram (std::string const &conf, char const *key);
+  ~colvarbias_histogram ();
+
+  cvm::real update ();
+
+private:
+
+  /// n-dim histogram
+  colvar_grid_count    *grid;
+  std::vector<int>  bin;
+  std::string	  out_name;
+
+  int		  output_freq;  
+  void		  write_grid ();
+  std::ofstream	  grid_os;  /// Stream for writing grid to disk
+
+  std::istream& read_restart  (std::istream&);
+  std::ostream& write_restart (std::ostream&);
+};
+
+#endif

lib/colvars/colvarbias_meta.h

@@ -0,0 +1,414 @@
+#ifndef COLVARBIAS_META_H
+#define COLVARBIAS_META_H
+
+#include <vector>
+#include <list>
+#include <sstream>
+#include <fstream>
+
+#include "colvarbias.h"
+#include "colvargrid.h"
+
+/// Metadynamics bias (implementation of \link colvarbias \endlink)
+class colvarbias_meta : public colvarbias {
+
+public:
+
+  /// Communication between different replicas
+  enum Communication {
+    /// One replica (default)
+    single_replica,
+    /// Hills added concurrently by several replicas
+    multiple_replicas
+  };
+
+  /// Communication between different replicas
+  Communication comm;
+
+  /// Constructor
+  colvarbias_meta (std::string const &conf, char const *key);
+
+  /// Default constructor
+  colvarbias_meta();
+
+  /// Destructor
+  virtual ~colvarbias_meta();
+  
+  virtual cvm::real update();
+
+  virtual std::istream & read_restart (std::istream &is);
+
+  virtual std::ostream & write_restart (std::ostream &os);
+
+  virtual void write_pmf();
+
+  class hill;
+  typedef std::list<hill>::iterator hill_iter;
+
+protected:
+
+  /// \brief width of a hill
+  ///
+  /// The local width of each collective variable, multiplied by this
+  /// number, provides the hill width along that direction
+  cvm::real  hill_width;
+
+  /// \brief Number of simulation steps between two hills
+  size_t     new_hill_freq;
+
+  /// Write the hill logfile
+  bool           b_hills_traj;
+  /// Logfile of hill management (creation and deletion)
+  std::ofstream  hills_traj_os;
+
+  /// \brief List of hills used on this bias (total); if a grid is
+  /// employed, these don't need to be updated at every time step
+  std::list<hill> hills;
+
+  /// \brief Iterator to the first of the "newest" hills (when using
+  /// grids, those who haven't been mapped yet)
+  hill_iter new_hills_begin;
+
+  /// \brief List of hills used on this bias that are on the boundary
+  /// edges; these are updated regardless of whether hills are used
+  std::list<hill> hills_off_grid;
+
+  /// \brief Same as new_hills_begin, but for the off-grid ones
+  hill_iter new_hills_off_grid_begin;
+
+  /// Regenerate the hills_off_grid list
+  void recount_hills_off_grid (hill_iter h_first, hill_iter h_last,
+                               colvar_grid_scalar *ge);
+
+  /// Read a hill from a file
+  std::istream & read_hill (std::istream &is);
+
+  /// \brief step present in a state file
+  /// 
+  /// When using grids and reading state files containing them
+  /// (multiple replicas), this is used to check whether a hill is
+  /// newer or older than the grids
+  size_t                   state_file_step;
+
+  /// \brief Add a new hill; if a .hills trajectory is written,
+  /// write it there; if there is more than one replica, communicate
+  /// it to the others
+  virtual std::list<hill>::const_iterator create_hill (hill const &h);
+
+  /// \brief Remove a previously saved hill (returns an iterator for
+  /// the next hill in the list)
+  virtual std::list<hill>::const_iterator delete_hill (hill_iter &h);
+
+  /// \brief Calculate the values of the hills, incrementing
+  /// bias_energy
+  virtual void calc_hills (hill_iter  h_first,
+                           hill_iter  h_last,
+                           cvm::real &energy,
+                           std::vector<colvarvalue> const &values = std::vector<colvarvalue> (0));
+
+  /// \brief Calculate the forces acting on the i-th colvar,
+  /// incrementing colvar_forces[i]; must be called after calc_hills
+  /// each time the values of the colvars are changed
+  virtual void calc_hills_force (size_t const &i,
+                                 hill_iter h_first,
+                                 hill_iter h_last,
+                                 std::vector<colvarvalue> &forces,
+                                 std::vector<colvarvalue> const &values = std::vector<colvarvalue> (0));
+
+
+  /// Height of new hills
+  cvm::real  hill_weight;
+
+  /// \brief Bin the hills on grids of energy and forces, and use them
+  /// to force the colvars (as opposed to deriving the hills analytically)
+  bool       use_grids;
+
+  /// \brief Rebin the hills upon restarting
+  bool       rebin_grids;
+
+  /// \brief Should the grids be expanded if necessary?
+  bool       expand_grids;
+
+  /// \brief How often the hills should be projected onto the grids
+  size_t     grids_freq;
+
+  /// \brief Whether to keep the hills in the restart file (e.g. to do
+  /// meaningful accurate rebinning afterwards)
+  bool       keep_hills;
+
+  /// \brief Dump the free energy surface (.pmf file) every restartFrequency
+  bool       dump_fes;
+
+  /// \brief Keep the free energy surface files at different
+  /// iterations, appending a step number to each
+  bool       dump_fes_save;
+
+  /// \brief Try to read the restart information by allocating new
+  /// grids before replacing the current ones (used e.g. in
+  /// multiple_replicas)
+  bool       safely_read_restart;
+
+  /// Hill energy, cached on a grid
+  colvar_grid_scalar    *hills_energy;
+
+  /// Hill forces, cached on a grid
+  colvar_grid_gradient  *hills_energy_gradients;
+
+  /// \brief Project the selected hills onto grids
+  void project_hills (hill_iter h_first, hill_iter h_last,
+                      colvar_grid_scalar *ge, colvar_grid_gradient *gf);
+
+
+  // Multiple Replicas variables and functions
+
+  /// \brief Identifier for this replica
+  std::string            replica_id;
+
+  /// \brief File containing the paths to the output files from this replica
+  std::string            replica_file_name;
+
+  /// \brief Read the existing replicas on registry
+  virtual void update_replicas_registry();
+
+  /// \brief Read new data from replicas' files
+  virtual void read_replica_files();
+
+  /// \brief Write data to other replicas
+  virtual void write_replica_state_file();
+
+  /// \brief Additional, "mirror" metadynamics biases, to collect info
+  /// from the other replicas
+  ///
+  /// These are supposed to be synchronized by reading data from the
+  /// other replicas, and not be modified by the "local" replica
+  std::vector<colvarbias_meta *> replicas;
+
+  /// \brief Frequency at which data the "mirror" biases are updated
+  size_t                 replica_update_freq;
+
+  /// List of replicas (and their output list files): contents are
+  /// copied into replicas_registry for convenience
+  std::string            replicas_registry_file;
+  /// List of replicas (and their output list files)
+  std::string            replicas_registry;
+  /// List of files written by this replica
+  std::string            replica_list_file;
+
+  /// Hills energy and gradients written specifically for other
+  /// replica (in addition to its own restart file)
+  std::string            replica_state_file;
+  /// Whether a mirror bias has read the latest version of its state file
+  bool                   replica_state_file_in_sync;
+
+  /// If there was a failure reading one of the files (because they
+  /// are not complete), this counter is incremented
+  size_t                 update_status;
+
+  /// Explicit hills communicated between replicas
+  ///
+  /// This file becomes empty after replica_state_file is rewritten
+  std::string            replica_hills_file;
+
+  /// \brief Output stream corresponding to replica_hills_file
+  std::ofstream          replica_hills_os;
+
+  /// Position within replica_hills_file (when reading it)
+  size_t                 replica_hills_file_pos;
+
+};
+
+
+
+
+/// \brief A hill for the metadynamics bias
+class colvarbias_meta::hill {
+
+protected:
+
+  /// Value of the hill function (ranges between 0 and 1)
+  cvm::real hill_value;
+
+  /// Scale factor, which could be modified at runtime (default: 1)
+  cvm::real sW;
+
+  /// Maximum height in energy of the hill
+  cvm::real W;
+
+  /// Center of the hill in the collective variable space
+  std::vector<colvarvalue>  centers;
+
+  /// Widths of the hill in the collective variable space
+  std::vector<cvm::real>    widths;
+
+public:
+
+  friend class colvarbias_meta;
+
+  /// Time step at which this hill was added
+  size_t      it;
+
+  /// Identity of the replica who added this hill (only in multiple replica simulations)
+  std::string replica;
+
+  /// \brief Runtime constructor: data are read directly from
+  /// collective variables \param weight Weight of the hill \param
+  /// cv Pointer to the array of collective variables involved \param
+  /// replica (optional) Identity of the replica which creates the
+  /// hill
+  inline hill (cvm::real             const &W_in,
+               std::vector<colvar *>       &cv,
+               cvm::real             const &hill_width,
+               std::string           const &replica_in = "")
+    : sW (1.0),
+      W (W_in),
+      centers (cv.size()),
+      widths (cv.size()),
+      it (cvm::it),
+      replica (replica_in)
+  {
+    for (size_t i = 0; i < cv.size(); i++) {
+      centers[i].type (cv[i]->type());
+      centers[i] = cv[i]->value();
+      widths[i] = cv[i]->width * hill_width;
+    }
+    if (cvm::debug()) 
+      cvm::log ("New hill, applied to "+cvm::to_str (cv.size())+
+                " collective variables, with centers "+
+                cvm::to_str (centers)+", widths "+
+                cvm::to_str (widths)+" and weight "+
+                cvm::to_str (W)+".\n");
+  }
+
+  /// \brief General constructor: all data are explicitly passed as
+  /// arguments (used for instance when reading hills saved on a
+  /// file) \param it Time step of creation of the hill \param
+  /// weight Weight of the hill \param centers Center of the hill
+  /// \param widths Width of the hill around centers \param replica
+  /// (optional) Identity of the replica which creates the hill
+  inline hill (size_t                    const &it_in,
+               cvm::real                 const &W_in,
+               std::vector<colvarvalue>  const &centers_in,
+               std::vector<cvm::real>    const &widths_in,
+               std::string               const &replica_in = "")
+    : sW (1.0),
+      W (W_in),
+      centers (centers_in),
+      widths (widths_in),
+      it (it_in),
+      replica (replica_in)
+  {}
+
+  /// Copy constructor
+  inline hill (colvarbias_meta::hill const &h)
+    : sW (1.0),
+      W (h.W),
+      centers (h.centers),
+      widths (h.widths),
+      it (h.it),
+      replica (h.replica)
+  {}
+
+  /// Destructor
+  inline ~hill()
+  {}
+  
+  /// Get the energy
+  inline cvm::real energy()
+  {
+    return W * sW * hill_value;
+  }
+
+  /// Get the energy using another hill weight
+  inline cvm::real energy (cvm::real const &new_weight)
+  {
+    return new_weight * sW * hill_value;
+  }
+
+  /// Get the current hill value
+  inline cvm::real const &value()
+  {
+    return hill_value;
+  }
+
+  /// Set the hill value as specified
+  inline void value (cvm::real const &new_value)
+  {
+    hill_value = new_value;
+  }
+
+  /// Get the weight
+  inline cvm::real weight()
+  {
+    return W * sW;
+  }
+
+  /// Scale the weight with this factor (by default 1.0 is used)
+  inline void scale (cvm::real const &new_scale_fac)
+  {
+    sW = new_scale_fac;
+  }
+
+  /// Get the center of the hill
+  inline std::vector<colvarvalue> & center()
+  {
+    return centers;
+  }
+
+  /// Get the i-th component of the center
+  inline colvarvalue & center (size_t const &i)
+  {
+    return centers[i];
+  }
+
+  /// Comparison operator
+  inline friend bool operator < (hill const &h1, hill const &h2)
+  {
+    if (h1.it < h2.it) return true;
+    else return false;
+  }
+
+  /// Comparison operator
+  inline friend bool operator <= (hill const &h1, hill const &h2)
+  {
+    if (h1.it <= h2.it) return true;
+    else return false;
+  }
+
+  /// Comparison operator
+  inline friend bool operator > (hill const &h1, hill const &h2)
+  {
+    if (h1.it > h2.it) return true;
+    else return false;
+  }
+
+  /// Comparison operator
+  inline friend bool operator >= (hill const &h1, hill const &h2)
+  {
+    if (h1.it >= h2.it) return true;
+    else return false;
+  }
+
+  /// Comparison operator
+  inline friend bool operator == (hill const &h1, hill const &h2)
+  {
+    if ( (h1.it >= h2.it) && (h1.replica == h2.replica) ) return true;
+    else return false;
+  }
+
+  /// Represent the hill ina string suitable for a trajectory file
+  std::string output_traj();
+
+  /// Write the hill to an output stream
+  inline friend std::ostream & operator << (std::ostream &os,
+                                            hill const &h);
+
+};
+
+
+#endif
+
+
+// Emacs
+// Local Variables:
+// mode: C++
+// End:

lib/colvars/colvarcomp.C

@@ -0,0 +1,91 @@
+#include "colvarmodule.h"
+#include "colvarvalue.h"
+#include "colvar.h"
+#include "colvarcomp.h"
+
+
+
+colvar::cvc::cvc()
+  : sup_coeff (1.0), sup_np (1),
+    b_periodic (false),
+    b_debug_gradients (false),
+    b_inverse_gradients (false),
+    b_Jacobian_derivative (false)
+{}
+
+
+colvar::cvc::cvc (std::string const &conf)
+  : sup_coeff (1.0), sup_np (1),
+    b_periodic (false),
+    b_debug_gradients (false),
+    b_inverse_gradients (false),
+    b_Jacobian_derivative (false)
+{
+  if (cvm::debug())
+    cvm::log ("Initializing cvc base object.\n");
+
+  get_keyval (conf, "name", this->name, std::string (""), parse_silent);
+
+  get_keyval (conf, "componentCoeff", sup_coeff, 1.0);
+  get_keyval (conf, "componentExp", sup_np, 1);
+
+  get_keyval (conf, "period", period, 0.0);
+  get_keyval (conf, "wrapAround", wrap_center, 0.0);
+
+  get_keyval (conf, "debugGradients", b_debug_gradients, false, parse_silent);
+
+  if (cvm::debug())
+    cvm::log ("Done initializing cvc base object.\n");
+}
+
+
+void colvar::cvc::parse_group (std::string const &conf, 
+                               char const *group_key,
+                               cvm::atom_group &group,
+                               bool optional)
+{
+  if (key_lookup (conf, group_key)) {
+    group.parse (conf, group_key);
+  } else {
+    if (! optional) {
+      cvm::fatal_error ("Error: definition for atom group \""+
+                      std::string (group_key)+"\" not found.\n");
+    }
+  }
+}
+
+
+colvar::cvc::~cvc()
+{}
+
+
+void colvar::cvc::calc_force_invgrads()
+{
+  cvm::fatal_error ("Error: calculation of inverse gradients is not implemented "
+                    "for colvar components of type \""+function_type+"\".\n");
+}
+
+
+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");
+}
+
+
+colvarvalue colvar::cvc::fdiff_change (cvm::atom_group &group)
+{  
+  colvarvalue change (x.type());
+
+  if (group.old_pos.size()) {
+    for (size_t i = 0; i < group.size(); i++) {
+      cvm::rvector const &pold = group.old_pos[i];
+      cvm::rvector const &p = group[i].pos;
+      change += group[i].grad * (p - pold);
+    }
+  }
+
+  // save for next step
+  group.old_pos = group.positions();
+  return change;
+}

lib/colvars/colvarcomp_angles.C

@@ -0,0 +1,377 @@
+
+#include "colvarmodule.h"
+#include "colvar.h"
+#include "colvarcomp.h"
+
+#include <cmath>
+
+
+
+colvar::angle::angle (std::string const &conf)
+  : cvc (conf)
+{
+  function_type = "angle";
+  b_inverse_gradients = true;
+  b_Jacobian_derivative = true;
+  parse_group (conf, "group1", group1);
+  parse_group (conf, "group2", group2);
+  parse_group (conf, "group3", group3);
+  atom_groups.push_back (&group1);
+  atom_groups.push_back (&group2);
+  atom_groups.push_back (&group3);
+  if (get_keyval (conf, "oneSiteSystemForce", b_1site_force, false)) {
+    cvm::log ("Computing system force on group 1 only");
+  }
+  x.type (colvarvalue::type_scalar);
+}
+
+
+colvar::angle::angle (cvm::atom const &a1,
+                      cvm::atom const &a2,
+                      cvm::atom const &a3)
+  : group1 (std::vector<cvm::atom> (1, a1)),
+    group2 (std::vector<cvm::atom> (1, a2)),
+    group3 (std::vector<cvm::atom> (1, a3))
+{
+  function_type = "angle";
+  b_inverse_gradients = true;
+  b_Jacobian_derivative = true;
+  b_1site_force = false;
+  atom_groups.push_back (&group1);
+  atom_groups.push_back (&group2);
+  atom_groups.push_back (&group3);
+
+  x.type (colvarvalue::type_scalar);
+}
+
+
+colvar::angle::angle()
+{
+  function_type = "angle";
+  x.type (colvarvalue::type_scalar);
+}
+
+
+void colvar::angle::calc_value()
+{
+  group1.read_positions();
+  group2.read_positions();
+  group3.read_positions();
+
+  cvm::atom_pos const g1_pos = group1.center_of_mass();
+  cvm::atom_pos const g2_pos = group2.center_of_mass();
+  cvm::atom_pos const g3_pos = group3.center_of_mass();
+
+  r21  = cvm::position_distance (g2_pos, g1_pos);
+  r21l = r21.norm();
+  r23  = cvm::position_distance (g2_pos, g3_pos);
+  r23l = r23.norm();
+
+  cvm::real     const cos_theta = (r21*r23)/(r21l*r23l);
+
+  x.real_value = (180.0/PI) * std::acos (cos_theta);
+}
+
+
+void colvar::angle::calc_gradients()
+{
+  cvm::real const cos_theta = (r21*r23)/(r21l*r23l);
+  cvm::real const dxdcos = -1.0 / std::sqrt (1.0 - cos_theta*cos_theta);
+    
+  dxdr1 = (180.0/PI) * dxdcos *
+    (1.0/r21l) * ( r23/r23l + (-1.0) * cos_theta * r21/r21l );
+
+  dxdr3 = (180.0/PI) * dxdcos *
+    (1.0/r23l) * ( r21/r21l + (-1.0) * cos_theta * r23/r23l );
+
+  for (size_t i = 0; i < group1.size(); i++) {
+    group1[i].grad = (group1[i].mass/group1.total_mass) *
+      (dxdr1);
+  }
+
+  for (size_t i = 0; i < group2.size(); i++) {
+    group2[i].grad = (group2[i].mass/group2.total_mass) *
+      (dxdr1 + dxdr3) * (-1.0);
+  }
+
+  for (size_t i = 0; i < group3.size(); i++) {
+    group3[i].grad = (group3[i].mass/group3.total_mass) *
+      (dxdr3);
+  }
+}
+
+void colvar::angle::calc_force_invgrads()
+{
+  // This uses a force measurement on groups 1 and 3 only
+  // to keep in line with the implicit variable change used to
+  // evaluate the Jacobian term (essentially polar coordinates
+  // centered on group2, which means group2 is kept fixed
+  // when propagating changes in the angle)
+
+  if (b_1site_force) {
+    group1.read_system_forces();
+    cvm::real norm_fact = 1.0 / dxdr1.norm2();
+    ft.real_value = norm_fact * dxdr1 * group1.system_force();
+  } else {
+    group1.read_system_forces();
+    group3.read_system_forces();
+    cvm::real norm_fact = 1.0 / (dxdr1.norm2() + dxdr3.norm2());
+    ft.real_value = norm_fact * ( dxdr1 * group1.system_force()
+                                + dxdr3 * group3.system_force());
+  }
+  return;
+}
+
+void colvar::angle::calc_Jacobian_derivative()
+{
+  // det(J) = (2 pi) r^2 * sin(theta)
+  // hence Jd = cot(theta)
+  const cvm::real theta = x.real_value * PI / 180.0;
+  jd = PI / 180.0 * (theta != 0.0 ? std::cos(theta) / std::sin(theta) : 0.0);
+}
+
+
+void colvar::angle::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);
+
+  if (!group3.noforce)
+    group3.apply_colvar_force (force.real_value);
+}
+
+
+
+
+colvar::dihedral::dihedral (std::string const &conf)
+  : cvc (conf)
+{
+  function_type = "dihedral";
+  period = 360.0;
+  b_periodic = true;
+  b_inverse_gradients = true;
+  b_Jacobian_derivative = true;
+  if (get_keyval (conf, "oneSiteSystemForce", b_1site_force, false)) {
+    cvm::log ("Computing system force on group 1 only");
+  }
+  parse_group (conf, "group1", group1);
+  parse_group (conf, "group2", group2);
+  parse_group (conf, "group3", group3);
+  parse_group (conf, "group4", group4);
+  atom_groups.push_back (&group1);
+  atom_groups.push_back (&group2);
+  atom_groups.push_back (&group3);
+  atom_groups.push_back (&group4);
+
+  x.type (colvarvalue::type_scalar);
+}
+
+
+colvar::dihedral::dihedral (cvm::atom const &a1,
+                            cvm::atom const &a2,
+                            cvm::atom const &a3,
+                            cvm::atom const &a4)
+  : group1 (std::vector<cvm::atom> (1, a1)),
+    group2 (std::vector<cvm::atom> (1, a2)),
+    group3 (std::vector<cvm::atom> (1, a3)),
+    group4 (std::vector<cvm::atom> (1, a4))
+{
+  if (cvm::debug())
+    cvm::log ("Initializing dihedral object from atom groups.\n");
+
+  function_type = "dihedral";
+  period = 360.0;
+  b_periodic = true;
+  b_inverse_gradients = true;
+  b_Jacobian_derivative = true;
+  b_1site_force = false;
+
+  atom_groups.push_back (&group1);
+  atom_groups.push_back (&group2);
+  atom_groups.push_back (&group3);
+  atom_groups.push_back (&group4);
+
+  x.type (colvarvalue::type_scalar);
+
+  if (cvm::debug())
+    cvm::log ("Done initializing dihedral object from atom groups.\n");
+}
+
+
+colvar::dihedral::dihedral()
+{
+  function_type = "dihedral";
+  period = 360.0;
+  b_periodic = true;
+  b_inverse_gradients = true;
+  b_Jacobian_derivative = true;
+  x.type (colvarvalue::type_scalar);
+}
+
+
+void colvar::dihedral::calc_value()
+{
+  group1.read_positions();
+  group2.read_positions();
+  group3.read_positions();
+  group4.read_positions();
+
+  cvm::atom_pos const g1_pos = group1.center_of_mass();
+  cvm::atom_pos const g2_pos = group2.center_of_mass();
+  cvm::atom_pos const g3_pos = group3.center_of_mass();
+  cvm::atom_pos const g4_pos = group4.center_of_mass();
+
+  // Usual sign convention: r12 = r2 - r1
+  r12 = cvm::position_distance (g1_pos, g2_pos);
+  r23 = cvm::position_distance (g2_pos, g3_pos);
+  r34 = cvm::position_distance (g3_pos, g4_pos);
+
+  cvm::rvector const n1 = cvm::rvector::outer (r12, r23);
+  cvm::rvector const n2 = cvm::rvector::outer (r23, r34);
+
+  cvm::real const cos_phi = n1 * n2;
+  cvm::real const sin_phi = n1 * r34 * r23.norm();
+
+  x.real_value = (180.0/PI) * std::atan2 (sin_phi, cos_phi);
+  this->wrap (x);
+}
+
+
+void colvar::dihedral::calc_gradients()
+{
+  cvm::rvector A = cvm::rvector::outer (r12, r23);
+  cvm::real   rA = A.norm();
+  cvm::rvector B = cvm::rvector::outer (r23, r34);
+  cvm::real   rB = B.norm();
+  cvm::rvector C = cvm::rvector::outer (r23, A); 
+  cvm::real   rC = C.norm();
+
+  cvm::real const cos_phi = (A*B)/(rA*rB);
+  cvm::real const sin_phi = (C*B)/(rC*rB);
+
+  cvm::rvector f1, f2, f3;
+
+  rB = 1.0/rB;
+  B *= rB;
+
+  if (std::fabs (sin_phi) > 0.1) {
+    rA = 1.0/rA;
+    A *= rA;
+    cvm::rvector const dcosdA = rA*(cos_phi*A-B);
+    cvm::rvector const dcosdB = rB*(cos_phi*B-A);
+    rA = 1.0;
+
+    cvm::real const K = (1.0/sin_phi) * (180.0/PI);
+
+	f1 = K * cvm::rvector::outer (r23, dcosdA);
+	f3 = K * cvm::rvector::outer (dcosdB, r23);
+	f2 = K * (cvm::rvector::outer (dcosdA, r12)
+		   +  cvm::rvector::outer (r34, dcosdB));
+  }
+  else {
+    rC = 1.0/rC;
+    C *= rC;
+    cvm::rvector const dsindC = rC*(sin_phi*C-B);
+    cvm::rvector const dsindB = rB*(sin_phi*B-C);
+    rC = 1.0;
+
+    cvm::real    const K = (-1.0/cos_phi) * (180.0/PI);
+
+    f1.x = K*((r23.y*r23.y + r23.z*r23.z)*dsindC.x
+              - r23.x*r23.y*dsindC.y
+              - r23.x*r23.z*dsindC.z);
+    f1.y = K*((r23.z*r23.z + r23.x*r23.x)*dsindC.y
+              - r23.y*r23.z*dsindC.z
+              - r23.y*r23.x*dsindC.x);
+    f1.z = K*((r23.x*r23.x + r23.y*r23.y)*dsindC.z
+              - r23.z*r23.x*dsindC.x
+              - r23.z*r23.y*dsindC.y);
+
+    f3 = cvm::rvector::outer (dsindB, r23);
+    f3 *= K;
+
+    f2.x = K*(-(r23.y*r12.y + r23.z*r12.z)*dsindC.x
+              +(2.0*r23.x*r12.y - r12.x*r23.y)*dsindC.y
+              +(2.0*r23.x*r12.z - r12.x*r23.z)*dsindC.z
+              +dsindB.z*r34.y - dsindB.y*r34.z);
+    f2.y = K*(-(r23.z*r12.z + r23.x*r12.x)*dsindC.y
+              +(2.0*r23.y*r12.z - r12.y*r23.z)*dsindC.z
+              +(2.0*r23.y*r12.x - r12.y*r23.x)*dsindC.x
+              +dsindB.x*r34.z - dsindB.z*r34.x);
+    f2.z = K*(-(r23.x*r12.x + r23.y*r12.y)*dsindC.z
+              +(2.0*r23.z*r12.x - r12.z*r23.x)*dsindC.x
+              +(2.0*r23.z*r12.y - r12.z*r23.y)*dsindC.y
+              +dsindB.y*r34.x - dsindB.x*r34.y);
+  }
+
+  for (size_t i = 0; i < group1.size(); i++)
+    group1[i].grad = (group1[i].mass/group1.total_mass) * (-f1);
+
+  for (size_t i = 0; i < group2.size(); i++)
+    group2[i].grad = (group2[i].mass/group2.total_mass) * (-f2 + f1);
+
+  for (size_t i = 0; i < group3.size(); i++)
+	group3[i].grad = (group3[i].mass/group3.total_mass) * (-f3 + f2);
+
+  for (size_t i = 0; i < group4.size(); i++)
+    group4[i].grad = (group4[i].mass/group4.total_mass) * (f3);
+}
+
+
+void colvar::dihedral::calc_force_invgrads()
+{
+  cvm::rvector const u12 = r12.unit();
+  cvm::rvector const u23 = r23.unit();
+  cvm::rvector const u34 = r34.unit();
+
+  cvm::real const d12 = r12.norm();
+  cvm::real const d34 = r34.norm();
+
+  cvm::rvector const cross1 = (cvm::rvector::outer (u23, u12)).unit();
+  cvm::rvector const cross4 = (cvm::rvector::outer (u23, u34)).unit();
+
+  cvm::real const dot1 = u23 * u12;
+  cvm::real const dot4 = u23 * u34;
+
+  cvm::real const fact1 = d12 * std::sqrt (1.0 - dot1 * dot1);
+  cvm::real const fact4 = d34 * std::sqrt (1.0 - dot4 * dot4);
+
+  group1.read_system_forces();
+  if ( b_1site_force ) {
+    // This is only measuring the force on group 1
+    ft.real_value = PI/180.0 * fact1 * (cross1 * group1.system_force());
+  } else {
+    // Default case: use groups 1 and 4
+    group4.read_system_forces();
+    ft.real_value = PI/180.0 * 0.5 * (fact1 * (cross1 * group1.system_force())
+				      + fact4 * (cross4 * group4.system_force()));
+  }
+}
+
+
+void colvar::dihedral::calc_Jacobian_derivative()
+{
+  // With this choice of inverse gradient ("internal coordinates"), Jacobian correction is 0
+  jd = 0.0;
+}
+
+
+void colvar::dihedral::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);
+
+  if (!group3.noforce)
+    group3.apply_colvar_force (force.real_value);
+
+  if (!group4.noforce)
+    group4.apply_colvar_force (force.real_value);
+}
+
+

lib/colvars/colvarcomp_coordnums.C

@@ -0,0 +1,404 @@
+#include <cmath>
+
+#include "colvarmodule.h"
+#include "colvarparse.h"
+#include "colvaratoms.h"
+#include "colvarvalue.h"
+#include "colvar.h"
+#include "colvarcomp.h"
+
+
+
+template<bool calculate_gradients>
+cvm::real colvar::coordnum::switching_function (cvm::real const &r0,
+                                                int const &en,
+                                                int const &ed,
+                                                cvm::atom &A1,
+                                                cvm::atom &A2)
+{
+  cvm::rvector const diff = cvm::position_distance (A1.pos, A2.pos);
+  cvm::real const l2 = diff.norm2()/(r0*r0);
+
+  // Assume en and ed are even integers, and avoid sqrt in the following
+  int const en2 = en/2;
+  int const ed2 = ed/2;
+
+  cvm::real const xn = std::pow (l2, en2);
+  cvm::real const xd = std::pow (l2, ed2);
+  cvm::real const func = (1.0-xn)/(1.0-xd);
+
+  if (calculate_gradients) {
+    cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) - func*ed2*(xd/l2))*(-1.0);
+    cvm::rvector const dl2dx = (2.0/(r0*r0))*diff;
+    A1.grad += (-1.0)*dFdl2*dl2dx;
+    A2.grad +=        dFdl2*dl2dx;
+  }
+
+  return func;
+}
+
+
+template<bool calculate_gradients>
+cvm::real colvar::coordnum::switching_function (cvm::rvector const &r0_vec,
+                                                int const &en,
+                                                int const &ed,
+                                                cvm::atom &A1,
+                                                cvm::atom &A2)
+{
+  cvm::rvector const diff = cvm::position_distance (A1.pos, A2.pos);
+  cvm::rvector const scal_diff (diff.x/r0_vec.x, diff.y/r0_vec.y, diff.z/r0_vec.z);
+  cvm::real const l2 = scal_diff.norm2();
+
+  // Assume en and ed are even integers, and avoid sqrt in the following
+  int const en2 = en/2;
+  int const ed2 = ed/2;
+
+  cvm::real const xn = std::pow (l2, en2);
+  cvm::real const xd = std::pow (l2, ed2);
+  cvm::real const func = (1.0-xn)/(1.0-xd);
+
+  if (calculate_gradients) {
+    cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) - func*ed2*(xd/l2))*(-1.0);
+    cvm::rvector const dl2dx ((2.0/(r0_vec.x*r0_vec.x))*diff.x,
+                              (2.0/(r0_vec.y*r0_vec.y))*diff.y,
+                              (2.0/(r0_vec.z*r0_vec.z))*diff.z);
+    A1.grad += (-1.0)*dFdl2*dl2dx;
+    A2.grad +=        dFdl2*dl2dx;
+  }
+  return func;
+}
+
+
+
+colvar::coordnum::coordnum (std::string const &conf)
+  : distance (conf), b_anisotropic (false), b_group2_center_only (false)
+{ 
+  function_type = "coordnum";
+  x.type (colvarvalue::type_scalar);
+
+  // group1 and group2 are already initialized by distance()
+
+  // need to specify this explicitly because the distance() constructor
+  // has set it to true
+  b_inverse_gradients = false;
+
+  bool const b_scale = get_keyval (conf, "cutoff", r0,
+                                   cvm::real (4.0 * cvm::unit_angstrom()));
+
+  if (get_keyval (conf, "cutoff3", r0_vec,
+                  cvm::rvector (4.0, 4.0, 4.0), parse_silent)) {
+
+    if (b_scale)
+      cvm::fatal_error ("Error: cannot specify \"scale\" and "
+                        "\"scale3\" at the same time.\n");
+    b_anisotropic = true;
+    // remove meaningless negative signs
+    if (r0_vec.x < 0.0) r0_vec.x *= -1.0;
+    if (r0_vec.y < 0.0) r0_vec.y *= -1.0;
+    if (r0_vec.z < 0.0) r0_vec.z *= -1.0;
+  }
+
+  get_keyval (conf, "expNumer", en, int (6) );
+  get_keyval (conf, "expDenom", ed, int (12));
+
+  if ( (en%2) || (ed%2) ) {
+    cvm::fatal_error ("Error: odd exponents provided, can only use even ones.\n");
+  }
+
+  get_keyval (conf, "group2CenterOnly", b_group2_center_only, false);
+}
+
+
+colvar::coordnum::coordnum()
+  : b_anisotropic (false), b_group2_center_only (false)
+{
+  function_type = "coordnum";
+  x.type (colvarvalue::type_scalar);
+}
+
+
+void colvar::coordnum::calc_value()
+{
+  x.real_value = 0.0;
+
+  // these are necessary: for each atom, gradients are summed together
+  // by multiple calls to switching_function()
+  group1.reset_atoms_data();
+  group2.reset_atoms_data();
+
+  group1.read_positions();
+  group2.read_positions();
+
+  if (b_group2_center_only) {
+
+    // create a fake atom to hold the group2 com coordinates
+    cvm::atom group2_com_atom (group2[0]);
+    group2_com_atom.pos = group2.center_of_mass();
+
+    if (b_anisotropic) {
+      for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)
+        x.real_value += switching_function<false> (r0_vec, en, ed, *ai1, group2_com_atom);
+    } else {
+      for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)
+        x.real_value += switching_function<false> (r0, en, ed, *ai1, group2_com_atom);
+    }
+
+  } else {
+
+    if (b_anisotropic) {
+      for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)
+        for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
+          x.real_value += switching_function<false> (r0_vec, en, ed, *ai1, *ai2);
+        }
+    } else {
+      for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)
+        for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
+          x.real_value += switching_function<false> (r0, en, ed, *ai1, *ai2);
+        }
+    }
+  }
+}
+
+
+void colvar::coordnum::calc_gradients()
+{
+  if (b_group2_center_only) {
+
+    // create a fake atom to hold the group2 com coordinates
+    cvm::atom group2_com_atom (group2[0]);
+    group2_com_atom.pos = group2.center_of_mass();
+
+    if (b_anisotropic) {
+      for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)
+        switching_function<true> (r0_vec, en, ed, *ai1, group2_com_atom);
+    } else {
+      for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)
+        switching_function<true> (r0, en, ed, *ai1, group2_com_atom);
+    }
+
+    group2.set_weighted_gradient (group2_com_atom.grad);
+
+  } else {
+
+    if (b_anisotropic) {
+      for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)
+        for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
+          switching_function<true> (r0_vec, en, ed, *ai1, *ai2);
+        }
+    } else {
+      for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)
+        for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
+          switching_function<true> (r0, en, ed, *ai1, *ai2);
+        }
+    }
+  }
+
+  //   if (cvm::debug()) {
+  //     for (size_t i = 0; i < group1.size(); i++) {
+  //       cvm::log ("atom["+cvm::to_str (group1[i].id+1)+"] gradient: "+
+  //                 cvm::to_str (group1[i].grad)+"\n");
+  //     }
+
+  //     for (size_t i = 0; i < group2.size(); i++) {
+  //       cvm::log ("atom["+cvm::to_str (group2[i].id+1)+"] gradient: "+
+  //                 cvm::to_str (group2[i].grad)+"\n");
+  //     }
+  //   }
+}
+
+void colvar::coordnum::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);
+}
+
+
+
+// h_bond member functions
+
+colvar::h_bond::h_bond (std::string const &conf)
+  : cvc (conf)
+{
+  if (cvm::debug())
+    cvm::log ("Initializing h_bond object.\n");
+
+  function_type = "h_bond";
+  x.type (colvarvalue::type_scalar);
+
+  int a_num, d_num;
+  get_keyval (conf, "acceptor", a_num, -1);
+  get_keyval (conf, "donor",    d_num, -1);
+
+  if ( (a_num == -1) || (d_num == -1) ) {
+    cvm::fatal_error ("Error: either acceptor or donor undefined.\n");
+  }
+
+  acceptor = cvm::atom (a_num);
+  donor    = cvm::atom (d_num);
+  atom_groups.push_back (new cvm::atom_group);
+  atom_groups[0]->add_atom (acceptor);
+  atom_groups[0]->add_atom (donor);
+
+  get_keyval (conf, "cutoff",   r0, (3.3 * cvm::unit_angstrom()));
+  get_keyval (conf, "expNumer", en, 6);
+  get_keyval (conf, "expDenom", ed, 8);
+
+  if ( (en%2) || (ed%2) ) {
+    cvm::fatal_error ("Error: odd exponents provided, can only use even ones.\n");
+  }
+
+  if (cvm::debug())
+    cvm::log ("Done initializing h_bond object.\n");
+}
+
+
+colvar::h_bond::h_bond (cvm::atom const &acceptor_i,
+                        cvm::atom const &donor_i,
+                        cvm::real r0_i, int en_i, int ed_i)
+  : acceptor (acceptor_i),
+    donor (donor_i),
+    r0 (r0_i), en (en_i), ed (ed_i)
+{
+  function_type = "h_bond";
+  x.type (colvarvalue::type_scalar);
+
+  atom_groups.push_back (new cvm::atom_group);
+  atom_groups[0]->add_atom (acceptor);
+  atom_groups[0]->add_atom (donor);
+}
+
+colvar::h_bond::h_bond()
+  : cvc ()
+{
+  function_type = "h_bond";
+  x.type (colvarvalue::type_scalar);
+}
+
+
+colvar::h_bond::~h_bond()
+{
+  for (int i=0; i<atom_groups.size(); i++) {
+    delete atom_groups[i];
+  }
+}
+
+
+void colvar::h_bond::calc_value()
+{
+  // this is necessary, because switching_function() will sum the new
+  // gradient to the current one
+  acceptor.reset_data();
+  donor.reset_data();
+
+  acceptor.read_position();
+  donor.read_position();
+
+  x.real_value = colvar::coordnum::switching_function<false> (r0, en, ed, acceptor, donor);
+}
+
+
+void colvar::h_bond::calc_gradients()
+{
+  colvar::coordnum::switching_function<true> (r0, en, ed, acceptor, donor);
+  (*atom_groups[0])[0].grad = acceptor.grad;
+  (*atom_groups[0])[1].grad = donor.grad;
+}
+
+
+void colvar::h_bond::apply_force (colvarvalue const &force)
+{
+  acceptor.apply_force (force.real_value * acceptor.grad);
+  donor.apply_force    (force.real_value * donor.grad);
+}
+
+
+// Self-coordination number for a group
+
+template<bool calculate_gradients>
+cvm::real colvar::selfcoordnum::switching_function (cvm::real const &r0,
+                                                int const &en,
+                                                int const &ed,
+                                                cvm::atom &A1,
+                                                cvm::atom &A2)
+{
+  cvm::rvector const diff = cvm::position_distance (A1.pos, A2.pos);
+  cvm::real const l2 = diff.norm2()/(r0*r0);
+
+  // Assume en and ed are even integers, and avoid sqrt in the following
+  int const en2 = en/2;
+  int const ed2 = ed/2;
+
+  cvm::real const xn = std::pow (l2, en2);
+  cvm::real const xd = std::pow (l2, ed2);
+  cvm::real const func = (1.0-xn)/(1.0-xd);
+
+  if (calculate_gradients) {
+    cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) - func*ed2*(xd/l2))*(-1.0);
+    cvm::rvector const dl2dx = (2.0/(r0*r0))*diff;
+    A1.grad += (-1.0)*dFdl2*dl2dx;
+    A2.grad +=        dFdl2*dl2dx;
+  }
+
+  return func;
+}
+
+
+colvar::selfcoordnum::selfcoordnum (std::string const &conf)
+ : distance (conf, false)
+{ 
+  function_type = "selfcoordnum";
+  x.type (colvarvalue::type_scalar);
+
+  // group1 is already initialized by distance()
+
+  // need to specify this explicitly because the distance() constructor
+  // has set it to true
+  b_inverse_gradients = false;
+
+  get_keyval (conf, "cutoff", r0, cvm::real (4.0 * cvm::unit_angstrom()));
+  get_keyval (conf, "expNumer", en, int (6) );
+  get_keyval (conf, "expDenom", ed, int (12));
+
+  if ( (en%2) || (ed%2) ) {
+    cvm::fatal_error ("Error: odd exponents provided, can only use even ones.\n");
+  }
+}
+
+
+colvar::selfcoordnum::selfcoordnum()
+{
+  function_type = "selfcoordnum";
+  x.type (colvarvalue::type_scalar);
+}
+
+
+void colvar::selfcoordnum::calc_value()
+{
+  x.real_value = 0.0;
+
+  // for each atom, gradients are summed by multiple calls to switching_function()
+  group1.reset_atoms_data();
+  group1.read_positions();
+
+  for (size_t i = 0; i < group1.size() - 1; i++)
+    for (size_t j = i + 1; j < group1.size(); j++)
+      x.real_value += switching_function<false> (r0, en, ed, group1[i], group1[j]);
+}
+
+
+void colvar::selfcoordnum::calc_gradients()
+{
+  for (size_t i = 0; i < group1.size() - 1; i++)
+    for (size_t j = i + 1; j < group1.size(); j++)
+      switching_function<true> (r0, en, ed, group1[i], group1[j]);
+}
+
+void colvar::selfcoordnum::apply_force (colvarvalue const &force)
+{
+  if (!group1.noforce)
+    group1.apply_colvar_force (force.real_value);
+}
+

lib/colvars/colvarcomp_protein.C

@@ -0,0 +1,385 @@
+#include <cmath>
+
+#include "colvarmodule.h"
+#include "colvarvalue.h"
+#include "colvarparse.h"
+#include "colvar.h"
+#include "colvarcomp.h"
+
+
+//////////////////////////////////////////////////////////////////////
+// alpha component
+//////////////////////////////////////////////////////////////////////
+
+colvar::alpha_angles::alpha_angles (std::string const &conf)
+  : cvc (conf)
+{
+  if (cvm::debug())
+    cvm::log ("Initializing alpha_angles object.\n");
+
+  function_type = "alpha_angles";
+  x.type (colvarvalue::type_scalar);
+
+  std::string segment_id;
+  get_keyval (conf, "psfSegID", segment_id, std::string ("MAIN"));
+
+  std::vector<int> residues;
+  {
+    std::string residues_conf = "";
+    key_lookup (conf, "residueRange", residues_conf);
+    if (residues_conf.size()) {
+      std::istringstream is (residues_conf);
+      int initial, final;
+      char dash;
+      if ( (is >> initial) && (initial > 0) &&
+           (is >> dash) && (dash == '-') &&
+           (is >> final) && (final > 0) ) {
+        for (int rnum = initial; rnum <= final; rnum++) {
+          residues.push_back (rnum);
+        }
+      }
+    } else {
+      cvm::fatal_error ("Error: no residues defined in \"residueRange\".\n");
+    }
+  }
+
+  if (residues.size() < 5) {
+    cvm::fatal_error ("Error: not enough residues defined in \"residueRange\".\n");
+  }
+
+  std::string const &sid    = segment_id;
+  std::vector<int> const &r = residues;
+
+
+  get_keyval (conf, "hBondCoeff", hb_coeff, 0.5);
+  if ( (hb_coeff < 0.0) || (hb_coeff > 1.0) ) {
+    cvm::fatal_error ("Error: hBondCoeff must be defined between 0 and 1.\n");
+  }
+
+
+  get_keyval (conf, "angleRef", theta_ref, 88.0);
+  get_keyval (conf, "angleTol", theta_tol, 15.0);
+
+  if (hb_coeff < 1.0) {
+
+    for (size_t i = 0; i < residues.size()-2; i++) {
+      theta.push_back (new colvar::angle (cvm::atom (r[i  ], "CA", sid),
+                                          cvm::atom (r[i+1], "CA", sid),
+                                          cvm::atom (r[i+2], "CA", sid)));
+    }
+
+  } else {
+    cvm::log ("The hBondCoeff specified will disable the Calpha-Calpha-Calpha angle terms.\n");
+  }
+
+  {
+    cvm::real r0;
+    size_t en, ed;
+    get_keyval (conf, "hBondCutoff",   r0, (3.3 * cvm::unit_angstrom()));
+    get_keyval (conf, "hBondExpNumer", en, 6);
+    get_keyval (conf, "hBondExpDenom", ed, 8);
+
+    if (hb_coeff > 0.0) {
+
+      for (size_t i = 0; i < residues.size()-4; i++) {
+        hb.push_back (new colvar::h_bond (cvm::atom (r[i  ], "O",  sid),
+                                          cvm::atom (r[i+4], "N",  sid),
+                                          r0, en, ed));
+      }
+
+    } else {
+      cvm::log ("The hBondCoeff specified will disable the hydrogen bond terms.\n");
+    }
+  }
+
+  if (cvm::debug())
+    cvm::log ("Done initializing alpha_angles object.\n");
+}
+
+
+colvar::alpha_angles::alpha_angles()
+  : cvc ()
+{
+  function_type = "alpha_angles";
+  x.type (colvarvalue::type_scalar);
+}
+
+colvar::alpha_angles::~alpha_angles()
+{
+  while (theta.size() != 0) {
+    delete theta.back();
+    theta.pop_back();
+  }
+  while (hb.size() != 0) {
+    delete hb.back();
+    hb.pop_back();
+  }
+}
+
+void colvar::alpha_angles::calc_value()
+{
+  x.real_value = 0.0;
+
+  if (theta.size()) {
+
+    cvm::real const theta_norm = 
+      (1.0-hb_coeff) / cvm::real (theta.size());
+
+    for (size_t i = 0; i < theta.size(); i++) {
+
+      (theta[i])->calc_value();
+
+      cvm::real const t = ((theta[i])->value().real_value-theta_ref)/theta_tol;
+      cvm::real const f = ( (1.0 - std::pow (t, (int) 2)) /
+                            (1.0 - std::pow (t, (int) 4)) );
+
+      x.real_value += theta_norm * f;
+
+      if (cvm::debug())
+        cvm::log ("Calpha-Calpha angle no. "+cvm::to_str (i+1)+" in \""+
+                  this->name+"\" has a value of "+
+                  (cvm::to_str ((theta[i])->value().real_value))+
+                  " degrees, f = "+cvm::to_str (f)+".\n");
+    }
+  }
+
+  if (hb.size()) {
+
+    cvm::real const hb_norm =
+      hb_coeff / cvm::real (hb.size());
+
+    for (size_t i = 0; i < hb.size(); i++) {
+      (hb[i])->calc_value();
+      x.real_value += hb_norm * (hb[i])->value().real_value;
+      if (cvm::debug())
+        cvm::log ("Hydrogen bond no. "+cvm::to_str (i+1)+" in \""+
+                  this->name+"\" has a value of "+
+                  (cvm::to_str ((hb[i])->value().real_value))+".\n");
+    }
+  }
+}
+
+
+void colvar::alpha_angles::calc_gradients()
+{
+  for (size_t i = 0; i < theta.size(); i++) 
+    (theta[i])->calc_gradients();
+
+  for (size_t i = 0; i < hb.size(); i++)
+    (hb[i])->calc_gradients();
+}
+
+
+void colvar::alpha_angles::apply_force (colvarvalue const &force)
+{
+
+  if (theta.size()) {
+
+    cvm::real const theta_norm = 
+      (1.0-hb_coeff) / cvm::real (theta.size());
+    
+    for (size_t i = 0; i < theta.size(); i++) {
+
+      cvm::real const t = ((theta[i])->value().real_value-theta_ref)/theta_tol;
+      cvm::real const f = ( (1.0 - std::pow (t, (int) 2)) /
+                            (1.0 - std::pow (t, (int) 4)) );
+
+      cvm::real const dfdt =
+        1.0/(1.0 - std::pow (t, (int) 4)) * 
+        ( (-2.0 * t) + (-1.0*f)*(-4.0 * std::pow (t, (int) 3)) );
+
+      (theta[i])->apply_force (theta_norm * 
+                               dfdt * (1.0/theta_tol) *
+                               force.real_value );
+    }
+  }
+
+  if (hb.size()) {
+
+    cvm::real const hb_norm =
+      hb_coeff / cvm::real (hb.size());
+
+    for (size_t i = 0; i < hb.size(); i++) {
+      (hb[i])->apply_force (0.5 * hb_norm * force.real_value);
+    }
+  }
+}
+
+
+//////////////////////////////////////////////////////////////////////
+// dihedral principal component
+//////////////////////////////////////////////////////////////////////
+
+    // FIXME: this will not make collect_gradients work
+    // because gradients in individual atom groups
+    // are those of the sub-cvcs (angle, hb), not those
+    // of this cvc (alpha)
+    // This is true of all cvcs with sub-cvcs, and those
+    // that do not calculate explicit gradients
+    // SO: we need a flag giving the availability of 
+    // atomic gradients
+colvar::dihedPC::dihedPC (std::string const &conf)
+  : cvc (conf)
+{
+  if (cvm::debug())
+    cvm::log ("Initializing dihedral PC object.\n");
+
+  function_type = "dihedPC";
+  x.type (colvarvalue::type_scalar);
+
+  std::string segment_id;
+  get_keyval (conf, "psfSegID", segment_id, std::string ("MAIN"));
+
+  std::vector<int> residues;
+  {
+    std::string residues_conf = "";
+    key_lookup (conf, "residueRange", residues_conf);
+    if (residues_conf.size()) {
+      std::istringstream is (residues_conf);
+      int initial, final;
+      char dash;
+      if ( (is >> initial) && (initial > 0) &&
+           (is >> dash) && (dash == '-') &&
+           (is >> final) && (final > 0) ) {
+        for (int rnum = initial; rnum <= final; rnum++) {
+          residues.push_back (rnum);
+        }
+      }
+    } else {
+      cvm::fatal_error ("Error: no residues defined in \"residueRange\".\n");
+    }
+  }
+
+  if (residues.size() < 2) {
+    cvm::fatal_error ("Error: dihedralPC requires at least two residues.\n");
+  }
+
+  std::string const &sid    = segment_id;
+  std::vector<int> const &r = residues;
+
+  std::string vecFileName;
+  int         vecNumber;
+  if (get_keyval (conf, "vectorFile", vecFileName, vecFileName)) {
+    get_keyval (conf, "vectorNumber", vecNumber, 0); 
+    if (vecNumber < 1)
+      cvm::fatal_error ("A positive value of vectorNumber is required.");
+
+    std::ifstream vecFile;
+    vecFile.open (vecFileName.c_str());
+    if (!vecFile.good())
+      cvm::fatal_error ("Error opening dihedral PCA vector file " + vecFileName + " for reading");
+
+    // TODO: adapt to different formats by setting this flag
+    bool eigenvectors_as_columns = true;
+
+    if (eigenvectors_as_columns) {
+      // Carma-style dPCA file
+      std::string line;
+      cvm::real c;
+      while (vecFile.good()) {
+        getline(vecFile, line);
+        if (line.length() < 2) break;
+        std::istringstream ls (line);
+        for (int i=0; i<vecNumber; i++) ls >> c;
+        coeffs.push_back(c);
+      }
+    }
+/*  TODO Uncomment this when different formats are recognized
+    else {
+      // Eigenvectors as lines
+      // Skip to the right line
+      for (int i = 1; i<vecNumber; i++)
+        vecFile.ignore(999999, '\n');
+
+      if (!vecFile.good())
+        cvm::fatal_error ("Error reading dihedral PCA vector file " + vecFileName);
+
+      std::string line;
+      getline(vecFile, line);
+      std::istringstream ls (line);
+      cvm::real c;
+      while (ls.good()) {
+        ls >> c;
+        coeffs.push_back(c);
+      }
+    }
+ */
+    vecFile.close();
+
+  } else {
+    get_keyval (conf, "vector", coeffs, coeffs);
+  }
+
+  if ( coeffs.size() != 4 * (residues.size() - 1)) {
+    cvm::fatal_error ("Error: wrong number of coefficients: " +
+        cvm::to_str(coeffs.size()) + ". Expected " +
+        cvm::to_str(4 * (residues.size() - 1)) +
+        " (4 coeffs per residue, minus one residue).\n");
+  }
+
+  for (size_t i = 0; i < residues.size()-1; i++) {
+    // Psi
+    theta.push_back (new colvar::dihedral (cvm::atom (r[i  ], "N", sid),
+                                           cvm::atom (r[i  ], "CA", sid),
+                                           cvm::atom (r[i  ], "C", sid),
+                                           cvm::atom (r[i+1], "N", sid)));
+    // Phi (next res)
+    theta.push_back (new colvar::dihedral (cvm::atom (r[i  ], "C", sid),
+                                           cvm::atom (r[i+1], "N", sid),
+                                           cvm::atom (r[i+1], "CA", sid),
+                                           cvm::atom (r[i+1], "C", sid)));
+  }
+
+  if (cvm::debug())
+    cvm::log ("Done initializing dihedPC object.\n");
+}
+
+
+colvar::dihedPC::dihedPC()
+  : cvc ()
+{
+  function_type = "dihedPC";
+  x.type (colvarvalue::type_scalar);
+}
+
+
+colvar::dihedPC::~dihedPC()
+{
+  while (theta.size() != 0) {
+    delete theta.back();
+    theta.pop_back();
+  }
+}
+
+
+void colvar::dihedPC::calc_value()
+{
+  x.real_value = 0.0;
+  for (size_t i = 0; i < theta.size(); i++) {
+    theta[i]->calc_value();
+    cvm::real const t = (PI / 180.) * theta[i]->value().real_value;
+    x.real_value += coeffs[2*i  ] * std::cos(t)
+                  + coeffs[2*i+1] * std::sin(t);
+  }
+}
+
+
+void colvar::dihedPC::calc_gradients()
+{
+  for (size_t i = 0; i < theta.size(); i++) {
+    theta[i]->calc_gradients();
+  }
+}
+
+
+void colvar::dihedPC::apply_force (colvarvalue const &force)
+{
+  for (size_t i = 0; i < theta.size(); i++) {
+    cvm::real const t = (PI / 180.) * theta[i]->value().real_value;
+    cvm::real const dcosdt = - (PI / 180.) * std::sin(t);
+    cvm::real const dsindt =   (PI / 180.) * std::cos(t);
+
+    theta[i]->apply_force((coeffs[2*i  ] * dcosdt +
+                           coeffs[2*i+1] * dsindt) * force);
+  }
+}

lib/colvars/colvarcomp_rotations.C

@@ -0,0 +1,371 @@
+#include <cmath>
+
+#include "colvarmodule.h"
+#include "colvarvalue.h"
+#include "colvarparse.h"
+#include "colvar.h"
+#include "colvarcomp.h"
+
+
+
+colvar::orientation::orientation (std::string const &conf)
+  : cvc (conf)
+{
+  function_type = "orientation";
+  parse_group (conf, "atoms", atoms);
+  atom_groups.push_back (&atoms);
+  x.type (colvarvalue::type_quaternion);
+
+  ref_pos.reserve (atoms.size());
+
+  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");
+    }
+  }
+
+  {
+    std::string file_name;
+    if (get_keyval (conf, "refPositionsFile", file_name)) {
+
+      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);
+    }
+  }
+
+  if (!ref_pos.size()) {
+    cvm::fatal_error ("Error: must define a set of "
+                      "reference coordinates.\n");
+  }
+
+
+  cvm::log ("Centering the reference coordinates: it is "
+            "assumed that each atom is the closest "
+            "periodic image to the center of geometry.\n");
+  cvm::rvector cog (0.0, 0.0, 0.0);
+  for (size_t i = 0; i < ref_pos.size(); i++) {
+    cog += ref_pos[i];
+  }
+  cog /= cvm::real (ref_pos.size());
+  for (size_t i = 0; i < ref_pos.size(); i++) {
+    ref_pos[i] -= cog;
+  }
+
+  get_keyval (conf, "closestToQuaternion", ref_quat, cvm::quaternion (1.0, 0.0, 0.0, 0.0));
+
+  // initialize rot member data
+  if (!atoms.noforce) {
+    rot.request_group2_gradients (atoms.size());
+  }
+
+  rot.b_debug_gradients = this->b_debug_gradients;
+}
+
+  
+colvar::orientation::orientation()
+  : cvc ()
+{
+  function_type = "orientation";
+  x.type (colvarvalue::type_quaternion);
+}
+
+
+void colvar::orientation::calc_value()
+{
+  atoms.reset_atoms_data();
+  atoms.read_positions();
+
+  atoms_cog = atoms.center_of_geometry();
+
+  rot.calc_optimal_rotation (ref_pos, atoms.positions_shifted (-1.0 * atoms_cog));
+
+  if ((rot.q).inner (ref_quat) >= 0.0) {
+    x.quaternion_value = rot.q;
+  } else {
+    x.quaternion_value = -1.0 * rot.q;
+  }
+}
+
+
+void colvar::orientation::calc_gradients()
+{
+  // gradients have already been calculated and stored within the
+  // member object "rot"; we're not using the "grad" member of each
+  // atom object, because it only can represent the gradient of a
+  // scalar colvar
+}
+
+
+void colvar::orientation::apply_force (colvarvalue const &force)
+{
+  cvm::quaternion const &FQ = force.quaternion_value;
+
+  if (!atoms.noforce) {
+    for (size_t ia = 0; ia < atoms.size(); ia++) {
+      for (size_t i = 0; i < 4; i++) {
+        atoms[ia].apply_force (FQ[i] * rot.dQ0_2[ia][i]);
+      }
+    }
+  }
+}
+
+
+
+colvar::orientation_angle::orientation_angle (std::string const &conf)
+  : orientation (conf)
+{
+  function_type = "orientation_angle";
+  x.type (colvarvalue::type_scalar);
+}
+
+
+colvar::orientation_angle::orientation_angle()
+  : orientation()
+{
+  function_type = "orientation_angle";
+  x.type (colvarvalue::type_scalar);
+}
+
+
+void colvar::orientation_angle::calc_value()
+{
+  atoms.reset_atoms_data();
+  atoms.read_positions();
+
+  atoms_cog = atoms.center_of_geometry();
+
+  rot.calc_optimal_rotation (ref_pos, atoms.positions_shifted (-1.0 * atoms_cog));
+
+  if ((rot.q).q0 >= 0.0) {
+    x.real_value = (180.0/PI) * 2.0 * std::acos ((rot.q).q0);
+  } else {
+    x.real_value = (180.0/PI) * 2.0 * std::acos (-1.0 * (rot.q).q0);
+  }
+}
+
+
+void colvar::orientation_angle::calc_gradients()
+{
+  cvm::real const dxdq0 =
+    ( ((rot.q).q0 * (rot.q).q0 < 1.0) ? 
+      ((180.0 / PI) * (-2.0) / std::sqrt (1.0 - ((rot.q).q0 * (rot.q).q0))) :
+      0.0 );
+
+  for (size_t ia = 0; ia < atoms.size(); ia++) {
+    atoms[ia].grad = (dxdq0 * (rot.dQ0_2[ia])[0]);
+  }
+}
+
+
+void colvar::orientation_angle::apply_force (colvarvalue const &force)
+{
+  cvm::real const &fw = force.real_value;
+
+  if (!atoms.noforce) {
+    atoms.apply_colvar_force (fw);
+  }
+}
+
+
+colvar::tilt::tilt (std::string const &conf)
+  : orientation (conf)
+{
+  function_type = "tilt";
+
+  get_keyval (conf, "axis", axis, cvm::rvector (0.0, 0.0, 1.0));
+
+  if (axis.norm2() != 1.0) {
+    axis /= axis.norm();
+    cvm::log ("Normalizing rotation axis to "+cvm::to_str (axis)+".\n");
+  }
+
+  x.type (colvarvalue::type_scalar);
+}
+
+
+colvar::tilt::tilt()
+  : orientation()
+{
+  function_type = "tilt";
+  x.type (colvarvalue::type_scalar);
+}
+
+
+void colvar::tilt::calc_value()
+{
+  atoms.reset_atoms_data();
+  atoms.read_positions();
+
+  atoms_cog = atoms.center_of_geometry();
+
+  rot.calc_optimal_rotation (ref_pos, atoms.positions_shifted (-1.0 * atoms_cog));
+
+  x.real_value = rot.cos_theta (axis);
+}
+
+
+void colvar::tilt::calc_gradients()
+{
+  cvm::quaternion const dxdq = rot.dcos_theta_dq (axis);
+
+  for (size_t ia = 0; ia < atoms.size(); ia++) {
+    atoms[ia].grad = cvm::rvector (0.0, 0.0, 0.0);
+    for (size_t iq = 0; iq < 4; iq++) {
+      atoms[ia].grad += (dxdq[iq] * (rot.dQ0_2[ia])[iq]);
+    }
+  }
+
+  if (cvm::debug()) {
+
+    std::vector<cvm::atom_pos> pos_test (atoms.positions_shifted (-1.0 * atoms_cog));
+
+    for (size_t comp = 0; comp < 3; comp++) {
+      
+      // correct this cartesian component for the "new" center of geometry
+      for (size_t ia = 0; ia < atoms.size(); ia++) {
+        pos_test[ia][comp] -=
+          cvm::debug_gradients_step_size / cvm::real (atoms.size());
+      }
+
+      for (size_t ia = 0; ia < atoms.size(); ia++) {
+
+        pos_test[ia][comp] += cvm::debug_gradients_step_size;
+        rot.calc_optimal_rotation (ref_pos, pos_test);
+        pos_test[ia][comp] -= cvm::debug_gradients_step_size;
+
+        cvm::log ("|dx(real) - dx(pred)|/dx(real) = "+
+                  cvm::to_str (std::fabs (rot.cos_theta (axis) - x.real_value -
+                                       cvm::debug_gradients_step_size * atoms[ia].grad[comp]) /
+                               std::fabs (rot.cos_theta (axis) - x.real_value),
+                               cvm::cv_width, cvm::cv_prec)+".\n");
+      }
+
+      for (size_t ia = 0; ia < atoms.size(); ia++) {
+        pos_test[ia][comp] +=
+          cvm::debug_gradients_step_size / cvm::real (atoms.size());
+      }
+    }
+  }
+}
+
+
+void colvar::tilt::apply_force (colvarvalue const &force)
+{
+  cvm::real const &fw = force.real_value;
+
+  if (!atoms.noforce) {
+    atoms.apply_colvar_force (fw);
+  }
+}
+
+
+
+colvar::spin_angle::spin_angle (std::string const &conf)
+  : orientation (conf)
+{
+  function_type = "spin_angle";
+
+  get_keyval (conf, "axis", axis, cvm::rvector (0.0, 0.0, 1.0));
+
+  if (axis.norm2() != 1.0) {
+    axis /= axis.norm();
+    cvm::log ("Normalizing rotation axis to "+cvm::to_str (axis)+".\n");
+  }
+
+  period = 360.0;
+  b_periodic = true;
+  x.type (colvarvalue::type_scalar);
+}
+
+
+colvar::spin_angle::spin_angle()
+  : orientation()
+{
+  function_type = "spin_angle";
+  period = 360.0;
+  b_periodic = true;
+  x.type (colvarvalue::type_scalar);
+}
+
+
+void colvar::spin_angle::calc_value()
+{
+  atoms.reset_atoms_data();
+  atoms.read_positions();
+
+  atoms_cog = atoms.center_of_geometry();
+
+  rot.calc_optimal_rotation (ref_pos, atoms.positions_shifted (-1.0 * atoms_cog));
+
+  x.real_value = rot.spin_angle (axis);
+  this->wrap (x);
+}
+
+
+void colvar::spin_angle::calc_gradients()
+{
+  cvm::quaternion const dxdq = rot.dspin_angle_dq (axis);
+
+  for (size_t ia = 0; ia < atoms.size(); ia++) {
+    atoms[ia].grad = cvm::rvector (0.0, 0.0, 0.0);
+    for (size_t iq = 0; iq < 4; iq++) {
+      atoms[ia].grad += (dxdq[iq] * (rot.dQ0_2[ia])[iq]);
+    }
+  }
+
+  if (cvm::debug()) {
+
+    std::vector<cvm::atom_pos> pos_test (atoms.positions_shifted (-1.0 * atoms_cog));
+
+    for (size_t comp = 0; comp < 3; comp++) {
+      
+      // correct this cartesian component for the "new" center of geometry
+      for (size_t ia = 0; ia < atoms.size(); ia++) {
+        pos_test[ia][comp] -=
+          cvm::debug_gradients_step_size / cvm::real (atoms.size());
+      }
+
+      for (size_t ia = 0; ia < atoms.size(); ia++) {
+
+        pos_test[ia][comp] += cvm::debug_gradients_step_size;
+        rot.calc_optimal_rotation (ref_pos, pos_test);
+        pos_test[ia][comp] -= cvm::debug_gradients_step_size;
+
+        cvm::log ("|dx(real) - dx(pred)|/dx(real) = "+
+                  cvm::to_str (std::fabs (rot.spin_angle (axis) - x.real_value -
+                                       cvm::debug_gradients_step_size * atoms[ia].grad[comp]) /
+                               std::fabs (rot.spin_angle (axis) - x.real_value),
+                               cvm::cv_width, cvm::cv_prec)+".\n");
+      }
+
+      for (size_t ia = 0; ia < atoms.size(); ia++) {
+        pos_test[ia][comp] +=
+          cvm::debug_gradients_step_size / cvm::real (atoms.size());
+      }
+    }
+  }
+}
+
+
+void colvar::spin_angle::apply_force (colvarvalue const &force)
+{
+  cvm::real const &fw = force.real_value;
+
+  if (!atoms.noforce) {
+    atoms.apply_colvar_force (fw);
+  }
+}

lib/colvars/colvargrid.C

@@ -0,0 +1,217 @@
+// -*- c++ -*-
+
+#include "colvarmodule.h"
+#include "colvarvalue.h"
+#include "colvarparse.h"
+#include "colvar.h"
+#include "colvarcomp.h"
+#include "colvargrid.h"
+
+
+colvar_grid_count::colvar_grid_count()
+  : colvar_grid<size_t>()
+{}
+
+colvar_grid_count::colvar_grid_count (std::vector<int> const &nx_i,
+                                      size_t const           &def_count)
+  : colvar_grid<size_t> (nx_i, def_count)
+{}
+
+colvar_grid_count::colvar_grid_count (std::vector<colvar *>  &colvars,
+                                      size_t const           &def_count)
+  : colvar_grid<size_t> (colvars, def_count)
+{}
+
+std::istream & colvar_grid_count::read_restart (std::istream &is)
+{
+  size_t const start_pos = is.tellg();
+  std::string key, conf;
+  if ((is >> key) && (key == std::string ("grid_parameters"))) {
+    is.seekg (start_pos, std::ios::beg);
+    is >> colvarparse::read_block ("grid_parameters", conf);
+    parse_params (conf);
+  } else {
+    cvm::log ("Grid parameters are missing in the restart file, using those from the configuration.\n");
+    is.seekg (start_pos, std::ios::beg);
+  }
+  read_raw (is);
+  return is;
+}
+
+std::ostream & colvar_grid_count::write_restart (std::ostream &os)
+{
+  write_params (os);
+  write_raw (os);
+  return os;
+}
+
+
+
+colvar_grid_scalar::colvar_grid_scalar()
+  : colvar_grid<cvm::real>(), samples (NULL), grad (NULL)
+{}
+
+colvar_grid_scalar::colvar_grid_scalar (colvar_grid_scalar const &g)
+  : colvar_grid<cvm::real> (g), samples (NULL), grad (NULL)
+{
+  grad = new cvm::real[nd];
+}
+
+colvar_grid_scalar::colvar_grid_scalar (std::vector<int> const &nx_i)
+  : colvar_grid<cvm::real> (nx_i, 0.0, 1), samples (NULL)
+{
+  grad = new cvm::real[nd];
+}
+
+colvar_grid_scalar::colvar_grid_scalar (std::vector<colvar *> &colvars, bool margin)
+  : colvar_grid<cvm::real> (colvars, 0.0, 1, margin), samples (NULL)
+{
+  grad = new cvm::real[nd];
+}
+
+colvar_grid_scalar::~colvar_grid_scalar()
+{
+  if (grad) {
+    delete [] grad;
+    grad = NULL;
+  }
+}
+
+std::istream & colvar_grid_scalar::read_restart (std::istream &is)
+{
+  size_t const start_pos = is.tellg();
+  std::string key, conf;
+  if ((is >> key) && (key == std::string ("grid_parameters"))) {
+    is.seekg (start_pos, std::ios::beg);
+    is >> colvarparse::read_block ("grid_parameters", conf);
+    parse_params (conf);
+  } else {
+    cvm::log ("Grid parameters are missing in the restart file, using those from the configuration.\n");
+    is.seekg (start_pos, std::ios::beg);
+  }
+  read_raw (is);
+  return is;
+}
+
+std::ostream & colvar_grid_scalar::write_restart (std::ostream &os)
+{
+  write_params (os);
+  write_raw (os);
+  return os;
+}
+
+
+
+colvar_grid_gradient::colvar_grid_gradient()
+  : colvar_grid<cvm::real>(), samples (NULL)
+{}
+
+colvar_grid_gradient::colvar_grid_gradient (std::vector<int> const &nx_i)
+  : colvar_grid<cvm::real> (nx_i, 0.0, nx_i.size()), samples (NULL)
+{}
+
+colvar_grid_gradient::colvar_grid_gradient (std::vector<colvar *> &colvars)
+  : colvar_grid<cvm::real> (colvars, 0.0, colvars.size()), samples (NULL)
+{}
+
+std::istream & colvar_grid_gradient::read_restart (std::istream &is)
+{
+  size_t const start_pos = is.tellg();
+  std::string key, conf;
+  if ((is >> key) && (key == std::string ("grid_parameters"))) {
+    is.seekg (start_pos, std::ios::beg);
+    is >> colvarparse::read_block ("grid_parameters", conf);
+    parse_params (conf);
+  } else {
+    cvm::log ("Grid parameters are missing in the restart file, using those from the configuration.\n");
+    is.seekg (start_pos, std::ios::beg);
+  }
+  read_raw (is);
+  return is;
+}
+
+std::ostream & colvar_grid_gradient::write_restart (std::ostream &os)
+{
+  write_params (os);
+  write_raw (os);
+  return os;
+}
+
+void colvar_grid_gradient::write_1D_integral (std::ostream &os)
+{
+  cvm::real bin, min, integral;
+  std::vector<cvm::real> int_vals;
+  
+  os << "#       xi            A(xi)\n";
+
+  if ( cv.size() != 1 ) {
+    cvm::fatal_error ("Cannot write integral for multi-dimensional gradient grids.");
+  }
+
+  integral = 0.0;
+  int_vals.push_back ( 0.0 );
+  bin = 0.0;
+  min = 0.0;
+
+  // correction for periodic colvars, so that the PMF is periodic
+  cvm::real corr;
+  if ( periodic[0] ) {
+    corr = average();
+  } else {
+    corr = 0.0;
+  }
+
+  for (std::vector<int> ix = new_index(); index_ok (ix); incr (ix), bin += 1.0 ) {
+	  
+    if (samples) {
+      size_t const samples_here = samples->value (ix);
+      if (samples_here)
+        integral += (value (ix) / cvm::real (samples_here) - corr) * cv[0]->width;
+    } else {
+      integral += (value (ix) - corr) * cv[0]->width;
+    }
+
+    if ( integral < min ) min = integral;
+    int_vals.push_back ( integral );
+  }
+
+  bin = 0.0;
+  for ( int i = 0; i < nx[0]; i++, bin += 1.0 ) {
+    os << std::setw (10) << cv[0]->lower_boundary.real_value + cv[0]->width * bin << " "
+       << std::setw (cvm::cv_width)
+       << std::setprecision (cvm::cv_prec)
+       << int_vals[i] - min << "\n";
+  }
+
+  os << std::setw (10) << cv[0]->lower_boundary.real_value + cv[0]->width * bin << " "
+     << std::setw (cvm::cv_width)
+     << std::setprecision (cvm::cv_prec)
+     << int_vals[nx[0]] - min << "\n";
+
+  return;
+}
+
+
+
+
+// quaternion_grid::quaternion_grid (std::vector<colvar *>      const &cv_i,
+//                                   std::vector<std::string>   const &grid_str)
+// {
+//   cv = cv_i;
+
+//   std::istringstream is (grid_str[0]);
+//   is >> grid_size;
+
+//   min.assign (3, -1.0);
+//   max.assign (3,  1.0);
+//   np.assign  (3, grid_size);
+//   dx.assign  (3, 2.0/(cvm::real (grid_size)));
+
+//   // assumes a uniform grid in the three directions; change
+//   // get_value() if you want to use different sizes
+//   cvm::log ("Allocating quaternion grid ("+cvm::to_str (np.size())+" dimensional)...");
+//   data.create (np, 0.0);
+//   cvm::log ("done.\n");
+//   if (cvm::debug()) cvm::log ("Grid size = "+data.size());
+// }
+

lib/colvars/colvarmodule.h

@@ -0,0 +1,483 @@
+#ifndef COLVARMODULE_H
+#define COLVARMODULE_H
+
+#ifndef COLVARS_VERSION
+#define COLVARS_VERSION "2012-03-03"
+#endif
+
+#ifndef COLVARS_DEBUG
+#define COLVARS_DEBUG false
+#endif
+
+/// \file colvarmodule.h 
+/// \brief Collective variables main module
+/// 
+/// This file declares the main class for defining and manipulating
+/// collective variables: there should be only one instance of this
+/// class, because several variables are made static (i.e. they are
+/// shared between all object instances) to be accessed from other
+/// objects.
+
+
+#include <iostream>
+#include <iomanip>
+#include <string>
+#include <sstream>
+#include <fstream>
+#include <cmath>
+#include <vector>
+#include <list>
+
+
+class colvarparse;
+class colvar;
+class colvarbias;
+class colvarproxy;
+
+
+/// \brief Collective variables module (main class)
+///
+/// Class to control the collective variables calculation.  An object
+/// (usually one) of this class is spawned from the MD program,
+/// containing all i/o routines and general interface.
+/// 
+/// At initialization, the colvarmodule object creates a proxy object
+/// to provide a transparent interface between the MD program and the
+/// child objects
+class colvarmodule {
+
+private:
+
+  /// Impossible to initialize the main object without arguments
+  colvarmodule();
+
+public:
+
+  friend class colvarproxy;
+  
+  /// Defining an abstract real number allows to switch precision
+  typedef  double    real;
+  /// Residue identifier
+  typedef  int       residue_id;
+
+  class rvector;
+  template <class T,
+            size_t const length> class vector1d;
+  template <class T,
+            size_t const outer_length,
+            size_t const inner_length> class matrix2d;
+  class quaternion;
+  class rotation;
+
+  /// \brief Atom position (different type name from rvector, to make
+  /// possible future PBC-transparent implementations)
+  typedef rvector atom_pos;
+
+  /// \brief 3x3 matrix of real numbers
+  class rmatrix;
+
+  // allow these classes to access protected data
+  class atom;
+  class atom_group;
+  friend class atom;
+  friend class atom_group;
+  typedef std::vector<atom>::iterator       atom_iter;
+  typedef std::vector<atom>::const_iterator atom_const_iter;
+
+  /// Current step number
+  static size_t it;
+  /// Starting step number for this run
+  static size_t it_restart;
+
+  /// Return the current step number from the beginning of this run
+  static inline size_t step_relative()
+  {
+    return it - it_restart;
+  }
+
+  /// Return the current step number from the beginning of the whole
+  /// calculation
+  static inline size_t step_absolute()
+  {
+    return it;
+  }
+
+  /// \brief Finite difference step size (if there is no dynamics, or
+  /// if gradients need to be tested independently from the size of
+  /// dt)
+  static real debug_gradients_step_size;
+
+  /// Prefix for all output files for this run
+  static std::string output_prefix;
+
+  /// Prefix for files from a previous run (including restart/output)
+  static std::string input_prefix;
+
+  /// input restart file name (determined from input_prefix)
+  static std::string restart_in_name;
+
+
+  /// Array of collective variables
+  static std::vector<colvar *>     colvars;
+
+  /* TODO: implement named CVCs
+  /// Array of named (reusable) collective variable components
+  static std::vector<cvc *>     cvcs;
+  /// Named cvcs register themselves at initialization time
+  inline void register_cvc (cvc *p) {
+    cvcs.push_back(p);
+  }
+  */
+
+  /// Array of collective variable biases
+  static std::vector<colvarbias *> biases;
+  /// \brief Number of ABF biases initialized (in normal conditions
+  /// should be 1)
+  static size_t n_abf_biases;
+  /// \brief Number of metadynamics biases initialized (in normal
+  /// conditions should be 1)
+  static size_t n_meta_biases;
+  /// \brief Number of harmonic biases initialized (no limit on the
+  /// number)
+  static size_t n_harm_biases;
+  /// \brief Number of histograms initialized (no limit on the
+  /// number)
+  static size_t n_histo_biases;
+   
+  /// \brief Whether debug output should be enabled (compile-time option)
+  static inline bool debug()
+  {
+    return COLVARS_DEBUG;
+  }
+
+
+  /// \brief Constructor \param config_name Configuration file name
+  /// \param restart_name (optional) Restart file name
+  colvarmodule (char const *config_name, 
+                colvarproxy *proxy_in);
+  /// Destructor
+  ~colvarmodule();
+
+  /// Initialize collective variables
+  void init_colvars (std::string const &conf);
+
+  /// Initialize collective variable biases
+  void init_biases (std::string const &conf);
+
+  /// Calculate collective variables and biases
+  void calc();
+  /// Read the input restart file
+  std::istream & read_restart (std::istream &is);
+  /// Write the output restart file
+  std::ostream & write_restart (std::ostream &os);
+  /// Write all output files (called by the proxy)
+  void write_output_files();
+  /// \brief Call colvarproxy::backup_file()
+  static void backup_file (char const *filename);
+
+  /// Perform analysis
+  void analyze();
+  /// \brief Read a collective variable trajectory (post-processing
+  /// only, not called at runtime)
+  bool read_traj (char const *traj_filename,
+                  size_t      traj_read_begin,
+                  size_t      traj_read_end);
+ 
+  /// Get the pointer of a colvar from its name (returns NULL if not found)
+  static colvar * colvar_p (std::string const &name);
+
+  /// Quick conversion of an object to a string
+  template<typename T> static std::string to_str (T const &x, 
+                                                  size_t const &width = 0,
+                                                  size_t const &prec = 0);
+  /// Quick conversion of a vector of objects to a string
+  template<typename T> static std::string to_str (std::vector<T> const &x, 
+                                                  size_t const &width = 0,
+                                                  size_t const &prec = 0);
+
+  /// Reduce the number of characters in a string
+  static inline std::string wrap_string (std::string const &s,
+                                         size_t const &nchars)
+  {
+    if (!s.size())
+      return std::string (nchars, ' ');
+    else
+      return ( (s.size() <= size_t (nchars)) ?
+               (s+std::string (nchars-s.size(), ' ')) :
+               (std::string (s, 0, nchars)) );
+  }
+
+  /// Number of characters to represent a time step
+  static size_t const it_width;
+  /// Number of digits to represent a collective variables value(s)
+  static size_t const cv_prec;
+  /// Number of characters to represent a collective variables value(s)
+  static size_t const cv_width;
+  /// Number of digits to represent the collective variables energy
+  static size_t const en_prec;
+  /// Number of characters to represent the collective variables energy
+  static size_t const en_width;
+  /// Line separator in the log output
+  static std::string const line_marker;
+
+
+  // proxy functions
+
+  /// \brief Value of the unit for atomic coordinates with respect to
+  /// angstroms (used by some variables for hard-coded default values)
+  static real unit_angstrom();
+
+  /// \brief Boltmann constant
+  static real boltzmann();
+
+  /// \brief Temperature of the simulation (K)
+  static real temperature();
+
+  /// \brief Time step of MD integrator (fs)
+  static real dt();
+  
+  /// Request calculation of system force from MD engine
+  static void request_system_force();
+  
+  /// Print a message to the main log
+  static void log (std::string const &message);
+
+  /// Print a message to the main log and exit with error code
+  static void fatal_error (std::string const &message);
+
+  /// Print a message to the main log and exit normally
+  static void exit (std::string const &message);
+
+
+  /// \brief Get the distance between two atomic positions with pbcs handled
+  /// correctly
+  static rvector position_distance (atom_pos const &pos1,
+                                    atom_pos const &pos2);
+
+
+  /// \brief Get the square distance between two positions (with
+  /// periodic boundary conditions handled transparently)
+  ///
+  /// Note: in the case of periodic boundary conditions, this provides
+  /// an analytical square distance (while taking the square of
+  /// position_distance() would produce leads to a cusp)
+  static real position_dist2 (atom_pos const &pos1,
+                              atom_pos const &pos2);
+
+  /// \brief Get the closest periodic image to a reference position
+  /// \param pos The position to look for the closest periodic image
+  /// \param ref_pos (optional) The reference position 
+  static void select_closest_image (atom_pos &pos,
+                                    atom_pos const &ref_pos);
+
+  /// \brief Perform select_closest_image() on a set of atomic positions
+  ///
+  /// After that, distance vectors can then be calculated directly,
+  /// without using position_distance()
+  static void select_closest_images (std::vector<atom_pos> &pos,
+                                     atom_pos const &ref_pos);
+
+
+  /// \brief Create atoms from a file \param filename name of the file
+  /// (usually a PDB) \param atoms array of the atoms to be allocated
+  /// \param pdb_field (optiona) if "filename" is a PDB file, use this
+  /// field to determine which are the atoms to be set
+  static void load_atoms (char const *filename,
+                          std::vector<atom> &atoms,
+                          std::string const &pdb_field,
+                          double const pdb_field_value = 0.0);
+
+  /// \brief Load the coordinates for a group of atoms from a file
+  /// (usually a PDB); the number of atoms in "filename" must match
+  /// the number of elements in "pos"
+  static void load_coords (char const *filename,
+                           std::vector<atom_pos> &pos,
+                           const std::vector<int> &indices,
+                           std::string const &pdb_field,
+                           double const pdb_field_value = 0.0);
+
+
+  /// Frequency for collective variables trajectory output
+  static size_t cv_traj_freq;
+
+  /// \brief True if only analysis is performed and not a run
+  static bool   b_analysis;
+
+  /// Frequency for saving output restarts
+  static size_t restart_out_freq;
+  /// Output restart file name
+  std::string   restart_out_name;
+
+  /// Pseudo-random number with Gaussian distribution
+  static real rand_gaussian (void);
+protected:
+
+  /// Configuration file
+  std::ifstream config_s;
+
+  /// Configuration file parser object
+  colvarparse *parse;
+
+  /// Name of the trajectory file
+  std::string   cv_traj_name;
+
+  /// Collective variables output trajectory file
+  std::ofstream cv_traj_os;
+
+  /// Appending to the existing trajectory file?
+  bool          cv_traj_append;
+
+  /// Output restart file
+  std::ofstream restart_out_os;
+
+  /// \brief Pointer to the proxy object, used to retrieve atomic data
+  /// from the hosting program; it is static in order to be accessible
+  /// from static functions in the colvarmodule class
+  static colvarproxy *proxy;
+
+  /// \brief Counter for the current depth in the object hierarchy (useg e.g. in outpu
+  static size_t depth;
+
+public:
+
+  /// Increase the depth (number of indentations in the output)
+  static void increase_depth();
+
+  /// Decrease the depth (number of indentations in the output)
+  static void decrease_depth();
+};
+
+
+/// Shorthand for the frequently used type prefix
+typedef colvarmodule cvm;
+
+
+#include "colvartypes.h"
+
+
+std::ostream & operator << (std::ostream &os, cvm::rvector const &v);
+std::istream & operator >> (std::istream &is, cvm::rvector &v);
+
+
+template<typename T> std::string cvm::to_str (T const &x, 
+                                              size_t const &width,
+                                              size_t const &prec) {
+  std::ostringstream os;
+  if (width) os.width (width);
+  if (prec) {
+    os.setf (std::ios::scientific, std::ios::floatfield);
+    os.precision (prec);
+  }
+  os << x;
+  return os.str();
+}
+
+template<typename T> std::string cvm::to_str (std::vector<T> const &x, 
+                                              size_t const &width,
+                                              size_t const &prec) {
+  if (!x.size()) return std::string ("");
+  std::ostringstream os;
+  if (prec) {
+    os.setf (std::ios::scientific, std::ios::floatfield);
+  }
+  os << "{ ";
+  if (width) os.width (width);
+  if (prec) os.precision (prec);
+  os << x[0];
+  for (size_t i = 1; i < x.size(); i++) {
+    os << ", ";
+    if (width) os.width (width);
+    if (prec) os.precision (prec);
+    os << x[i];
+  }
+  os << " }";
+  return os.str();
+}
+
+
+#include "colvarproxy.h"
+
+
+inline cvm::real cvm::unit_angstrom()
+{
+  return proxy->unit_angstrom();
+}
+
+inline cvm::real cvm::boltzmann()
+{
+  return proxy->boltzmann();
+}
+
+inline cvm::real cvm::temperature()
+{
+  return proxy->temperature();
+}
+
+inline cvm::real cvm::dt()
+{
+  return proxy->dt();
+}
+
+inline void cvm::request_system_force()
+{
+  proxy->request_system_force (true);
+}
+  
+inline void cvm::select_closest_image (atom_pos &pos,
+                                       atom_pos const &ref_pos)
+{
+  proxy->select_closest_image (pos, ref_pos);
+}
+
+inline void cvm::select_closest_images (std::vector<atom_pos> &pos,
+                                        atom_pos const &ref_pos)
+{
+  proxy->select_closest_images (pos, ref_pos);
+}
+
+inline cvm::rvector cvm::position_distance (atom_pos const &pos1,
+                                            atom_pos const &pos2)
+{
+  return proxy->position_distance (pos1, pos2);
+}
+
+inline cvm::real cvm::position_dist2 (cvm::atom_pos const &pos1,
+                                      cvm::atom_pos const &pos2)
+{
+  return proxy->position_dist2 (pos1, pos2);
+}
+
+inline void cvm::load_atoms (char const *file_name,
+                             std::vector<cvm::atom> &atoms,
+                             std::string const &pdb_field,
+                             double const pdb_field_value)
+{
+  proxy->load_atoms (file_name, atoms, pdb_field, pdb_field_value);
+}
+
+inline void cvm::load_coords (char const *file_name,
+                              std::vector<cvm::atom_pos> &pos,
+                              const std::vector<int> &indices,
+                              std::string const &pdb_field,
+                              double const pdb_field_value)
+{
+  proxy->load_coords (file_name, pos, indices, pdb_field, pdb_field_value);
+}
+
+inline void cvm::backup_file (char const *filename)
+{
+  proxy->backup_file (filename);
+}
+
+inline cvm::real cvm::rand_gaussian (void)
+{
+  return proxy->rand_gaussian();
+}
+
+#endif
+
+
+// Emacs
+// Local Variables:
+// mode: C++
+// End:

lib/colvars/colvarparse.C

@@ -0,0 +1,636 @@
+#include <sstream>
+#include <iostream>
+
+
+
+#include "colvarmodule.h"
+#include "colvarvalue.h"
+#include "colvarparse.h"
+
+
+
+// space & tab
+std::string const colvarparse::white_space = " \t";
+
+std::string colvarparse::dummy_string = "";
+size_t      colvarparse::dummy_pos = 0;
+
+
+// definition of single-value keyword parsers
+
+#define _get_keyval_scalar_string_(TYPE)                                \
+                                                                        \
+  bool colvarparse::get_keyval (std::string const &conf,                \
+                                char const *key,                        \
+                                TYPE &value,                            \
+                                TYPE const &def_value,                  \
+                                Parse_Mode const parse_mode)            \
+  {                                                                     \
+    std::string data;                                                   \
+    bool b_found = false, b_found_any = false;                          \
+    size_t save_pos = 0, found_count = 0;                               \
+                                                                        \
+    do {                                                                \
+      std::string data_this = "";                                       \
+      b_found = key_lookup (conf, key, data_this, save_pos);            \
+      if (b_found) {                                                    \
+        if (!b_found_any)                                               \
+          b_found_any = true;                                           \
+        found_count++;                                                  \
+        data = data_this;                                               \
+      }                                                                 \
+    } while (b_found);                                                  \
+                                                                        \
+    if (found_count > 1)                                                \
+      cvm::log ("Warning: found more than one instance of \""+          \
+                std::string (key)+"\".\n");                             \
+                                                                        \
+    if (data.size()) {                                                  \
+      std::istringstream is (data);                                     \
+      TYPE x (def_value);                                               \
+      if (is >> x)                                                      \
+        value = x;                                                      \
+      else                                                              \
+        cvm::fatal_error ("Error: in parsing \""+                       \
+                          std::string (key)+"\".\n");                   \
+      if (parse_mode != parse_silent) {                                 \
+        cvm::log ("# "+std::string (key)+" = "+                         \
+                  cvm::to_str (value)+"\n");                            \
+      }                                                                 \
+    } else {                                                            \
+                                                                        \
+      if (b_found_any)                                                  \
+        cvm::fatal_error ("Error: improper or missing value "           \
+                          "for \""+std::string (key)+"\".\n");          \
+      value = def_value;                                                \
+      if (parse_mode != parse_silent) {                                 \
+        cvm::log ("# "+std::string (key)+" = \""+                       \
+                  cvm::to_str (def_value)+"\" [default]\n");            \
+      }                                                                 \
+    }                                                                   \
+                                                                        \
+    return b_found_any;                                                 \
+  }
+
+
+#define _get_keyval_scalar_(TYPE)                                       \
+                                                                        \
+  bool colvarparse::get_keyval (std::string const &conf,                \
+                                char const *key,                        \
+                                TYPE &value,                            \
+                                TYPE const &def_value,                  \
+                                Parse_Mode const parse_mode)            \
+  {                                                                     \
+    std::string data;                                                   \
+    bool b_found = false, b_found_any = false;                          \
+    size_t save_pos = 0, found_count = 0;                               \
+                                                                        \
+    do {                                                                \
+      std::string data_this = "";                                       \
+      b_found = key_lookup (conf, key, data_this, save_pos);            \
+      if (b_found) {                                                    \
+        if (!b_found_any)                                               \
+          b_found_any = true;                                           \
+        found_count++;                                                  \
+        data = data_this;                                               \
+      }                                                                 \
+    } while (b_found);                                                  \
+                                                                        \
+    if (found_count > 1)                                                \
+      cvm::log ("Warning: found more than one instance of \""+          \
+                std::string (key)+"\".\n");                             \
+                                                                        \
+    if (data.size()) {                                                  \
+      std::istringstream is (data);                                     \
+      TYPE x (def_value);                                               \
+      if (is >> x)                                                      \
+        value = x;                                                      \
+      else                                                              \
+        cvm::fatal_error ("Error: in parsing \""+                       \
+                          std::string (key)+"\".\n");                   \
+      if (parse_mode != parse_silent) {                                 \
+        cvm::log ("# "+std::string (key)+" = "+                         \
+                  cvm::to_str (value)+"\n");                            \
+      }                                                                 \
+    } else {                                                            \
+                                                                        \
+      if (b_found_any)                                                  \
+        cvm::fatal_error ("Error: improper or missing value "           \
+                          "for \""+std::string (key)+"\".\n");          \
+      value = def_value;                                                \
+      if (parse_mode != parse_silent) {                                 \
+        cvm::log ("# "+std::string (key)+" = "+                         \
+                  cvm::to_str (def_value)+" [default]\n");              \
+      }                                                                 \
+    }                                                                   \
+                                                                        \
+    return b_found_any;                                                 \
+  }
+
+
+// definition of multiple-value keyword parsers
+
+#define _get_keyval_vector_(TYPE)                                       \
+                                                                        \
+  bool colvarparse::get_keyval (std::string const &conf,                \
+                                char const *key,                        \
+                                std::vector<TYPE> &values,              \
+                                std::vector<TYPE> const &def_values,    \
+                                Parse_Mode const parse_mode)            \
+  {                                                                     \
+    std::string data;                                                   \
+    bool b_found = false, b_found_any = false;                          \
+    size_t save_pos = 0, found_count = 0;                               \
+                                                                        \
+    do {                                                                \
+      std::string data_this = "";                                       \
+      b_found = key_lookup (conf, key, data_this, save_pos);            \
+      if (b_found) {                                                    \
+        if (!b_found_any)                                               \
+          b_found_any = true;                                           \
+        found_count++;                                                  \
+        data = data_this;                                               \
+      }                                                                 \
+    } while (b_found);                                                  \
+                                                                        \
+    if (found_count > 1)                                                \
+      cvm::log ("Warning: found more than one instance of \""+          \
+                std::string (key)+"\".\n");                             \
+                                                                        \
+    if (data.size()) {                                                  \
+      std::istringstream is (data);                                     \
+                                                                        \
+      if (values.size() == 0) {                                         \
+                                                                        \
+        std::vector<TYPE> x;                                            \
+        if (def_values.size())                                          \
+          x = def_values;                                               \
+        else                                                            \
+          x.assign (1, TYPE());                                         \
+                                                                        \
+        for (size_t i = 0;                                              \
+             ( is >> x[ ((i<x.size()) ? i : x.size()-1) ] );            \
+             i++) {                                                     \
+          values.push_back (x[ ((i<x.size()) ? i : x.size()-1) ]);      \
+        }                                                               \
+                                                                        \
+      } else {                                                          \
+                                                                        \
+        size_t i = 0;                                                   \
+        for ( ; i < values.size(); i++) {                               \
+          TYPE x (values[i]);                                           \
+          if (is >> x)                                                  \
+            values[i] = x;                                              \
+          else                                                          \
+            cvm::fatal_error ("Error: in parsing \""+                   \
+                              std::string (key)+"\".\n");               \
+        }                                                               \
+      }                                                                 \
+                                                                        \
+      if (parse_mode != parse_silent) {                                 \
+        cvm::log ("# "+std::string (key)+" = "+                         \
+                  cvm::to_str (values)+"\n");                           \
+      }                                                                 \
+                                                                        \
+    } else {                                                            \
+                                                                        \
+      if (b_found_any)                                                  \
+        cvm::fatal_error ("Error: improper or missing values for \""+   \
+                          std::string (key)+"\".\n");                   \
+                                                                        \
+      for (size_t i = 0; i < values.size(); i++)                        \
+        values[i] = def_values[ (i > def_values.size()) ? 0 : i ];      \
+                                                                        \
+      if (parse_mode != parse_silent) {                                 \
+        cvm::log ("# "+std::string (key)+" = "+                         \
+                  cvm::to_str (def_values)+" [default]\n");             \
+      }                                                                 \
+    }                                                                   \
+                                                                        \
+    return b_found_any;                                                 \
+  }
+
+
+// single-value keyword parsers
+
+_get_keyval_scalar_ (int);
+_get_keyval_scalar_ (size_t);
+_get_keyval_scalar_string_ (std::string);
+_get_keyval_scalar_ (cvm::real);
+_get_keyval_scalar_ (cvm::rvector);
+_get_keyval_scalar_ (cvm::quaternion);
+_get_keyval_scalar_ (colvarvalue);
+
+
+// multiple-value keyword parsers
+
+_get_keyval_vector_ (int);
+_get_keyval_vector_ (size_t);
+_get_keyval_vector_ (std::string);
+_get_keyval_vector_ (cvm::real);
+_get_keyval_vector_ (cvm::rvector);
+_get_keyval_vector_ (cvm::quaternion);
+_get_keyval_vector_ (colvarvalue);
+
+
+
+bool colvarparse::get_keyval (std::string const &conf,        
+                              char const *key,                
+                              bool &value,                    
+                              bool const &def_value,          
+                              Parse_Mode const parse_mode)    
+{
+  std::string data;
+  bool b_found = false, b_found_any = false;
+  size_t save_pos = 0, found_count = 0;
+
+  do {
+    std::string data_this = "";
+    b_found = key_lookup (conf, key, data_this, save_pos);
+    if (b_found) {
+      if (!b_found_any)
+        b_found_any = true;
+      found_count++;
+      data = data_this;
+    }
+  } while (b_found);
+
+  if (found_count > 1)
+    cvm::log ("Warning: found more than one instance of \""+
+              std::string (key)+"\".\n");
+
+  if (data.size()) {
+    std::istringstream is (data);
+    if ( (data == std::string ("on")) ||
+         (data == std::string ("yes")) ||
+         (data == std::string ("true")) ) {
+      value = true;
+    } else if ( (data == std::string ("off")) ||
+              (data == std::string ("no")) ||
+              (data == std::string ("false")) ) {
+      value = false;
+    } else
+      cvm::fatal_error ("Error: boolean values only are allowed "  
+                        "for \""+std::string (key)+"\".\n");
+    if (parse_mode != parse_silent) {
+      cvm::log ("# "+std::string (key)+" = "+
+                (value ? "on" : "off")+"\n");
+    }
+  } else {
+
+    if (b_found_any) {
+      if (parse_mode != parse_silent) {
+        cvm::log ("# "+std::string (key)+" = on\n");
+      }
+      value = true;
+    } else {
+      value = def_value;
+      if (parse_mode != parse_silent) {
+        cvm::log ("# "+std::string (key)+" = "+
+                  (def_value ? "on" : "off")+" [default]\n");
+      }
+    }
+  }
+
+  return b_found_any;
+}
+
+
+void colvarparse::add_keyword (char const *key)
+{
+  for (std::list<std::string>::iterator ki = allowed_keywords.begin();
+       ki != allowed_keywords.end(); ki++) {
+    if (to_lower_cppstr (std::string (key)) == *ki)
+      return;
+  }
+  // not found in the list
+  //   if (cvm::debug())
+  //     cvm::log ("Registering a new keyword, \""+std::string (key)+"\".\n");
+  allowed_keywords.push_back (to_lower_cppstr (std::string (key)));
+}
+
+
+void colvarparse::strip_values (std::string &conf)
+{
+  size_t offset = 0;
+  data_begin_pos.sort();
+  data_end_pos.sort();
+
+  std::list<size_t>::iterator data_begin = data_begin_pos.begin();
+  std::list<size_t>::iterator data_end   = data_end_pos.begin();
+
+  for ( ; (data_begin != data_begin_pos.end()) &&
+          (data_end   != data_end_pos.end()) ;
+        data_begin++, data_end++) {
+
+    //     std::cerr << "data_begin, data_end "
+    //               << *data_begin << ", " << *data_end
+    //               << "\n";
+
+    size_t const nchars = *data_end-*data_begin;
+
+    //     std::cerr << "conf[data_begin:data_end] = \""
+    //               << std::string (conf, *data_begin - offset, nchars)
+    //               << "\"\n";
+
+    conf.erase (*data_begin - offset, nchars);
+    offset += nchars;
+
+    //     std::cerr << ("Stripped config = \"\n"+conf+"\"\n");
+
+  }
+}
+
+
+void colvarparse::check_keywords (std::string &conf, char const *key)
+{
+  if (cvm::debug())
+    cvm::log ("Configuration string for \""+std::string (key)+
+              "\": \"\n"+conf+"\".\n");
+
+  strip_values (conf);
+  // after stripping, the config string has either empty lines, or
+  // lines beginning with a keyword
+
+  std::string line;
+  std::istringstream is (conf);
+  while (std::getline (is, line)) {
+    if (line.size() == 0)
+      continue;
+    if (line.find_first_not_of (white_space) ==
+        std::string::npos)
+      continue;
+
+    std::string uk;
+    std::istringstream line_is (line);
+    line_is >> uk;
+    if (cvm::debug())
+      cvm::log ("Checking the validity of \""+uk+"\" from line:\n" + line);
+    uk = to_lower_cppstr (uk);
+
+    bool found_keyword = false;
+    for (std::list<std::string>::iterator ki = allowed_keywords.begin();
+         ki != allowed_keywords.end(); ki++) {
+      if (uk == *ki) {
+        found_keyword = true;
+        break;
+      }
+    }
+    if (!found_keyword)
+      cvm::fatal_error ("Error: keyword \""+uk+"\" is not supported, "
+                        "or not recognized in this context.\n");
+  }
+}
+
+
+std::istream & colvarparse::getline_nocomments (std::istream &is,
+                                                std::string &line,
+                                                char const delim)
+{
+  std::getline (is, line, delim);
+  size_t const comment = line.find ('#');
+  if (comment != std::string::npos) {
+    line.erase (comment);
+  }
+  return is;
+}
+
+
+bool colvarparse::key_lookup (std::string const &conf,
+                              char const *key_in,
+                              std::string &data,
+                              size_t &save_pos)
+{
+  // add this keyword to the register (in its camelCase version)
+  add_keyword (key_in);
+
+  // use the lowercase version from now on
+  std::string const key (to_lower_cppstr (key_in));
+
+  // "conf_lower" is only used to lookup the keyword, but its value
+  // will be read from "conf", in order not to mess up file names
+  std::string const conf_lower (to_lower_cppstr (conf));
+
+  // by default, there is no value, unless we found one
+  data = "";
+
+  // when the function is invoked without save_pos, ensure that we
+  // start from zero
+  colvarparse::dummy_pos = 0;
+
+  // start from the first occurrence of key
+  size_t pos = conf_lower.find (key, save_pos); 
+  
+  // iterate over all instances until it finds the isolated keyword
+  while (true) {
+
+    if (pos == std::string::npos) {
+      // no valid instance of the keyword has been found
+      return false;
+    }
+
+    bool b_isolated_left = true, b_isolated_right = true;
+
+    if (pos > 0) {
+      if ( std::string ("\n"+white_space+
+                        "}").find (conf[pos-1]) ==
+           std::string::npos ) {
+        // none of the valid delimiting characters is on the left of key 
+        b_isolated_left = false;
+      }
+    }
+
+    if (pos < conf.size()-key.size()-1) {
+      if ( std::string ("\n"+white_space+
+                        "{").find (conf[pos+key.size()]) ==
+           std::string::npos ) {
+        // none of the valid delimiting characters is on the right of key 
+        b_isolated_right = false;
+      }
+    }
+
+    // check that there are matching braces between here and the end of conf
+    bool const b_not_within_block = brace_check (conf, pos);
+
+    bool const b_isolated = (b_isolated_left && b_isolated_right &&
+                             b_not_within_block);
+    
+    if (b_isolated) {
+      // found it
+      break;
+    } else {
+      // try the next occurrence of key
+      pos = conf_lower.find (key, pos+key.size());
+    }
+  }
+
+  // check it is not between quotes
+  //   if ( (conf.find_last_of  ("\"",
+  //                             conf.find_last_of  (white_space, pos)) !=
+  //         std::string::npos) &&
+  //        (conf.find_first_of ("\"",
+  //                             conf.find_first_of (white_space, pos)) !=
+  //         std::string::npos) )
+  //     return false;
+
+
+  // save the pointer for a future call (when iterating over multiple
+  // valid instances of the same keyword)
+  save_pos = pos + key.size();
+
+  // get the remainder of the line
+  size_t pl = conf.rfind ("\n", pos);
+  size_t line_begin = (pl == std::string::npos) ? 0 : pos;
+  size_t nl = conf.find  ("\n", pos);
+  size_t line_end = (nl == std::string::npos) ? conf.size() : nl;
+  std::string line (conf, line_begin, (line_end-line_begin)); 
+
+  size_t data_begin = (to_lower_cppstr (line)).find (key) + key.size();
+  data_begin = line.find_first_not_of (white_space, data_begin+1);
+
+  //   size_t data_begin_absolute = data_begin + line_begin;
+  //   size_t data_end_absolute   = data_begin;
+
+  if (data_begin != std::string::npos) {
+
+    size_t data_end = line.find_last_not_of (white_space) + 1;
+    data_end = (data_end == std::string::npos) ? line.size() : data_end;
+    //     data_end_absolute = data_end + line_begin;
+
+    if (line.find ('{', data_begin) != std::string::npos) {
+
+      size_t brace_count = 1;
+      size_t brace = line.find ('{', data_begin);  // start from the first opening brace
+
+      while (brace_count > 0) {
+
+        // find the matching closing brace
+        brace = line.find_first_of ("{}", brace+1);
+        while (brace == std::string::npos) {
+          // add a new line
+          if (line_end >= conf.size()) {
+            cvm::fatal_error ("Parse error: reached the end while "
+                              "looking for closing brace; until now "
+                              "the following was parsed: \"\n"+
+                              line+"\".\n");
+            return false;
+          }
+          size_t const old_end = line.size();
+          //           data_end_absolute += old_end+1;
+
+          line_begin = line_end;
+          nl = conf.find ('\n', line_begin+1);
+          if (nl == std::string::npos) 
+            line_end = conf.size();
+          else 
+            line_end = nl;
+          line.append (conf, line_begin, (line_end-line_begin));
+
+          brace = line.find_first_of ("{}", old_end);
+        }
+
+        if (line[brace] == '{') brace_count++;
+        if (line[brace] == '}') brace_count--;
+      }
+
+      // set data_begin afterthe opening brace
+      data_begin = line.find_first_of ('{', data_begin) + 1;
+      data_begin = line.find_first_not_of (white_space,
+                                           data_begin);
+      // set data_end before the closing brace
+      data_end = brace;
+      data_end = line.find_last_not_of (white_space+"}",
+                                        data_end) + 1;
+      //       data_end_absolute = line_end;
+
+      if (data_end > line.size())
+        data_end = line.size();
+    }
+
+    data.append (line, data_begin, (data_end-data_begin));
+
+    if (data.size() && save_delimiters) {
+      data_begin_pos.push_back (conf.find (data, pos+key.size()));
+      data_end_pos.push_back   (data_begin_pos.back()+data.size());
+      //       std::cerr << "key = " << key << ", data = \""
+      //                 << data << "\", data_begin, data_end = "
+      //                 << data_begin_pos.back() << ", " << data_end_pos.back()
+      //                 << "\n";
+    }
+  }
+  
+  save_pos = line_end;
+
+  return true;
+}
+
+
+std::istream & operator>> (std::istream &is, colvarparse::read_block const &rb)
+{
+  size_t start_pos = is.tellg();
+  std::string read_key, next;
+
+  if ( !(is >> read_key) || !(read_key == rb.key) ||
+       !(is >> next) ) {
+    // the requested keyword has not been found, or it is not possible
+    // to read data after it
+    is.clear();
+    is.seekg (start_pos, std::ios::beg);
+    is.setstate (std::ios::failbit);
+    return is;
+  }
+
+  if (next != "{") {
+    (*rb.data) = next;
+    return is;
+  }
+
+  size_t brace_count = 1;
+  std::string line;
+  while (colvarparse::getline_nocomments (is, line)) {
+    size_t br = 0, br_old;
+    while ( (br = line.find_first_of ("{}", br)) != std::string::npos) {
+      if (line[br] == '{') brace_count++;
+      if (line[br] == '}') brace_count--;
+      br_old = br;
+      br++;
+    }
+    if (brace_count) (*rb.data).append (line + "\n");
+    else {
+      (*rb.data).append (line, 0, br_old);
+      break;
+    }
+  }
+  if (brace_count)  {
+    // end-of-file reached
+    // restore initial position
+    is.clear();
+    is.seekg (start_pos, std::ios::beg);
+    is.setstate (std::ios::failbit);
+  }
+  return is;
+}
+
+
+bool colvarparse::brace_check (std::string const &conf,
+                               size_t const start_pos)
+{
+  size_t brace_count = 0;
+  size_t brace = start_pos;
+  while ( (brace = conf.find_first_of ("{}", brace)) != std::string::npos) {
+    if (conf[brace] == '{') brace_count++;
+    if (conf[brace] == '}') brace_count--;
+    brace++;
+  }
+
+  if (brace_count != 0)
+    return false;
+  else
+    return true;
+}
+
+
+// Emacs
+// Local Variables:
+// mode: C++
+// End:

lib/colvars/colvarparse.h

@@ -0,0 +1,201 @@
+#ifndef COLVARPARSE_H
+#define COLVARPARSE_H
+
+#include <cstring>
+#include <string>
+
+#include "colvarmodule.h"
+#include "colvarvalue.h"
+
+
+/// \file colvarparse.h Parsing functions for collective variables
+
+
+/// \brief Base class containing parsing functions; all objects which
+/// need to parse input inherit from this
+class colvarparse {
+
+protected:
+
+  /// \brief List of legal keywords for this object: this is updated
+  /// by each call to colvarparse::get_keyval() or
+  /// colvarparse::key_lookup()
+  std::list<std::string> allowed_keywords;
+
+  /// \brief List of delimiters for the values of each keyword in the
+  /// configuration string; all keywords will be stripped of their
+  /// values before the keyword check is performed
+  std::list<size_t>      data_begin_pos;
+
+  /// \brief List of delimiters for the values of each keyword in the
+  /// configuration string; all keywords will be stripped of their
+  /// values before the keyword check is performed
+  std::list<size_t>      data_end_pos;
+
+  /// \brief Whether or not to accumulate data_begin_pos and
+  /// data_end_pos in key_lookup(); it may be useful to disable
+  /// this after the constructor is called, because other files may be
+  /// read (e.g. restart) that would mess up with the registry; in any
+  /// case, nothing serious happens until check_keywords() is invoked
+  /// (which should happen only right after construction)
+  bool save_delimiters;
+
+  /// \brief Add a new valid keyword to the list
+  void add_keyword (char const *key);
+
+  /// \brief Remove all the values from the config string
+  void strip_values (std::string &conf);
+
+public:
+
+  inline colvarparse()
+    : save_delimiters (true)
+  {}
+
+  /// How a keyword is parsed in a string
+  enum Parse_Mode { 
+    /// \brief (default) Read the first instance of a keyword (if
+    /// any), report its value, and print a warning when there is more
+    /// than one
+    parse_normal, 
+    /// \brief Like parse_normal, but don't send any message to the log
+    /// (useful e.g. in restart files when such messages are very
+    /// numerous and redundant)
+    parse_silent 
+  };
+
+  /// \fn get_keyval bool const get_keyval (std::string const &conf,
+  /// char const *key, _type_ &value, _type_ const &def_value,
+  /// Parse_Mode const parse_mode) \brief Helper function to parse
+  /// keywords in the configuration and get their values
+  ///
+  /// In normal circumstances, you should use either version the
+  /// get_keyval function.  Both of them look for the C string "key"
+  /// in the C++ string "conf", and assign the corresponding value (if
+  /// available) to the variable "value" (first version), or assign as
+  /// many values as found to the vector "values" (second version).
+  ///
+  /// If "key" is found but no value is associated to it, the default
+  /// value is provided (either one copy or as many copies as the
+  /// current length of the vector "values" specifies).  A message
+  /// will print, unless parse_mode is equal to parse_silent.  The
+  /// return value of both forms of get_keyval is true if "key" is
+  /// found (with or without value), and false when "key" is absent in
+  /// the string "conf".  If there is more than one instance of the
+  /// keyword, a warning will be raised; instead, to loop over
+  /// multiple instances key_lookup() should be invoked directly.
+  ///
+  /// If you introduce a new data type, add two new instances of this
+  /// functions, or insert this type in the \link colvarvalue \endlink
+  /// wrapper class (colvarvalue.h).
+
+#define _get_keyval_scalar_proto_(_type_,_def_value_)           \
+  bool get_keyval (std::string const &conf,                     \
+                   char const *key,                             \
+                   _type_ &value,                               \
+                   _type_ const &def_value = _def_value_,       \
+                   Parse_Mode const parse_mode = parse_normal)
+
+    _get_keyval_scalar_proto_ (int, (int)0);
+    _get_keyval_scalar_proto_ (size_t, (size_t)0);
+    _get_keyval_scalar_proto_ (std::string, std::string (""));
+    _get_keyval_scalar_proto_ (cvm::real, (cvm::real)0.0);
+    _get_keyval_scalar_proto_ (cvm::rvector, cvm::rvector());
+    _get_keyval_scalar_proto_ (cvm::quaternion, cvm::quaternion());
+    _get_keyval_scalar_proto_ (colvarvalue, colvarvalue (colvarvalue::type_notset));
+    _get_keyval_scalar_proto_ (bool, false);
+
+#define _get_keyval_vector_proto_(_type_,_def_value_)                   \
+  bool get_keyval (std::string const &conf,                             \
+                   char const *key,                                     \
+                   std::vector<_type_> &values,                         \
+                   std::vector<_type_> const &def_values =              \
+                   std::vector<_type_> (0, static_cast<_type_>(_def_value_)),                \
+                   Parse_Mode const parse_mode = parse_normal)
+  
+    _get_keyval_vector_proto_ (int, 0);
+    _get_keyval_vector_proto_ (size_t, 0);
+    _get_keyval_vector_proto_ (std::string, std::string (""));
+    _get_keyval_vector_proto_ (cvm::real, 0.0);
+    _get_keyval_vector_proto_ (cvm::rvector, cvm::rvector());
+    _get_keyval_vector_proto_ (cvm::quaternion, cvm::quaternion());
+    _get_keyval_vector_proto_ (colvarvalue, colvarvalue (colvarvalue::type_notset));
+
+
+  /// \brief Check that all the keywords within "conf" are in the list
+  /// of allowed keywords; this will invoke strip_values() first and
+  /// then loop over all words
+  void check_keywords (std::string &conf, char const *key);
+
+
+  /// \brief Return a lowercased copy of the string
+  static inline std::string to_lower_cppstr (std::string const &in)
+  {
+    std::string out = "";
+    for (size_t i = 0; i < in.size(); i++) {
+      out.append (1, (char) ::tolower (in[i]) );
+    }
+    return out;
+  }
+
+  /// \brief Helper class to read a block of the type "key { ... }"
+  /// from a stream and store it in a string
+  ///
+  /// Useful on restarts, where the file is too big to be loaded in a
+  /// string by key_lookup; it can only check that the keyword is
+  /// correct and the block is properly delimited by braces, not
+  /// skipping other blocks
+  class read_block {
+
+    std::string const   key;
+    std::string * const data;
+
+  public:
+    inline read_block (std::string const &key_in, std::string &data_in)
+      : key (key_in), data (&data_in)
+    {}
+    inline ~read_block() {}
+    friend std::istream & operator >> (std::istream &is, read_block const &rb);
+  };
+
+
+  /// Accepted white space delimiters, used in key_lookup()
+  static std::string const white_space;
+
+  /// \brief Low-level function for parsing configuration strings;
+  /// automatically adds the requested keywords to the list of valid
+  /// ones.  \param conf the content of the configuration file or one
+  /// of its blocks \param key the keyword to search in "conf" \param
+  /// data (optional) holds the string provided after "key", if any
+  /// \param save_pos (optional) stores the position of the keyword
+  /// within "conf", useful when doing multiple calls \param
+  /// save_delimiters (optional) 
+  bool key_lookup (std::string const &conf,
+                   char const *key,
+                   std::string &data = dummy_string,
+                   size_t &save_pos = dummy_pos);
+
+  /// Used as a default argument by key_lookup
+  static std::string dummy_string;
+  /// Used as a default argument by key_lookup
+  static size_t dummy_pos;
+
+  /// \brief Works as std::getline() but also removes everything
+  /// between a comment character and the following newline
+  static std::istream & getline_nocomments (std::istream &is,
+                                            std::string &s,
+                                            char const delim = '\n');
+
+  /// Check if the content of the file has matching braces
+  bool brace_check (std::string const &conf,
+                    size_t const start_pos = 0);
+
+};
+
+
+#endif
+
+// Emacs
+// Local Variables:
+// mode: C++
+// End:

lib/colvars/colvarproxy.h

@@ -0,0 +1,137 @@
+#ifndef COLVARPROXY_H
+#define COLVARPROXY_H
+
+
+
+#include "colvarmodule.h"
+
+
+/// \brief Interface class between the collective variables module and
+/// the simulation program
+
+class colvarproxy {
+
+public:
+
+  /// Pointer to the instance of colvarmodule
+  colvarmodule *colvars;
+
+  /// \brief Value of the unit for atomic coordinates with respect to
+  /// angstroms (used by some variables for hard-coded default values)
+  virtual cvm::real unit_angstrom() = 0;
+
+  /// \brief Boltzmann constant
+  virtual cvm::real boltzmann() = 0;
+
+  /// \brief Temperature of the simulation (K)
+  virtual cvm::real temperature() = 0;
+
+  /// \brief Time step of the simulation (fs)
+  virtual cvm::real dt() = 0;
+
+  /// Pass restraint energy value for current timestep to MD engine
+  virtual void add_energy (cvm::real energy) = 0;
+
+  /// Tell the proxy whether system forces are needed
+  virtual void request_system_force (bool yesno) = 0;
+
+  /// Print a message to the main log
+  virtual void log (std::string const &message) = 0;
+
+  /// Print a message to the main log and exit with error code
+  virtual void fatal_error (std::string const &message) = 0;
+
+  /// Print a message to the main log and exit normally
+  virtual void exit (std::string const &message) = 0;
+
+
+  /// \brief Prefix to be used for input files (restarts, not
+  /// configuration)
+  virtual std::string input_prefix() = 0;
+
+  /// \brief Prefix to be used for output restart files
+  virtual std::string restart_output_prefix() = 0;
+
+  /// \brief Prefix to be used for output files (final system
+  /// configuration)
+  virtual std::string output_prefix() = 0;
+
+  /// \brief Restarts will be fritten each time this number of steps has passed
+  virtual size_t restart_frequency() = 0;
+
+  // **************** PERIODIC BOUNDARY CONDITIONS ****************
+
+  /// \brief Get the simple distance vector between two positions
+  /// (with periodic boundary conditions handled transparently)
+  virtual cvm::rvector position_distance (cvm::atom_pos const &pos1,
+                                          cvm::atom_pos const &pos2) = 0;
+
+  /// \brief Get the square distance between two positions (with
+  /// periodic boundary conditions handled transparently)
+  ///
+  /// Note: in the case of periodic boundary conditions, this provides
+  /// an analytical square distance (while taking the square of
+  /// position_distance() would produce leads to a cusp)
+  virtual cvm::real position_dist2 (cvm::atom_pos const &pos1,
+                                    cvm::atom_pos const &pos2) = 0;
+
+  /// \brief Get the closest periodic image to a reference position
+  /// \param pos The position to look for the closest periodic image
+  /// \param ref_pos The reference position 
+  virtual void select_closest_image (cvm::atom_pos &pos,
+                                     cvm::atom_pos const &ref_pos) = 0;
+
+  /// \brief Perform select_closest_image() on a set of atomic positions
+  ///
+  /// After that, distance vectors can then be calculated directly,
+  /// without using position_distance()
+  void select_closest_images (std::vector<cvm::atom_pos> &pos,
+                              cvm::atom_pos const &ref_pos);
+
+
+  /// \brief Read atom identifiers from a file \param filename name of
+  /// the file (usually a PDB) \param atoms array to which atoms read
+  /// from "filename" will be appended \param pdb_field (optiona) if
+  /// "filename" is a PDB file, use this field to determine which are
+  /// the atoms to be set
+  virtual void load_atoms (char const *filename,
+                           std::vector<cvm::atom> &atoms,
+                           std::string const pdb_field,
+                           double const pdb_field_value = 0.0) = 0;
+
+  /// \brief Load the coordinates for a group of atoms from a file
+  /// (usually a PDB); if "pos" is already allocated, the number of its
+  /// elements must match the number of atoms in "filename"
+  virtual void load_coords (char const *filename,
+                            std::vector<cvm::atom_pos> &pos,
+                            const std::vector<int> &indices,
+                            std::string const pdb_field,
+                            double const pdb_field_value = 0.0) = 0;
+
+  /// \brief Rename the given file, under the convention provided by
+  /// the MD program
+  virtual void backup_file (char const *filename) = 0;
+
+  /// \brief Pseudo-random number with Gaussian distribution
+  virtual cvm::real rand_gaussian (void) = 0;
+
+  virtual inline ~colvarproxy() {}
+};
+
+
+inline void colvarproxy::select_closest_images (std::vector<cvm::atom_pos> &pos,
+                                                cvm::atom_pos const &ref_pos)
+{
+  for (std::vector<cvm::atom_pos>::iterator pi = pos.begin();
+       pi != pos.end(); pi++) {
+    select_closest_image (*pi, ref_pos);
+  }
+}
+
+#endif
+
+
+// Emacs
+// Local Variables:
+// mode: C++
+// End:

lib/colvars/colvarproxy_namd.C

@@ -0,0 +1,613 @@
+#include "common.h"
+#include "BackEnd.h"
+#include "InfoStream.h"
+#include "Node.h"
+#include "Molecule.h"
+#include "PDB.h"
+#include "PDBData.h"
+#include "ReductionMgr.h"
+
+#include "colvarmodule.h"
+#include "colvaratoms.h"
+#include "colvarproxy.h"
+#include "colvarproxy_namd.h"
+
+
+colvarproxy_namd::colvarproxy_namd()
+{
+  first_timestep = true;
+  system_force_requested = false;
+
+  // initialize pointers to NAMD configuration data
+  simparams = Node::Object()->simParameters;
+  lattice = &(simparams->lattice);
+
+  if (cvm::debug())
+    iout << "Info: initializing the colvars proxy object.\n" << endi;
+
+  // find the configuration file
+  StringList *config = Node::Object()->configList->find ("colvarsConfig");
+  if (!config)
+    NAMD_die ("No configuration file for collective variables: exiting.\n");
+
+  // find the input state file
+  StringList *input_restart = Node::Object()->configList->find ("colvarsInput");
+  input_prefix_str = std::string (input_restart ? input_restart->data : "");
+  if (input_prefix_str.rfind (".colvars.state") != std::string::npos) {
+    // strip the extension, if present
+    input_prefix_str.erase (input_prefix_str.rfind (".colvars.state"),
+                            std::string (".colvars.state").size());
+  }
+
+  // get the thermostat temperature
+  if (simparams->rescaleFreq > 0)
+    thermostat_temperature = simparams->rescaleTemp;
+  else if (simparams->reassignFreq > 0)
+    thermostat_temperature = simparams->reassignTemp;
+  else if (simparams->langevinOn)
+    thermostat_temperature = simparams->langevinTemp;
+  else if (simparams->tCoupleOn)
+    thermostat_temperature = simparams->tCoupleTemp;
+  //else if (simparams->loweAndersenOn)
+  //  thermostat_temperature = simparams->loweAndersenTemp;
+  else 
+    thermostat_temperature = 0.0;
+
+  random = Random (simparams->randomSeed);
+
+  // take the output prefixes from the namd input
+  output_prefix_str = std::string (simparams->outputFilename);
+  restart_output_prefix_str = std::string (simparams->restartFilename);
+  restart_frequency_s = simparams->restartFrequency;
+
+  // initiate the colvarmodule, this object will be the communication
+  // proxy
+  colvars = new colvarmodule (config->data, this);
+
+  if (cvm::debug()) {
+    cvm::log ("colvars_atoms = "+cvm::to_str (colvars_atoms)+"\n");
+    cvm::log ("colvars_atoms_ncopies = "+cvm::to_str (colvars_atoms_ncopies)+"\n");
+    cvm::log ("positions = "+cvm::to_str (positions)+"\n");
+    cvm::log ("total_forces = "+cvm::to_str (total_forces)+"\n");
+    cvm::log ("applied_forces = "+cvm::to_str (applied_forces)+"\n");
+    cvm::log (cvm::line_marker);
+  }
+
+  // Initialize reduction object to submit restraint energy as MISC
+  reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
+
+  if (cvm::debug())
+    iout << "Info: done initializing the colvars proxy object.\n" << endi;
+}
+
+
+colvarproxy_namd::~colvarproxy_namd()
+{
+  delete reduction;
+  if (colvars != NULL) {
+    delete colvars;
+    colvars = NULL;
+  }
+}
+
+// Reimplemented function from GlobalMaster
+void colvarproxy_namd::calculate()
+{
+
+  if (first_timestep) {
+    first_timestep = false;
+  } else {
+    // Use the time step number inherited from GlobalMaster
+    if ( step - previous_NAMD_step == 1 ) {
+      colvars->it++;
+    }
+    // Other cases could mean:
+    // - run 0
+    // - beginning of a new run statement
+    // then the internal counter should not be incremented
+  }
+  previous_NAMD_step = step;
+
+  if (cvm::debug()) {
+    cvm::log (cvm::line_marker+
+              "colvarproxy_namd, step no. "+cvm::to_str (colvars->it)+"\n"+
+              "Updating internal data.\n");
+  }
+
+  // must delete the forces applied at the previous step: they have
+  // already been used and copied to other memory locations
+  modifyForcedAtoms().resize (0);
+  modifyAppliedForces().resize (0);
+
+  // prepare the local arrays to contain the sorted copies of the NAMD
+  // arrays
+  for (size_t i = 0; i < colvars_atoms.size(); i++) {
+    positions[i] = cvm::rvector (0.0, 0.0, 0.0);
+    total_forces[i] = cvm::rvector (0.0, 0.0, 0.0);
+    applied_forces[i] = cvm::rvector (0.0, 0.0, 0.0);
+  }
+
+  // sort the positions array
+  for (size_t i = 0; i < colvars_atoms.size(); i++) {
+    bool found_position = false;
+    AtomIDList::const_iterator a_i = this->getAtomIdBegin();
+    AtomIDList::const_iterator a_e = this->getAtomIdEnd();
+    PositionList::const_iterator p_i = this->getAtomPositionBegin();
+    for ( ; a_i != a_e; ++a_i, ++p_i ) {
+      if ( *a_i == colvars_atoms[i] ) {
+        found_position = true;
+        Position const &namd_pos = *p_i;
+        positions[i] = cvm::rvector (namd_pos.x, namd_pos.y, namd_pos.z);
+        break;
+      }
+    }
+    if (!found_position)
+      cvm::fatal_error ("Error: cannot find the position of atom "+
+                        cvm::to_str (colvars_atoms[i]+1)+"\n");
+  }
+
+
+  if (system_force_requested && cvm::step_relative() > 0) {
+
+    // sort the array of total forces from the previous step (but only
+    // do it if there *is* a previous step!)
+    for (size_t i = 0; i < colvars_atoms.size(); i++) {
+      bool found_total_force = false;
+      //found_total_force = false;
+      AtomIDList::const_iterator a_i = this->getForceIdBegin();
+      AtomIDList::const_iterator a_e = this->getForceIdEnd();
+      PositionList::const_iterator f_i = this->getTotalForce();
+      for ( ; a_i != a_e; ++a_i, ++f_i ) {
+        if ( *a_i == colvars_atoms[i] ) {
+          found_total_force = true;
+          Vector const &namd_force = *f_i;
+          total_forces[i] = cvm::rvector (namd_force.x, namd_force.y, namd_force.z);
+          //           if (cvm::debug()) 
+          //             cvm::log ("Found the total force of atom "+
+          //                       cvm::to_str (colvars_atoms[i]+1)+", which is "+
+          //                       cvm::to_str (total_forces[i])+".\n");
+          break;
+        }
+      }
+      if (!found_total_force)
+        cvm::fatal_error ("Error: system forces were requested, but total force on atom "+
+              cvm::to_str (colvars_atoms[i]+1) + " was not\n"
+              "found. The most probable cause is combination of energy minimization with a\n"
+              "biasing method that requires MD (e.g. ABF). Always run minimization\n"
+              "and ABF separately.");
+    }
+
+    // do the same for applied forces
+    for (size_t i = 0; i < colvars_atoms.size(); i++) {
+      AtomIDList::const_iterator a_i = this->getLastAtomsForcedBegin();
+      AtomIDList::const_iterator a_e = this->getLastAtomsForcedEnd();
+      PositionList::const_iterator f_i = this->getLastForcesBegin();
+      for ( ; a_i != a_e; ++a_i, ++f_i ) {
+        if ( *a_i == colvars_atoms[i] ) {
+          Vector const &namd_force = *f_i;
+          if (cvm::debug())
+            cvm::log ("Found a force applied to atom "+
+                      cvm::to_str (colvars_atoms[i]+1)+": "+
+                      cvm::to_str (cvm::rvector (namd_force.x, namd_force.y, namd_force.z))+
+                      "; current total is "+
+                      cvm::to_str (applied_forces[i])+".\n");
+          applied_forces[i] += cvm::rvector (namd_force.x, namd_force.y, namd_force.z);
+        }
+      }
+    }
+  }
+
+  // call the collective variable module
+  colvars->calc();
+  // send MISC energy
+  reduction->submit();
+
+  // NAMD does not destruct GlobalMaster objects, so we must remember
+  // to write all output files at the end of the run
+  if (step == simparams->N) {
+    colvars->write_output_files();
+  }
+}
+
+
+void colvarproxy_namd::add_energy (cvm::real energy)
+{
+  reduction->item(REDUCTION_MISC_ENERGY) += energy;
+}
+
+void colvarproxy_namd::request_system_force (bool yesno)
+{
+  system_force_requested = yesno;
+}
+
+void colvarproxy_namd::log (std::string const &message)
+{
+  std::istringstream is (message);
+  std::string line;
+  while (std::getline (is, line))
+    iout << "colvars: " << line << "\n";
+  iout << endi;
+}
+
+
+void colvarproxy_namd::fatal_error (std::string const &message)
+{
+  cvm::log (message);
+  if (!cvm::debug())
+    cvm::log ("If this error message is unclear, "
+              "try recompiling with -DCOLVARS_DEBUG.\n");
+  NAMD_die ("Error in the collective variables module: exiting.\n");
+}
+
+
+void colvarproxy_namd::exit (std::string const &message)
+{
+  cvm::log (message);
+  BackEnd::exit();
+}
+
+
+enum e_pdb_field {
+  e_pdb_none,
+  e_pdb_occ,
+  e_pdb_beta,
+  e_pdb_x,
+  e_pdb_y,
+  e_pdb_z,
+  e_pdb_ntot
+};
+
+
+e_pdb_field pdb_field_str2enum (std::string const &pdb_field_str)
+{
+  e_pdb_field pdb_field = e_pdb_none;
+
+  if (colvarparse::to_lower_cppstr (pdb_field_str) ==
+      colvarparse::to_lower_cppstr ("O")) {
+    pdb_field = e_pdb_occ;
+  }
+
+  if (colvarparse::to_lower_cppstr (pdb_field_str) ==
+      colvarparse::to_lower_cppstr ("B")) {
+    pdb_field = e_pdb_beta;
+  }
+
+  if (colvarparse::to_lower_cppstr (pdb_field_str) ==
+      colvarparse::to_lower_cppstr ("X")) {
+    pdb_field = e_pdb_x;
+  }
+  
+  if (colvarparse::to_lower_cppstr (pdb_field_str) ==
+      colvarparse::to_lower_cppstr ("Y")) {
+    pdb_field = e_pdb_y;
+  }
+
+  if (colvarparse::to_lower_cppstr (pdb_field_str) ==
+      colvarparse::to_lower_cppstr ("Z")) {
+    pdb_field = e_pdb_z;
+  }
+
+  if (pdb_field == e_pdb_none) {
+    cvm::fatal_error ("Error: unsupported PDB field, \""+
+                      pdb_field_str+"\".\n");
+  }
+
+  return pdb_field;
+}
+
+
+void colvarproxy_namd::load_coords (char const *pdb_filename,
+                                    std::vector<cvm::atom_pos> &pos,
+                                    const std::vector<int> &indices,
+                                    std::string const pdb_field_str,
+                                    double const pdb_field_value)
+{
+  if (pdb_field_str.size() == 0 && indices.size() == 0) {
+    cvm::fatal_error ("Bug alert: either PDB field should be defined or list of "
+                      "atom IDs should be available when loading atom coordinates!\n");
+  }
+
+  e_pdb_field pdb_field_index;
+  bool const use_pdb_field = (pdb_field_str.size() > 0);
+  if (use_pdb_field) {
+    pdb_field_index = pdb_field_str2enum (pdb_field_str);
+  }
+
+  // next index to be looked up in PDB file (if list is supplied)
+  std::vector<int>::const_iterator current_index = indices.begin();
+
+  PDB *pdb = new PDB (pdb_filename);
+  size_t const pdb_natoms = pdb->num_atoms();
+  
+  if (pos.size() != pdb_natoms) {
+
+    bool const pos_allocated = (pos.size() > 0);
+
+    size_t ipos = 0, ipdb = 0;
+    for ( ; ipdb < pdb_natoms; ipdb++) {
+
+      if (use_pdb_field) {
+        // PDB field mode: skip atoms with wrong value in PDB field
+        double atom_pdb_field_value = 0.0;
+
+        switch (pdb_field_index) {
+        case e_pdb_occ:
+          atom_pdb_field_value = (pdb->atom (ipdb))->occupancy();
+          break;
+        case e_pdb_beta:
+          atom_pdb_field_value = (pdb->atom (ipdb))->temperaturefactor();
+          break;
+        case e_pdb_x:
+          atom_pdb_field_value = (pdb->atom (ipdb))->xcoor();
+          break;
+        case e_pdb_y:
+          atom_pdb_field_value = (pdb->atom (ipdb))->ycoor();
+          break;
+        case e_pdb_z:
+          atom_pdb_field_value = (pdb->atom (ipdb))->zcoor();
+          break;
+        default:
+          break;
+        }
+
+        if ( (pdb_field_value) &&
+             (atom_pdb_field_value != pdb_field_value) ) {
+          continue;
+        } else if (atom_pdb_field_value == 0.0) {
+          continue;
+        }
+
+      } else {
+        // Atom ID mode: use predefined atom IDs from the atom group
+        if (ipdb != *current_index) {
+          // Skip atoms not in the list
+          continue;
+        } else {
+          current_index++;
+        }
+      }
+      
+      if (!pos_allocated) {
+        pos.push_back (cvm::atom_pos (0.0, 0.0, 0.0));
+      } else if (ipos >= pos.size()) {
+        cvm::fatal_error ("Error: the PDB file \""+
+                          std::string (pdb_filename)+
+                          "\" contains coordinates for "
+                          "more atoms than needed.\n");
+      }
+
+      pos[ipos] = cvm::atom_pos ((pdb->atom (ipdb))->xcoor(),
+                                 (pdb->atom (ipdb))->ycoor(),
+                                 (pdb->atom (ipdb))->zcoor());
+      ipos++;
+      if (!use_pdb_field && current_index == indices.end())
+        break;
+    }
+
+    if (ipos < pos.size())
+      cvm::fatal_error ("Error: the PDB file \""+
+                        std::string (pdb_filename)+
+                        "\" contains coordinates for only "+
+                        cvm::to_str (ipos)+
+                        " atoms, but "+cvm::to_str (pos.size())+
+                        " are needed.\n");
+
+    if (current_index != indices.end())
+      cvm::fatal_error ("Error: not all atoms found in PDB file.\n");
+
+  } else {
+
+    // when the PDB contains exactly the number of atoms of the array,
+    // ignore the fields and just read coordinates
+    for (size_t ia = 0; ia < pos.size(); ia++) {
+      pos[ia] = cvm::atom_pos ((pdb->atom (ia))->xcoor(),
+                               (pdb->atom (ia))->ycoor(),
+                               (pdb->atom (ia))->zcoor());
+    }
+  }
+
+  delete pdb;
+}
+
+
+void colvarproxy_namd::load_atoms (char const *pdb_filename,
+                                   std::vector<cvm::atom> &atoms,
+                                   std::string const pdb_field_str,
+                                   double const pdb_field_value)
+{
+  if (pdb_field_str.size() == 0)
+    cvm::fatal_error ("Error: must define which PDB field to use "
+                      "in order to define atoms from a PDB file.\n");
+
+  PDB *pdb = new PDB (pdb_filename);
+  size_t const pdb_natoms = pdb->num_atoms();
+
+  e_pdb_field pdb_field_index = pdb_field_str2enum (pdb_field_str);
+
+  for (size_t ipdb = 0; ipdb < pdb_natoms; ipdb++) {
+
+    double atom_pdb_field_value = 0.0;
+
+    switch (pdb_field_index) {
+    case e_pdb_occ:
+      atom_pdb_field_value = (pdb->atom (ipdb))->occupancy();
+      break;
+    case e_pdb_beta:
+      atom_pdb_field_value = (pdb->atom (ipdb))->temperaturefactor();
+      break;
+    case e_pdb_x:
+      atom_pdb_field_value = (pdb->atom (ipdb))->xcoor();
+      break;
+    case e_pdb_y:
+      atom_pdb_field_value = (pdb->atom (ipdb))->ycoor();
+      break;
+    case e_pdb_z:
+      atom_pdb_field_value = (pdb->atom (ipdb))->zcoor();
+      break;
+    default:
+      break;
+    }
+
+    if ( (pdb_field_value) &&
+         (atom_pdb_field_value != pdb_field_value) ) {
+      continue;
+    } else if (atom_pdb_field_value == 0.0) {
+      continue;
+    }
+     
+    atoms.push_back (cvm::atom (ipdb+1));
+  }
+
+  delete pdb;
+}
+
+
+void colvarproxy_namd::backup_file (char const *filename)
+{
+  if (std::string (filename).rfind (std::string (".colvars.state")) != std::string::npos) {
+    NAMD_backup_file (filename, ".old");
+  } else {
+    NAMD_backup_file (filename, ".BAK");
+  }
+}
+
+
+size_t colvarproxy_namd::init_namd_atom (AtomID const &aid)
+{
+  modifyRequestedAtoms().add (aid);
+  for (size_t i = 0; i < colvars_atoms.size(); i++) {
+    if (colvars_atoms[i] == aid) {
+      // this atom id was already recorded
+      colvars_atoms_ncopies[i] += 1;
+      return i;
+    }
+  }
+
+  // allocate a new slot for this atom
+  colvars_atoms_ncopies.push_back (1);
+  colvars_atoms.push_back (aid);
+  positions.push_back (cvm::rvector());
+  total_forces.push_back (cvm::rvector());
+  applied_forces.push_back (cvm::rvector());
+
+  return (colvars_atoms.size()-1);
+}
+
+// atom member functions, NAMD specific implementations
+
+cvm::atom::atom (int const &atom_number)
+{
+  // NAMD internal numbering starts from zero
+  AtomID const aid (atom_number-1);
+
+  if (cvm::debug())
+    cvm::log ("Adding atom "+cvm::to_str (aid+1)+
+              " for collective variables calculation.\n");
+
+  if ( (aid < 0) || (aid >= Node::Object()->molecule->numAtoms) ) 
+    cvm::fatal_error ("Error: invalid atom number specified, "+
+                      cvm::to_str (atom_number)+"\n");
+  this->index = ((colvarproxy_namd *) cvm::proxy)->init_namd_atom (aid);
+  if (cvm::debug())
+    cvm::log ("The index of this atom in the colvarproxy_namd arrays is "+
+              cvm::to_str (this->index)+".\n");
+  this->id = aid;
+  this->mass = Node::Object()->molecule->atommass (aid);
+  this->reset_data();
+}
+
+
+/// For AMBER topologies, the segment id is automatically set to
+/// "MAIN" (the segment id assigned by NAMD's AMBER topology parser),
+/// and is therefore optional when an AMBER topology is used
+cvm::atom::atom (cvm::residue_id const &residue,
+                 std::string const     &atom_name,
+                 std::string const     &segment_id)
+{
+  AtomID const aid =
+    (segment_id.size() ? 
+       Node::Object()->molecule->get_atom_from_name (segment_id.c_str(),
+                                                     residue,
+                                                     atom_name.c_str()) :
+     Node::Object()->molecule->get_atom_from_name ("MAIN",
+                                                   residue,
+                                                   atom_name.c_str()));
+    
+
+  if (cvm::debug())
+    cvm::log ("Adding atom \""+
+              atom_name+"\" in residue "+
+              cvm::to_str (residue)+
+              " (index "+cvm::to_str (aid)+
+              ") for collective variables calculation.\n");
+
+  if (aid < 0) {
+    // get_atom_from_name() has returned an error value
+    cvm::fatal_error ("Error: could not find atom \""+
+                      atom_name+"\" in residue "+
+                      cvm::to_str (residue)+
+                      ( (segment_id != "MAIN") ?
+                        (", segment \""+segment_id+"\"") :
+                        ("") )+
+                      "\n");
+  }
+
+  this->index = ((colvarproxy_namd *) cvm::proxy)->init_namd_atom (aid);
+  if (cvm::debug())
+    cvm::log ("The index of this atom in the colvarproxy_namd arrays is "+
+              cvm::to_str (this->index)+".\n");
+  this->id = aid;
+  this->mass = Node::Object()->molecule->atommass (aid);
+  this->reset_data();
+}
+
+
+// copy constructor
+cvm::atom::atom (cvm::atom const &a)
+  : index (a.index), id (a.id), mass (a.mass)
+{
+  // init_namd_atom() has already been called by a's constructor, no
+  // need to call it again
+
+  // need to increment the counter anyway
+  colvarproxy_namd *gm = (colvarproxy_namd *) cvm::proxy;
+  gm->colvars_atoms_ncopies[this->index] += 1;
+}
+
+
+cvm::atom::~atom() 
+{
+  colvarproxy_namd *gm = (colvarproxy_namd *) cvm::proxy;
+  if (gm->colvars_atoms_ncopies[this->index] > 0)
+    gm->colvars_atoms_ncopies[this->index] -= 1;
+}
+
+
+void cvm::atom::read_position()
+{
+  colvarproxy_namd const * const gm = (colvarproxy_namd *) cvm::proxy;
+  this->pos = gm->positions[this->index];
+}
+
+
+void cvm::atom::read_velocity()
+{
+  cvm::fatal_error ("Error: NAMD does not have yet a way to communicate "
+                    "atom velocities to the colvars.\n");
+}
+
+
+void cvm::atom::read_system_force()
+{
+  colvarproxy_namd const * const gm = (colvarproxy_namd *) cvm::proxy;
+  this->system_force = gm->total_forces[this->index] - gm->applied_forces[this->index];
+}
+
+
+void cvm::atom::apply_force (cvm::rvector const &new_force)
+{
+  colvarproxy_namd *gm = (colvarproxy_namd *) cvm::proxy;
+  gm->modifyForcedAtoms().add (this->id);
+  gm->modifyAppliedForces().add (Vector (new_force.x, new_force.y, new_force.z));
+}
+

lib/colvars/colvarproxy_namd.h

@@ -0,0 +1,172 @@
+#ifndef COLVARPROXY_NAMD_H
+#define COLVARPROXY_NAMD_H
+
+#include "Vector.h"
+#include "ResizeArray.h"
+#include "NamdTypes.h"
+#include "SimParameters.h"
+#include "Lattice.h"
+#include "GlobalMaster.h"
+#include "Random.h"
+
+#include "colvarmodule.h"
+#include "colvarproxy.h"
+
+
+/// \brief Communication between colvars and NAMD (implementation of
+/// \link colvarproxy \endlink)
+class colvarproxy_namd : public colvarproxy, public GlobalMaster {
+
+protected:
+
+  /// Pointer to the NAMD simulation input object
+  SimParameters const *simparams;
+  /// Pointer to the NAMD boundary conditions object
+  Lattice const       *lattice;
+
+  BigReal thermostat_temperature;
+
+  /// NAMD-style PRNG object
+  Random random;
+
+  std::string input_prefix_str, output_prefix_str, restart_output_prefix_str;
+  size_t      restart_frequency_s;
+  size_t      previous_NAMD_step;
+  bool        first_timestep;
+  bool        system_force_requested;
+
+  std::vector<int>          colvars_atoms;
+  std::vector<size_t>       colvars_atoms_ncopies;
+  std::vector<cvm::rvector> positions;
+  std::vector<cvm::rvector> total_forces;
+  std::vector<cvm::rvector> applied_forces;
+
+  size_t init_namd_atom (AtomID const &aid);
+
+  SubmitReduction *reduction;
+public:
+
+  friend class cvm::atom;
+
+  colvarproxy_namd();
+  ~colvarproxy_namd();
+
+  /// \brief Reimplements GlobalMaster member function, to be called
+  /// at every step
+  void calculate();
+
+  void add_energy (cvm::real energy);
+  void request_system_force (bool yesno);
+  void log (std::string const &message);
+  void fatal_error (std::string const &message);
+  void exit (std::string const &message);
+
+  inline cvm::real unit_angstrom()
+  {
+    return 1.0;
+  }
+
+  cvm::real boltzmann()
+  {
+    return 0.001987191;
+  }
+
+  cvm::real temperature()
+  {
+    return thermostat_temperature;
+  }
+
+  cvm::real dt()
+  {
+    return simparams->dt;
+  }
+
+  inline std::string input_prefix()
+  {
+    return input_prefix_str;
+  }
+  inline std::string restart_output_prefix()
+  {
+    return restart_output_prefix_str;
+  }
+  inline std::string output_prefix()
+  {
+    return output_prefix_str;
+  }
+  inline size_t restart_frequency()
+  {
+    return restart_frequency_s;
+  }
+
+  cvm::rvector position_distance (cvm::atom_pos const &pos1,
+                                  cvm::atom_pos const &pos2);
+  cvm::real position_dist2 (cvm::atom_pos const &pos1,
+                            cvm::atom_pos const &pos2);
+
+  void select_closest_image (cvm::atom_pos &pos,
+                             cvm::atom_pos const &ref_pos);
+
+
+  void load_atoms (char const *filename,
+                   std::vector<cvm::atom> &atoms,
+                   std::string const pdb_field,
+                   double const pdb_field_value = 0.0);
+
+  void load_coords (char const *filename,
+                    std::vector<cvm::atom_pos> &pos,
+                    const std::vector<int> &indices,
+                    std::string const pdb_field,
+                    double const pdb_field_value = 0.0);
+
+  void backup_file (char const *filename);
+
+  cvm::real rand_gaussian (void)
+  {
+    return random.gaussian();
+  }
+};
+
+
+inline cvm::rvector colvarproxy_namd::position_distance (cvm::atom_pos const &pos1,
+                                                         cvm::atom_pos const &pos2)
+{
+  Position const p1 (pos1.x, pos1.y, pos1.z);
+  Position const p2 (pos2.x, pos2.y, pos2.z);
+  // return p2 - p1
+  Vector const d = this->lattice->delta (p2, p1);
+  return cvm::rvector (d.x, d.y, d.z);
+}
+
+
+inline void colvarproxy_namd::select_closest_image (cvm::atom_pos &pos,
+                                                    cvm::atom_pos const &ref_pos)
+{
+  Position const p (pos.x, pos.y, pos.z);
+  Position const rp (ref_pos.x, ref_pos.y, ref_pos.z);
+  ScaledPosition const srp = this->lattice->scale (rp);
+  Position const np = this->lattice->nearest (p, srp);
+  pos.x = np.x;
+  pos.y = np.y;
+  pos.z = np.z;
+}
+
+
+inline cvm::real colvarproxy_namd::position_dist2 (cvm::atom_pos const &pos1,
+                                                   cvm::atom_pos const &pos2)
+{
+  Lattice const *l = this->lattice;
+  Vector const p1 (pos1.x, pos1.y, pos1.z);
+  Vector const p2 (pos2.x, pos2.y, pos2.z);
+  Vector const d = l->delta (p1, p2);
+  return cvm::real (d.x*d.x + d.y*d.y + d.z*d.z);
+}
+
+
+
+#endif
+
+
+// Emacs
+// Local Variables:
+// mode: C++
+// End:

lib/colvars/colvarvalue.C

@@ -0,0 +1,261 @@
+#include <vector>
+
+#include "colvarmodule.h"
+#include "colvarvalue.h"
+
+
+
+std::string const colvarvalue::type_desc[colvarvalue::type_quaternion+1] =
+  { "not_set",
+    "scalar number",
+    "3-dimensional vector",
+    "3-dimensional unit vector",
+    "4-dimensional unit vector" };
+
+size_t const      colvarvalue::dof_num[  colvarvalue::type_quaternion+1] =
+  { 0, 1, 3, 2, 3 };
+
+
+void colvarvalue::undef_op() const
+{
+  cvm::fatal_error ("Error: Undefined operation on a colvar of type \""+
+                    colvarvalue::type_desc[this->value_type]+"\".\n");
+}
+
+void colvarvalue::error_rside
+(colvarvalue::Type const &vt) const
+{
+  cvm::fatal_error ("Trying to assign a colvar value with type \""+
+                    type_desc[this->value_type]+"\" to one with type \""+
+                    type_desc[vt]+"\".\n");
+}    
+
+void colvarvalue::error_lside
+(colvarvalue::Type const &vt) const
+{
+  cvm::fatal_error ("Trying to use a colvar value with type \""+
+                    type_desc[vt]+"\" as one of type \""+
+                    type_desc[this->value_type]+"\".\n");
+}    
+
+
+
+void colvarvalue::inner_opt (colvarvalue                        const &x,
+                             std::vector<colvarvalue>::iterator       &xv,
+                             std::vector<colvarvalue>::iterator const &xv_end,
+                             std::vector<cvm::real>::iterator         &inner)
+{
+  // doing type check only once, here
+  colvarvalue::check_types (x, *xv);
+
+  std::vector<colvarvalue>::iterator &xvi = xv;
+  std::vector<cvm::real>::iterator    &ii = inner;
+
+  switch (x.value_type) {
+  case colvarvalue::type_scalar:
+    while (xvi != xv_end) {
+      *(ii++) += (xvi++)->real_value * x.real_value;
+    }
+    break;
+  case colvarvalue::type_vector:
+  case colvarvalue::type_unitvector:
+    while (xvi != xv_end) {
+      *(ii++) += (xvi++)->rvector_value * x.rvector_value;
+    }
+    break;
+  case colvarvalue::type_quaternion:
+    while (xvi != xv_end) {
+      *(ii++) += ((xvi++)->quaternion_value).cosine (x.quaternion_value);
+    }
+    break;
+  default:
+    x.undef_op();
+  };
+}
+
+void colvarvalue::inner_opt (colvarvalue const                      &x,
+                             std::list<colvarvalue>::iterator       &xv,
+                             std::list<colvarvalue>::iterator const &xv_end,
+                             std::vector<cvm::real>::iterator       &inner)
+{
+  // doing type check only once, here
+  colvarvalue::check_types (x, *xv);
+
+  std::list<colvarvalue>::iterator &xvi = xv;
+  std::vector<cvm::real>::iterator  &ii = inner;
+
+  switch (x.value_type) {
+  case colvarvalue::type_scalar:
+    while (xvi != xv_end) {
+      *(ii++) += (xvi++)->real_value * x.real_value;
+    }
+    break;
+  case colvarvalue::type_vector:
+  case colvarvalue::type_unitvector:
+    while (xvi != xv_end) {
+      *(ii++) += (xvi++)->rvector_value * x.rvector_value;
+    }
+    break;
+  case colvarvalue::type_quaternion:
+    while (xvi != xv_end) {
+      *(ii++) += ((xvi++)->quaternion_value).cosine (x.quaternion_value);
+    }
+    break;
+  default:
+    x.undef_op();
+  };
+}
+
+
+void colvarvalue::p2leg_opt (colvarvalue const                        &x,
+                             std::vector<colvarvalue>::iterator       &xv,
+                             std::vector<colvarvalue>::iterator const &xv_end,
+                             std::vector<cvm::real>::iterator         &inner)
+{
+  // doing type check only once, here
+  colvarvalue::check_types (x, *xv);
+
+  std::vector<colvarvalue>::iterator &xvi = xv;
+  std::vector<cvm::real>::iterator    &ii = inner;
+
+  switch (x.value_type) {
+  case colvarvalue::type_scalar:
+    cvm::fatal_error ("Error: cannot calculate Legendre polynomials "
+                      "for scalar variables.\n");
+    break;
+  case colvarvalue::type_vector:
+    while (xvi != xv_end) {
+      cvm::real const cosine = 
+        ((xvi)->rvector_value * x.rvector_value) /
+        ((xvi)->rvector_value.norm() * x.rvector_value.norm());
+      xvi++;
+      *(ii++) += 1.5*cosine*cosine - 0.5;
+    }
+    break;
+  case colvarvalue::type_unitvector:
+    while (xvi != xv_end) {
+      cvm::real const cosine = (xvi++)->rvector_value * x.rvector_value;
+      *(ii++) += 1.5*cosine*cosine - 0.5;
+    }
+    break;
+  case colvarvalue::type_quaternion:
+    while (xvi != xv_end) {
+      cvm::real const cosine = (xvi++)->quaternion_value.cos (x.quaternion_value);
+      *(ii++) += 1.5*cosine*cosine - 0.5;
+    }
+    break;
+  default:
+    x.undef_op();
+  };
+}
+
+void colvarvalue::p2leg_opt (colvarvalue const                        &x,
+                             std::list<colvarvalue>::iterator         &xv,
+                             std::list<colvarvalue>::iterator const   &xv_end,
+                             std::vector<cvm::real>::iterator         &inner)
+{
+  // doing type check only once, here
+  colvarvalue::check_types (x, *xv);
+
+  std::list<colvarvalue>::iterator &xvi = xv;
+  std::vector<cvm::real>::iterator  &ii = inner;
+
+  switch (x.value_type) {
+  case colvarvalue::type_scalar:
+    cvm::fatal_error ("Error: cannot calculate Legendre polynomials "
+                      "for scalar variables.\n");
+    break;
+  case colvarvalue::type_vector:
+    while (xvi != xv_end) {
+      cvm::real const cosine = 
+        ((xvi)->rvector_value * x.rvector_value) /
+        ((xvi)->rvector_value.norm() * x.rvector_value.norm());
+      xvi++;
+      *(ii++) += 1.5*cosine*cosine - 0.5;
+    }
+    break;
+  case colvarvalue::type_unitvector:
+    while (xvi != xv_end) {
+      cvm::real const cosine = (xvi++)->rvector_value * x.rvector_value;
+      *(ii++) += 1.5*cosine*cosine - 0.5;
+    }
+    break;
+  case colvarvalue::type_quaternion:
+    while (xvi != xv_end) {
+      cvm::real const cosine = (xvi++)->quaternion_value.cosine (x.quaternion_value);
+      *(ii++) += 1.5*cosine*cosine - 0.5;
+    }
+    break;
+  default:
+    x.undef_op();
+  };
+}
+
+
+std::ostream & operator << (std::ostream &os, colvarvalue const &x)
+{
+  switch (x.type()) {
+  case colvarvalue::type_scalar:
+    os << x.real_value; break;
+  case colvarvalue::type_vector:
+  case colvarvalue::type_unitvector:
+    os << x.rvector_value; break;
+  case colvarvalue::type_quaternion:
+    os << x.quaternion_value; break;
+  case colvarvalue::type_notset:
+    os << "not set"; break;
+  }
+  return os;
+}
+
+
+std::ostream & operator << (std::ostream &os, std::vector<colvarvalue> const &v)
+{
+  for (size_t i = 0; i < v.size(); i++) os << v[i];
+  return os;
+}
+
+
+std::istream & operator >> (std::istream &is, colvarvalue &x)
+{
+  if (x.type() == colvarvalue::type_notset) {
+    cvm::fatal_error ("Trying to read from a stream a colvarvalue, "
+                      "which has not yet been assigned a data type.\n");
+  }
+
+  switch (x.type()) {
+  case colvarvalue::type_scalar:
+    is >> x.real_value; 
+    break;
+  case colvarvalue::type_vector:
+  case colvarvalue::type_unitvector:
+    is >> x.rvector_value; 
+    break;
+  case colvarvalue::type_quaternion:
+    is >> x.quaternion_value;
+    break;
+  default:
+    x.undef_op();
+  }
+  x.apply_constraints();
+  return is;
+}
+
+
+size_t colvarvalue::output_width (size_t const &real_width)
+{
+  switch (this->value_type) {
+  case colvarvalue::type_scalar:
+    return real_width;
+  case colvarvalue::type_vector:
+  case colvarvalue::type_unitvector:
+    return cvm::rvector::output_width (real_width);
+  case colvarvalue::type_quaternion:
+    return cvm::quaternion::output_width (real_width);
+  case colvarvalue::type_notset:
+  default:
+    return 0;
+  }
+}
+
+

lib/colvars/colvarvalue.h

@@ -0,0 +1,727 @@
+#ifndef COLVARVALUE_H
+#define COLVARVALUE_H
+
+#include "colvarmodule.h"
+
+
+/// \brief Value of a collective variable: this is a metatype which
+/// can be set at runtime.  By default it is set to be a scalar
+/// number, and can be treated like that in all operations (this is
+/// done by most \link cvc \endlink implementations).
+///
+/// \link colvarvalue \endlink allows \link colvar \endlink to be
+/// treat different data types.  By default, a \link colvarvalue
+/// \endlink variable is a scalar number.  If you want to use it as
+/// another type, you should declare and initialize a variable as
+/// \code colvarvalue x (colvarvalue::type_xxx); \endcode where
+/// type_xxx is a value within the \link Type \endlink enum.
+/// Alternatively, initialize x with \link x.type
+/// (colvarvalue::type_xxx) \endlink at a later stage.
+///
+/// Given a colvarvalue variable x which is not yet assigned (and
+/// thus has not yet a type) it is also possible to correctly assign
+/// the type with \code x = y; \endcode if y is correctly set.
+/// Otherwise, an error will be raised if the \link Type \endlink of x
+/// is different from the \link Type \endlink of y.
+///
+/// Also, every operator (either unary or binary) on a \link
+/// colvarvalue \endlink object performs one or more checks on the
+/// \link Type \endlink to avoid errors such as initializing a
+/// three-dimensional vector with a scalar number (legal otherwise).
+///
+/// \b Special \b case: when reading from a stream, there is no way to
+/// detect the \link Type \endlink and safely handle errors at the
+/// same time.  Hence, when using \code is >> x; \endcode x \b MUST
+/// already have a type correcly set up for properly parsing the
+/// stream.  An error will be raised otherwise.  Usually this is not
+/// the problem, because \link colvarvalue \endlink objects are first
+/// initialized in the configuration, and the stream operation will be
+/// performed only when reading restart files.
+/// 
+/// No problem of course with the output streams: \code os << x;
+/// \endcode will print a different output according to the value of
+/// colvarvalue::value_type, and the width of such output is returned
+/// by colvarvalue::output_width()
+///
+/// \em Note \em on \em performance: at every operation between two
+/// \link colvarvalue \endlink objects, their two \link Type \endlink
+/// flags will be checked for a match.  In a long array of \link
+/// colvarvalue \endlink objects this is time consuming: a few static
+/// functions are defined ("xxx_opt" functions) within the \link
+/// colvarvalue \endlink class, which only check the matching once for
+/// a large array, and execute different loops according to the type.
+/// You should do the same for every time consuming loop involving
+/// operations on \link colvarvalue \endlink objects if you want
+/// e.g. to optimize your colvar bias.
+
+class colvarvalue {
+
+public:
+
+  /// \brief Possible types of value
+  ///
+  /// These three cover most possibilities of data type one can
+  /// devise.  If you need to implement a new colvar with a very
+  /// complex data type, it's better to put an allocatable class here
+  enum Type {
+    /// Undefined type
+    type_notset,
+    /// Scalar number, implemented as \link colvarmodule::real \endlink (default)
+    type_scalar,
+    /// 3-dimensional vector, implemented as \link colvarmodule::rvector \endlink
+    type_vector,
+    /// 3-dimensional unit vector, implemented as \link colvarmodule::rvector \endlink
+    type_unitvector,
+    /// 4-dimensional unit vector representing a rotation, implemented as \link colvarmodule::quaternion \endlink
+    type_quaternion
+  };
+
+  /// Runtime description of value types
+  std::string static const type_desc[colvarvalue::type_quaternion+1];
+
+  /// Number of degrees of freedom for each type
+  size_t static const      dof_num[  colvarvalue::type_quaternion+1];
+
+  /// \brief Real data member
+  cvm::real       real_value;
+
+  /// \brief Vector data member
+  cvm::rvector    rvector_value;
+
+  /// \brief Quaternion data member
+  cvm::quaternion quaternion_value;
+
+  /// Current type of this colvarvalue object
+  Type  value_type; 
+
+  static inline bool type_checking()
+  {
+    return true;
+  }
+
+  /// \brief Default constructor: this class defaults to a scalar
+  /// number and always behaves like it unless you change its type
+  inline colvarvalue()
+    : real_value (0.0), value_type (type_scalar)
+  {}
+
+  /// Constructor from a type specification
+  inline colvarvalue (Type const &vti)
+    : value_type (vti)
+  {
+    reset();
+  }
+
+  /// Copy constructor from real base type
+  inline colvarvalue (cvm::real const &x)
+    : real_value (x), value_type (type_scalar)
+  {}
+
+  /// \brief Copy constructor from rvector base type (Note: this sets
+  /// automatically a type \link type_vector \endlink , if you want a
+  /// \link type_unitvector \endlink you must set it explicitly)
+  inline colvarvalue (cvm::rvector const &v)
+    : rvector_value (v), value_type (type_vector)
+  {}
+
+  /// \brief Copy constructor from rvector base type (additional
+  /// argument to make possible to choose a \link type_unitvector
+  /// \endlink
+  inline colvarvalue (cvm::rvector const &v, Type const &vti)
+    : rvector_value (v), value_type (vti)
+  {}
+
+  /// \brief Copy constructor from quaternion base type
+  inline colvarvalue (cvm::quaternion const &q)
+    : quaternion_value (q), value_type (type_quaternion)
+  {}
+
+  /// Copy constructor from another \link colvarvalue \endlink
+  inline colvarvalue (colvarvalue const &x)
+    : value_type (x.value_type)
+  {
+    reset();
+
+    switch (x.value_type) {
+    case type_scalar:
+      real_value = x.real_value;
+      break;
+    case type_vector: 
+    case type_unitvector:
+      rvector_value = x.rvector_value;
+      break;
+    case type_quaternion:
+      quaternion_value = x.quaternion_value;
+      break;
+    case type_notset:
+    default:
+      break;
+    }
+  }
+
+
+
+  /// Set to the null value for the data type currently defined
+  void reset();
+
+
+  /// \brief If the variable has constraints (e.g. unitvector or
+  /// quaternion), transform it to satisfy them; use it when the \link
+  /// colvarvalue \endlink is not calculated from \link cvc
+  /// \endlink objects, but manipulated by you
+  void apply_constraints();
+
+
+  /// Get the current type
+  inline Type type() const
+  {
+    return value_type;
+  }
+
+  /// Set the type explicitly
+  inline void type (Type const &vti)
+  {
+    reset();
+    value_type = vti;
+    reset();
+  }
+
+  /// Set the type after another \link colvarvalue \endlink
+  inline void type (colvarvalue const &x)
+  {
+    reset();
+    value_type = x.value_type;
+    reset();
+  }
+
+  /// Square norm of this colvarvalue
+  cvm::real norm2() const;
+
+  /// Norm of this colvarvalue
+  inline cvm::real norm() const
+  {
+    return std::sqrt (this->norm2());
+  }
+
+  /// \brief Return the value whose scalar product with this value is
+  /// 1
+  inline colvarvalue inverse() const;
+
+  /// Square distance between this \link colvarvalue \endlink and another
+  cvm::real dist2 (colvarvalue const &x2) const;
+
+  /// Derivative with respect to this \link colvarvalue \endlink of the square distance
+  colvarvalue dist2_grad (colvarvalue const &x2) const;
+
+  /// Assignment operator (type of x is checked)
+  colvarvalue & operator = (colvarvalue const &x);
+
+  void operator += (colvarvalue const &x);
+  void operator -= (colvarvalue const &x);
+  void operator *= (cvm::real const &a);
+  void operator /= (cvm::real const &a);
+
+
+  // Casting operator
+  inline operator cvm::real() const
+  {
+    if (value_type != type_scalar) error_rside (type_scalar);
+    return real_value;
+  }
+
+  // Casting operator
+  inline operator cvm::rvector() const
+  {
+    if ( (value_type != type_vector) && (value_type != type_unitvector))
+      error_rside (type_vector);
+    return rvector_value;
+  }
+
+  // Casting operator
+  inline operator cvm::quaternion() const
+  {
+    if (value_type != type_quaternion) error_rside (type_quaternion);
+    return quaternion_value;
+  }
+
+  /// Special case when the variable is a real number, and all operations are defined
+  inline bool is_scalar()
+  {
+    return (value_type == type_scalar);
+  }
+
+
+  /// Ensure that the two types are the same within a binary operator
+  void static check_types (colvarvalue const &x1, colvarvalue const &x2);
+
+  /// Undefined operation 
+  void undef_op() const;
+
+  /// Trying to assign this \link colvarvalue \endlink object to
+  /// another object set with a different type
+  void error_lside (Type const &vt) const;
+
+  /// Trying to assign another \link colvarvalue \endlink object set
+  /// with a different type to this object
+  void error_rside (Type const &vt) const;
+
+  ///<2F>Give the number of characters required to output this
+  /// colvarvalue, given the current type assigned and the number of
+  /// characters for a real number
+  size_t output_width (size_t const &real_width);
+
+
+  // optimized routines for operations with an array; xv and inner as
+  // vectors are assumed to have the same number of elements (i.e. the
+  // end for inner comes together with xv_end in the loop)
+
+  /// \brief Optimized routine for the inner product of one collective
+  /// variable with an array
+  void static inner_opt (colvarvalue const                        &x,
+                         std::vector<colvarvalue>::iterator       &xv,
+                         std::vector<colvarvalue>::iterator const &xv_end,
+                         std::vector<cvm::real>::iterator         &inner);
+
+  /// \brief Optimized routine for the inner product of one collective
+  /// variable with an array
+  void static inner_opt (colvarvalue const                        &x,
+                         std::list<colvarvalue>::iterator         &xv,
+                         std::list<colvarvalue>::iterator const   &xv_end,
+                         std::vector<cvm::real>::iterator         &inner);
+
+  /// \brief Optimized routine for the second order Legendre
+  /// polynomial, (3cos^2(w)-1)/2, of one collective variable with an
+  /// array
+  void static p2leg_opt (colvarvalue const                        &x,
+                         std::vector<colvarvalue>::iterator       &xv,
+                         std::vector<colvarvalue>::iterator const &xv_end,
+                         std::vector<cvm::real>::iterator         &inner);
+
+  /// \brief Optimized routine for the second order Legendre
+  /// polynomial of one collective variable with an array
+  void static p2leg_opt (colvarvalue const                        &x,
+                         std::list<colvarvalue>::iterator         &xv,
+                         std::list<colvarvalue>::iterator const   &xv_end,
+                         std::vector<cvm::real>::iterator         &inner);
+
+};
+
+
+
+inline void colvarvalue::reset()
+{
+  switch (value_type) {
+  case colvarvalue::type_scalar:
+    real_value = cvm::real (0.0);
+    break;
+  case colvarvalue::type_vector:
+  case colvarvalue::type_unitvector:
+    rvector_value = cvm::rvector (0.0, 0.0, 0.0);
+    break;
+  case colvarvalue::type_quaternion:
+    quaternion_value = cvm::quaternion (0.0, 0.0, 0.0, 0.0);
+    break;
+  case colvarvalue::type_notset:
+  default:
+    break;
+  }  
+}
+
+
+inline void colvarvalue::apply_constraints()
+{
+  switch (value_type) {
+  case colvarvalue::type_scalar:
+    break;
+  case colvarvalue::type_vector:
+    break;
+  case colvarvalue::type_unitvector:
+    rvector_value /= std::sqrt (rvector_value.norm2());
+    break;
+  case colvarvalue::type_quaternion:
+    quaternion_value /= std::sqrt (quaternion_value.norm2());
+    break;
+  case colvarvalue::type_notset:
+  default:
+    break;
+  }  
+}
+
+
+
+inline cvm::real colvarvalue::norm2() const
+{
+  switch (value_type) {
+  case colvarvalue::type_scalar:
+    return (this->real_value)*(this->real_value);
+  case colvarvalue::type_vector:
+  case colvarvalue::type_unitvector:
+    return (this->rvector_value).norm2();
+  case colvarvalue::type_quaternion:
+    return (this->quaternion_value).norm2();
+  case colvarvalue::type_notset:
+  default:
+    return 0.0;
+  }  
+}
+
+
+inline colvarvalue colvarvalue::inverse() const
+{
+  switch (value_type) {
+  case colvarvalue::type_scalar:
+    return colvarvalue (1.0/real_value);
+  case colvarvalue::type_vector:
+    return colvarvalue (cvm::rvector (1.0/rvector_value.x,
+                                      1.0/rvector_value.y,
+                                      1.0/rvector_value.z));
+  case colvarvalue::type_quaternion:
+    return colvarvalue (cvm::quaternion (1.0/quaternion_value.q0,
+                                         1.0/quaternion_value.q1,
+                                         1.0/quaternion_value.q2,
+                                         1.0/quaternion_value.q3));
+  case colvarvalue::type_notset:
+  default:
+    undef_op();
+  }  
+  return colvarvalue();
+}
+
+
+inline colvarvalue & colvarvalue::operator = (colvarvalue const &x)
+{
+  if (this->value_type != type_notset)
+    if (this->value_type != x.value_type) 
+      error_lside (x.value_type);
+
+  this->value_type = x.value_type;
+  
+  switch (this->value_type) {
+  case colvarvalue::type_scalar:
+    this->real_value = x.real_value;
+    break;
+  case colvarvalue::type_vector:
+  case colvarvalue::type_unitvector:
+    this->rvector_value = x.rvector_value;
+    break;
+  case colvarvalue::type_quaternion:
+    this->quaternion_value = x.quaternion_value;
+    break;
+  case colvarvalue::type_notset:
+  default:
+    undef_op();
+    break;
+  }
+  return *this;
+}
+
+inline void colvarvalue::operator += (colvarvalue const &x)
+{
+  if (colvarvalue::type_checking())
+    colvarvalue::check_types (*this, x);
+
+  switch (this->value_type) {
+  case colvarvalue::type_scalar:
+    this->real_value += x.real_value;
+    break;
+  case colvarvalue::type_vector:
+  case colvarvalue::type_unitvector:
+    this->rvector_value += x.rvector_value;
+    break;
+  case colvarvalue::type_quaternion:
+    this->quaternion_value += x.quaternion_value;
+    break;
+  case colvarvalue::type_notset:
+  default:
+    undef_op();
+  }
+}
+
+inline void colvarvalue::operator -= (colvarvalue const &x)
+{
+  if (colvarvalue::type_checking())
+    colvarvalue::check_types (*this, x);
+
+  switch (value_type) {
+  case colvarvalue::type_scalar:
+    real_value -= x.real_value; break;
+  case colvarvalue::type_vector:
+  case colvarvalue::type_unitvector:
+    rvector_value -= x.rvector_value; break;
+  case colvarvalue::type_quaternion:
+    quaternion_value -= x.quaternion_value; break;
+  case colvarvalue::type_notset:
+  default:
+    undef_op();
+  }
+}
+
+inline void colvarvalue::operator *= (cvm::real const &a)
+{
+  switch (value_type) {
+  case colvarvalue::type_scalar:
+    real_value *= a;
+    break;
+  case colvarvalue::type_vector:
+  case colvarvalue::type_unitvector:
+    rvector_value *= a;
+    break;
+  case colvarvalue::type_quaternion:
+    quaternion_value *= a;
+    break;
+  case colvarvalue::type_notset:
+  default:
+    undef_op();
+  }
+}
+
+inline void colvarvalue::operator /= (cvm::real const &a)
+{
+  switch (value_type) {
+  case colvarvalue::type_scalar:
+    real_value /= a; break;
+  case colvarvalue::type_vector:
+  case colvarvalue::type_unitvector:
+    rvector_value /= a; break;
+  case colvarvalue::type_quaternion:
+    quaternion_value /= a; break;
+  case colvarvalue::type_notset:
+  default:
+    undef_op();
+  }
+}
+
+
+// binary operations between two colvarvalues
+
+inline colvarvalue operator + (colvarvalue const &x1,
+                               colvarvalue const &x2)
+{
+  if (colvarvalue::type_checking())
+    colvarvalue::check_types (x1, x2);
+
+  switch (x1.value_type) {
+  case colvarvalue::type_scalar:
+    return colvarvalue (x1.real_value + x2.real_value);
+  case colvarvalue::type_vector:
+    return colvarvalue (x1.rvector_value + x2.rvector_value);
+  case colvarvalue::type_unitvector:
+    return colvarvalue (x1.rvector_value + x2.rvector_value,
+                        colvarvalue::type_unitvector);
+  case colvarvalue::type_quaternion:
+    return colvarvalue (x1.quaternion_value + x2.quaternion_value);
+  case colvarvalue::type_notset:
+  default:
+    x1.undef_op();
+    return colvarvalue (colvarvalue::type_notset);
+  };
+}
+
+inline colvarvalue operator - (colvarvalue const &x1,
+                               colvarvalue const &x2)
+{
+  if (colvarvalue::type_checking())
+    colvarvalue::check_types (x1, x2);
+
+  switch (x1.value_type) {
+  case colvarvalue::type_scalar:
+    return colvarvalue (x1.real_value - x2.real_value);
+  case colvarvalue::type_vector:
+    return colvarvalue (x1.rvector_value - x2.rvector_value);
+  case colvarvalue::type_unitvector:
+    return colvarvalue (x1.rvector_value - x2.rvector_value,
+                        colvarvalue::type_unitvector);
+  case colvarvalue::type_quaternion:
+    return colvarvalue (x1.quaternion_value - x2.quaternion_value);
+  default:
+    x1.undef_op();
+    return colvarvalue (colvarvalue::type_notset);
+  };
+}
+
+
+// binary operations with real numbers
+
+inline colvarvalue operator * (cvm::real const &a,
+                               colvarvalue const &x)
+{
+  switch (x.value_type) {
+  case colvarvalue::type_scalar:
+    return colvarvalue (a * x.real_value);
+  case colvarvalue::type_vector:
+    return colvarvalue (a * x.rvector_value);
+  case colvarvalue::type_unitvector:
+    return colvarvalue (a * x.rvector_value,
+                        colvarvalue::type_unitvector);
+  case colvarvalue::type_quaternion:
+    return colvarvalue (a * x.quaternion_value);
+  case colvarvalue::type_notset:
+  default:
+    x.undef_op();
+    return colvarvalue (colvarvalue::type_notset);
+  }
+}
+
+inline colvarvalue operator * (colvarvalue const &x,
+                               cvm::real const &a)
+{
+  switch (x.value_type) {
+  case colvarvalue::type_scalar:
+    return colvarvalue (x.real_value * a);
+  case colvarvalue::type_vector:
+    return colvarvalue (x.rvector_value * a);
+  case colvarvalue::type_unitvector:
+    return colvarvalue (x.rvector_value * a,
+                        colvarvalue::type_unitvector);
+  case colvarvalue::type_quaternion:
+    return colvarvalue (x.quaternion_value * a);
+  case colvarvalue::type_notset:
+  default:
+    x.undef_op();
+    return colvarvalue (colvarvalue::type_notset);
+  }
+}
+
+inline colvarvalue operator / (colvarvalue const &x,
+                               cvm::real const &a)
+{
+  switch (x.value_type) {
+  case colvarvalue::type_scalar:
+    return colvarvalue (x.real_value / a);
+  case colvarvalue::type_vector:
+    return colvarvalue (x.rvector_value / a);
+  case colvarvalue::type_unitvector:
+    return colvarvalue (x.rvector_value / a,
+                        colvarvalue::type_unitvector);
+  case colvarvalue::type_quaternion:
+    return colvarvalue (x.quaternion_value / a);
+  case colvarvalue::type_notset:
+  default:
+    x.undef_op();
+    return colvarvalue (colvarvalue::type_notset);
+  }
+}
+
+
+// inner product between two colvarvalues
+
+inline cvm::real operator * (colvarvalue const &x1,
+                             colvarvalue const &x2)
+{
+  if (colvarvalue::type_checking())
+    colvarvalue::check_types (x1, x2);
+
+  switch (x1.value_type) {
+  case colvarvalue::type_scalar:
+    return (x1.real_value * x2.real_value);
+  case colvarvalue::type_vector:
+  case colvarvalue::type_unitvector:
+    return (x1.rvector_value * x2.rvector_value);
+  case colvarvalue::type_quaternion:
+    // the "*" product is the quaternion product, here the inner
+    // member function is used instead
+    return (x1.quaternion_value.inner (x2.quaternion_value));
+  case colvarvalue::type_notset:
+  default:
+    x1.undef_op();
+    return 0.0;
+  };
+}
+
+
+
+inline cvm::real colvarvalue::dist2 (colvarvalue const &x2) const
+{
+  if (colvarvalue::type_checking())
+    colvarvalue::check_types (*this, x2);
+
+  switch (this->value_type) {
+  case colvarvalue::type_scalar:
+    return std::pow (this->real_value - x2.real_value, int (2));
+  case colvarvalue::type_vector:
+    return (this->rvector_value - x2.rvector_value).norm2();
+  case colvarvalue::type_unitvector:
+    // angle between (*this) and x2 is the distance
+    return std::pow (std::acos (this->rvector_value * x2.rvector_value), int (2));
+  case colvarvalue::type_quaternion:
+    // angle between (*this) and x2 is the distance, the quaternion
+    // object has it implemented internally
+    return this->quaternion_value.dist2 (x2.quaternion_value);
+  case colvarvalue::type_notset:
+  default:
+    this->undef_op();
+    return 0.0;
+  };
+}
+
+
+inline colvarvalue colvarvalue::dist2_grad (colvarvalue const &x2) const
+{
+  if (colvarvalue::type_checking())
+    colvarvalue::check_types (*this, x2);
+
+  switch (this->value_type) {
+  case colvarvalue::type_scalar:
+    return 2.0 * (this->real_value - x2.real_value);
+  case colvarvalue::type_vector:
+    return 2.0 * (this->rvector_value - x2.rvector_value);
+  case colvarvalue::type_unitvector:
+    {
+      cvm::rvector const &v1 = this->rvector_value;
+      cvm::rvector const &v2 = x2.rvector_value;
+      cvm::real const cos_t = v1 * v2;
+      cvm::real const sin_t = std::sqrt (1.0 - cos_t*cos_t);
+      return colvarvalue ( 2.0 * sin_t *
+                           cvm::rvector ( (-1.0) * sin_t * v2.x +
+                                          cos_t/sin_t * (v1.x - cos_t*v2.x),
+                                          (-1.0) * sin_t * v2.y +
+                                          cos_t/sin_t * (v1.y - cos_t*v2.y),
+                                          (-1.0) * sin_t * v2.z +
+                                          cos_t/sin_t * (v1.z - cos_t*v2.z)
+                                          ),
+                           colvarvalue::type_unitvector );
+    }
+  case colvarvalue::type_quaternion:
+    return this->quaternion_value.dist2_grad (x2.quaternion_value);
+  case colvarvalue::type_notset:
+  default:
+    this->undef_op();
+    return colvarvalue (colvarvalue::type_notset);
+  };
+}
+
+
+inline void colvarvalue::check_types (colvarvalue const &x1,
+                                      colvarvalue const &x2)
+{
+  if (x1.value_type != x2.value_type) {
+    cvm::log ("x1 type = "+cvm::to_str (x1.value_type)+
+              ", x2 type = "+cvm::to_str (x2.value_type)+"\n");
+    cvm::fatal_error ("Performing an operation between two colvar "
+                      "values with different types, \""+
+                      colvarvalue::type_desc[x1.value_type]+
+                      "\" and \""+
+                      colvarvalue::type_desc[x2.value_type]+
+                      "\".\n");
+  }
+}
+
+
+
+
+std::ostream & operator << (std::ostream &os, colvarvalue const &x);
+std::ostream & operator << (std::ostream &os, std::vector<colvarvalue> const &v);
+
+std::istream & operator >> (std::istream &is, colvarvalue &x);
+
+
+
+
+
+#endif
+
+
+// Emacs
+// Local Variables:
+// mode: C++
+// End: