Update Colvars library to version 2016-12-22

Significant code cleanup and several fixes (walls + extended Lagrangian)

New harmonicWalls bias to apply confining boundaries with time-dependent force
constant & integration

doc/src/fix_colvars.txt

@@ -10,7 +10,7 @@ fix colvars command :h3
 
 [Syntax:]
 
-fix ID group-ID colvars configfile keyword values ... :pre
+fix ID group-ID Colvars configfile keyword values ... :pre
 
 ID, group-ID are documented in "fix"_fix.html command :ulb,l
 colvars = style name of this fix command :l
@@ -31,21 +31,19 @@ fix abf all colvars colvars.inp tstat 1 :pre
 
 [Description:]
 
-This fix interfaces LAMMPS to a "collective variables" or "colvars"
+This fix interfaces LAMMPS to the collective variables "Colvars"
-module library which allows to calculate potentials of mean force
+library, which allows to calculate potentials of mean force
 (PMFs) for any set of colvars, using different sampling methods:
 currently implemented are the Adaptive Biasing Force (ABF) method,
 metadynamics, Steered Molecular Dynamics (SMD) and Umbrella Sampling
-(US) via a flexible harmonic restraint bias. The colvars library is
+(US) via a flexible harmonic restraint bias.
-hosted at "http://colvars.github.io/"_http://colvars.github.io/
 
 This documentation describes only the fix colvars command itself and
 LAMMPS specific parts of the code.  The full documentation of the
 colvars library is available as "this supplementary PDF document"_PDF/colvars-refman-lammps.pdf
 
-A detailed discussion of the implementation of the portable collective
+The Colvars library is developed at "https://github.com/colvars/colvars"_https://github.com/colvars/colvars
-variable library is in "(Fiorin)"_#Fiorin. Additional information can
+A detailed discussion of its implementation is in "(Fiorin)"_#Fiorin.
-be found in "(Henin)"_#Henin.
 
 There are some example scripts for using this package with LAMMPS in the
 examples/USER/colvars directory.
@@ -129,8 +127,3 @@ and tstat = NULL.
 
 :link(Fiorin)
 [(Fiorin)] Fiorin , Klein, Henin, Mol. Phys., DOI:10.1080/00268976.2013.813594
-
-:link(Henin)
-[(Henin)] Henin, Fiorin, Chipot, Klein, J. Chem. Theory Comput., 6,
-35-47 (2010)
-

lib/colvars/colvar.cpp

@@ -1,5 +1,12 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #include "colvarmodule.h"
 #include "colvarvalue.h"
 #include "colvarparse.h"
@@ -163,38 +170,38 @@ colvar::colvar(std::string const &conf)
     }
     feature_states[f_cv_homogeneous]->enabled = homogeneous;
   }
-  // Colvar is deemed periodic iff:
+
+  // Colvar is deemed periodic if:
   // - it is homogeneous
   // - all cvcs are periodic
   // - all cvcs have the same period
-
+  if (cvcs[0]->b_periodic) { // TODO make this a CVC feature
-  b_periodic = cvcs[0]->b_periodic && is_enabled(f_cv_homogeneous);
+    bool b_periodic = true;
-  period = cvcs[0]->period;
+    period = cvcs[0]->period;
-  for (i = 1; i < cvcs.size(); i++) {
+    for (i = 1; i < cvcs.size(); i++) {
-    if (!cvcs[i]->b_periodic || cvcs[i]->period != period) {
+      if (!cvcs[i]->b_periodic || cvcs[i]->period != period) {
-      b_periodic = false;
+        b_periodic = false;
-      period = 0.0;
+        period = 0.0;
+        cvm::log("Warning: although one component is periodic, this colvar will "
+                 "not be treated as periodic, either because the exponent is not "
+                 "1, or because components of different periodicity are defined.  "
+                 "Make sure that you know what you are doing!");
+      }
     }
+    feature_states[f_cv_periodic]->enabled = b_periodic;
   }
-  feature_states[f_cv_periodic]->enabled = b_periodic;
 
   // check that cvcs are compatible
 
   for (i = 0; i < cvcs.size(); i++) {
-    if ((cvcs[i])->b_periodic && !b_periodic) {
-        cvm::log("Warning: although this component is periodic, the colvar will "
-                  "not be treated as periodic, either because the exponent is not "
-                  "1, or because multiple components are present. Make sure that "
-                  "you know what you are doing!");
-    }
 
     // components may have different types only for scripted functions
     if (!is_enabled(f_cv_scripted) && (colvarvalue::check_types(cvcs[i]->value(),
-                                                           cvcs[0]->value())) ) {
+                                                                cvcs[0]->value())) ) {
       cvm::error("ERROR: you are definining this collective variable "
                  "by using components of different types. "
                  "You must use the same type in order to "
-                 " sum them together.\n", INPUT_ERROR);
+                 "sum them together.\n", INPUT_ERROR);
       return;
     }
   }
@@ -207,16 +214,15 @@ colvar::colvar(std::string const &conf)
   // at this point, the colvar's type is defined
   f.type(value());
   f_accumulated.type(value());
-  fb.type(value());
+
+  reset_bias_force();
 
   get_keyval(conf, "width", width, 1.0);
   if (width <= 0.0) {
     cvm::error("Error: \"width\" must be positive.\n", INPUT_ERROR);
+    return;
   }
 
-  // NOTE: not porting wall stuff to new deps, as this will change to a separate bias
-  // the grid functions will wait a little as well
-
   lower_boundary.type(value());
   lower_wall.type(value());
 
@@ -400,6 +406,9 @@ colvar::colvar(std::string const &conf)
   f_old.type(value());
   f_old.reset();
 
+  x_restart.type(value());
+  after_restart = false;
+
   if (cvm::b_analysis)
     parse_analysis(conf);
 
@@ -761,9 +770,13 @@ int colvar::calc_cvcs(int first_cvc, size_t num_cvcs)
     return error_code;
   }
 
+  if (cvm::step_relative() > 0) {
+    // Total force depends on Jacobian derivative from previous timestep
+    error_code |= calc_cvc_total_force(first_cvc, num_cvcs);
+  }
+  // atom coordinates are updated by the next line
   error_code |= calc_cvc_values(first_cvc, num_cvcs);
   error_code |= calc_cvc_gradients(first_cvc, num_cvcs);
-  error_code |= calc_cvc_total_force(first_cvc, num_cvcs);
   error_code |= calc_cvc_Jacobians(first_cvc, num_cvcs);
 
   if (cvm::debug())
@@ -780,9 +793,12 @@ int colvar::collect_cvc_data()
 
   int error_code = COLVARS_OK;
 
+  if (cvm::step_relative() > 0) {
+    // Total force depends on Jacobian derivative from previous timestep
+    error_code |= collect_cvc_total_forces();
+  }
   error_code |= collect_cvc_values();
   error_code |= collect_cvc_gradients();
-  error_code |= collect_cvc_total_forces();
   error_code |= collect_cvc_Jacobians();
   error_code |= calc_colvar_properties();
 
@@ -872,6 +888,22 @@ int colvar::collect_cvc_values()
     cvm::log("Colvar \""+this->name+"\" has value "+
               cvm::to_str(x, cvm::cv_width, cvm::cv_prec)+".\n");
 
+  if (after_restart) {
+    after_restart = false;
+    if (cvm::proxy->simulation_running()) {
+      cvm::real const jump2 = dist2(x, x_restart) / (width*width);
+      if (jump2 > 0.25) {
+        cvm::error("Error: the calculated value of colvar \""+name+
+                   "\":\n"+cvm::to_str(x)+"\n differs greatly from the value "
+                   "last read from the state file:\n"+cvm::to_str(x_restart)+
+                   "\nPossible causes are changes in configuration, "
+                   "wrong state file, or how PBC wrapping is handled.\n",
+                   INPUT_ERROR);
+        return INPUT_ERROR;
+      }
+    }
+  }
+
   return COLVARS_OK;
 }
 
@@ -979,22 +1011,17 @@ int colvar::calc_cvc_total_force(int first_cvc, size_t num_cvcs)
     if (cvm::debug())
       cvm::log("Calculating total force of colvar \""+this->name+"\".\n");
 
-    // if (!tasks[task_extended_lagrangian] && (cvm::step_relative() > 0)) {
+    cvm::increase_depth();
-   // Disabled check to allow for explicit total force calculation
-    // even with extended Lagrangian
 
-    if (cvm::step_relative() > 0) {
+    for (i = first_cvc, cvc_count = 0;
-      cvm::increase_depth();
+        (i < cvcs.size()) && (cvc_count < cvc_max_count);
-      // get from the cvcs the total forces from the PREVIOUS step
+        i++) {
-      for (i = first_cvc, cvc_count = 0;
+      if (!cvcs[i]->is_enabled()) continue;
-          (i < cvcs.size()) && (cvc_count < cvc_max_count);
+      cvc_count++;
-          i++) {
+      (cvcs[i])->calc_force_invgrads();
-        if (!cvcs[i]->is_enabled()) continue;
-        cvc_count++;
-        (cvcs[i])->calc_force_invgrads();
-      }
-      cvm::decrease_depth();
     }
+    cvm::decrease_depth();
+
 
     if (cvm::debug())
       cvm::log("Done calculating total force of colvar \""+this->name+"\".\n");
@@ -1013,6 +1040,11 @@ int colvar::collect_cvc_total_forces()
       // get from the cvcs the total forces from the PREVIOUS step
       for (size_t i = 0; i < cvcs.size();  i++) {
         if (!cvcs[i]->is_enabled()) continue;
+            if (cvm::debug())
+            cvm::log("Colvar component no. "+cvm::to_str(i+1)+
+                " within colvar \""+this->name+"\" has total force "+
+                cvm::to_str((cvcs[i])->total_force(),
+                cvm::cv_width, cvm::cv_prec)+".\n");
         // linear combination is assumed
         ft += (cvcs[i])->total_force() * (cvcs[i])->sup_coeff / active_cvc_square_norm;
       }
@@ -1056,6 +1088,11 @@ int colvar::collect_cvc_Jacobians()
     fj.reset();
     for (size_t i = 0; i < cvcs.size(); i++) {
       if (!cvcs[i]->is_enabled()) continue;
+        if (cvm::debug())
+          cvm::log("Colvar component no. "+cvm::to_str(i+1)+
+            " within colvar \""+this->name+"\" has Jacobian derivative"+
+            cvm::to_str((cvcs[i])->Jacobian_derivative(),
+            cvm::cv_width, cvm::cv_prec)+".\n");
       // linear combination is assumed
       fj += (cvcs[i])->Jacobian_derivative() * (cvcs[i])->sup_coeff / active_cvc_square_norm;
     }
@@ -1128,15 +1165,16 @@ cvm::real colvar::update_forces_energy()
   if (is_enabled(f_cv_Jacobian)) {
     // the instantaneous Jacobian force was not included in the reported total force;
     // instead, it is subtracted from the applied force (silent Jacobian correction)
+    // This requires the Jacobian term for the *current* timestep
     if (is_enabled(f_cv_hide_Jacobian))
       f -= fj;
   }
 
-  if (is_enabled(f_cv_lower_wall) || is_enabled(f_cv_upper_wall)) {
+  // Wall force
+  colvarvalue fw(x);
+  fw.reset();
 
-    // Wall force
+  if (is_enabled(f_cv_lower_wall) || is_enabled(f_cv_upper_wall)) {
-    colvarvalue fw(x);
-    fw.reset();
 
     if (cvm::debug())
       cvm::log("Calculating wall forces for colvar \""+this->name+"\".\n");
@@ -1144,12 +1182,11 @@ cvm::real colvar::update_forces_energy()
     // For a periodic colvar, both walls may be applicable at the same time
     // in which case we pick the closer one
     if ( (!is_enabled(f_cv_upper_wall)) ||
-         (this->dist2(x_reported, lower_wall) < this->dist2(x_reported, upper_wall)) ) {
+         (this->dist2(x, lower_wall) < this->dist2(x, upper_wall)) ) {
 
-      cvm::real const grad = this->dist2_lgrad(x_reported, lower_wall);
+      cvm::real const grad = this->dist2_lgrad(x, lower_wall);
       if (grad < 0.0) {
         fw = -0.5 * lower_wall_k * grad;
-        f += fw;
         if (cvm::debug())
           cvm::log("Applying a lower wall force("+
                     cvm::to_str(fw)+") to \""+this->name+"\".\n");
@@ -1157,10 +1194,9 @@ cvm::real colvar::update_forces_energy()
 
     } else {
 
-      cvm::real const grad = this->dist2_lgrad(x_reported, upper_wall);
+      cvm::real const grad = this->dist2_lgrad(x, upper_wall);
       if (grad > 0.0) {
         fw = -0.5 * upper_wall_k * grad;
-        f += fw;
         if (cvm::debug())
           cvm::log("Applying an upper wall force("+
                     cvm::to_str(fw)+") to \""+this->name+"\".\n");
@@ -1168,6 +1204,9 @@ cvm::real colvar::update_forces_energy()
     }
   }
 
+  // At this point f is the force f from external biases that will be applied to the
+  // extended variable if there is one
+
   if (is_enabled(f_cv_extended_Lagrangian)) {
 
     cvm::real dt = cvm::dt();
@@ -1175,11 +1214,12 @@ cvm::real colvar::update_forces_energy()
     f_ext.reset();
 
     // the total force is applied to the fictitious mass, while the
-    // atoms only feel the harmonic force
+    // atoms only feel the harmonic force + wall force
     // fr: bias force on extended variable (without harmonic spring), for output in trajectory
     // f_ext: total force on extended variable (including harmonic spring)
-    // f: - initially, external biasing force (including wall forces)
+    // f: - initially, external biasing force
-    //    - after this code block, colvar force to be applied to atomic coordinates, ie. spring force
+    //    - after this code block, colvar force to be applied to atomic coordinates
+    //      ie. spring force + wall force
     fr    = f;
     f_ext = f + (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x);
     f     =     (-0.5 * ext_force_k) * this->dist2_rgrad(xr, x);
@@ -1209,9 +1249,16 @@ cvm::real colvar::update_forces_energy()
     vr  += (0.5 * dt) * f_ext / ext_mass;
     xr  += dt * vr;
     xr.apply_constraints();
-    if (this->b_periodic) this->wrap(xr);
+    if (this->is_enabled(f_cv_periodic)) this->wrap(xr);
   }
 
+  // TODO remove the wall force
+  f += fw;
+  // Now adding the force on the actual colvar (for those biases who
+  // bypass the extended Lagrangian mass)
+  f += fb_actual;
+
+  // Store force to be applied, possibly summed over several timesteps
   f_accumulated += f;
 
   if (is_enabled(f_cv_fdiff_velocity)) {
@@ -1428,14 +1475,15 @@ std::istream & colvar::read_restart(std::istream &is)
     }
   }
 
-  if ( !(get_keyval(conf, "x", x,
+  if ( !(get_keyval(conf, "x", x, x, colvarparse::parse_silent)) ) {
-                    colvarvalue(x.type()), colvarparse::parse_silent)) ) {
     cvm::log("Error: restart file does not contain "
              "the value of the colvar \""+
              name+"\" .\n");
   } else {
     cvm::log("Restarting collective variable \""+name+"\" from value: "+
              cvm::to_str(x)+"\n");
+    x_restart = x;
+    after_restart = true;
   }
 
   if (is_enabled(f_cv_extended_Lagrangian)) {

lib/colvars/colvar.h

@@ -1,5 +1,12 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #ifndef COLVAR_H
 #define COLVAR_H
 
@@ -170,6 +177,9 @@ public:
   /// the biases are updated
   colvarvalue fb;
 
+  /// \brief Bias force to the actual value (only useful with extended Lagrangian)
+  colvarvalue fb_actual;
+
   /// \brief Total \em applied force; fr (if extended_lagrangian
   /// is defined), fb (if biases are applied) and the walls' forces
   /// (if defined) contribute to it
@@ -183,13 +193,9 @@ public:
   colvarvalue ft;
 
 
-  /// Period, if it is a constant
+  /// Period, if this variable is periodic
   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
@@ -290,6 +296,9 @@ public:
   /// Add to the total force from biases
   void add_bias_force(colvarvalue const &force);
 
+  /// Apply a force to the actual value (only meaningful with extended Lagrangian)
+  void add_bias_force_actual_value(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()
@@ -386,6 +395,12 @@ protected:
   /// Previous value (to calculate velocities during analysis)
   colvarvalue            x_old;
 
+  /// Value read from the most recent state file (if any)
+  colvarvalue            x_restart;
+
+  /// True if a state file was just read
+  bool                   after_restart;
+
   /// Time series of values and velocities used in correlation
   /// functions
   std::list< std::list<colvarvalue> > acf_x_history, acf_v_history;
@@ -577,9 +592,20 @@ inline void colvar::add_bias_force(colvarvalue const &force)
 }
 
 
+inline void colvar::add_bias_force_actual_value(colvarvalue const &force)
+{
+  if (cvm::debug()) {
+    cvm::log("Adding biasing force "+cvm::to_str(force)+" to colvar \""+name+"\".\n");
+  }
+  fb_actual += force;
+}
+
+
 inline void colvar::reset_bias_force() {
   fb.type(value());
   fb.reset();
+  fb_actual.type(value());
+  fb_actual.reset();
 }
 
 #endif

lib/colvars/colvaratoms.cpp

@@ -1,5 +1,12 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #include "colvarmodule.h"
 #include "colvarparse.h"
 #include "colvaratoms.h"
@@ -171,7 +178,10 @@ int cvm::atom_group::remove_atom(cvm::atom_iter ai)
 
 int cvm::atom_group::init()
 {
-  if (!key.size()) key = "atoms";
+  if (!key.size()) key = "unnamed";
+  description = "atom group " + key;
+  // These will be overwritten by parse(), if initializing from a config string
+
   atoms.clear();
 
   // TODO: check with proxy whether atom forces etc are available
@@ -179,6 +189,7 @@ int cvm::atom_group::init()
 
   index = -1;
 
+  b_dummy = false;
   b_center = false;
   b_rotate = false;
   b_user_defined_fit = false;
@@ -440,6 +451,7 @@ int cvm::atom_group::parse(std::string const &conf)
 
   if (b_print_atom_ids) {
     cvm::log("Internal definition of the atom group:\n");
+    cvm::log(print_atom_ids());
   }
 
   cvm::decrease_depth();

lib/colvars/colvaratoms.h

@@ -1,5 +1,12 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #ifndef COLVARATOMS_H
 #define COLVARATOMS_H
 
@@ -64,7 +71,7 @@ public:
 
   /// \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)
+  /// the system topology (1-based)
   atom(int atom_number);
 
   /// \brief Initialize an atom for collective variable calculation
@@ -453,6 +460,8 @@ public:
   /// 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.
+  /// This function will be phased out eventually, in favor of
+  /// apply_colvar_force() once that is implemented for non-scalar values
   void apply_force(cvm::rvector const &force);
 
 };

lib/colvars/colvarbias.cpp

@@ -1,5 +1,12 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #include "colvarmodule.h"
 #include "colvarvalue.h"
 #include "colvarbias.h"
@@ -29,6 +36,7 @@ colvarbias::colvarbias(char const *key)
   has_data = false;
   b_output_energy = false;
   reset();
+  state_file_step = 0;
 
   // Start in active state by default
   enable(f_cvb_active);
@@ -141,20 +149,25 @@ int colvarbias::clear()
 int colvarbias::add_colvar(std::string const &cv_name)
 {
   if (colvar *cv = cvm::colvar_by_name(cv_name)) {
-    // Removed this as nor all biases apply forces eg histogram
+
-    // cv->enable(colvar::task_gradients);
     if (cvm::debug()) {
       cvm::log("Applying this bias to collective variable \""+
                cv->name+"\".\n");
     }
+
     colvars.push_back(cv);
 
     colvar_forces.push_back(colvarvalue());
     colvar_forces.back().type(cv->value()); // make sure each force is initialized to zero
+    colvar_forces.back().is_derivative(); // colvar constraints are not applied to the force
     colvar_forces.back().reset();
 
     cv->biases.push_back(this); // add back-reference to this bias to colvar
 
+    if (is_enabled(f_cvb_apply_force)) {
+      cv->enable(f_cv_gradient);
+    }
+
     // Add dependency link.
     // All biases need at least the value of each colvar
     // although possibly not at all timesteps
@@ -165,6 +178,7 @@ int colvarbias::add_colvar(std::string const &cv_name)
                cv_name+"\".\n", INPUT_ERROR);
     return INPUT_ERROR;
   }
+
   return COLVARS_OK;
 }
 
@@ -190,16 +204,19 @@ void colvarbias::communicate_forces()
 }
 
 
-void colvarbias::change_configuration(std::string const &conf)
+int colvarbias::change_configuration(std::string const &conf)
 {
-  cvm::error("Error: change_configuration() not implemented.\n");
+  cvm::error("Error: change_configuration() not implemented.\n",
+             COLVARS_NOT_IMPLEMENTED);
+  return COLVARS_NOT_IMPLEMENTED;
 }
 
 
 cvm::real colvarbias::energy_difference(std::string const &conf)
 {
-  cvm::error("Error: energy_difference() not implemented.\n");
+  cvm::error("Error: energy_difference() not implemented.\n",
-  return 0.;
+             COLVARS_NOT_IMPLEMENTED);
+  return 0.0;
 }
 
 
@@ -225,6 +242,118 @@ int colvarbias::replica_share()
   return COLVARS_NOT_IMPLEMENTED;
 }
 
+
+std::string const colvarbias::get_state_params() const
+{
+  std::ostringstream os;
+  os << "step " << cvm::step_absolute() << "\n"
+     << "name " << this->name << "\n";
+  return os.str();
+}
+
+
+int colvarbias::set_state_params(std::string const &conf)
+{
+  std::string new_name = "";
+  if (colvarparse::get_keyval(conf, "name", new_name,
+                              std::string(""), colvarparse::parse_silent) &&
+      (new_name != this->name)) {
+    cvm::error("Error: in the state file, the "
+               "\""+bias_type+"\" block has a different name, \""+new_name+
+               "\": different system?\n", INPUT_ERROR);
+  }
+
+  if (name.size() == 0) {
+    cvm::error("Error: \""+bias_type+"\" block within the restart file "
+               "has no identifiers.\n", INPUT_ERROR);
+  }
+
+  colvarparse::get_keyval(conf, "step", state_file_step,
+                          cvm::step_absolute(), colvarparse::parse_silent);
+
+  return COLVARS_OK;
+}
+
+
+std::ostream & colvarbias::write_state(std::ostream &os)
+{
+  if (cvm::debug()) {
+    cvm::log("Writing state file for bias \""+name+"\"\n");
+  }
+  os.setf(std::ios::scientific, std::ios::floatfield);
+  os.precision(cvm::cv_prec);
+  os << bias_type << " {\n"
+     << "  configuration {\n";
+  std::istringstream is(get_state_params());
+  std::string line;
+  while (std::getline(is, line)) {
+    os << "    " << line << "\n";
+  }
+  os << "  }\n";
+  write_state_data(os);
+  os << "}\n\n";
+  return os;
+}
+
+
+std::istream & colvarbias::read_state(std::istream &is)
+{
+  size_t const start_pos = is.tellg();
+
+  std::string key, brace, conf;
+  if ( !(is >> key)   || !(key == bias_type) ||
+       !(is >> brace) || !(brace == "{") ||
+       !(is >> colvarparse::read_block("configuration", conf)) ||
+       (set_state_params(conf) != COLVARS_OK) ) {
+    cvm::error("Error: in reading state configuration for \""+bias_type+"\" bias \""+
+               this->name+"\" at position "+
+               cvm::to_str(is.tellg())+" in stream.\n", INPUT_ERROR);
+    is.clear();
+    is.seekg(start_pos, std::ios::beg);
+    is.setstate(std::ios::failbit);
+    return is;
+  }
+
+  if (!read_state_data(is)) {
+    cvm::error("Error: in reading state data for \""+bias_type+"\" bias \""+
+               this->name+"\" at position "+
+               cvm::to_str(is.tellg())+" in stream.\n", INPUT_ERROR);
+    is.clear();
+    is.seekg(start_pos, std::ios::beg);
+    is.setstate(std::ios::failbit);
+  }
+
+  is >> brace;
+  if (brace != "}") {
+    cvm::error("Error: corrupt restart information for \""+bias_type+"\" bias \""+
+               this->name+"\": no matching brace at position "+
+               cvm::to_str(is.tellg())+" in stream.\n");
+    is.setstate(std::ios::failbit);
+  }
+
+  return is;
+}
+
+
+std::istream & colvarbias::read_state_data_key(std::istream &is, char const *key)
+{
+  size_t const start_pos = is.tellg();
+  std::string key_in;
+  if ( !(is >> key_in) ||
+       !(key_in == to_lower_cppstr(std::string(key))) ) {
+    cvm::error("Error: in reading restart configuration for "+
+               bias_type+" bias \""+this->name+"\" at position "+
+               cvm::to_str(is.tellg())+" in stream.\n", INPUT_ERROR);
+    is.clear();
+    is.seekg(start_pos, std::ios::beg);
+    is.setstate(std::ios::failbit);
+    return is;
+  }
+  return is;
+}
+
+
+
 std::ostream & colvarbias::write_traj_label(std::ostream &os)
 {
   os << " ";

lib/colvars/colvarbias.h

@@ -1,5 +1,12 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #ifndef COLVARBIAS_H
 #define COLVARBIAS_H
 
@@ -9,7 +16,8 @@
 
 
 /// \brief Collective variable bias, base class
-class colvarbias : public colvarparse, public colvardeps {
+class colvarbias
+  : public virtual colvarparse, public virtual colvardeps {
 public:
 
   /// Name of this bias
@@ -24,6 +32,12 @@ public:
   /// Add a new collective variable to this bias
   int add_colvar(std::string const &cv_name);
 
+  /// Add a new collective variable to this bias
+  size_t number_of_colvars() const
+  {
+    return colvars.size();
+  }
+
   /// Retrieve colvar values and calculate their biasing forces
   /// Return bias energy
   virtual int update();
@@ -31,7 +45,7 @@ public:
   // TODO: move update_bias here (share with metadynamics)
 
   /// Load new configuration - force constant and/or centers only
-  virtual void change_configuration(std::string const &conf);
+  virtual int change_configuration(std::string const &conf);
 
   /// Calculate change in energy from using alternate configuration
   virtual cvm::real energy_difference(std::string const &conf);
@@ -49,7 +63,7 @@ public:
   virtual void analyze() {}
 
   /// Send forces to the collective variables
-  void communicate_forces();
+  virtual void communicate_forces();
 
   /// \brief Constructor
   colvarbias(char const *key);
@@ -60,13 +74,11 @@ public:
   /// \brief Set to zero all mutable data
   virtual int reset();
 
-protected:
+private:
 
   /// Default constructor
   colvarbias();
 
-private:
-
   /// Copy constructor
   colvarbias(colvarbias &);
 
@@ -78,28 +90,59 @@ public:
   /// Destructor
   virtual ~colvarbias();
 
-  /// Read the bias configuration from a restart file
+  /// Write the values of specific mutable properties to a string
-  virtual std::istream & read_restart(std::istream &is) = 0;
+  virtual std::string const get_state_params() const;
 
-  /// Write the bias configuration to a restart file
+  /// Read the values of specific mutable properties from a string
-  virtual std::ostream & write_restart(std::ostream &os) = 0;
+  virtual int set_state_params(std::string const &state_conf);
+
+  /// Write all mutable data not already written by get_state_params()
+  virtual std::ostream & write_state_data(std::ostream &os)
+  {
+    return os;
+  }
+
+  /// Read all mutable data not already set by set_state_params()
+  virtual std::istream & read_state_data(std::istream &is)
+  {
+    return is;
+  }
+
+  /// Read a keyword from the state data (typically a header)
+  std::istream & read_state_data_key(std::istream &is, char const *key);
+
+  /// Write the bias configuration to a restart file or other stream
+  virtual std::ostream & write_state(std::ostream &os);
+
+  /// Read the bias configuration from a restart file or other stream
+  virtual std::istream & read_state(std::istream &is);
 
   /// Write a label to the trajectory file (comment line)
   virtual std::ostream & write_traj_label(std::ostream &os);
 
-  /// (Re)initialize the output files (does not write them yet)
-  virtual int setup_output() { return COLVARS_OK; }
-
   /// Output quantities such as the bias energy to the trajectory file
   virtual std::ostream & write_traj(std::ostream &os);
 
-  /// Write output files (if defined, e.g. in analysis mode)
+  /// (Re)initialize the output files (does not write them yet)
+  virtual int setup_output()
+  {
+    return COLVARS_OK;
+  }
+
+  /// Write any output files that this bias may have (e.g. PMF files)
   virtual int write_output_files()
   {
     return COLVARS_OK;
   }
 
-  inline cvm::real get_energy() {
+  /// If this bias is communicating with other replicas through files, send it to them
+  virtual int write_state_to_replicas()
+  {
+    return COLVARS_OK;
+  }
+
+  inline cvm::real get_energy()
+  {
     return bias_energy;
   }
 
@@ -107,9 +150,11 @@ public:
   static std::vector<feature *> cvb_features;
 
   /// \brief Implementation of the feature list accessor for colvarbias
-  virtual std::vector<feature *> &features() {
+  virtual std::vector<feature *> &features()
+  {
     return cvb_features;
   }
+
 protected:
 
   /// \brief Pointers to collective variables to which the bias is
@@ -130,6 +175,9 @@ protected:
   /// (for history-dependent biases)
   bool                     has_data;
 
+  /// \brief Step number read from the last state file
+  size_t                   state_file_step;
+
 };
 
 #endif

lib/colvars/colvarbias_abf.cpp

@@ -1,5 +1,12 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #include "colvarmodule.h"
 #include "colvar.h"
 #include "colvarbias_abf.h"
@@ -23,6 +30,8 @@ int colvarbias_abf::init(std::string const &conf)
 {
   colvarbias::init(conf);
 
+  provide(f_cvb_scalar_variables);
+  enable(f_cvb_scalar_variables);
   provide(f_cvb_history_dependent);
 
   // TODO relax this in case of VMD plugin
@@ -590,16 +599,10 @@ void colvarbias_abf::read_gradients_samples()
 }
 
 
-std::ostream & colvarbias_abf::write_restart(std::ostream& os)
+std::ostream & colvarbias_abf::write_state_data(std::ostream& os)
 {
-
   std::ios::fmtflags flags(os.flags());
 
-  os << "abf {\n"
-     << "  configuration {\n"
-     << "    name " << this->name << "\n";
-  os << "  }\n";
-
   os.setf(std::ios::fmtflags(0), std::ios::floatfield); // default floating-point format
   os << "\nsamples\n";
   samples->write_raw(os, 8);
@@ -617,117 +620,47 @@ std::ostream & colvarbias_abf::write_restart(std::ostream& os)
     z_gradients->write_raw(os, 8);
   }
 
-  os << "}\n\n";
-
   os.flags(flags);
   return os;
 }
 
 
-std::istream & colvarbias_abf::read_restart(std::istream& is)
+std::istream & colvarbias_abf::read_state_data(std::istream& is)
 {
   if ( input_prefix.size() > 0 ) {
     cvm::error("ERROR: cannot provide both inputPrefix and a colvars state file.\n", INPUT_ERROR);
   }
 
-  size_t const start_pos = is.tellg();
+  if (! read_state_data_key(is, "samples")) {
-
-  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::error("Error: in the restart file, the "
-                      "\"abf\" block has wrong name(" + name + ")\n");
-  if ( name == "" ) {
-    cvm::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)) {
-    is.clear();
-    is.seekg(start_pos, std::ios::beg);
-    is.setstate(std::ios::failbit);
     return is;
   }
 
-  if ( !(is >> key)   || !(key == "gradient")) {
+  if (! read_state_data_key(is, "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)) {
-    is.clear();
-    is.seekg(start_pos, std::ios::beg);
-    is.setstate(std::ios::failbit);
     return is;
   }
 
   if (z_gradients) {
-    if ( !(is >> key) || !(key == "z_samples")) {
+
-      cvm::log("Error: in reading restart configuration for ABF bias \""+
+    if (! read_state_data_key(is, "z_samples")) {
-                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 (! z_samples->read_raw(is)) {
-      is.clear();
-      is.seekg(start_pos, std::ios::beg);
-      is.setstate(std::ios::failbit);
       return is;
     }
 
-    if ( !(is >> key)   || !(key == "z_gradient")) {
+    if (! read_state_data_key(is, "z_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 (! z_gradients->read_raw(is)) {
-      is.clear();
-      is.seekg(start_pos, std::ios::beg);
-      is.setstate(std::ios::failbit);
       return is;
     }
   }
-  is >> brace;
+
-  if (brace != "}") {
-    cvm::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;
 }

lib/colvars/colvarbias_abf.h

@@ -1,5 +1,12 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #ifndef COLVARBIAS_ABF_H
 #define COLVARBIAS_ABF_H
 
@@ -107,8 +114,8 @@ private:
   /// Read human-readable FE gradients and sample count (if not using restart)
   void		  read_gradients_samples();
 
-  std::istream& read_restart(std::istream&);
+  std::istream& read_state_data(std::istream&);
-  std::ostream& write_restart(std::ostream&);
+  std::ostream& write_state_data(std::ostream&);
 };
 
 #endif

lib/colvars/colvarbias_alb.cpp

@@ -1,5 +1,12 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
@@ -33,6 +40,9 @@ int colvarbias_alb::init(std::string const &conf)
 {
   colvarbias::init(conf);
 
+  provide(f_cvb_scalar_variables);
+  enable(f_cvb_scalar_variables);
+
   provide(f_cvb_history_dependent);
 
   size_t i;
@@ -239,37 +249,8 @@ int colvarbias_alb::update()
 }
 
 
-std::istream & colvarbias_alb::read_restart(std::istream &is)
+int colvarbias_alb::set_state_params(std::string const &conf)
 {
-  size_t const start_pos = is.tellg();
-
-  cvm::log("Restarting adaptive linear bias \""+
-            this->name+"\".\n");
-
-  std::string key, brace, conf;
-  if ( !(is >> key)   || !(key == "ALB") ||
-       !(is >> brace) || !(brace == "{") ||
-       !(is >> colvarparse::read_block("configuration", conf)) ) {
-
-    cvm::log("Error: in reading restart configuration for restraint 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 "
-                      "\"ALB\" block has a wrong name\n");
-  if (name.size() == 0) {
-    cvm::fatal_error("Error: \"ALB\" block in the restart file "
-                      "has no identifiers.\n");
-  }
-
   if (!get_keyval(conf, "setCoupling", set_coupling))
     cvm::fatal_error("Error: current setCoupling  is missing from the restart.\n");
 
@@ -299,23 +280,13 @@ std::istream & colvarbias_alb::read_restart(std::istream &is)
   if (!get_keyval(conf, "b_equilibration", b_equilibration))
     cvm::fatal_error("Error: current updateCalls is missing from the restart.\n");
 
-  is >> brace;
+  return COLVARS_OK;
-  if (brace != "}") {
-    cvm::fatal_error("Error: corrupt restart information for adaptive linear 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_alb::write_restart(std::ostream &os)
+std::string const colvarbias_alb::get_state_params() const
 {
-  os << "ALB {\n"
+  std::ostringstream os;
-     << "  configuration {\n"
-     << "    name " << this->name << "\n";
   os << "    setCoupling ";
   size_t i;
   for (i = 0; i < colvars.size(); i++) {
@@ -358,10 +329,7 @@ std::ostream & colvarbias_alb::write_restart(std::ostream &os)
   else
     os << "    b_equilibration no\n";
 
-  os << "  }\n"
+  return os.str();
-     << "}\n\n";
-
-  return os;
 }
 
 

lib/colvars/colvarbias_alb.h

@@ -1,10 +1,18 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #ifndef COLVARBIAS_ALB_H
 #define COLVARBIAS_ALB_H
 
 #include "colvar.h"
-#include "colvarbias_restraint.h"
+#include "colvarbias.h"
+
 
 class colvarbias_alb : public colvarbias {
 
@@ -15,8 +23,8 @@ public:
   virtual int init(std::string const &conf);
   virtual int update();
 
-  virtual std::istream & read_restart(std::istream &is);
+  virtual std::string const get_state_params() const;
-  virtual std::ostream & write_restart(std::ostream &os);
+  virtual int set_state_params(std::string const &conf);
   virtual std::ostream & write_traj_label(std::ostream &os);
   virtual std::ostream & write_traj(std::ostream &os);
 

lib/colvars/colvarbias_histogram.cpp

@@ -1,5 +1,12 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #include "colvarmodule.h"
 #include "colvar.h"
 #include "colvarbias_histogram.h"
@@ -17,6 +24,9 @@ int colvarbias_histogram::init(std::string const &conf)
 {
   colvarbias::init(conf);
 
+  provide(f_cvb_scalar_variables);
+  enable(f_cvb_scalar_variables);
+
   provide(f_cvb_history_dependent);
   enable(f_cvb_history_dependent);
 
@@ -196,78 +206,25 @@ int colvarbias_histogram::write_output_files()
 }
 
 
-std::istream & colvarbias_histogram::read_restart(std::istream& is)
+std::istream & colvarbias_histogram::read_state_data(std::istream& is)
 {
-  size_t const start_pos = is.tellg();
+  if (! read_state_data_key(is, "grid")) {
-
-  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::error("Error: in the restart file, the "
-                      "\"histogram\" block has a wrong name: different system?\n");
-  if ( (id == -1) && (name == "") ) {
-    cvm::error("Error: \"histogram\" block in the restart file "
-                      "has no name.\n");
-  }
-
-  if ( !(is >> key)   || !(key == "grid")) {
-    cvm::error("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)) {
-    is.clear();
-    is.seekg(start_pos, std::ios::beg);
-    is.setstate(std::ios::failbit);
     return is;
   }
 
-  is >> brace;
-  if (brace != "}") {
-    cvm::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)
+
+std::ostream & colvarbias_histogram::write_state_data(std::ostream& os)
 {
   std::ios::fmtflags flags(os.flags());
   os.setf(std::ios::fmtflags(0), std::ios::floatfield);
-
-  os << "histogram {\n"
-     << "  configuration {\n"
-     << "    name " << this->name << "\n";
-  os << "  }\n";
-
   os << "grid\n";
   grid->write_raw(os, 8);
-
-  os << "}\n\n";
-
   os.flags(flags);
   return os;
 }

lib/colvars/colvarbias_histogram.h

@@ -1,5 +1,12 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #ifndef COLVARBIAS_HISTOGRAM_H
 #define COLVARBIAS_HISTOGRAM_H
 
@@ -35,8 +42,8 @@ protected:
   /// If colvar_array_size is larger than 1, weigh each one by this number before accumulating the histogram
   std::vector<cvm::real> weights;
 
-  virtual std::istream& read_restart(std::istream&);
+  virtual std::istream & read_state_data(std::istream &is);
-  virtual std::ostream& write_restart(std::ostream&);
+  virtual std::ostream & write_state_data(std::ostream &os);
 };
 
 #endif

lib/colvars/colvarbias_meta.cpp

@@ -1,5 +1,12 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #include <iostream>
 #include <sstream>
 #include <fstream>
@@ -26,10 +33,9 @@
 
 
 colvarbias_meta::colvarbias_meta(char const *key)
-  : colvarbias(key),
+  : colvarbias(key)
-    new_hills_begin(hills.end()),
-    state_file_step(0)
 {
+  new_hills_begin = hills.end();
 }
 
 
@@ -164,11 +170,32 @@ int colvarbias_meta::init(std::string const &conf)
   target_dist = NULL;
   get_keyval(conf, "ebMeta", ebmeta, false);
   if(ebmeta){
+    if (use_grids && expand_grids) {
+      cvm::fatal_error("Error: expandBoundaries is not supported with "
+                       "ebMeta please allocate wide enough boundaries for "
+                       "each colvar ahead of time and set targetdistfile "
+                       "accordingly. \n");
+    }
     target_dist = new colvar_grid_scalar();
     target_dist->init_from_colvars(colvars);
     get_keyval(conf, "targetdistfile", target_dist_file);
     std::ifstream targetdiststream(target_dist_file.c_str());
     target_dist->read_multicol(targetdiststream);
+    cvm::real min_val = target_dist->minimum_value();
+    if(min_val<0){
+      cvm::error("Error: Target distribution of ebMeta "
+                 "has negative values!.\n", INPUT_ERROR);
+    }
+    cvm::real min_pos_val = target_dist->minimum_pos_value();
+    if(min_pos_val<=0){
+      cvm::error("Error: Target distribution of ebMeta has negative "
+                 "or zero minimum positive value!.\n", INPUT_ERROR);
+    }
+    if(min_val==0){
+      cvm::log("WARNING: Target distribution has zero values.\n");
+      cvm::log("Zeros will be converted to the minimum positive value.\n");
+      target_dist->remove_zeros(min_pos_val);
+    }
     // normalize target distribution and multiply by effective volume = exp(differential entropy)
     target_dist->multiply_constant(1.0/target_dist->integral());
     cvm::real volume = std::exp(target_dist->entropy());
@@ -405,7 +432,9 @@ int colvarbias_meta::update()
     if (ebmeta) {
        hills_scale *= 1.0/target_dist->value(target_dist->get_colvars_index());
        if(cvm::step_absolute() <= long(ebmeta_equil_steps)) {
-         cvm::real const hills_lambda=(cvm::real(ebmeta_equil_steps - cvm::step_absolute()))/(cvm::real(ebmeta_equil_steps));
+         cvm::real const hills_lambda =
+           (cvm::real(long(ebmeta_equil_steps) - cvm::step_absolute())) /
+           (cvm::real(ebmeta_equil_steps));
          hills_scale = hills_lambda + (1-hills_lambda)*hills_scale;
        }
     }
@@ -972,7 +1001,7 @@ void colvarbias_meta::read_replica_files()
                (replicas[ir])->replica_state_file+"\".\n");
 
       std::ifstream is((replicas[ir])->replica_state_file.c_str());
-      if (! (replicas[ir])->read_restart(is)) {
+      if (! (replicas[ir])->read_state(is)) {
         cvm::log("Reading from file \""+(replicas[ir])->replica_state_file+
                  "\" failed or incomplete: will try again in "+
                  cvm::to_str(replica_update_freq)+" steps.\n");
@@ -1068,63 +1097,24 @@ void colvarbias_meta::read_replica_files()
 }
 
 
-// **********************************************************************
+int colvarbias_meta::set_state_params(std::string const &state_conf)
-// input functions
-// **********************************************************************
-
-
-std::istream & colvarbias_meta::read_restart(std::istream& is)
 {
-  size_t const start_pos = is.tellg();
+  std::string new_replica = "";
-
+  if (colvarparse::get_keyval(state_conf, "replicaID", new_replica,
-  if (comm == single_replica) {
+                              std::string(""), colvarparse::parse_silent) &&
-    // if using a multiple replicas scheme, output messages
+      (new_replica != this->replica_id)) {
-    // are printed before and after calling this function
+    cvm::error("Error: in the state file, the "
-    cvm::log("Restarting metadynamics bias \""+this->name+"\""+
+               "\"metadynamics\" block has a different replicaID: different system?\n",
-             ".\n");
+               INPUT_ERROR);
-  }
+    return INPUT_ERROR;
-  std::string key, brace, conf;
-  if ( !(is >> key)   || !(key == "metadynamics") ||
-       !(is >> brace) || !(brace == "{") ||
-       !(is >> colvarparse::read_block("configuration", conf)) ) {
-
-    if (comm == single_replica)
-      cvm::log("Error: in reading restart configuration for metadynamics bias \""+
-               this->name+"\""+
-               ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
-               (replica_state_file_in_sync ? ("at position "+
-                                              cvm::to_str(start_pos)+
-                                              " in the state file") : "")+".\n");
-    is.clear();
-    is.seekg(start_pos, std::ios::beg);
-    is.setstate(std::ios::failbit);
-    return is;
   }
 
-  std::string name = "";
+  return COLVARS_OK;
-  if ( colvarparse::get_keyval(conf, "name", name,
+}
-                               std::string(""), colvarparse::parse_silent) &&
-       (name != this->name) )
-    cvm::fatal_error("Error: in the restart file, the "
-                     "\"metadynamics\" block has a different name: different system?\n");
 
-  if (name.size() == 0) {
-    cvm::fatal_error("Error: \"metadynamics\" block within the restart file "
-                     "has no identifiers.\n");
-  }
-
-  if (comm != single_replica) {
-    std::string replica = "";
-    if (colvarparse::get_keyval(conf, "replicaID", replica,
-                                std::string(""), colvarparse::parse_silent) &&
-        (replica != this->replica_id))
-      cvm::fatal_error("Error: in the restart file, the "
-                       "\"metadynamics\" block has a different replicaID: different system?\n");
-
-    colvarparse::get_keyval(conf, "step", state_file_step,
-                            cvm::step_absolute(), colvarparse::parse_silent);
-  }
 
+std::istream & colvarbias_meta::read_state_data(std::istream& is)
+{
   bool grids_from_restart_file = use_grids;
 
   if (use_grids) {
@@ -1155,6 +1145,7 @@ std::istream & colvarbias_meta::read_restart(std::istream& is)
     }
 
     size_t const hills_energy_pos = is.tellg();
+    std::string key;
     if (!(is >> key)) {
       if (hills_energy_backup != NULL) {
         delete hills_energy;
@@ -1316,17 +1307,6 @@ std::istream & colvarbias_meta::read_restart(std::istream& is)
     }
   }
 
-  is >> brace;
-  if (brace != "}") {
-    cvm::log("Incomplete restart information for metadynamics bias \""+
-             this->name+"\""+
-             ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
-             ": no closing brace at position "+
-             cvm::to_str(is.tellg())+" in the file.\n");
-    is.setstate(std::ios::failbit);
-    return is;
-  }
-
   if (cvm::debug())
     cvm::log("colvarbias_meta::read_restart() done\n");
 
@@ -1424,13 +1404,6 @@ std::istream & colvarbias_meta::read_hill(std::istream &is)
 }
 
 
-
-
-// **********************************************************************
-// output functions
-// **********************************************************************
-
-
 int colvarbias_meta::setup_output()
 {
 
@@ -1537,16 +1510,17 @@ int colvarbias_meta::setup_output()
 }
 
 
-std::ostream & colvarbias_meta::write_restart(std::ostream& os)
+std::string const colvarbias_meta::get_state_params() const
 {
-  os << "metadynamics {\n"
+  std::ostringstream os;
-     << "  configuration {\n"
-     << "    step " << cvm::step_absolute() << "\n"
-     << "    name " << this->name << "\n";
   if (this->comm != single_replica)
-    os << "    replicaID " << this->replica_id << "\n";
+    os << "replicaID " << this->replica_id << "\n";
-  os << "  }\n\n";
+  return (colvarbias::get_state_params() + os.str());
+}
 
+
+std::ostream & colvarbias_meta::write_state_data(std::ostream& os)
+{
   if (use_grids) {
 
     // this is a very good time to project hills, if you haven't done
@@ -1578,8 +1552,12 @@ std::ostream & colvarbias_meta::write_restart(std::ostream& os)
     }
   }
 
-  os << "}\n\n";
+  return os;
+}
 
+
+int colvarbias_meta::write_state_to_replicas()
+{
   if (comm != single_replica) {
     write_replica_state_file();
     // schedule to reread the state files of the other replicas (they
@@ -1588,12 +1566,16 @@ std::ostream & colvarbias_meta::write_restart(std::ostream& os)
       (replicas[ir])->replica_state_file_in_sync = false;
     }
   }
+  return COLVARS_OK;
+}
 
+
+int colvarbias_meta::write_output_files()
+{
   if (dump_fes) {
     write_pmf();
   }
-
+  return COLVARS_OK;
-  return os;
 }
 
 
@@ -1665,53 +1647,67 @@ void colvarbias_meta::write_pmf()
 
 
 
-void colvarbias_meta::write_replica_state_file()
+int colvarbias_meta::write_replica_state_file()
 {
-  // write down also the restart for the other replicas: TODO: this
+  if (cvm::debug()) {
-  // is duplicated code, that could be cleaned up later
+    cvm::log("Writing replica state file for bias \""+name+"\"\n");
+  }
+  // write down also the restart for the other replicas
   cvm::backup_file(replica_state_file.c_str());
-  cvm::ofstream rep_state_os(replica_state_file.c_str());
+  std::ostream *rep_state_os = cvm::proxy->output_stream(replica_state_file);
-  if (!rep_state_os.is_open())
+  if (rep_state_os == NULL) {
-    cvm::fatal_error("Error: in opening file \""+
+    cvm::error("Error: in opening file \""+
-                     replica_state_file+"\" for writing.\n");
+               replica_state_file+"\" for writing.\n", FILE_ERROR);
-
+    return FILE_ERROR;
-  rep_state_os.setf(std::ios::scientific, std::ios::floatfield);
-  rep_state_os << "\n"
-               << "metadynamics {\n"
-               << "  configuration {\n"
-               << "    name " << this->name << "\n"
-               << "    step " << cvm::step_absolute() << "\n";
-  if (this->comm != single_replica) {
-    rep_state_os << "    replicaID " << this->replica_id << "\n";
-  }
-  rep_state_os << "  }\n\n";
-  rep_state_os << "  hills_energy\n";
-  rep_state_os << std::setprecision(cvm::cv_prec)
-               << std::setw(cvm::cv_width);
-  hills_energy->write_restart(rep_state_os);
-  rep_state_os << "  hills_energy_gradients\n";
-  rep_state_os << std::setprecision(cvm::cv_prec)
-               << std::setw(cvm::cv_width);
-  hills_energy_gradients->write_restart(rep_state_os);
-
-  if ( (!use_grids) || keep_hills ) {
-    // write all hills currently in memory
-    for (std::list<hill>::const_iterator h = this->hills.begin();
-         h != this->hills.end();
-         h++) {
-      rep_state_os << *h;
-    }
-  } else {
-    // write just those that are near the grid boundaries
-    for (std::list<hill>::const_iterator h = this->hills_off_grid.begin();
-         h != this->hills_off_grid.end();
-         h++) {
-      rep_state_os << *h;
-    }
   }
 
-  rep_state_os << "}\n\n";
+  rep_state_os->setf(std::ios::scientific, std::ios::floatfield);
-  rep_state_os.close();
+
+  if (!write_state(*rep_state_os)) {
+    cvm::error("Error: in writing to file \""+
+               replica_state_file+"\".\n", FILE_ERROR);
+    cvm::proxy->close_output_stream(replica_state_file);
+    return FILE_ERROR;
+  }
+
+  cvm::proxy->close_output_stream(replica_state_file);
+
+  // rep_state_os.setf(std::ios::scientific, std::ios::floatfield);
+  // rep_state_os << "\n"
+  //              << "metadynamics {\n"
+  //              << "  configuration {\n"
+  //              << "    name " << this->name << "\n"
+  //              << "    step " << cvm::step_absolute() << "\n";
+  // if (this->comm != single_replica) {
+  //   rep_state_os << "    replicaID " << this->replica_id << "\n";
+  // }
+  // rep_state_os << "  }\n\n";
+  // rep_state_os << "  hills_energy\n";
+  // rep_state_os << std::setprecision(cvm::cv_prec)
+  //              << std::setw(cvm::cv_width);
+  // hills_energy->write_restart(rep_state_os);
+  // rep_state_os << "  hills_energy_gradients\n";
+  // rep_state_os << std::setprecision(cvm::cv_prec)
+  //              << std::setw(cvm::cv_width);
+  // hills_energy_gradients->write_restart(rep_state_os);
+
+  // if ( (!use_grids) || keep_hills ) {
+  //   // write all hills currently in memory
+  //   for (std::list<hill>::const_iterator h = this->hills.begin();
+  //        h != this->hills.end();
+  //        h++) {
+  //     rep_state_os << *h;
+  //   }
+  // } else {
+  //   // write just those that are near the grid boundaries
+  //   for (std::list<hill>::const_iterator h = this->hills_off_grid.begin();
+  //        h != this->hills_off_grid.end();
+  //        h++) {
+  //     rep_state_os << *h;
+  //   }
+  // }
+  // rep_state_os << "}\n\n";
+  // rep_state_os.close();
 
   // reopen the hills file
   replica_hills_os.close();
@@ -1721,8 +1717,11 @@ void colvarbias_meta::write_replica_state_file()
     cvm::fatal_error("Error: in opening file \""+
                      replica_hills_file+"\" for writing.\n");
   replica_hills_os.setf(std::ios::scientific, std::ios::floatfield);
+
+  return COLVARS_OK;
 }
 
+
 std::string colvarbias_meta::hill::output_traj()
 {
   std::ostringstream os;

lib/colvars/colvarbias_meta.h

@@ -1,5 +1,12 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #ifndef COLVARBIAS_META_H
 #define COLVARBIAS_META_H
 
@@ -31,10 +38,16 @@ public:
   virtual int init(std::string const &conf);
   virtual ~colvarbias_meta();
   virtual int update();
-  virtual std::istream & read_restart(std::istream &is);
+
-  virtual std::ostream & write_restart(std::ostream &os);
+  virtual std::string const get_state_params() const;
+  virtual int set_state_params(std::string const &state_conf);
+  virtual std::ostream & write_state_data(std::ostream &os);
+  virtual std::istream & read_state_data(std::istream &os);
+
   virtual int setup_output();
+  virtual int write_output_files();
   virtual void write_pmf();
+  virtual int write_state_to_replicas();
 
   class hill;
   typedef std::list<hill>::iterator hill_iter;
@@ -77,13 +90,6 @@ protected:
   /// 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
@@ -187,7 +193,7 @@ protected:
   virtual void read_replica_files();
 
   /// \brief Write data to other replicas
-  virtual void write_replica_state_file();
+  virtual int write_replica_state_file();
 
   /// \brief Additional, "mirror" metadynamics biases, to collect info
   /// from the other replicas

lib/colvars/colvarbias_restraint.h

@@ -1,36 +1,43 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #ifndef COLVARBIAS_RESTRAINT_H
 #define COLVARBIAS_RESTRAINT_H
 
 #include "colvarbias.h"
 
-/// \brief Bias restraint, optionally moving towards a target
+/// \brief Most general definition of a colvar restraint:
+/// see derived classes for specific types
 /// (implementation of \link colvarbias \endlink)
-class colvarbias_restraint : public colvarbias {
+class colvarbias_restraint
+  : public virtual colvarbias
+{
 
 public:
 
   /// Retrieve colvar values and calculate their biasing forces
   virtual int update();
 
-  // TODO the following can be supplanted by a new call to init()
   /// Load new configuration - force constant and/or centers only
-  virtual void change_configuration(std::string const &conf);
+  virtual int change_configuration(std::string const &conf) { return COLVARS_NOT_IMPLEMENTED; }
 
   /// Calculate change in energy from using alternate configuration
-  virtual cvm::real energy_difference(std::string const &conf);
+  virtual cvm::real energy_difference(std::string const &conf) { return 0.0; }
 
-  /// Read the bias configuration from a restart file
+  virtual std::string const get_state_params() const;
-  virtual std::istream & read_restart(std::istream &is);
+  virtual int set_state_params(std::string const &conf);
+  // virtual std::ostream & write_state_data(std::ostream &os);
+  // virtual std::istream & read_state_data(std::istream &os);
+  virtual std::ostream & write_state(std::ostream &os);
+  virtual std::istream & read_state(std::istream &is);
 
-  /// Write the bias configuration to a restart file
-  virtual std::ostream & write_restart(std::ostream &os);
-
-  /// Write a label to the trajectory file (comment line)
   virtual std::ostream & write_traj_label(std::ostream &os);
-
-  /// Output quantities such as the bias energy to the trajectory file
   virtual std::ostream & write_traj(std::ostream &os);
 
   /// \brief Constructor
@@ -42,26 +49,110 @@ public:
 
 protected:
 
-  /// \brief Potential function
+  /// \brief Potential function for the i-th colvar
-  virtual cvm::real restraint_potential(cvm::real k, colvar const *x,
+  virtual cvm::real restraint_potential(size_t i) const = 0;
-                                        colvarvalue const &xcenter) const = 0;
 
-  /// \brief Force function
+  /// \brief Force function for the i-th colvar
-  virtual colvarvalue restraint_force(cvm::real k, colvar const *x,
+  virtual colvarvalue const restraint_force(size_t i) const = 0;
-                                      colvarvalue const &xcenter) const = 0;
 
-  ///\brief Unit scaling
+  /// \brief Derivative of the potential function with respect to the force constant
-  virtual cvm::real restraint_convert_k(cvm::real k, cvm::real dist_measure) const = 0;
+  virtual cvm::real d_restraint_potential_dk(size_t i) const = 0;
+};
+
+
+/// Definition and parsing of the restraint centers
+class colvarbias_restraint_centers
+  : public virtual colvarbias_restraint
+{
+public:
+
+  colvarbias_restraint_centers(char const *key);
+  virtual int init(std::string const &conf);
+  virtual int change_configuration(std::string const &conf);
+
+protected:
 
   /// \brief Restraint centers
   std::vector<colvarvalue> colvar_centers;
 
-  /// \brief Restraint centers without wrapping or constraints applied
+  /// \brief Restraint centers outside the domain of the colvars (no wrapping or constraints applied)
   std::vector<colvarvalue> colvar_centers_raw;
+};
+
+
+/// Definition and parsing of the force constant
+class colvarbias_restraint_k
+  : public virtual colvarbias_restraint
+{
+public:
+
+  colvarbias_restraint_k(char const *key);
+  virtual int init(std::string const &conf);
+  virtual int change_configuration(std::string const &conf);
+
+protected:
+  /// \brief Restraint force constant
+  cvm::real force_k;
+};
+
+
+/// Options to change the restraint configuration over time (shared between centers and k moving)
+class colvarbias_restraint_moving
+  : public virtual colvarparse {
+public:
+
+  colvarbias_restraint_moving(char const *key);
+  // Note: despite the diamond inheritance, most of this function gets only executed once
+  virtual int init(std::string const &conf);
+  virtual int update() { return COLVARS_OK; }
+  virtual int change_configuration(std::string const &conf) { return COLVARS_NOT_IMPLEMENTED; }
+
+  virtual std::string const get_state_params() const;
+  virtual int set_state_params(std::string const &conf);
+
+protected:
 
   /// \brief Moving target?
   bool b_chg_centers;
 
+  /// \brief Changing force constant?
+  bool b_chg_force_k;
+
+  /// \brief Number of stages over which to perform the change
+  /// If zero, perform a continuous change
+  int target_nstages;
+
+  /// \brief Number of current stage of the perturbation
+  int stage;
+
+    /// \brief Lambda-schedule for custom varying force constant
+  std::vector<cvm::real> lambda_schedule;
+
+  /// \brief Number of steps required to reach the target force constant
+  /// or restraint centers
+  long target_nsteps;
+};
+
+
+/// Options to change the restraint centers over time
+class colvarbias_restraint_centers_moving
+  : public virtual colvarbias_restraint_centers,
+    public virtual colvarbias_restraint_moving
+{
+public:
+
+  colvarbias_restraint_centers_moving(char const *key);
+  virtual int init(std::string const &conf);
+  virtual int update();
+  virtual int change_configuration(std::string const &conf) { return COLVARS_NOT_IMPLEMENTED; }
+
+  virtual std::string const get_state_params() const;
+  virtual int set_state_params(std::string const &conf);
+  virtual std::ostream & write_traj_label(std::ostream &os);
+  virtual std::ostream & write_traj(std::ostream &os);
+
+protected:
+
   /// \brief New restraint centers
   std::vector<colvarvalue> target_centers;
 
@@ -78,11 +169,29 @@ protected:
   /// \brief Accumulated work
   cvm::real acc_work;
 
-  /// \brief Restraint force constant
+  /// Update the accumulated work
-  cvm::real force_k;
+  int update_acc_work();
+};
 
-  /// \brief Changing force constant?
+
-  bool b_chg_force_k;
+/// Options to change the restraint force constant over time
+class colvarbias_restraint_k_moving
+  : public virtual colvarbias_restraint_k,
+    public virtual colvarbias_restraint_moving
+{
+public:
+
+  colvarbias_restraint_k_moving(char const *key);
+  virtual int init(std::string const &conf);
+  virtual int update();
+  virtual int change_configuration(std::string const &conf) { return COLVARS_NOT_IMPLEMENTED; }
+
+  virtual std::string const get_state_params() const;
+  virtual int set_state_params(std::string const &conf);
+  virtual std::ostream & write_traj_label(std::ostream &os);
+  virtual std::ostream & write_traj(std::ostream &os);
+
+protected:
 
   /// \brief Restraint force constant (target value)
   cvm::real target_force_k;
@@ -90,9 +199,6 @@ protected:
   /// \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 Exponent for varying the force constant
   cvm::real force_k_exp;
 
@@ -100,71 +206,91 @@ protected:
   /// (in TI, would be the accumulating FE derivative)
   cvm::real restraint_FE;
 
-
   /// \brief Equilibration steps for restraint FE calculation through TI
   cvm::real target_equil_steps;
-
-  /// \brief Number of stages over which to perform the change
-  /// If zero, perform a continuous change
-  int target_nstages;
-
-  /// \brief Number of current stage of the perturbation
-  int stage;
-
-  /// \brief Number of steps required to reach the target force constant
-  /// or restraint centers
-  long target_nsteps;
 };
 
 
 /// \brief Harmonic bias restraint
 /// (implementation of \link colvarbias_restraint \endlink)
-class colvarbias_restraint_harmonic : public colvarbias_restraint {
+class colvarbias_restraint_harmonic
-
+  : public colvarbias_restraint_centers_moving,
+    public colvarbias_restraint_k_moving
+{
 public:
-
   colvarbias_restraint_harmonic(char const *key);
   virtual int init(std::string const &conf);
-  // no additional members, destructor not needed
+  virtual int update();
+  virtual std::string const get_state_params() const;
+  virtual int set_state_params(std::string const &conf);
+  virtual std::ostream & write_traj_label(std::ostream &os);
+  virtual std::ostream & write_traj(std::ostream &os);
+  virtual int change_configuration(std::string const &conf);
+  virtual cvm::real energy_difference(std::string const &conf);
 
 protected:
 
-  /// \brief Potential function
+  virtual cvm::real restraint_potential(size_t i) const;
-  virtual cvm::real restraint_potential(cvm::real k, colvar const *x,
+  virtual colvarvalue const restraint_force(size_t i) const;
-                                        colvarvalue const &xcenter) const;
+  virtual cvm::real d_restraint_potential_dk(size_t i) const;
+};
 
-  /// \brief Force function
-  virtual colvarvalue restraint_force(cvm::real k, colvar const *x,
-                                      colvarvalue const &xcenter) const;
 
-  ///\brief Unit scaling
+/// \brief Wall restraint
-  virtual cvm::real restraint_convert_k(cvm::real k, cvm::real dist_measure) const;
+/// (implementation of \link colvarbias_restraint \endlink)
+class colvarbias_restraint_harmonic_walls
+  : public colvarbias_restraint_k_moving
+{
+public:
 
+  colvarbias_restraint_harmonic_walls(char const *key);
+  virtual int init(std::string const &conf);
+  virtual int update();
+  virtual void communicate_forces();
+  virtual std::string const get_state_params() const;
+  virtual int set_state_params(std::string const &conf);
+  virtual std::ostream & write_traj_label(std::ostream &os);
+  virtual std::ostream & write_traj(std::ostream &os);
+
+protected:
+
+  /// \brief Location of the lower walls
+  std::vector<colvarvalue> lower_walls;
+
+  /// \brief Location of the upper walls
+  std::vector<colvarvalue> upper_walls;
+
+  virtual cvm::real colvar_distance(size_t i) const;
+  virtual cvm::real restraint_potential(size_t i) const;
+  virtual colvarvalue const restraint_force(size_t i) const;
+  virtual cvm::real d_restraint_potential_dk(size_t i) const;
 };
 
 
 /// \brief Linear bias restraint
 /// (implementation of \link colvarbias_restraint \endlink)
-class colvarbias_restraint_linear : public colvarbias_restraint {
+class colvarbias_restraint_linear
+  : public colvarbias_restraint_centers_moving,
+    public colvarbias_restraint_k_moving
+{
 
 public:
   colvarbias_restraint_linear(char const *key);
   virtual int init(std::string const &conf);
-  // no additional members, destructor not needed
+  virtual int update();
+  virtual int change_configuration(std::string const &conf);
+  virtual cvm::real energy_difference(std::string const &conf);
+
+  virtual std::string const get_state_params() const;
+  virtual int set_state_params(std::string const &conf);
+  virtual std::ostream & write_traj_label(std::ostream &os);
+  virtual std::ostream & write_traj(std::ostream &os);
 
 protected:
 
-  /// \brief Potential function
+  virtual cvm::real restraint_potential(size_t i) const;
-  virtual cvm::real restraint_potential(cvm::real k, colvar const *x,
+  virtual colvarvalue const restraint_force(size_t i) const;
-                                        colvarvalue const &xcenter) const;
+  virtual cvm::real d_restraint_potential_dk(size_t i) const;
-
-  /// \brief Force function
-  virtual colvarvalue restraint_force(cvm::real k, colvar const *x,
-                                      colvarvalue const &xcenter) const;
-
-  ///\brief Unit scaling
-  virtual cvm::real restraint_convert_k(cvm::real k, cvm::real dist_measure) const;
-
 };
 
 

lib/colvars/colvarcomp.cpp

@@ -1,5 +1,12 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #include "colvarmodule.h"
 #include "colvarvalue.h"
 #include "colvar.h"
@@ -154,15 +161,17 @@ void colvar::cvc::read_data()
 
 void colvar::cvc::calc_force_invgrads()
 {
-  cvm::fatal_error("Error: calculation of inverse gradients is not implemented "
+  cvm::error("Error: calculation of inverse gradients is not implemented "
-                    "for colvar components of type \""+function_type+"\".\n");
+             "for colvar components of type \""+function_type+"\".\n",
+             COLVARS_NOT_IMPLEMENTED);
 }
 
 
 void colvar::cvc::calc_Jacobian_derivative()
 {
-  cvm::fatal_error("Error: calculation of inverse gradients is not implemented "
+  cvm::error("Error: calculation of inverse gradients is not implemented "
-                    "for colvar components of type \""+function_type+"\".\n");
+             "for colvar components of type \""+function_type+"\".\n",
+             COLVARS_NOT_IMPLEMENTED);
 }
 
 
@@ -281,6 +290,33 @@ void colvar::cvc::debug_gradients(cvm::atom_group *group)
 }
 
 
+cvm::real colvar::cvc::dist2(colvarvalue const &x1,
+                             colvarvalue const &x2) const
+{
+  return x1.dist2(x2);
+}
+
+
+colvarvalue colvar::cvc::dist2_lgrad(colvarvalue const &x1,
+                                     colvarvalue const &x2) const
+{
+  return x1.dist2_grad(x2);
+}
+
+
+colvarvalue colvar::cvc::dist2_rgrad(colvarvalue const &x1,
+                                     colvarvalue const &x2) const
+{
+  return x2.dist2_grad(x1);
+}
+
+
+void colvar::cvc::wrap(colvarvalue &x) const
+{
+  return;
+}
+
+
 // Static members
 
 std::vector<colvardeps::feature *> colvar::cvc::cvc_features;

lib/colvars/colvarcomp.h

@@ -1,5 +1,12 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #ifndef COLVARCOMP_H
 #define COLVARCOMP_H
 
@@ -105,8 +112,8 @@ public:
   /// options from the provided configuration string
   /// Returns reference to new group
   cvm::atom_group *parse_group(std::string const &conf,
-                   char const *group_key,
+                               char const *group_key,
-                   bool optional = false);
+                               bool optional = false);
 
   /// \brief Parse options pertaining to total force calculation
   virtual int init_total_force_params(std::string const &conf);
@@ -130,7 +137,7 @@ public:
   }
 
   /// \brief Obtain data needed for the calculation for the backend
-  void read_data();
+  virtual void read_data();
 
   /// \brief Calculate the variable
   virtual void calc_value() = 0;
@@ -151,17 +158,14 @@ public:
 
 
   /// \brief Return the previously calculated value
-  virtual colvarvalue const & value() const;
+  colvarvalue const & value() const;
-
-  // /// \brief Return const pointer to the previously calculated value
-  // virtual const colvarvalue *p_value() const;
 
   /// \brief Return the previously calculated total force
-  virtual colvarvalue const & total_force() const;
+  colvarvalue const & total_force() const;
 
   /// \brief Return the previously calculated divergence of the
   /// inverse atomic gradients
-  virtual colvarvalue const & Jacobian_derivative() const;
+  colvarvalue const & Jacobian_derivative() const;
 
   /// \brief Apply the collective variable force, by communicating the
   /// atomic forces to the simulation program (\b Note: the \link ft
@@ -247,52 +251,24 @@ protected:
 };
 
 
-
-
 inline colvarvalue const & colvar::cvc::value() const
 {
   return x;
 }
 
-// inline const colvarvalue * colvar::cvc::p_value() const
-// {
-//   return &x;
-// }
 
 inline colvarvalue const & colvar::cvc::total_force() const
 {
   return ft;
 }
 
+
 inline colvarvalue const & colvar::cvc::Jacobian_derivative() const
 {
   return jd;
 }
 
 
-inline cvm::real colvar::cvc::dist2(colvarvalue const &x1,
-                                    colvarvalue const &x2) const
-{
-  return x1.dist2(x2);
-}
-
-inline colvarvalue colvar::cvc::dist2_lgrad(colvarvalue const &x1,
-                                            colvarvalue const &x2) const
-{
-  return x1.dist2_grad(x2);
-}
-
-inline colvarvalue colvar::cvc::dist2_rgrad(colvarvalue const &x1,
-                                            colvarvalue const &x2) const
-{
-  return x2.dist2_grad(x1);
-}
-
-inline void colvar::cvc::wrap(colvarvalue &x) const
-{
-  return;
-}
-
 
 /// \brief Colvar component: distance between the centers of mass of
 /// two groups (colvarvalue::type_scalar type, range [0:*))
@@ -312,7 +288,7 @@ protected:
 public:
   distance(std::string const &conf);
   distance();
-  virtual inline ~distance() {}
+  virtual ~distance() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void calc_force_invgrads();
@@ -327,6 +303,7 @@ public:
 };
 
 
+
 // \brief Colvar component: distance vector between centers of mass
 // of two groups (\link colvarvalue::type_3vector \endlink type,
 // range (-*:*)x(-*:*)x(-*:*))
@@ -336,7 +313,7 @@ class colvar::distance_vec
 public:
   distance_vec(std::string const &conf);
   distance_vec();
-  virtual inline ~distance_vec() {}
+  virtual ~distance_vec() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force(colvarvalue const &force);
@@ -352,6 +329,7 @@ public:
 };
 
 
+
 /// \brief Colvar component: distance unit vector (direction) between
 /// centers of mass of two groups (colvarvalue::type_unit3vector type,
 /// range [-1:1]x[-1:1]x[-1:1])
@@ -361,19 +339,14 @@ class colvar::distance_dir
 public:
   distance_dir(std::string const &conf);
   distance_dir();
-  virtual inline ~distance_dir() {}
+  virtual ~distance_dir() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force(colvarvalue const &force);
-  virtual cvm::real dist2(colvarvalue const &x1,
-                          colvarvalue const &x2) const;
-  virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
-                                  colvarvalue const &x2) const;
-  virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
-                                  colvarvalue const &x2) const;
 };
 
 
+
 /// \brief Colvar component: projection of the distance vector along
 /// an axis(colvarvalue::type_scalar type, range (-*:*))
 class colvar::distance_z
@@ -399,7 +372,7 @@ protected:
 public:
   distance_z(std::string const &conf);
   distance_z();
-  virtual inline ~distance_z() {}
+  virtual ~distance_z() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void calc_force_invgrads();
@@ -416,6 +389,7 @@ public:
 };
 
 
+
 /// \brief Colvar component: projection of the distance vector on a
 /// plane (colvarvalue::type_scalar type, range [0:*))
 class colvar::distance_xy
@@ -429,7 +403,7 @@ protected:
 public:
   distance_xy(std::string const &conf);
   distance_xy();
-  virtual inline ~distance_xy() {}
+  virtual ~distance_xy() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void calc_force_invgrads();
@@ -444,6 +418,7 @@ public:
 };
 
 
+
 /// \brief Colvar component: average distance between two groups of atoms, weighted as the sixth power,
 /// as in NMR refinements(colvarvalue::type_scalar type, range (0:*))
 class colvar::distance_inv
@@ -455,7 +430,7 @@ protected:
 public:
   distance_inv(std::string const &conf);
   distance_inv();
-  virtual inline ~distance_inv() {}
+  virtual ~distance_inv() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force(colvarvalue const &force);
@@ -468,6 +443,7 @@ public:
 };
 
 
+
 /// \brief Colvar component: N1xN2 vector of pairwise distances
 /// (colvarvalue::type_vector type, range (0:*) for each component)
 class colvar::distance_pairs
@@ -483,16 +459,10 @@ protected:
 public:
   distance_pairs(std::string const &conf);
   distance_pairs();
-  virtual inline ~distance_pairs() {}
+  virtual ~distance_pairs() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force(colvarvalue const &force);
-  virtual cvm::real dist2(colvarvalue const &x1,
-                          colvarvalue const &x2) const;
-  virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
-                                  colvarvalue const &x2) const;
-  virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
-                                  colvarvalue const &x2) const;
 };
 
 
@@ -509,7 +479,7 @@ public:
   /// Constructor
   gyration(std::string const &conf);
   gyration();
-  virtual inline ~gyration() {}
+  virtual ~gyration() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void calc_force_invgrads();
@@ -524,6 +494,7 @@ public:
 };
 
 
+
 /// \brief Colvar component: moment of inertia of an atom group
 /// (colvarvalue::type_scalar type, range [0:*))
 class colvar::inertia
@@ -533,7 +504,7 @@ public:
   /// Constructor
   inertia(std::string const &conf);
   inertia();
-  virtual inline ~inertia() {}
+  virtual ~inertia() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force(colvarvalue const &force);
@@ -546,6 +517,7 @@ public:
 };
 
 
+
 /// \brief Colvar component: moment of inertia of an atom group
 /// around a user-defined axis (colvarvalue::type_scalar type, range [0:*))
 class colvar::inertia_z
@@ -558,7 +530,7 @@ public:
   /// Constructor
   inertia_z(std::string const &conf);
   inertia_z();
-  virtual inline ~inertia_z() {}
+  virtual ~inertia_z() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force(colvarvalue const &force);
@@ -571,6 +543,7 @@ public:
 };
 
 
+
 /// \brief Colvar component: projection of 3N coordinates onto an
 /// eigenvector(colvarvalue::type_scalar type, range (-*:*))
 class colvar::eigenvector
@@ -597,7 +570,7 @@ public:
 
   /// Constructor
   eigenvector(std::string const &conf);
-  virtual inline ~eigenvector() {}
+  virtual ~eigenvector() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void calc_force_invgrads();
@@ -645,7 +618,7 @@ public:
   /// \brief Initialize the three groups after three atoms
   angle(cvm::atom const &a1, cvm::atom const &a2, cvm::atom const &a3);
   angle();
-  virtual inline ~angle() {}
+  virtual ~angle() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void calc_force_invgrads();
@@ -659,6 +632,8 @@ public:
                                   colvarvalue const &x2) const;
 };
 
+
+
 /// \brief Colvar component: angle between the dipole of a molecule and an axis
 /// formed by two groups of atoms(colvarvalue::type_scalar type, range [0:PI])
 class colvar::dipole_angle
@@ -691,7 +666,7 @@ public:
   /// \brief Initialize the three groups after three atoms
   dipole_angle (cvm::atom const &a1, cvm::atom const &a2, cvm::atom const &a3);
   dipole_angle();
-  virtual inline ~dipole_angle() {}
+  virtual ~dipole_angle() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force (colvarvalue const &force);
@@ -703,6 +678,8 @@ public:
                                    colvarvalue const &x2) const;
 };
 
+
+
 /// \brief Colvar component: dihedral between the centers of mass of
 /// four groups (colvarvalue::type_scalar type, range [-PI:PI])
 class colvar::dihedral
@@ -732,7 +709,7 @@ public:
   /// \brief Initialize the four groups after four atoms
   dihedral(cvm::atom const &a1, cvm::atom const &a2, cvm::atom const &a3, cvm::atom const &a4);
   dihedral();
-  virtual inline ~dihedral() {}
+  virtual ~dihedral() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void calc_force_invgrads();
@@ -753,6 +730,7 @@ public:
 };
 
 
+
 /// \brief Colvar component: coordination number between two groups
 /// (colvarvalue::type_scalar type, range [0:N1*N2])
 class colvar::coordnum
@@ -781,7 +759,7 @@ public:
   /// Constructor
   coordnum(std::string const &conf);
   coordnum();
-  virtual inline ~coordnum() {}
+  virtual ~coordnum() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force(colvarvalue const &force);
@@ -812,6 +790,8 @@ public:
                                   colvarvalue const &x2) const;
 };
 
+
+
 /// \brief Colvar component: self-coordination number within a group
 /// (colvarvalue::type_scalar type, range [0:N*(N-1)/2])
 class colvar::selfcoordnum
@@ -830,7 +810,7 @@ public:
   /// Constructor
   selfcoordnum(std::string const &conf);
   selfcoordnum();
-  virtual inline ~selfcoordnum() {}
+  virtual ~selfcoordnum() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force(colvarvalue const &force);
@@ -852,6 +832,7 @@ public:
 };
 
 
+
 /// \brief Colvar component: coordination number between two groups
 /// (colvarvalue::type_scalar type, range [0:N1*N2])
 class colvar::groupcoordnum
@@ -873,7 +854,7 @@ public:
   /// Constructor
   groupcoordnum(std::string const &conf);
   groupcoordnum();
-  virtual inline ~groupcoordnum() {}
+  virtual ~groupcoordnum() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force(colvarvalue const &force);
@@ -896,6 +877,7 @@ public:
   static cvm::real switching_function(cvm::rvector const &r0_vec,
                                       int const &exp_num, int const &exp_den,
                                       cvm::atom &A1, cvm::atom &A2);
+  */
 
   virtual cvm::real dist2(colvarvalue const &x1,
                           colvarvalue const &x2) const;
@@ -903,10 +885,10 @@ public:
                                   colvarvalue const &x2) const;
   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
                                   colvarvalue const &x2) const;
-  */
 };
 
 
+
 /// \brief Colvar component: hydrogen bond, defined as the product of
 /// a colvar::coordnum and 1/2*(1-cos((180-ang)/ang_tol))
 /// (colvarvalue::type_scalar type, range [0:1])
@@ -941,6 +923,7 @@ public:
 };
 
 
+
 /// \brief Colvar component: alpha helix content of a contiguous
 /// segment of 5 or more residues, implemented as a sum of phi/psi
 /// dihedral angles and hydrogen bonds (colvarvalue::type_scalar type,
@@ -969,7 +952,7 @@ public:
 
 //   alpha_dihedrals (std::string const &conf);
 //   alpha_dihedrals();
-//   virtual inline ~alpha_dihedrals() {}
+//   virtual ~alpha_dihedrals() {}
 //   virtual void calc_value();
 //   virtual void calc_gradients();
 //   virtual void apply_force (colvarvalue const &force);
@@ -982,6 +965,7 @@ public:
 // };
 
 
+
 /// \brief Colvar component: alpha helix content of a contiguous
 /// segment of 5 or more residues, implemented as a sum of Ca-Ca-Ca
 /// angles and hydrogen bonds (colvarvalue::type_scalar type, range
@@ -1022,6 +1006,8 @@ public:
                                   colvarvalue const &x2) const;
 };
 
+
+
 /// \brief Colvar component: dihedPC
 /// Projection of the config onto a dihedral principal component
 /// See e.g. Altis et al., J. Chem. Phys 126, 244111 (2007)
@@ -1050,6 +1036,8 @@ public:
                                   colvarvalue const &x2) const;
 };
 
+
+
 /// \brief Colvar component: orientation in space of an atom group,
 /// with respect to a set of reference coordinates
 /// (colvarvalue::type_quaternion type, range
@@ -1078,7 +1066,7 @@ public:
 
   orientation(std::string const &conf);
   orientation();
-  virtual inline ~orientation() {}
+  virtual ~orientation() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force(colvarvalue const &force);
@@ -1091,6 +1079,7 @@ public:
 };
 
 
+
 /// \brief Colvar component: angle of rotation with respect to a set
 /// of reference coordinates (colvarvalue::type_scalar type, range
 /// [0:PI))
@@ -1101,7 +1090,7 @@ public:
 
   orientation_angle(std::string const &conf);
   orientation_angle();
-  virtual inline ~orientation_angle() {}
+  virtual ~orientation_angle() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force(colvarvalue const &force);
@@ -1114,6 +1103,7 @@ public:
 };
 
 
+
 /// \brief Colvar component: cosine of the angle of rotation with respect to a set
 /// of reference coordinates (colvarvalue::type_scalar type, range
 /// [-1:1])
@@ -1124,7 +1114,7 @@ public:
 
   orientation_proj(std::string const &conf);
   orientation_proj();
-  virtual inline ~orientation_proj() {}
+  virtual ~orientation_proj() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force(colvarvalue const &force);
@@ -1137,6 +1127,7 @@ public:
 };
 
 
+
 /// \brief Colvar component: projection of the orientation vector onto
 /// a predefined axis (colvarvalue::type_scalar type, range [-1:1])
 class colvar::tilt
@@ -1150,7 +1141,7 @@ public:
 
   tilt(std::string const &conf);
   tilt();
-  virtual inline ~tilt() {}
+  virtual ~tilt() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force(colvarvalue const &force);
@@ -1177,7 +1168,7 @@ public:
 
   spin_angle(std::string const &conf);
   spin_angle();
-  virtual inline ~spin_angle() {}
+  virtual ~spin_angle() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force(colvarvalue const &force);
@@ -1215,7 +1206,7 @@ public:
 
   /// Constructor
   rmsd(std::string const &conf);
-  virtual inline ~rmsd() {}
+  virtual ~rmsd() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void calc_force_invgrads();
@@ -1230,6 +1221,7 @@ public:
 };
 
 
+
 // \brief Colvar component: flat vector of Cartesian coordinates
 // Mostly useful to compute scripted colvar values
 class colvar::cartesian
@@ -1243,7 +1235,7 @@ protected:
 public:
   cartesian(std::string const &conf);
   cartesian();
-  virtual inline ~cartesian() {}
+  virtual ~cartesian() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force(colvarvalue const &force);
@@ -1260,255 +1252,26 @@ public:
 
 #define simple_scalar_dist_functions(TYPE)                              \
                                                                         \
-  inline cvm::real colvar::TYPE::dist2(colvarvalue const &x1,           \
+                                                                        \
-                                       colvarvalue const &x2) const     \
+  cvm::real colvar::TYPE::dist2(colvarvalue const &x1,                  \
+                                colvarvalue const &x2) const            \
   {                                                                     \
     return (x1.real_value - x2.real_value)*(x1.real_value - x2.real_value); \
   }                                                                     \
                                                                         \
-  inline colvarvalue colvar::TYPE::dist2_lgrad(colvarvalue const &x1,   \
+                                                                        \
-                                               colvarvalue const &x2) const \
+  colvarvalue colvar::TYPE::dist2_lgrad(colvarvalue const &x1,          \
+                                        colvarvalue const &x2) const    \
   {                                                                     \
     return 2.0 * (x1.real_value - x2.real_value);                       \
   }                                                                     \
                                                                         \
-  inline colvarvalue colvar::TYPE::dist2_rgrad(colvarvalue const &x1,   \
+                                                                        \
-                                               colvarvalue const &x2) const \
+  colvarvalue colvar::TYPE::dist2_rgrad(colvarvalue const &x1,          \
+                                        colvarvalue const &x2) const    \
   {                                                                     \
     return this->dist2_lgrad(x2, x1);                                   \
   }                                                                     \
                                                                         \
 
-simple_scalar_dist_functions(distance)
-// NOTE: distance_z has explicit functions, see below
-simple_scalar_dist_functions(distance_xy)
-simple_scalar_dist_functions(distance_inv)
-simple_scalar_dist_functions(angle)
-simple_scalar_dist_functions(dipole_angle)
-simple_scalar_dist_functions(coordnum)
-simple_scalar_dist_functions(selfcoordnum)
-simple_scalar_dist_functions(h_bond)
-simple_scalar_dist_functions(gyration)
-simple_scalar_dist_functions(inertia)
-simple_scalar_dist_functions(inertia_z)
-simple_scalar_dist_functions(rmsd)
-simple_scalar_dist_functions(orientation_angle)
-simple_scalar_dist_functions(orientation_proj)
-simple_scalar_dist_functions(tilt)
-simple_scalar_dist_functions(eigenvector)
-//  simple_scalar_dist_functions (alpha_dihedrals)
-simple_scalar_dist_functions(alpha_angles)
-simple_scalar_dist_functions(dihedPC)
-
-
-// metrics functions for cvc implementations with a periodicity
-
-inline cvm::real colvar::dihedral::dist2(colvarvalue const &x1,
-                                           colvarvalue const &x2) const
-{
-cvm::real diff = x1.real_value - x2.real_value;
-diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
-return diff * diff;
-}
-
-inline colvarvalue colvar::dihedral::dist2_lgrad(colvarvalue const &x1,
-                                                   colvarvalue const &x2) const
-{
-cvm::real diff = x1.real_value - x2.real_value;
-diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
-return 2.0 * diff;
-}
-
-inline colvarvalue colvar::dihedral::dist2_rgrad(colvarvalue const &x1,
-                                                   colvarvalue const &x2) const
-{
-cvm::real diff = x1.real_value - x2.real_value;
-diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
-return (-2.0) * diff;
-}
-
-inline void colvar::dihedral::wrap(colvarvalue &x) const
-{
-if ((x.real_value - wrap_center) >= 180.0) {
-x.real_value -= 360.0;
-return;
-}
-
-if ((x.real_value - wrap_center) < -180.0) {
-x.real_value += 360.0;
-return;
-}
-
-return;
-}
-
-inline cvm::real colvar::spin_angle::dist2(colvarvalue const &x1,
-                                             colvarvalue const &x2) const
-{
-cvm::real diff = x1.real_value - x2.real_value;
-diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
-return diff * diff;
-}
-
-inline colvarvalue colvar::spin_angle::dist2_lgrad(colvarvalue const &x1,
-                                                     colvarvalue const &x2) const
-{
-cvm::real diff = x1.real_value - x2.real_value;
-diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
-return 2.0 * diff;
-}
-
-inline colvarvalue colvar::spin_angle::dist2_rgrad(colvarvalue const &x1,
-                                                     colvarvalue const &x2) const
-{
-cvm::real diff = x1.real_value - x2.real_value;
-diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
-return (-2.0) * diff;
-}
-
-inline void colvar::spin_angle::wrap(colvarvalue &x) const
-{
-if ((x.real_value - wrap_center) >= 180.0) {
-x.real_value -= 360.0;
-return;
-}
-
-if ((x.real_value - wrap_center) < -180.0) {
-x.real_value += 360.0;
-return;
-}
-
-return;
-}
-
-
-// Projected distance
-// Differences should always be wrapped around 0 (ignoring wrap_center)
-inline cvm::real colvar::distance_z::dist2(colvarvalue const &x1,
-                                             colvarvalue const &x2) const
-{
-  cvm::real diff = x1.real_value - x2.real_value;
-  if (period != 0.0) {
-    cvm::real shift = std::floor(diff/period + 0.5);
-    diff -= shift * period;
-  }
-  return diff * diff;
-}
-
-inline colvarvalue colvar::distance_z::dist2_lgrad(colvarvalue const &x1,
-                                                   colvarvalue const &x2) const
-{
-  cvm::real diff = x1.real_value - x2.real_value;
-  if (period != 0.0) {
-    cvm::real shift = std::floor(diff/period + 0.5);
-    diff -= shift * period;
-  }
-  return 2.0 * diff;
-}
-
-inline colvarvalue colvar::distance_z::dist2_rgrad(colvarvalue const &x1,
-                                                   colvarvalue const &x2) const
-{
-  cvm::real diff = x1.real_value - x2.real_value;
-  if (period != 0.0) {
-    cvm::real shift = std::floor(diff/period + 0.5);
-    diff -= shift * period;
-  }
-  return (-2.0) * diff;
-}
-
-inline void colvar::distance_z::wrap(colvarvalue &x) const
-{
-  if (! this->b_periodic) {
-    // don't wrap if the period has not been set
-    return;
-  }
-
-  cvm::real shift = std::floor((x.real_value - wrap_center) / period + 0.5);
-  x.real_value -= shift * period;
-  return;
-}
-
-
-// distance between three dimensional vectors
-//
-// TODO apply PBC to distance_vec
-// Note: differences should be centered around (0, 0, 0)!
-
-inline cvm::real colvar::distance_vec::dist2(colvarvalue const &x1,
-                                             colvarvalue const &x2) const
-{
-  return cvm::position_dist2(x1.rvector_value, x2.rvector_value);
-}
-
-inline colvarvalue colvar::distance_vec::dist2_lgrad(colvarvalue const &x1,
-                                                     colvarvalue const &x2) const
-{
-  return 2.0 * cvm::position_distance(x2.rvector_value, x1.rvector_value);
-}
-
-inline colvarvalue colvar::distance_vec::dist2_rgrad(colvarvalue const &x1,
-                                                     colvarvalue const &x2) const
-{
-  return 2.0 * cvm::position_distance(x2.rvector_value, x1.rvector_value);
-}
-
-inline cvm::real colvar::distance_dir::dist2(colvarvalue const &x1,
-                                             colvarvalue const &x2) const
-{
-  return (x1.rvector_value - x2.rvector_value).norm2();
-}
-
-inline colvarvalue colvar::distance_dir::dist2_lgrad(colvarvalue const &x1,
-                                                     colvarvalue const &x2) const
-{
-  return colvarvalue((x1.rvector_value - x2.rvector_value), colvarvalue::type_unit3vector);
-}
-
-inline colvarvalue colvar::distance_dir::dist2_rgrad(colvarvalue const &x1,
-                                                     colvarvalue const &x2) const
-{
-  return colvarvalue((x2.rvector_value - x1.rvector_value), colvarvalue::type_unit3vector);
-}
-
-inline cvm::real colvar::distance_pairs::dist2(colvarvalue const &x1,
-                                               colvarvalue const &x2) const
-{
-  return (x1.vector1d_value - x2.vector1d_value).norm2();
-}
-
-inline colvarvalue colvar::distance_pairs::dist2_lgrad(colvarvalue const &x1,
-                                                       colvarvalue const &x2) const
-{
-  return colvarvalue((x1.vector1d_value - x2.vector1d_value), colvarvalue::type_vector);
-}
-
-inline colvarvalue colvar::distance_pairs::dist2_rgrad(colvarvalue const &x1,
-                                                       colvarvalue const &x2) const
-{
-  return colvarvalue((x2.vector1d_value - x1.vector1d_value), colvarvalue::type_vector);
-}
-
-
-// distance between quaternions
-
-inline cvm::real colvar::orientation::dist2(colvarvalue const &x1,
-                                            colvarvalue const &x2) const
-{
-  return x1.quaternion_value.dist2(x2);
-}
-
-inline colvarvalue colvar::orientation::dist2_lgrad(colvarvalue const &x1,
-                                                    colvarvalue const &x2) const
-{
-  return x1.quaternion_value.dist2_grad(x2);
-}
-
-inline colvarvalue colvar::orientation::dist2_rgrad(colvarvalue const &x1,
-                                                    colvarvalue const &x2) const
-{
-  return x2.quaternion_value.dist2_grad(x1);
-}
-
-
 #endif

lib/colvars/colvarcomp_angles.cpp

@@ -1,5 +1,12 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #include "colvarmodule.h"
 #include "colvar.h"
 #include "colvarcomp.h"
@@ -7,6 +14,7 @@
 #include <cmath>
 
 
+
 colvar::angle::angle(std::string const &conf)
   : cvc(conf)
 {
@@ -85,6 +93,7 @@ void colvar::angle::calc_gradients()
   group3->set_weighted_gradient(dxdr3);
 }
 
+
 void colvar::angle::calc_force_invgrads()
 {
   // This uses a force measurement on groups 1 and 3 only
@@ -107,6 +116,7 @@ void colvar::angle::calc_force_invgrads()
   return;
 }
 
+
 void colvar::angle::calc_Jacobian_derivative()
 {
   // det(J) = (2 pi) r^2 * sin(theta)
@@ -129,6 +139,8 @@ void colvar::angle::apply_force(colvarvalue const &force)
 }
 
 
+simple_scalar_dist_functions(angle)
+
 
 
 colvar::dipole_angle::dipole_angle(std::string const &conf)
@@ -235,6 +247,8 @@ void colvar::dipole_angle::apply_force(colvarvalue const &force)
 }
 
 
+simple_scalar_dist_functions(dipole_angle)
+
 
 
 colvar::dihedral::dihedral(std::string const &conf)
@@ -453,3 +467,46 @@ void colvar::dihedral::apply_force(colvarvalue const &force)
 }
 
 
+// metrics functions for cvc implementations with a periodicity
+
+cvm::real colvar::dihedral::dist2(colvarvalue const &x1,
+                                  colvarvalue const &x2) const
+{
+  cvm::real diff = x1.real_value - x2.real_value;
+  diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
+  return diff * diff;
+}
+
+
+colvarvalue colvar::dihedral::dist2_lgrad(colvarvalue const &x1,
+                                          colvarvalue const &x2) const
+{
+  cvm::real diff = x1.real_value - x2.real_value;
+  diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
+  return 2.0 * diff;
+}
+
+
+colvarvalue colvar::dihedral::dist2_rgrad(colvarvalue const &x1,
+                                          colvarvalue const &x2) const
+{
+  cvm::real diff = x1.real_value - x2.real_value;
+  diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
+  return (-2.0) * diff;
+}
+
+
+void colvar::dihedral::wrap(colvarvalue &x) const
+{
+  if ((x.real_value - wrap_center) >= 180.0) {
+    x.real_value -= 360.0;
+    return;
+  }
+
+  if ((x.real_value - wrap_center) < -180.0) {
+    x.real_value += 360.0;
+    return;
+  }
+
+  return;
+}

lib/colvars/colvarcomp_coordnums.cpp

@@ -1,5 +1,12 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #include <cmath>
 
 #include "colvarmodule.h"
@@ -13,10 +20,10 @@
 
 template<bool calculate_gradients>
 cvm::real colvar::coordnum::switching_function(cvm::real const &r0,
-                                                int const &en,
+                                               int const &en,
-                                                int const &ed,
+                                               int const &ed,
-                                                cvm::atom &A1,
+                                               cvm::atom &A1,
-                                                cvm::atom &A2)
+                                               cvm::atom &A2)
 {
   cvm::rvector const diff = cvm::position_distance(A1.pos, A2.pos);
   cvm::real const l2 = diff.norm2()/(r0*r0);
@@ -42,10 +49,10 @@ cvm::real colvar::coordnum::switching_function(cvm::real const &r0,
 
 template<bool calculate_gradients>
 cvm::real colvar::coordnum::switching_function(cvm::rvector const &r0_vec,
-                                                int const &en,
+                                               int const &en,
-                                                int const &ed,
+                                               int const &ed,
-                                                cvm::atom &A1,
+                                               cvm::atom &A1,
-                                                cvm::atom &A2)
+                                               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);
@@ -190,6 +197,7 @@ void colvar::coordnum::calc_gradients()
   }
 }
 
+
 void colvar::coordnum::apply_force(colvarvalue const &force)
 {
   if (!group1->noforce)
@@ -200,6 +208,9 @@ void colvar::coordnum::apply_force(colvarvalue const &force)
 }
 
 
+simple_scalar_dist_functions(coordnum)
+
+
 
 // h_bond member functions
 
@@ -252,6 +263,7 @@ colvar::h_bond::h_bond(cvm::atom const &acceptor,
   atom_groups[0]->add_atom(donor);
 }
 
+
 colvar::h_bond::h_bond()
   : cvc()
 {
@@ -284,6 +296,8 @@ void colvar::h_bond::apply_force(colvarvalue const &force)
 }
 
 
+simple_scalar_dist_functions(h_bond)
+
 
 
 colvar::selfcoordnum::selfcoordnum(std::string const &conf)
@@ -339,6 +353,9 @@ void colvar::selfcoordnum::apply_force(colvarvalue const &force)
 }
 
 
+simple_scalar_dist_functions(selfcoordnum)
+
+
 // groupcoordnum member functions
 colvar::groupcoordnum::groupcoordnum(std::string const &conf)
   : distance(conf), b_anisotropic(false)
@@ -448,7 +465,6 @@ cvm::real colvar::groupcoordnum::switching_function(cvm::rvector const &r0_vec,
 #endif
 
 
-
 void colvar::groupcoordnum::calc_value()
 {
 
@@ -460,7 +476,6 @@ void colvar::groupcoordnum::calc_value()
 
   x.real_value = coordnum::switching_function<false>(r0, en, ed,
                                                      group1_com_atom, group2_com_atom);
-
 }
 
 
@@ -486,3 +501,6 @@ void colvar::groupcoordnum::apply_force(colvarvalue const &force)
   if (!group2->noforce)
     group2->apply_colvar_force(force.real_value);
 }
+
+
+simple_scalar_dist_functions(groupcoordnum)

lib/colvars/colvarcomp_distances.cpp

@@ -1,5 +1,12 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #include <cmath>
 
 #include "colvarmodule.h"
@@ -91,6 +98,9 @@ void colvar::distance::apply_force(colvarvalue const &force)
 }
 
 
+simple_scalar_dist_functions(distance)
+
+
 
 colvar::distance_vec::distance_vec(std::string const &conf)
   : distance(conf)
@@ -138,6 +148,27 @@ void colvar::distance_vec::apply_force(colvarvalue const &force)
 }
 
 
+cvm::real colvar::distance_vec::dist2(colvarvalue const &x1,
+                                      colvarvalue const &x2) const
+{
+  return cvm::position_dist2(x1.rvector_value, x2.rvector_value);
+}
+
+
+colvarvalue colvar::distance_vec::dist2_lgrad(colvarvalue const &x1,
+                                              colvarvalue const &x2) const
+{
+  return 2.0 * cvm::position_distance(x2.rvector_value, x1.rvector_value);
+}
+
+
+colvarvalue colvar::distance_vec::dist2_rgrad(colvarvalue const &x1,
+                                              colvarvalue const &x2) const
+{
+  return 2.0 * cvm::position_distance(x2.rvector_value, x1.rvector_value);
+}
+
+
 
 colvar::distance_z::distance_z(std::string const &conf)
   : cvc(conf)
@@ -191,6 +222,7 @@ colvar::distance_z::distance_z(std::string const &conf)
 
 }
 
+
 colvar::distance_z::distance_z()
 {
   function_type = "distance_z";
@@ -200,6 +232,7 @@ colvar::distance_z::distance_z()
   x.type(colvarvalue::type_scalar);
 }
 
+
 void colvar::distance_z::calc_value()
 {
   if (fixed_axis) {
@@ -227,6 +260,7 @@ void colvar::distance_z::calc_value()
   this->wrap(x);
 }
 
+
 void colvar::distance_z::calc_gradients()
 {
   main->set_weighted_gradient( axis );
@@ -248,6 +282,7 @@ void colvar::distance_z::calc_gradients()
   }
 }
 
+
 void colvar::distance_z::calc_force_invgrads()
 {
   main->read_total_forces();
@@ -260,11 +295,13 @@ void colvar::distance_z::calc_force_invgrads()
   }
 }
 
+
 void colvar::distance_z::calc_Jacobian_derivative()
 {
   jd.real_value = 0.0;
 }
 
+
 void colvar::distance_z::apply_force(colvarvalue const &force)
 {
   if (!ref1->noforce)
@@ -278,6 +315,56 @@ void colvar::distance_z::apply_force(colvarvalue const &force)
 }
 
 
+// Differences should always be wrapped around 0 (ignoring wrap_center)
+cvm::real colvar::distance_z::dist2(colvarvalue const &x1,
+                                    colvarvalue const &x2) const
+{
+  cvm::real diff = x1.real_value - x2.real_value;
+  if (b_periodic) {
+    cvm::real shift = std::floor(diff/period + 0.5);
+    diff -= shift * period;
+  }
+  return diff * diff;
+}
+
+
+colvarvalue colvar::distance_z::dist2_lgrad(colvarvalue const &x1,
+                                            colvarvalue const &x2) const
+{
+  cvm::real diff = x1.real_value - x2.real_value;
+  if (b_periodic) {
+    cvm::real shift = std::floor(diff/period + 0.5);
+    diff -= shift * period;
+  }
+  return 2.0 * diff;
+}
+
+
+colvarvalue colvar::distance_z::dist2_rgrad(colvarvalue const &x1,
+                                            colvarvalue const &x2) const
+{
+  cvm::real diff = x1.real_value - x2.real_value;
+  if (b_periodic) {
+    cvm::real shift = std::floor(diff/period + 0.5);
+    diff -= shift * period;
+  }
+  return (-2.0) * diff;
+}
+
+
+void colvar::distance_z::wrap(colvarvalue &x) const
+{
+  if (!b_periodic) {
+    // don't wrap if the period has not been set
+    return;
+  }
+
+  cvm::real shift = std::floor((x.real_value - wrap_center) / period + 0.5);
+  x.real_value -= shift * period;
+  return;
+}
+
+
 
 colvar::distance_xy::distance_xy(std::string const &conf)
   : distance_z(conf)
@@ -289,6 +376,7 @@ colvar::distance_xy::distance_xy(std::string const &conf)
   x.type(colvarvalue::type_scalar);
 }
 
+
 colvar::distance_xy::distance_xy()
   : distance_z()
 {
@@ -299,6 +387,7 @@ colvar::distance_xy::distance_xy()
   x.type(colvarvalue::type_scalar);
 }
 
+
 void colvar::distance_xy::calc_value()
 {
   if (b_no_PBC) {
@@ -321,6 +410,7 @@ void colvar::distance_xy::calc_value()
   x.real_value = dist_v_ortho.norm();
 }
 
+
 void colvar::distance_xy::calc_gradients()
 {
   // Intermediate quantity (r_P3 / r_12 where P is the projection
@@ -348,6 +438,7 @@ void colvar::distance_xy::calc_gradients()
   }
 }
 
+
 void colvar::distance_xy::calc_force_invgrads()
 {
   main->read_total_forces();
@@ -360,11 +451,13 @@ void colvar::distance_xy::calc_force_invgrads()
   }
 }
 
+
 void colvar::distance_xy::calc_Jacobian_derivative()
 {
   jd.real_value = x.real_value ? (1.0 / x.real_value) : 0.0;
 }
 
+
 void colvar::distance_xy::apply_force(colvarvalue const &force)
 {
   if (!ref1->noforce)
@@ -378,6 +471,9 @@ void colvar::distance_xy::apply_force(colvarvalue const &force)
 }
 
 
+simple_scalar_dist_functions(distance_xy)
+
+
 
 colvar::distance_dir::distance_dir(std::string const &conf)
   : distance(conf)
@@ -403,7 +499,7 @@ void colvar::distance_dir::calc_value()
     dist_v = group2->center_of_mass() - group1->center_of_mass();
   } else {
     dist_v = cvm::position_distance(group1->center_of_mass(),
-                                     group2->center_of_mass());
+                                    group2->center_of_mass());
   }
   x.rvector_value = dist_v.unit();
 }
@@ -460,6 +556,7 @@ colvar::distance_inv::distance_inv(std::string const &conf)
   x.type(colvarvalue::type_scalar);
 }
 
+
 colvar::distance_inv::distance_inv()
 {
   function_type = "distance_inv";
@@ -467,6 +564,7 @@ colvar::distance_inv::distance_inv()
   x.type(colvarvalue::type_scalar);
 }
 
+
 void colvar::distance_inv::calc_value()
 {
   x.real_value = 0.0;
@@ -504,6 +602,7 @@ void colvar::distance_inv::calc_value()
   x.real_value = std::pow(x.real_value, -1.0/(cvm::real(exponent)));
 }
 
+
 void colvar::distance_inv::calc_gradients()
 {
   cvm::real const dxdsum = (-1.0/(cvm::real(exponent))) * std::pow(x.real_value, exponent+1) / cvm::real(group1->size() * group2->size());
@@ -515,6 +614,7 @@ void colvar::distance_inv::calc_gradients()
   }
 }
 
+
 void colvar::distance_inv::apply_force(colvarvalue const &force)
 {
   if (!group1->noforce)
@@ -525,6 +625,9 @@ void colvar::distance_inv::apply_force(colvarvalue const &force)
 }
 
 
+simple_scalar_dist_functions(distance_inv)
+
+
 
 colvar::distance_pairs::distance_pairs(std::string const &conf)
   : cvc(conf)
@@ -579,11 +682,13 @@ void colvar::distance_pairs::calc_value()
   }
 }
 
+
 void colvar::distance_pairs::calc_gradients()
 {
   // will be calculated on the fly in apply_force()
 }
 
+
 void colvar::distance_pairs::apply_force(colvarvalue const &force)
 {
   if (b_no_PBC) {
@@ -608,6 +713,7 @@ void colvar::distance_pairs::apply_force(colvarvalue const &force)
 }
 
 
+
 colvar::gyration::gyration(std::string const &conf)
   : cvc(conf)
 {
@@ -682,6 +788,9 @@ void colvar::gyration::apply_force(colvarvalue const &force)
 }
 
 
+simple_scalar_dist_functions(gyration)
+
+
 
 colvar::inertia::inertia(std::string const &conf)
   : gyration(conf)
@@ -722,6 +831,10 @@ void colvar::inertia::apply_force(colvarvalue const &force)
 }
 
 
+simple_scalar_dist_functions(inertia_z)
+
+
+
 colvar::inertia_z::inertia_z(std::string const &conf)
   : inertia(conf)
 {
@@ -772,6 +885,10 @@ void colvar::inertia_z::apply_force(colvarvalue const &force)
 }
 
 
+simple_scalar_dist_functions(inertia)
+
+
+
 
 colvar::rmsd::rmsd(std::string const &conf)
   : cvc(conf)
@@ -971,6 +1088,8 @@ void colvar::rmsd::calc_Jacobian_derivative()
 }
 
 
+simple_scalar_dist_functions(rmsd)
+
 
 
 colvar::eigenvector::eigenvector(std::string const &conf)
@@ -1255,6 +1374,10 @@ void colvar::eigenvector::calc_Jacobian_derivative()
 }
 
 
+simple_scalar_dist_functions(eigenvector)
+
+
+
 colvar::cartesian::cartesian(std::string const &conf)
   : cvc(conf)
 {

lib/colvars/colvarcomp_protein.cpp

@@ -1,5 +1,12 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #include <cmath>
 
 #include "colvarmodule.h"
@@ -119,6 +126,7 @@ colvar::alpha_angles::alpha_angles()
   x.type(colvarvalue::type_scalar);
 }
 
+
 colvar::alpha_angles::~alpha_angles()
 {
   while (theta.size() != 0) {
@@ -131,6 +139,7 @@ colvar::alpha_angles::~alpha_angles()
   }
 }
 
+
 void colvar::alpha_angles::calc_value()
 {
   x.real_value = 0.0;
@@ -222,6 +231,10 @@ void colvar::alpha_angles::apply_force(colvarvalue const &force)
 }
 
 
+simple_scalar_dist_functions(alpha_angles)
+
+
+
 //////////////////////////////////////////////////////////////////////
 // dihedral principal component
 //////////////////////////////////////////////////////////////////////
@@ -337,14 +350,22 @@ colvar::dihedPC::dihedPC(std::string const &conf)
   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  ], "CA", sid),
-                                           cvm::atom(r[i  ], "C", sid),
+                                         cvm::atom(r[i  ], "C", sid),
-                                           cvm::atom(r[i+1], "N", sid)));
+                                         cvm::atom(r[i+1], "N", sid)));
+    atom_groups.push_back(theta.back()->atom_groups[0]);
+    atom_groups.push_back(theta.back()->atom_groups[1]);
+    atom_groups.push_back(theta.back()->atom_groups[2]);
+    atom_groups.push_back(theta.back()->atom_groups[3]);
     // 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], "N", sid),
-                                           cvm::atom(r[i+1], "CA", sid),
+                                         cvm::atom(r[i+1], "CA", sid),
-                                           cvm::atom(r[i+1], "C", sid)));
+                                         cvm::atom(r[i+1], "C", sid)));
+    atom_groups.push_back(theta.back()->atom_groups[0]);
+    atom_groups.push_back(theta.back()->atom_groups[1]);
+    atom_groups.push_back(theta.back()->atom_groups[2]);
+    atom_groups.push_back(theta.back()->atom_groups[3]);
   }
 
   if (cvm::debug())
@@ -400,3 +421,6 @@ void colvar::dihedPC::apply_force(colvarvalue const &force)
                            coeffs[2*i+1] * dsindt) * force);
   }
 }
+
+
+simple_scalar_dist_functions(dihedPC)

lib/colvars/colvarcomp_rotations.cpp

@@ -1,5 +1,12 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #include <cmath>
 
 #include "colvarmodule.h"
@@ -123,6 +130,27 @@ void colvar::orientation::apply_force(colvarvalue const &force)
 }
 
 
+cvm::real colvar::orientation::dist2(colvarvalue const &x1,
+                                     colvarvalue const &x2) const
+{
+  return x1.quaternion_value.dist2(x2);
+}
+
+
+colvarvalue colvar::orientation::dist2_lgrad(colvarvalue const &x1,
+                                             colvarvalue const &x2) const
+{
+  return x1.quaternion_value.dist2_grad(x2);
+}
+
+
+colvarvalue colvar::orientation::dist2_rgrad(colvarvalue const &x1,
+                                             colvarvalue const &x2) const
+{
+  return x2.quaternion_value.dist2_grad(x1);
+}
+
+
 
 colvar::orientation_angle::orientation_angle(std::string const &conf)
   : orientation(conf)
@@ -176,6 +204,9 @@ void colvar::orientation_angle::apply_force(colvarvalue const &force)
 }
 
 
+simple_scalar_dist_functions(orientation_angle)
+
+
 
 colvar::orientation_proj::orientation_proj(std::string const &conf)
   : orientation(conf)
@@ -220,6 +251,9 @@ void colvar::orientation_proj::apply_force(colvarvalue const &force)
 }
 
 
+simple_scalar_dist_functions(orientation_proj)
+
+
 
 colvar::tilt::tilt(std::string const &conf)
   : orientation(conf)
@@ -278,6 +312,9 @@ void colvar::tilt::apply_force(colvarvalue const &force)
 }
 
 
+simple_scalar_dist_functions(tilt)
+
+
 
 colvar::spin_angle::spin_angle(std::string const &conf)
   : orientation(conf)
@@ -339,3 +376,46 @@ void colvar::spin_angle::apply_force(colvarvalue const &force)
     atoms->apply_colvar_force(fw);
   }
 }
+
+
+cvm::real colvar::spin_angle::dist2(colvarvalue const &x1,
+                                    colvarvalue const &x2) const
+{
+  cvm::real diff = x1.real_value - x2.real_value;
+  diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
+  return diff * diff;
+}
+
+
+colvarvalue colvar::spin_angle::dist2_lgrad(colvarvalue const &x1,
+                                            colvarvalue const &x2) const
+{
+  cvm::real diff = x1.real_value - x2.real_value;
+  diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
+  return 2.0 * diff;
+}
+
+
+colvarvalue colvar::spin_angle::dist2_rgrad(colvarvalue const &x1,
+                                            colvarvalue const &x2) const
+{
+  cvm::real diff = x1.real_value - x2.real_value;
+  diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
+  return (-2.0) * diff;
+}
+
+
+void colvar::spin_angle::wrap(colvarvalue &x) const
+{
+  if ((x.real_value - wrap_center) >= 180.0) {
+    x.real_value -= 360.0;
+    return;
+  }
+
+  if ((x.real_value - wrap_center) < -180.0) {
+    x.real_value += 360.0;
+    return;
+  }
+
+  return;
+}

lib/colvars/colvardeps.cpp

@@ -1,3 +1,13 @@
+// -*- c++ -*-
+
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
+
 #include "colvardeps.h"
 
 
@@ -219,6 +229,9 @@ void colvardeps::init_cvb_requires() {
 
   f_description(f_cvb_history_dependent, "history-dependent");
 
+  f_description(f_cvb_scalar_variables, "require scalar variables");
+  f_req_children(f_cvb_scalar_variables, f_cv_scalar);
+
   // Initialize feature_states for each instance
   feature_states.reserve(f_cvb_ntot);
   for (i = 0; i < f_cvb_ntot; i++) {
@@ -229,6 +242,9 @@ void colvardeps::init_cvb_requires() {
 
   // some biases are not history-dependent
   feature_states[f_cvb_history_dependent]->available = false;
+
+  // by default, biases should work with vector variables, too
+  feature_states[f_cvb_scalar_variables]->available = false;
 }
 
 
@@ -321,6 +337,7 @@ void colvardeps::init_cv_requires() {
     // The features below are set programmatically
     f_description(f_cv_scripted, "scripted");
     f_description(f_cv_periodic, "periodic");
+    f_req_self(f_cv_periodic, f_cv_homogeneous);
     f_description(f_cv_scalar, "scalar");
     f_description(f_cv_linear, "linear");
     f_description(f_cv_homogeneous, "homogeneous");
@@ -407,6 +424,11 @@ void colvardeps::init_cvc_requires() {
   // Each cvc specifies what other features are available
   feature_states[f_cvc_active]->available = true;
   feature_states[f_cvc_gradient]->available = true;
+
+  // Features that are implemented by default if their requirements are
+  feature_states[f_cvc_one_site_total_force]->available = true;
+
+  // Features That are implemented only for certain simulation engine configurations
   feature_states[f_cvc_scalable_com]->available = (cvm::proxy->scalable_group_coms() == COLVARS_OK);
   feature_states[f_cvc_scalable]->available = feature_states[f_cvc_scalable_com]->available;
 }

lib/colvars/colvardeps.h

@@ -1,5 +1,12 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #ifndef COLVARDEPS_H
 #define COLVARDEPS_H
 
@@ -157,6 +164,7 @@ public:
     f_cvb_apply_force, // will apply forces
     f_cvb_get_total_force, // requires total forces
     f_cvb_history_dependent, // depends on simulation history
+    f_cvb_scalar_variables, // requires scalar colvars
     f_cvb_ntot
   };
 

lib/colvars/colvargrid.cpp

@@ -1,5 +1,12 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #include "colvarmodule.h"
 #include "colvarvalue.h"
 #include "colvarparse.h"
@@ -73,6 +80,21 @@ cvm::real colvar_grid_scalar::minimum_value() const
   return min;
 }
 
+cvm::real colvar_grid_scalar::minimum_pos_value() const
+{
+  cvm::real minpos = data[0];
+  size_t i;
+  for (i = 0; i < nt; i++) {
+    if(data[i] > 0) {
+      minpos = data[i];
+      break;
+    }
+  }
+  for (i = 0; i < nt; i++) {
+    if (data[i] > 0 && data[i] < minpos) minpos = data[i];
+  }
+  return minpos;
+}
 
 cvm::real colvar_grid_scalar::integral() const
 {

lib/colvars/colvargrid.h

@@ -1,5 +1,12 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #ifndef COLVARGRID_H
 #define COLVARGRID_H
 
@@ -378,6 +385,13 @@ public:
     return value_to_bin_scalar(actual_value[i] ? cv[i]->actual_value() : cv[i]->value(), i);
   }
 
+  /// \brief Report the bin corresponding to the current value of variable i
+  /// and assign first or last bin if out of boundaries
+  inline int current_bin_scalar_bound(int const i) const
+  {
+    return value_to_bin_scalar_bound(actual_value[i] ? cv[i]->actual_value() : cv[i]->value(), i);
+  }
+
   /// \brief Report the bin corresponding to the current value of item iv in variable i
   inline int current_bin_scalar(int const i, int const iv) const
   {
@@ -393,6 +407,16 @@ public:
     return (int) std::floor( (value.real_value - lower_boundaries[i].real_value) / widths[i] );
   }
 
+  /// \brief Use the lower boundary and the width to report which bin
+  /// the provided value is in and assign first or last bin if out of boundaries
+  inline int value_to_bin_scalar_bound(colvarvalue const &value, const int i) const
+  {
+    int bin_index = std::floor( (value.real_value - lower_boundaries[i].real_value) / widths[i] );
+    if (bin_index < 0) bin_index=0;
+    if (bin_index >=int(nx[i])) bin_index=int(nx[i])-1;
+    return (int) bin_index;
+  }
+
   /// \brief Same as the standard version, but uses another grid definition
   inline int value_to_bin_scalar(colvarvalue const &value,
                                  colvarvalue const &new_offset,
@@ -514,6 +538,13 @@ public:
       data[i] *= a;
   }
 
+  /// \brief Assign all zero elements a scalar constant (fast loop)
+  inline void remove_zeros(cvm::real const &a)
+  {
+    for (size_t i = 0; i < nt; i++)
+      if(data[i]==0) data[i] = a;
+  }
+
 
   /// \brief Get the bin indices corresponding to the provided values of
   /// the colvars
@@ -537,6 +568,17 @@ public:
     return index;
   }
 
+  /// \brief Get the bin indices corresponding to the provided values of
+  /// the colvars and assign first or last bin if out of boundaries
+  inline std::vector<int> const get_colvars_index_bound() const
+  {
+    std::vector<int> index = new_index();
+    for (size_t i = 0; i < nd; i++) {
+      index[i] = current_bin_scalar_bound(i);
+    }
+    return index;
+  }
+
   /// \brief Get the minimal distance (in number of bins) from the
   /// boundaries; a negative number is returned if the given point is
   /// off-grid
@@ -1169,42 +1211,46 @@ public:
   inline cvm::real log_gradient_finite_diff(const std::vector<int> &ix0,
                                             int n = 0)
   {
-    cvm::real A0, A1;
+    int A0, A1, A2;
-    std::vector<int> ix;
+    std::vector<int> ix = ix0;
-
-    // factor for mesh width, 2.0 for central finite difference
-    // but only 1.0 on edges for non-PBC coordinates
-    cvm::real factor;
 
     if (periodic[n]) {
-      factor = 2.;
-      ix = ix0;
       ix[n]--; wrap(ix);
       A0 = data[address(ix)];
       ix = ix0;
       ix[n]++; wrap(ix);
       A1 = data[address(ix)];
-    } else {
+      if (A0 * A1 == 0) {
-      factor = 0.;
+        return 0.; // can't handle empty bins
-      ix = ix0;
+      } else {
-      if (ix[n] > 0) { // not left edge
+        return (std::log((cvm::real)A1) - std::log((cvm::real)A0))
-        ix[n]--;
+          / (widths[n] * 2.);
-        factor += 1.;
       }
+    } else if (ix[n] > 0 && ix[n] < nx[n]-1) { // not an edge
+      ix[n]--;
       A0 = data[address(ix)];
       ix = ix0;
-      if (ix[n]+1 < nx[n]) { // not right edge
+      ix[n]++;
-        ix[n]++;
-        factor += 1.;
-      }
       A1 = data[address(ix)];
-    }
+      if (A0 * A1 == 0) {
-    if (A0 == 0 || A1 == 0) {
+        return 0.; // can't handle empty bins
-      // can't handle empty bins
+      } else {
-      return 0.;
+        return (std::log((cvm::real)A1) - std::log((cvm::real)A0))
+          / (widths[n] * 2.);
+      }
     } else {
-      return (std::log((cvm::real)A1) - std::log((cvm::real)A0))
+      // edge: use 2nd order derivative
-        / (widths[n] * factor);
+      int increment = (ix[n] == 0 ? 1 : -1);
+      // move right from left edge, or the other way around
+      A0 = data[address(ix)];
+      ix[n] += increment; A1 = data[address(ix)];
+      ix[n] += increment; A2 = data[address(ix)];
+      if (A0 * A1 * A2 == 0) {
+        return 0.; // can't handle empty bins
+      } else {
+        return (-1.5 * std::log((cvm::real)A0) + 2. * std::log((cvm::real)A1)
+          - 0.5 * std::log((cvm::real)A2)) * increment / widths[n];
+      }
     }
   }
 };
@@ -1322,6 +1368,9 @@ public:
   /// \brief Return the lowest value
   cvm::real minimum_value() const;
 
+  /// \brief Return the lowest positive value
+  cvm::real minimum_pos_value() const;
+
   /// \brief Calculates the integral of the map (uses widths if they are defined)
   cvm::real integral() const;
 

lib/colvars/colvarmodule.cpp

@@ -1,5 +1,12 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #include <sstream>
 #include <string.h>
 
@@ -296,6 +303,9 @@ int colvarmodule::parse_biases(std::string const &conf)
   /// initialize harmonic restraints
   parse_biases_type<colvarbias_restraint_harmonic>(conf, "harmonic", n_rest_biases);
 
+  /// initialize harmonic walls restraints
+  parse_biases_type<colvarbias_restraint_harmonic_walls>(conf, "harmonicWalls", n_rest_biases);
+
   /// initialize histograms
   parse_biases_type<colvarbias_histogram>(conf, "histogram", n_histo_biases);
 
@@ -562,7 +572,6 @@ int colvarmodule::calc_colvars()
     colvars_smp_items.reserve(colvars.size());
 
     // set up a vector containing all components
-    size_t num_colvar_items = 0;
     cvm::increase_depth();
     for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
 
@@ -576,8 +585,6 @@ int colvarmodule::calc_colvars()
         colvars_smp.push_back(*cvi);
         colvars_smp_items.push_back(icvc);
       }
-
-      num_colvar_items += num_items;
     }
     cvm::decrease_depth();
 
@@ -641,7 +648,7 @@ int colvarmodule::calc_biases()
     for (bi = biases.begin(); bi != biases.end(); bi++) {
       error_code |= (*bi)->update();
       if (cvm::get_error()) {
-        return COLVARS_ERROR;
+        return error_code;
       }
     }
     cvm::decrease_depth();
@@ -1007,7 +1014,7 @@ std::istream & colvarmodule::read_restart(std::istream &is)
   for (std::vector<colvarbias *>::iterator bi = biases.begin();
        bi != biases.end();
        bi++) {
-    if (!((*bi)->read_restart(is))) {
+    if (!((*bi)->read_state(is))) {
       cvm::error("Error: in reading restart configuration for bias \""+
                  (*bi)->name+"\".\n",
                  INPUT_ERROR);
@@ -1070,15 +1077,15 @@ continue the previous simulation.\n\n");
     cvm::log(cvm::line_marker);
 
     // update this ahead of time in this special case
-    output_prefix = proxy->output_prefix();
+    output_prefix = proxy->input_prefix();
     cvm::log("All output files will now be saved with the prefix \""+output_prefix+".tmp.*\".\n");
-    output_prefix = output_prefix+".tmp";
-    write_output_files();
     cvm::log(cvm::line_marker);
     cvm::log("Please review the important warning above. After that, you may rename:\n\
 \""+output_prefix+".tmp.colvars.state\"\n\
 to:\n\
-\""+output_prefix+".colvars.state\"\n");
+\""+ proxy->input_prefix()+".colvars.state\"\n");
+    output_prefix = output_prefix+".tmp";
+    write_output_files();
     cvm::error("Exiting with error until issue is addressed.\n", FATAL_ERROR);
   }
 
@@ -1120,6 +1127,7 @@ int colvarmodule::write_output_files()
        bi != biases.end();
        bi++) {
     (*bi)->write_output_files();
+    (*bi)->write_state_to_replicas();
   }
   cvm::decrease_depth();
 
@@ -1212,20 +1220,30 @@ std::ostream & colvarmodule::write_restart(std::ostream &os)
      << "  version " << std::string(COLVARS_VERSION) << "\n"
      << "}\n\n";
 
+  int error_code = COLVARS_OK;
+
   cvm::increase_depth();
   for (std::vector<colvar *>::iterator cvi = colvars.begin();
        cvi != colvars.end();
        cvi++) {
     (*cvi)->write_restart(os);
+    error_code |= (*cvi)->write_output_files();
   }
 
   for (std::vector<colvarbias *>::iterator bi = biases.begin();
        bi != biases.end();
        bi++) {
-    (*bi)->write_restart(os);
+    (*bi)->write_state(os);
+    error_code |= (*bi)->write_state_to_replicas();
+    error_code |= (*bi)->write_output_files();
   }
   cvm::decrease_depth();
 
+  if (error_code != COLVARS_OK) {
+    // TODO make this function return an int instead
+    os.setstate(std::ios::failbit);
+  }
+
   return os;
 }
 

lib/colvars/colvarmodule.h

@@ -1,10 +1,17 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #ifndef COLVARMODULE_H
 #define COLVARMODULE_H
 
 #ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2016-10-27"
+#define COLVARS_VERSION "2016-12-22"
 #endif
 
 #ifndef COLVARS_DEBUG

lib/colvars/colvarparse.cpp

@@ -1,3 +1,11 @@
+// -*- c++ -*-
+
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
 
 
 #include <sstream>

lib/colvars/colvarparse.h

@@ -1,5 +1,12 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #ifndef COLVARPARSE_H
 #define COLVARPARSE_H
 

lib/colvars/colvarproxy.h

@@ -1,5 +1,12 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #ifndef COLVARPROXY_H
 #define COLVARPROXY_H
 
@@ -24,11 +31,18 @@ public:
   /// Pointer to the main object
   colvarmodule *colvars;
 
-  /// Default constructor
+  /// Constructor
-  inline colvarproxy() : b_smp_active(true), script(NULL) {}
+  colvarproxy()
+  {
+    colvars = NULL;
+    b_simulation_running = true;
+    b_smp_active = true;
+    script = NULL;
+  }
 
-  /// Default destructor
+  /// Destructor
-  virtual ~colvarproxy() {}
+  virtual ~colvarproxy()
+  {}
 
   /// (Re)initialize required member data after construction
   virtual int setup()
@@ -104,6 +118,19 @@ public:
     return 0;
   }
 
+protected:
+
+  /// Whether a simulation is running (and try to prevent irrecovarable errors)
+  bool b_simulation_running;
+
+public:
+
+  /// Whether a simulation is running (and try to prevent irrecovarable errors)
+  virtual bool simulation_running() const
+  {
+    return b_simulation_running;
+  }
+
 protected:
 
   /// \brief Currently opened output files: by default, these are ofstream objects.

lib/colvars/colvarscript.cpp

@@ -1,5 +1,12 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #include <cstdlib>
 #include <stdlib.h>
 #include <string.h>

lib/colvars/colvarscript.h

@@ -1,5 +1,12 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #ifndef COLVARSCRIPT_H
 #define COLVARSCRIPT_H
 

lib/colvars/colvartypes.cpp

@@ -1,5 +1,12 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #include <stdlib.h>
 #include <string.h>
 

lib/colvars/colvartypes.h

@@ -1,5 +1,12 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #ifndef COLVARTYPES_H
 #define COLVARTYPES_H
 

lib/colvars/colvarvalue.cpp

@@ -1,5 +1,12 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #include <vector>
 #include <sstream>
 #include <iostream>
@@ -72,6 +79,7 @@ void colvarvalue::set_elem(int const icv, colvarvalue const &x)
 
 void colvarvalue::set_random()
 {
+  size_t ic;
   switch (this->type()) {
   case colvarvalue::type_scalar:
     this->real_value = cvm::rand_gaussian();
@@ -91,7 +99,7 @@ void colvarvalue::set_random()
     this->quaternion_value.q3 = cvm::rand_gaussian();
     break;
   case colvarvalue::type_vector:
-    for (size_t ic = 0; ic < this->vector1d_value.size(); ic++) {
+    for (ic = 0; ic < this->vector1d_value.size(); ic++) {
       this->vector1d_value[ic] = cvm::rand_gaussian();
     }
     break;
@@ -103,56 +111,6 @@ void colvarvalue::set_random()
 }
 
 
-colvarvalue colvarvalue::inverse() const
-{
-  switch (value_type) {
-  case colvarvalue::type_scalar:
-    return colvarvalue(1.0/real_value);
-    break;
-  case colvarvalue::type_3vector:
-  case colvarvalue::type_unit3vector:
-  case colvarvalue::type_unit3vectorderiv:
-    return colvarvalue(cvm::rvector(1.0/rvector_value.x,
-                                    1.0/rvector_value.y,
-                                    1.0/rvector_value.z));
-    break;
-  case colvarvalue::type_quaternion:
-  case colvarvalue::type_quaternionderiv:
-    return colvarvalue(cvm::quaternion(1.0/quaternion_value.q0,
-                                       1.0/quaternion_value.q1,
-                                       1.0/quaternion_value.q2,
-                                       1.0/quaternion_value.q3));
-    break;
-  case colvarvalue::type_vector:
-    {
-      cvm::vector1d<cvm::real> result(vector1d_value);
-      if (elem_types.size() > 0) {
-        // if we have information about non-scalar types, use it
-        size_t i;
-        for (i = 0; i < elem_types.size(); i++) {
-          result.sliceassign(elem_indices[i], elem_indices[i]+elem_sizes[i],
-                             cvm::vector1d<cvm::real>((this->get_elem(i)).inverse()));
-        }
-      } else {
-        size_t i;
-        for (i = 0; i < result.size(); i++) {
-          if (result[i] != 0.0) {
-            result = 1.0/result[i];
-          }
-        }
-      }
-      return colvarvalue(result, type_vector);
-    }
-    break;
-  case colvarvalue::type_notset:
-  default:
-    undef_op();
-    break;
-  }
-  return colvarvalue();
-}
-
-
 // binary operations between two colvarvalues
 
 colvarvalue operator + (colvarvalue const &x1,
@@ -321,7 +279,7 @@ colvarvalue colvarvalue::dist2_grad(colvarvalue const &x2) const
                                        (-1.0) * sin_t * v2.z +
                                        cos_t/sin_t * (v1.z - cos_t*v2.z)
                                        ),
-                          colvarvalue::type_unit3vector );
+                          colvarvalue::type_unit3vectorderiv );
     }
   case colvarvalue::type_quaternion:
   case colvarvalue::type_quaternionderiv:

lib/colvars/colvarvalue.h

@@ -1,5 +1,12 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 #ifndef COLVARVALUE_H
 #define COLVARVALUE_H
 
@@ -128,30 +135,23 @@ public:
   {}
 
   /// \brief Copy constructor from rvector base type (Note: this sets
-  /// automatically a type \link type_3vector \endlink , if you want a
+  /// by default a type \link type_3vector \endlink , if you want a
   /// \link type_unit3vector \endlink you must set it explicitly)
-  inline colvarvalue(cvm::rvector const &v)
+  inline colvarvalue(cvm::rvector const &v, Type vti = type_3vector)
-    : value_type(type_3vector), rvector_value(v)
-  {}
-
-  /// \brief Copy constructor from rvector base type (additional
-  /// argument to make possible to choose a \link type_unit3vector
-  /// \endlink
-  inline colvarvalue(cvm::rvector const &v, Type const &vti)
     : value_type(vti), rvector_value(v)
   {}
 
   /// \brief Copy constructor from quaternion base type
-  inline colvarvalue(cvm::quaternion const &q)
+  inline colvarvalue(cvm::quaternion const &q, Type vti = type_quaternion)
-    : value_type(type_quaternion), quaternion_value(q)
+    : value_type(vti), quaternion_value(q)
   {}
 
+  /// Copy constructor from vector1d base type
+  colvarvalue(cvm::vector1d<cvm::real> const &v, Type vti = type_vector);
+
   /// Copy constructor from another \link colvarvalue \endlink
   colvarvalue(colvarvalue const &x);
 
-  /// Copy constructor from vector1d base type
-  colvarvalue(cvm::vector1d<cvm::real> const &v, Type const &vti);
-
 
   /// Set to the null value for the data type currently defined
   void reset();
@@ -211,10 +211,6 @@ public:
     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;
 
@@ -536,7 +532,7 @@ inline colvarvalue::colvarvalue(colvarvalue const &x)
   }
 }
 
-inline colvarvalue::colvarvalue(cvm::vector1d<cvm::real> const &v, Type const &vti)
+inline colvarvalue::colvarvalue(cvm::vector1d<cvm::real> const &v, Type vti)
 {
   if ((vti != type_vector) && (v.size() != num_dimensions(vti))) {
     cvm::error("Error: trying to initialize a variable of type \""+type_desc(vti)+
@@ -620,11 +616,22 @@ inline int colvarvalue::check_types_assign(colvarvalue::Type const &vt1,
   }
 
   if (vt1 != type_notset) {
-    if (vt1 != vt2) {
+    if (((vt1 == type_unit3vector) &&
-      cvm::error("Trying to assign a colvar value with type \""+
+         (vt2 == type_unit3vectorderiv)) ||
-                 type_desc(vt2)+"\" to one with type \""+
+        ((vt2 == type_unit3vector) &&
-                 type_desc(vt1)+"\".\n");
+         (vt1 == type_unit3vectorderiv)) ||
-      return COLVARS_ERROR;
+        ((vt1 == type_quaternion) &&
+         (vt2 == type_quaternionderiv)) ||
+        ((vt2 == type_quaternion) &&
+         (vt1 == type_quaternionderiv))) {
+      return COLVARS_OK;
+    } else {
+      if (vt1 != vt2) {
+        cvm::error("Trying to assign a colvar value with type \""+
+                   type_desc(vt2)+"\" to one with type \""+
+                   type_desc(vt1)+"\".\n");
+        return COLVARS_ERROR;
+      }
     }
   }
   return COLVARS_OK;

src/USER-COLVARS/README

@@ -1,49 +1,46 @@
-This package implements the "fix colvars" command which can be used
+This package implements the "fix colvars" command, which can be used
 in a LAMMPS input script.
 
-This fix allows use of "collective variables" to implement Adaptive
+The fix allows use of enhanced sampling methods based on a reduced set
-Biasing Force, Metadynamics, Steered MD, Umbrella Sampling and
+of "collective variables", including free-energy estimators based on
-Restraints.
+thermodynamic forces, non-equilibrium work and probability
+distributions.
 
-This package uses an external library in lib/colvars which must be
+The package uses the "Colvars" library, whose source code is included
-compiled before making LAMMPS.  See the lib/colvars/README file and
+in the LAMMPS source code distribution and must be linked with LAMMPS.
-the LAMMPS manual for information on building LAMMPS with external
+See the lib/colvars/README file and the LAMMPS manual for information
-libraries.  The settings in the Makefile.lammps file in that directory
+on building LAMMPS with external libraries.  The settings in the
-must be correct for LAMMPS to build correctly with this package
+Makefile.lammps file in that directory must be correct for LAMMPS to
-installed.
+build correctly with this package installed.
 
-The external library is a portable collective variable module library
+The files in the USER-COLVARS package folder implement an interface
-written and maintained by Giacomo Fiorin (ICMS, Temple University,
+between LAMMPS and Colvars, originally written by Axel Kohlmeyer
-Philadelphia, PA, USA) and Jerome Henin (IBPC, CNRS, Paris, France).
+(akohlmey@gmail.com) and maintained by Giacomo Fiorin
+(giacomo.fiorin@gmail.com).
 
-More info about this library can be found at: http://colvars.github.io
+More info about the Colvars library can be found at:
 
-and in these publications:
+https://github.com/colvars/colvars
+
+and in the reference article:
 
 Using collective variables to drive molecular dynamics simulations,
-Giacomo Fiorin, Michael L. Klein & Jérôme Hénin: Molecular Physics,
+G. Fiorin, M. L. Klein, and J. Henin,
-111, 3345-3362 (2013)
+Molecular Physics 111, 3345 (2013)
+http://dx.doi.org/10.1080/00268976.2013.813594
 
-Exploring Multidimensional Free Energy Landscapes Using Time-Dependent
+A reference manual for the package and library is included with the
-Biases on Collective Variables, J. Hénin, G. Fiorin, C. Chipot, and
+LAMMPS doc pages:
-M. L. Klein, J. Chem. Theory Comput., 6, 35-47 (2010).
+doc/PDF/colvars-refman-lammps.pdf
+which also includes citations to the articles documenting the various
+methods that make use Colvars.
 
-The colvars fix implementes a thin interface layer, which exchanges
+There are also example scripts for using this package in the folder
-information between LAMMPS and the collective variable module.  This
+examples/USER/colvars, as well as the GitHub page for Colvars.
-interface was written and is maintained by Axel Kohlmeyer
-(akohlmey@gmail.com)
 
-See the doc page of fix colvars for more details.
+Please contact Giacomo Fiorin (giacomo.fiorin@gmail.com) for questions
-
+regarding this package.
-There is a reference manual for the package included with the LAMMPS
-doc pages: doc/PDF/colvars-refman-lammps.pdf
-
-There are example scripts for using this package in
-examples/USER/colvars.
-
-The person who created this package is Axel Kohlmeyer at Temple U
-(akohlmey at gmail.com).  Contact him directly if you have questions.
 
 ---------------------------------
 
-Version: 2015-01-08
+Version: 2016-12-22

src/USER-COLVARS/colvarproxy_lammps.cpp

@@ -1,3 +1,12 @@
+// -*- c++ -*-
+
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 
 #include <mpi.h>
 #include "lammps.h"

src/USER-COLVARS/colvarproxy_lammps.h

@@ -1,5 +1,13 @@
 // -*- c++ -*-
 
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
+
 #ifndef COLVARPROXY_LAMMPS_H
 #define COLVARPROXY_LAMMPS_H
 
@@ -21,7 +29,7 @@
 #endif
 
 #ifndef COLVARPROXY_VERSION
-#define COLVARPROXY_VERSION "2016-10-05"
+#define COLVARPROXY_VERSION "2016-12-27"
 #endif
 
 /* struct for packed data communication of coordinates and forces. */

src/USER-COLVARS/fix_colvars.cpp

@@ -1,3 +1,12 @@
+// -*- c++ -*-
+
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories

src/USER-COLVARS/fix_colvars.h

@@ -1,3 +1,12 @@
+// -*- c++ -*-
+
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
 /* -*- c++ -*- ----------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories

src/USER-COLVARS/group_ndx.cpp

@@ -1,3 +1,5 @@
+// -*- c++ -*-
+
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories

src/USER-COLVARS/group_ndx.h

@@ -1,4 +1,6 @@
-/* -*- c++ -*- ----------------------------------------------------------
+// -*- c++ -*-
+
+/* ----------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov

src/USER-COLVARS/ndx_group.cpp

@@ -1,3 +1,5 @@
+// -*- c++ -*-
+
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
@@ -10,9 +12,8 @@
 
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
-
 /* ----------------------------------------------------------------------
-   Contributing author: Axel Kohlmeyer (Temple U)
+   Contributing author:  Axel Kohlmeyer (Temple U)
 ------------------------------------------------------------------------- */
 
 #include "ndx_group.h"
@@ -35,7 +36,7 @@ static char *find_section(FILE *fp, const char *name)
 {
   char linebuf[BUFLEN];
   char *n,*p,*t,*r;
-  
+
   while ((p = fgets(linebuf,BUFLEN,fp))) {
     t = strtok(p," \t\n\r\f");
     if ((t != NULL) && *t == '[') {

src/USER-COLVARS/ndx_group.h

@@ -1,4 +1,6 @@
-/* -*- c++ -*- ----------------------------------------------------------
+// -*- c++ -*-
+
+/* ----------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov