Update Colvars to version 2020-02-27

This update contains several small new features or usability improvements.
    Descriptions and authorship information can be accessed from the pull
    requests listed below.

    Allow setting sigma parameter directly for metadynamics
    https://github.com/Colvars/colvars/pull/325

    Remove default values for lowerWall and upperWall legacy keywords
    https://github.com/Colvars/colvars/pull/324

    biasActualColvar option to bypass extended-Lagragian for a bias
    https://github.com/Colvars/colvars/pull/321

    Flexible restart-reading
    https://github.com/Colvars/colvars/pull/320

lib/colvars/Makefile.deps

@@ -27,46 +27,58 @@ $(COLVARS_OBJ_DIR)colvarbias_restraint.o: colvarbias_restraint.cpp \
 $(COLVARS_OBJ_DIR)colvarcomp_angles.o: colvarcomp_angles.cpp \
  colvarmodule.h colvars_version.h colvar.h colvarvalue.h colvartypes.h \
  colvarparse.h colvarparams.h colvardeps.h colvarcomp.h colvaratoms.h \
- colvarproxy.h colvar_arithmeticpath.h
+ colvarproxy.h colvar_arithmeticpath.h colvar_geometricpath.h
-$(COLVARS_OBJ_DIR)colvarcomp_apath.o: colvarcomp_apath.cpp
+$(COLVARS_OBJ_DIR)colvarcomp_apath.o: colvarcomp_apath.cpp colvarmodule.h \
+ colvars_version.h colvarvalue.h colvartypes.h colvarparse.h \
+ colvarparams.h colvar.h colvardeps.h colvarcomp.h colvaratoms.h \
+ colvarproxy.h colvar_arithmeticpath.h colvar_geometricpath.h
 $(COLVARS_OBJ_DIR)colvarcomp_coordnums.o: colvarcomp_coordnums.cpp \
  colvarmodule.h colvars_version.h colvarparse.h colvarvalue.h \
  colvartypes.h colvarparams.h colvaratoms.h colvarproxy.h colvardeps.h \
- colvar.h colvarcomp.h colvar_arithmeticpath.h
+ colvar.h colvarcomp.h colvar_arithmeticpath.h colvar_geometricpath.h
 $(COLVARS_OBJ_DIR)colvarcomp.o: colvarcomp.cpp colvarmodule.h \
  colvars_version.h colvarvalue.h colvartypes.h colvar.h colvarparse.h \
  colvarparams.h colvardeps.h colvarcomp.h colvaratoms.h colvarproxy.h \
- colvar_arithmeticpath.h
+ colvar_arithmeticpath.h colvar_geometricpath.h
 $(COLVARS_OBJ_DIR)colvarcomp_distances.o: colvarcomp_distances.cpp \
  colvarmodule.h colvars_version.h colvarvalue.h colvartypes.h \
  colvarparse.h colvarparams.h colvar.h colvardeps.h colvarcomp.h \
- colvaratoms.h colvarproxy.h colvar_arithmeticpath.h
+ colvaratoms.h colvarproxy.h colvar_arithmeticpath.h \
-$(COLVARS_OBJ_DIR)colvarcomp_gpath.o: colvarcomp_gpath.cpp
+ colvar_geometricpath.h
+$(COLVARS_OBJ_DIR)colvarcomp_gpath.o: colvarcomp_gpath.cpp colvarmodule.h \
+ colvars_version.h colvarvalue.h colvartypes.h colvarparse.h \
+ colvarparams.h colvar.h colvardeps.h colvarcomp.h colvaratoms.h \
+ colvarproxy.h colvar_arithmeticpath.h colvar_geometricpath.h
 $(COLVARS_OBJ_DIR)colvarcomp_protein.o: colvarcomp_protein.cpp \
  colvarmodule.h colvars_version.h colvarvalue.h colvartypes.h \
  colvarparse.h colvarparams.h colvar.h colvardeps.h colvarcomp.h \
- colvaratoms.h colvarproxy.h colvar_arithmeticpath.h
+ colvaratoms.h colvarproxy.h colvar_arithmeticpath.h \
+ colvar_geometricpath.h
 $(COLVARS_OBJ_DIR)colvarcomp_rotations.o: colvarcomp_rotations.cpp \
  colvarmodule.h colvars_version.h colvarvalue.h colvartypes.h \
  colvarparse.h colvarparams.h colvar.h colvardeps.h colvarcomp.h \
- colvaratoms.h colvarproxy.h colvar_arithmeticpath.h
+ colvaratoms.h colvarproxy.h colvar_arithmeticpath.h \
+ colvar_geometricpath.h
 $(COLVARS_OBJ_DIR)colvar.o: colvar.cpp colvarmodule.h colvars_version.h \
  colvarvalue.h colvartypes.h colvarparse.h colvarparams.h colvar.h \
  colvardeps.h colvarcomp.h colvaratoms.h colvarproxy.h \
- colvar_arithmeticpath.h colvarscript.h colvarbias.h
+ colvar_arithmeticpath.h colvar_geometricpath.h colvarscript.h \
+ colvarbias.h
 $(COLVARS_OBJ_DIR)colvardeps.o: colvardeps.cpp colvarmodule.h \
  colvars_version.h colvarproxy.h colvartypes.h colvarvalue.h colvardeps.h \
  colvarparse.h colvarparams.h
 $(COLVARS_OBJ_DIR)colvargrid.o: colvargrid.cpp colvarmodule.h \
  colvars_version.h colvarvalue.h colvartypes.h colvarparse.h \
  colvarparams.h colvar.h colvardeps.h colvarcomp.h colvaratoms.h \
- colvarproxy.h colvar_arithmeticpath.h colvargrid.h
+ colvarproxy.h colvar_arithmeticpath.h colvar_geometricpath.h \
+ colvargrid.h
 $(COLVARS_OBJ_DIR)colvarmodule.o: colvarmodule.cpp colvarmodule.h \
  colvars_version.h colvarparse.h colvarvalue.h colvartypes.h \
  colvarparams.h colvarproxy.h colvar.h colvardeps.h colvarbias.h \
  colvarbias_abf.h colvargrid.h colvar_UIestimator.h colvarbias_alb.h \
  colvarbias_histogram.h colvarbias_meta.h colvarbias_restraint.h \
- colvarscript.h colvaratoms.h colvarcomp.h colvar_arithmeticpath.h
+ colvarscript.h colvaratoms.h colvarcomp.h colvar_arithmeticpath.h \
+ colvar_geometricpath.h
 $(COLVARS_OBJ_DIR)colvarparams.o: colvarparams.cpp colvarmodule.h \
  colvars_version.h colvarvalue.h colvartypes.h colvarparams.h
 $(COLVARS_OBJ_DIR)colvarparse.o: colvarparse.cpp colvarmodule.h \

lib/colvars/colvar.cpp

@@ -33,6 +33,8 @@ colvar::colvar()
   dev_null = 0.0;
 #endif
 
+  expand_boundaries = false;
+
   description = "uninitialized colvar";
   init_dependencies();
 }
@@ -474,6 +476,8 @@ int colvar::init_custom_function(std::string const &conf)
 
 int colvar::init_grid_parameters(std::string const &conf)
 {
+  int error_code = COLVARS_OK;
+
   colvarmodule *cv = cvm::main();
 
   cvm::real default_width = width;
@@ -527,15 +531,18 @@ int colvar::init_grid_parameters(std::string const &conf)
       disable(f_cv_hard_upper_boundary);
     }
 
+    // Parse legacy wall options and set up a harmonicWalls bias if needed
+    cvm::real lower_wall_k = 0.0, upper_wall_k = 0.0;
+    cvm::real lower_wall = 0.0, upper_wall = 0.0;
     std::string lw_conf, uw_conf;
+
     if (get_keyval(conf, "lowerWallConstant", lower_wall_k, 0.0,
                    parse_silent)) {
       cvm::log("Reading legacy options lowerWall and lowerWallConstant: "
-               "consider using a harmonicWalls restraint\n(caution: force constant would then be scaled by width^2).\n");
+               "consider using a harmonicWalls restraint (caution: force constant would then be scaled by width^2).\n");
-      lower_wall.type(value());
+      if (!get_keyval(conf, "lowerWall", lower_wall)) {
-      if (!get_keyval(conf, "lowerWall", lower_wall, lower_boundary)) {
+        error_code != cvm::error("Error: the value of lowerWall must be set "
-        cvm::log("Warning: lowerWall will need to be "
+                                 "explicitly.\n", INPUT_ERROR);
-                 "defined explicitly in the next release.\n");
       }
       lw_conf = std::string("\n\
     lowerWallConstant "+cvm::to_str(lower_wall_k*width*width)+"\n\
@@ -545,11 +552,10 @@ int colvar::init_grid_parameters(std::string const &conf)
     if (get_keyval(conf, "upperWallConstant", upper_wall_k, 0.0,
                    parse_silent)) {
       cvm::log("Reading legacy options upperWall and upperWallConstant: "
-               "consider using a harmonicWalls restraint\n(caution: force constant would then be scaled by width^2).\n");
+               "consider using a harmonicWalls restraint (caution: force constant would then be scaled by width^2).\n");
-      upper_wall.type(value());
+      if (!get_keyval(conf, "upperWall", upper_wall)) {
-      if (!get_keyval(conf, "upperWall", upper_wall, upper_boundary)) {
+        error_code != cvm::error("Error: the value of upperWall must be set "
-        cvm::log("Warning: upperWall will need to be "
+                                 "explicitly.\n", INPUT_ERROR);
-                 "defined explicitly in the next release.\n");
       }
       uw_conf = std::string("\n\
     upperWallConstant "+cvm::to_str(upper_wall_k*width*width)+"\n\
@@ -558,12 +564,11 @@ int colvar::init_grid_parameters(std::string const &conf)
 
     if (lw_conf.size() && uw_conf.size()) {
       if (lower_wall >= upper_wall) {
-        cvm::error("Error: the upper wall, "+
+        error_code |= cvm::error("Error: the upper wall, "+
-                   cvm::to_str(upper_wall)+
+                                 cvm::to_str(upper_wall)+
-                   ", is not higher than the lower wall, "+
+                                 ", is not higher than the lower wall, "+
-                   cvm::to_str(lower_wall)+".\n",
+                                 cvm::to_str(lower_wall)+".\n",
-                   INPUT_ERROR);
+                                 INPUT_ERROR);
-        return INPUT_ERROR;
       }
     }
 
@@ -575,7 +580,7 @@ harmonicWalls {\n\
     colvars "+this->name+"\n"+lw_conf+uw_conf+"\
     timeStepFactor "+cvm::to_str(time_step_factor)+"\n"+
                              "}\n");
-      cv->append_new_config(walls_conf);
+      error_code |= cv->append_new_config(walls_conf);
     }
   }
 
@@ -588,29 +593,27 @@ harmonicWalls {\n\
   // consistency checks for boundaries and walls
   if (is_enabled(f_cv_lower_boundary) && is_enabled(f_cv_upper_boundary)) {
     if (lower_boundary >= upper_boundary) {
-      cvm::error("Error: the upper boundary, "+
+      error_code |= cvm::error("Error: the upper boundary, "+
-                        cvm::to_str(upper_boundary)+
+                               cvm::to_str(upper_boundary)+
-                        ", is not higher than the lower boundary, "+
+                               ", is not higher than the lower boundary, "+
-                        cvm::to_str(lower_boundary)+".\n",
+                               cvm::to_str(lower_boundary)+".\n",
-                INPUT_ERROR);
+                               INPUT_ERROR);
-      return INPUT_ERROR;
     }
   }
 
-  get_keyval(conf, "expandBoundaries", expand_boundaries, false);
+  get_keyval(conf, "expandBoundaries", expand_boundaries, expand_boundaries);
   if (expand_boundaries && periodic_boundaries()) {
-    cvm::error("Error: trying to expand boundaries that already "
+    error_code |= cvm::error("Error: trying to expand boundaries that already "
-               "cover a whole period of a periodic colvar.\n",
+                             "cover a whole period of a periodic colvar.\n",
-               INPUT_ERROR);
+                             INPUT_ERROR);
-    return INPUT_ERROR;
   }
+
   if (expand_boundaries && is_enabled(f_cv_hard_lower_boundary) &&
       is_enabled(f_cv_hard_upper_boundary)) {
-    cvm::error("Error: inconsistent configuration "
+    error_code |= cvm::error("Error: inconsistent configuration "
-               "(trying to expand boundaries with both "
+                             "(trying to expand boundaries, but both "
-               "hardLowerBoundary and hardUpperBoundary enabled).\n",
+                             "hardLowerBoundary and hardUpperBoundary "
-               INPUT_ERROR);
+                             "are enabled).\n", INPUT_ERROR);
-    return INPUT_ERROR;
   }
 
   return COLVARS_OK;
@@ -2077,12 +2080,12 @@ void colvar::wrap(colvarvalue &x_unwrapped) const
 
 // ******************** INPUT FUNCTIONS ********************
 
-std::istream & colvar::read_restart(std::istream &is)
+std::istream & colvar::read_state(std::istream &is)
 {
   size_t const start_pos = is.tellg();
 
   std::string conf;
-  if ( !(is >> colvarparse::read_block("colvar", conf)) ) {
+  if ( !(is >> colvarparse::read_block("colvar", &conf)) ) {
     // this is not a colvar block
     is.clear();
     is.seekg(start_pos, std::ios::beg);
@@ -2092,15 +2095,24 @@ std::istream & colvar::read_restart(std::istream &is)
 
   {
     std::string check_name = "";
-    if ( (get_keyval(conf, "name", check_name,
+    get_keyval(conf, "name", check_name,
-                     std::string(""), colvarparse::parse_silent)) &&
+               std::string(""), colvarparse::parse_silent);
-         (check_name != name) )  {
-      cvm::error("Error: the state file does not match the "
-                 "configuration file, at colvar \""+name+"\".\n");
-    }
     if (check_name.size() == 0) {
       cvm::error("Error: Collective variable in the "
-                 "restart file without any identifier.\n");
+                 "restart file without any identifier.\n", INPUT_ERROR);
+      is.clear();
+      is.seekg(start_pos, std::ios::beg);
+      is.setstate(std::ios::failbit);
+      return is;
+    }
+
+    if (check_name != name)  {
+      if (cvm::debug()) {
+        cvm::log("Ignoring state of colvar \""+check_name+
+                 "\": this colvar is named \""+name+"\".\n");
+      }
+      is.seekg(start_pos, std::ios::beg);
+      return is;
     }
   }
 
@@ -2199,7 +2211,7 @@ std::istream & colvar::read_traj(std::istream &is)
 
 // ******************** OUTPUT FUNCTIONS ********************
 
-std::ostream & colvar::write_restart(std::ostream &os) {
+std::ostream & colvar::write_state(std::ostream &os) {
 
   os << "colvar {\n"
      << "  name " << name << "\n"

lib/colvars/colvar.h

@@ -229,17 +229,9 @@ public:
 
   /// \brief Location of the lower boundary
   colvarvalue lower_boundary;
-  /// \brief Location of the lower wall
-  colvarvalue lower_wall;
-  /// \brief Force constant for the lower boundary potential (|x-xb|^2)
-  cvm::real   lower_wall_k;
 
   /// \brief Location of the upper boundary
   colvarvalue upper_boundary;
-  /// \brief Location of the upper wall
-  colvarvalue upper_wall;
-  /// \brief Force constant for the upper boundary potential (|x-xb|^2)
-  cvm::real   upper_wall_k;
 
   /// \brief Is the interval defined by the two boundaries periodic?
   bool periodic_boundaries() const;
@@ -453,9 +445,9 @@ public:
   std::ostream & write_traj_label(std::ostream &os);
 
   /// Read the collective variable from a restart file
-  std::istream & read_restart(std::istream &is);
+  std::istream & read_state(std::istream &is);
   /// Write the collective variable to a restart file
-  std::ostream & write_restart(std::ostream &os);
+  std::ostream & write_state(std::ostream &os);
 
   /// Write output files (if defined, e.g. in analysis mode)
   int write_output_files();

lib/colvars/colvarbias.cpp

@@ -15,8 +15,10 @@
 
 
 colvarbias::colvarbias(char const *key)
-  : bias_type(to_lower_cppstr(key))
 {
+  bias_type = to_lower_cppstr(key);
+  state_keyword = bias_type;
+
   description = "uninitialized " + cvm::to_str(key) + " bias";
   init_dependencies();
   rank = 1;
@@ -25,6 +27,7 @@ colvarbias::colvarbias(char const *key)
   b_output_energy = false;
   reset();
   state_file_step = 0L;
+  matching_state = false;
 }
 
 
@@ -80,10 +83,21 @@ int colvarbias::init(std::string const &conf)
     cvm::log("Reinitializing bias \""+name+"\".\n");
   }
 
+  colvar_values.resize(num_variables());
+
+  for (i = 0; i < num_variables(); i++) {
+    colvar_values[i].type(colvars[i]->value().type());
+    colvar_forces[i].type(colvar_values[i].type());
+    previous_colvar_forces[i].type(colvar_values[i].type());
+  }
+
   output_prefix = cvm::output_prefix();
 
   get_keyval(conf, "outputEnergy", b_output_energy, b_output_energy);
 
+  // Disabled by default in base class; default value can be overridden by derived class constructor
+  get_keyval_feature(this, conf, "bypassExtendedLagrangian", f_cvb_bypass_ext_lagrangian, is_enabled(f_cvb_bypass_ext_lagrangian), parse_silent);
+
   get_keyval(conf, "timeStepFactor", time_step_factor, 1);
   if (time_step_factor < 1) {
     cvm::error("Error: timeStepFactor must be 1 or greater.\n");
@@ -114,6 +128,11 @@ int colvarbias::init_dependencies() {
     init_feature(f_cvb_apply_force, "apply force", f_type_user);
     require_feature_children(f_cvb_apply_force, f_cv_gradient);
 
+    init_feature(f_cvb_bypass_ext_lagrangian, "bypass extended-Lagrangian coordinates", f_type_user);
+    // The exclusion below prevents the inconsistency where biasing forces are applied onto
+    // the actual colvar, while total forces are measured on the extended coordinate
+    exclude_feature_self(f_cvb_bypass_ext_lagrangian, f_cvb_get_total_force);
+
     init_feature(f_cvb_get_total_force, "obtain total force", f_type_dynamic);
     require_feature_children(f_cvb_get_total_force, f_cv_total_force);
 
@@ -158,6 +177,12 @@ int colvarbias::init_dependencies() {
   // only compute TI samples when deriving from colvarbias_ti
   feature_states[f_cvb_calc_ti_samples].available = false;
 
+  // The feature f_cvb_bypass_ext_lagrangian is only implemented by some derived classes
+  // (initially, harmonicWalls)
+  feature_states[f_cvb_bypass_ext_lagrangian].available = false;
+  // disabled by default; can be changed by derived classes that implement it
+  feature_states[f_cvb_bypass_ext_lagrangian].enabled = false;
+
   return COLVARS_OK;
 }
 
@@ -265,6 +290,11 @@ int colvarbias::update()
 
   has_data = true;
 
+  // Update the cached colvar values
+  for (size_t i = 0; i < num_variables(); i++) {
+    colvar_values[i] = colvars[i]->value();
+  }
+
   error_code |= calc_energy(NULL);
   error_code |= calc_forces(NULL);
 
@@ -304,9 +334,11 @@ void colvarbias::communicate_forces()
     // may send forces to the same colvar
     // which is why rescaling has to happen now: the colvar is not
     // aware of this bias' time_step_factor
-    variables(i)->add_bias_force(cvm::real(time_step_factor) * colvar_forces[i]);
+    if (is_enabled(f_cvb_bypass_ext_lagrangian)) {
-  }
+      variables(i)->add_bias_force_actual_value(cvm::real(time_step_factor) * colvar_forces[i]);
-  for (i = 0; i < num_variables(); i++) {
+    } else {
+      variables(i)->add_bias_force(cvm::real(time_step_factor) * colvar_forces[i]);
+    }
     previous_colvar_forces[i] = colvar_forces[i];
   }
 }
@@ -368,20 +400,27 @@ std::string const colvarbias::get_state_params() const
 
 int colvarbias::set_state_params(std::string const &conf)
 {
-  std::string new_name = "";
+  matching_state = false;
-  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) {
+  std::string check_name = "";
+  colvarparse::get_keyval(conf, "name", check_name,
+                          std::string(""), colvarparse::parse_silent);
+
+  if (check_name.size() == 0) {
     cvm::error("Error: \""+bias_type+"\" block within the restart file "
                "has no identifiers.\n", INPUT_ERROR);
   }
 
+  if (check_name != this->name) {
+    if (cvm::debug()) {
+      cvm::log("Ignoring state of bias \""+check_name+
+               "\": this bias is named \""+name+"\".\n");
+    }
+    return COLVARS_OK;
+  }
+
+  matching_state = true;
+
   colvarparse::get_keyval(conf, "step", state_file_step,
                           cvm::step_absolute(), colvarparse::parse_silent);
 
@@ -396,7 +435,7 @@ std::ostream & colvarbias::write_state(std::ostream &os)
   }
   os.setf(std::ios::scientific, std::ios::floatfield);
   os.precision(cvm::cv_prec);
-  os << bias_type << " {\n"
+  os << state_keyword << " {\n"
      << "  configuration {\n";
   std::istringstream is(get_state_params());
   std::string line;
@@ -415,13 +454,12 @@ 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) ||
+  if ( !(is >> key)   || !(key == state_keyword || key == bias_type) ||
        !(is >> brace) || !(brace == "{") ||
-       !(is >> colvarparse::read_block("configuration", conf)) ||
+       !(is >> colvarparse::read_block("configuration", &conf)) ||
        (set_state_params(conf) != COLVARS_OK) ) {
-    if (key != bias_type)
+    cvm::error("Error: in reading state configuration for \""+bias_type+
-      cvm::log("Found key \"" + key + "\" instead of \"" + bias_type + "\"\n");
+               "\" bias \""+
-    cvm::error("Error: in reading state configuration for \""+bias_type+"\" bias \""+
                this->name+"\" at position "+
                cvm::to_str(static_cast<size_t>(is.tellg()))+
                " in stream.\n", INPUT_ERROR);
@@ -431,6 +469,15 @@ std::istream & colvarbias::read_state(std::istream &is)
     return is;
   }
 
+  if (matching_state == false) {
+    // This state is not for this bias
+    is.seekg(start_pos, std::ios::beg);
+    return is;
+  }
+
+  cvm::log("Restarting "+bias_type+" bias \""+name+"\" from step number "+
+           cvm::to_str(state_file_step)+".\n");
+
   if (!read_state_data(is)) {
     cvm::error("Error: in reading state data for \""+bias_type+"\" bias \""+
                this->name+"\" at position "+

lib/colvars/colvarbias.h

@@ -23,9 +23,12 @@ public:
   /// Name of this bias
   std::string name;
 
-  /// Type of this bias
+  /// Keyword indicating the type of this bias
   std::string bias_type;
 
+  /// Keyword used in state files (== bias_type most of the time)
+  std::string state_keyword;
+
   /// If there is more than one bias of this type, record its rank
   int rank;
 
@@ -65,7 +68,7 @@ public:
   virtual int calc_forces(std::vector<colvarvalue> const *values);
 
   /// Send forces to the collective variables
-  virtual void communicate_forces();
+  void communicate_forces();
 
   /// Carry out operations needed before next step is run
   virtual int end_of_step();
@@ -144,10 +147,10 @@ public:
   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);
+  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);
+  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);
@@ -207,6 +210,9 @@ protected:
   /// through each colvar object
   std::vector<colvar *>    colvars;
 
+  /// \brief Up to date value of each colvar
+  std::vector<colvarvalue> colvar_values;
+
   /// \brief Current forces from this bias to the variables
   std::vector<colvarvalue> colvar_forces;
 
@@ -226,6 +232,9 @@ protected:
   /// \brief Step number read from the last state file
   cvm::step_number         state_file_step;
 
+  /// Flag used to tell if the state string being read is for this bias
+  bool matching_state;
+
 };
 
 

lib/colvars/colvarbias_abf.cpp

@@ -115,7 +115,7 @@ int colvarbias_abf::init(std::string const &conf)
     if (colvars[i]->value().type() != colvarvalue::type_scalar) {
       cvm::error("Error: ABF bias can only use scalar-type variables.\n");
     }
-    colvars[i]->enable(f_cv_grid);
+    colvars[i]->enable(f_cv_grid); // Could be a child dependency of a f_cvb_use_grids feature
     if (hide_Jacobian) {
       colvars[i]->enable(f_cv_hide_Jacobian);
     }

lib/colvars/colvarbias_alb.cpp

@@ -246,6 +246,12 @@ int colvarbias_alb::update()
 
 int colvarbias_alb::set_state_params(std::string const &conf)
 {
+  int error_code = colvarbias::set_state_params(conf);
+
+  if (error_code != COLVARS_OK) {
+    return error_code;
+  }
+
   if (!get_keyval(conf, "setCoupling", set_coupling))
     cvm::fatal_error("Error: current setCoupling  is missing from the restart.\n");
 

lib/colvars/colvarbias_histogram.cpp

@@ -17,6 +17,7 @@ colvarbias_histogram::colvarbias_histogram(char const *key)
   : colvarbias(key),
     grid(NULL), out_name("")
 {
+  provide(f_cvb_bypass_ext_lagrangian); // Allow histograms of actual cv for extended-Lagrangian
 }
 
 
@@ -84,6 +85,10 @@ int colvarbias_histogram::init(std::string const &conf)
   grid = new colvar_grid_scalar();
   grid->init_from_colvars(colvars);
 
+  if (is_enabled(f_cvb_bypass_ext_lagrangian)) {
+    grid->request_actual_value();
+  }
+
   {
     std::string grid_conf;
     if (key_lookup(conf, "histogramGrid", &grid_conf)) {

lib/colvars/colvarbias_meta.cpp

@@ -38,7 +38,12 @@ colvarbias_meta::colvarbias_meta(char const *key)
   new_hills_begin = hills.end();
   hills_traj_os = NULL;
 
+  hill_width = 0.0;
+
+  new_hill_freq = 1000;
+
   use_grids = true;
+  grids_freq = 0;
   rebin_grids = false;
   hills_energy = NULL;
   hills_energy_gradients = NULL;
@@ -48,6 +53,8 @@ colvarbias_meta::colvarbias_meta(char const *key)
   dump_fes_save = false;
   dump_replica_fes = false;
 
+  b_hills_traj = false;
+
   ebmeta_equil_steps = 0L;
 
   replica_update_freq = 0;
@@ -58,6 +65,7 @@ colvarbias_meta::colvarbias_meta(char const *key)
 int colvarbias_meta::init(std::string const &conf)
 {
   int error_code = COLVARS_OK;
+  size_t i = 0;
 
   error_code |= colvarbias::init(conf);
   error_code |= colvarbias_ti::init(conf);
@@ -71,16 +79,38 @@ int colvarbias_meta::init(std::string const &conf)
     cvm::error("Error: hillWeight must be provided, and a positive number.\n", INPUT_ERROR);
   }
 
-  get_keyval(conf, "newHillFrequency", new_hill_freq, 1000);
+  get_keyval(conf, "newHillFrequency", new_hill_freq, new_hill_freq);
   if (new_hill_freq > 0) {
     enable(f_cvb_history_dependent);
+    if (grids_freq == 0) {
+      grids_freq = new_hill_freq;
+    }
   }
 
-  get_keyval(conf, "hillWidth", hill_width, cvm::sqrt(2.0 * PI) / 2.0);
+  get_keyval(conf, "gaussianSigmas", colvar_sigmas, colvar_sigmas);
-  cvm::log("Half-widths of the Gaussian hills (sigma's):\n");
+
-  for (size_t i = 0; i < num_variables(); i++) {
+  get_keyval(conf, "hillWidth", hill_width, hill_width);
-    cvm::log(variables(i)->name+std::string(": ")+
+
-             cvm::to_str(0.5 * variables(i)->width * hill_width));
+  if ((colvar_sigmas.size() > 0) && (hill_width > 0.0)) {
+    error_code |= cvm::error("Error: hillWidth and gaussianSigmas are "
+                             "mutually exclusive.", INPUT_ERROR);
+  }
+
+  if (hill_width > 0.0) {
+    colvar_sigmas.resize(num_variables());
+    // Print the calculated sigma parameters
+    cvm::log("Half-widths of the Gaussian hills (sigma's):\n");
+    for (i = 0; i < num_variables(); i++) {
+      colvar_sigmas[i] = variables(i)->width * hill_width / 2.0;
+      cvm::log(variables(i)->name+std::string(": ")+
+               cvm::to_str(colvar_sigmas[i]));
+    }
+  }
+
+  if (colvar_sigmas.size() == 0) {
+    error_code |= cvm::error("Error: positive values are required for "
+                             "either hillWidth or gaussianSigmas.",
+                             INPUT_ERROR);
   }
 
   {
@@ -95,11 +125,18 @@ int colvarbias_meta::init(std::string const &conf)
   get_keyval(conf, "useGrids", use_grids, use_grids);
 
   if (use_grids) {
-    get_keyval(conf, "gridsUpdateFrequency", grids_freq, new_hill_freq);
+
+    for (i = 0; i < num_variables(); i++) {
+      if (2.0*colvar_sigmas[i] < variables(i)->width) {
+        cvm::log("Warning: gaussianSigmas is too narrow for the grid "
+                 "spacing along "+variables(i)->name+".");
+      }
+    }
+
+    get_keyval(conf, "gridsUpdateFrequency", grids_freq, grids_freq);
     get_keyval(conf, "rebinGrids", rebin_grids, rebin_grids);
 
     expand_grids = false;
-    size_t i;
     for (i = 0; i < num_variables(); i++) {
       variables(i)->enable(f_cv_grid); // Could be a child dependency of a f_cvb_use_grids feature
       if (variables(i)->expand_boundaries) {
@@ -116,15 +153,17 @@ int colvarbias_meta::init(std::string const &conf)
     get_keyval(conf, "keepHills", keep_hills, keep_hills);
     get_keyval(conf, "keepFreeEnergyFiles", dump_fes_save, dump_fes_save);
 
-    hills_energy           = new colvar_grid_scalar(colvars);
+    if (hills_energy == NULL) {
-    hills_energy_gradients = new colvar_grid_gradient(colvars);
+      hills_energy           = new colvar_grid_scalar(colvars);
+      hills_energy_gradients = new colvar_grid_gradient(colvars);
+    }
 
   } else {
 
     dump_fes = false;
   }
 
-  get_keyval(conf, "writeHillsTrajectory", b_hills_traj, false);
+  get_keyval(conf, "writeHillsTrajectory", b_hills_traj, b_hills_traj);
 
   error_code |= init_replicas_params(conf);
   error_code |= init_well_tempered_params(conf);
@@ -323,7 +362,7 @@ int colvarbias_meta::clear_state_data()
 // **********************************************************************
 
 std::list<colvarbias_meta::hill>::const_iterator
-colvarbias_meta::create_hill(colvarbias_meta::hill const &h)
+colvarbias_meta::add_hill(colvarbias_meta::hill const &h)
 {
   hill_iter const hills_end = hills.end();
   hills.push_back(h);
@@ -563,12 +602,14 @@ int colvarbias_meta::update_bias()
 
     case single_replica:
 
-      create_hill(hill(hill_weight*hills_scale, colvars, hill_width));
+      add_hill(hill(cvm::step_absolute(), hill_weight*hills_scale,
+                    colvar_values, colvar_sigmas));
 
       break;
 
     case multiple_replicas:
-      create_hill(hill(hill_weight*hills_scale, colvars, hill_width, replica_id));
+      add_hill(hill(cvm::step_absolute(), hill_weight*hills_scale,
+                    colvar_values, colvar_sigmas, replica_id));
       std::ostream *replica_hills_os =
         cvm::proxy->get_output_stream(replica_hills_file);
       if (replica_hills_os) {
@@ -730,35 +771,21 @@ void colvarbias_meta::calc_hills(colvarbias_meta::hill_iter      h_first,
                                  std::vector<colvarvalue> const *values)
 {
   size_t i = 0;
-  std::vector<colvarvalue> curr_values(num_variables());
-  for (i = 0; i < num_variables(); i++) {
-    curr_values[i].type(variables(i)->value());
-  }
-
-  if (values) {
-    for (i = 0; i < num_variables(); i++) {
-      curr_values[i] = (*values)[i];
-    }
-  } else {
-    for (i = 0; i < num_variables(); i++) {
-      curr_values[i] = variables(i)->value();
-    }
-  }
 
   for (hill_iter h = h_first; h != h_last; h++) {
 
     // compute the gaussian exponent
     cvm::real cv_sqdev = 0.0;
     for (i = 0; i < num_variables(); i++) {
-      colvarvalue const &x  = curr_values[i];
+      colvarvalue const &x  = values ? (*values)[i] : colvar_values[i];
       colvarvalue const &center = h->centers[i];
-      cvm::real const    half_width = 0.5 * h->widths[i];
+      cvm::real const sigma = h->sigmas[i];
-      cv_sqdev += (variables(i)->dist2(x, center)) / (half_width*half_width);
+      cv_sqdev += (variables(i)->dist2(x, center)) / (sigma*sigma);
     }
 
     // compute the gaussian
     if (cv_sqdev > 23.0) {
-      // set it to zero if the exponent is more negative than log(1.0E-05)
+      // set it to zero if the exponent is more negative than log(1.0E-06)
       h->value(0.0);
     } else {
       h->value(cvm::exp(-0.5*cv_sqdev));
@@ -775,22 +802,22 @@ void colvarbias_meta::calc_hills_force(size_t const &i,
                                        std::vector<colvarvalue> const *values)
 {
   // Retrieve the value of the colvar
-  colvarvalue const x(values ? (*values)[i] : variables(i)->value());
+  colvarvalue const x(values ? (*values)[i] : colvar_values[i]);
 
   // do the type check only once (all colvarvalues in the hills series
   // were already saved with their types matching those in the
   // colvars)
 
   hill_iter h;
-  switch (variables(i)->value().type()) {
+  switch (x.type()) {
 
   case colvarvalue::type_scalar:
     for (h = h_first; h != h_last; h++) {
       if (h->value() == 0.0) continue;
       colvarvalue const &center = h->centers[i];
-      cvm::real const    half_width = 0.5 * h->widths[i];
+      cvm::real const sigma = h->sigmas[i];
       forces[i].real_value +=
-        ( h->weight() * h->value() * (0.5 / (half_width*half_width)) *
+        ( h->weight() * h->value() * (0.5 / (sigma*sigma)) *
           (variables(i)->dist2_lgrad(x, center)).real_value );
     }
     break;
@@ -801,9 +828,9 @@ void colvarbias_meta::calc_hills_force(size_t const &i,
     for (h = h_first; h != h_last; h++) {
       if (h->value() == 0.0) continue;
       colvarvalue const &center = h->centers[i];
-      cvm::real const    half_width = 0.5 * h->widths[i];
+      cvm::real const sigma = h->sigmas[i];
       forces[i].rvector_value +=
-        ( h->weight() * h->value() * (0.5 / (half_width*half_width)) *
+        ( h->weight() * h->value() * (0.5 / (sigma*sigma)) *
           (variables(i)->dist2_lgrad(x, center)).rvector_value );
     }
     break;
@@ -813,9 +840,9 @@ void colvarbias_meta::calc_hills_force(size_t const &i,
     for (h = h_first; h != h_last; h++) {
       if (h->value() == 0.0) continue;
       colvarvalue const &center = h->centers[i];
-      cvm::real const    half_width = 0.5 * h->widths[i];
+      cvm::real const sigma = h->sigmas[i];
       forces[i].quaternion_value +=
-        ( h->weight() * h->value() * (0.5 / (half_width*half_width)) *
+        ( h->weight() * h->value() * (0.5 / (sigma*sigma)) *
           (variables(i)->dist2_lgrad(x, center)).quaternion_value );
     }
     break;
@@ -824,9 +851,9 @@ void colvarbias_meta::calc_hills_force(size_t const &i,
     for (h = h_first; h != h_last; h++) {
       if (h->value() == 0.0) continue;
       colvarvalue const &center = h->centers[i];
-      cvm::real const    half_width = 0.5 * h->widths[i];
+      cvm::real const sigma = h->sigmas[i];
       forces[i].vector1d_value +=
-        ( h->weight() * h->value() * (0.5 / (half_width*half_width)) *
+        ( h->weight() * h->value() * (0.5 / (sigma*sigma)) *
           (variables(i)->dist2_lgrad(x, center)).vector1d_value );
     }
     break;
@@ -1246,14 +1273,20 @@ void colvarbias_meta::read_replica_files()
 
 int colvarbias_meta::set_state_params(std::string const &state_conf)
 {
-  std::string new_replica = "";
+  int error_code = colvarbias::set_state_params(state_conf);
-  if (colvarparse::get_keyval(state_conf, "replicaID", new_replica,
+
+  if (error_code != COLVARS_OK) {
+    return error_code;
+  }
+
+  std::string check_replica = "";
+  if (colvarparse::get_keyval(state_conf, "replicaID", check_replica,
                               std::string(""), colvarparse::parse_silent) &&
-      (new_replica != this->replica_id)) {
+      (check_replica != this->replica_id)) {
-    cvm::error("Error: in the state file, the "
+    return cvm::error("Error: in the state file , the "
-               "\"metadynamics\" block has a different replicaID: different system?\n",
+                      "\"metadynamics\" block has a different replicaID ("+
-               INPUT_ERROR);
+                      check_replica+" instead of "+replica_id+").\n",
-    return INPUT_ERROR;
+                      INPUT_ERROR);
   }
 
   return COLVARS_OK;
@@ -1479,17 +1512,18 @@ std::istream & colvarbias_meta::read_hill(std::istream &is)
   if (!is) return is; // do nothing if failbit is set
 
   size_t const start_pos = is.tellg();
+  size_t i = 0;
 
   std::string data;
-  if ( !(is >> read_block("hill", data)) ) {
+  if ( !(is >> read_block("hill", &data)) ) {
     is.clear();
     is.seekg(start_pos, std::ios::beg);
     is.setstate(std::ios::failbit);
     return is;
   }
 
-  cvm::step_number h_it;
+  cvm::step_number h_it = 0L;
-  get_keyval(data, "step", h_it, 0L, parse_silent);
+  get_keyval(data, "step", h_it, h_it, parse_restart);
   if (h_it <= state_file_step) {
     if (cvm::debug())
       cvm::log("Skipping a hill older than the state file for metadynamics bias \""+
@@ -1499,31 +1533,24 @@ std::istream & colvarbias_meta::read_hill(std::istream &is)
   }
 
   cvm::real h_weight;
-  get_keyval(data, "weight", h_weight, hill_weight, parse_silent);
+  get_keyval(data, "weight", h_weight, hill_weight, parse_restart);
 
   std::vector<colvarvalue> h_centers(num_variables());
-  for (size_t i = 0; i < num_variables(); i++) {
+  for (i = 0; i < num_variables(); i++) {
     h_centers[i].type(variables(i)->value());
   }
-  {
+  get_keyval(data, "centers", h_centers, h_centers, parse_restart);
-    // it is safer to read colvarvalue objects one at a time;
-    // TODO: change this it later
-    std::string centers_input;
-    key_lookup(data, "centers", &centers_input);
-    std::istringstream centers_is(centers_input);
-    for (size_t i = 0; i < num_variables(); i++) {
-      centers_is >> h_centers[i];
-    }
-  }
 
-  std::vector<cvm::real> h_widths(num_variables());
+  std::vector<cvm::real> h_sigmas(num_variables());
-  get_keyval(data, "widths", h_widths,
+  get_keyval(data, "widths", h_sigmas, h_sigmas, parse_restart);
-             std::vector<cvm::real>(num_variables(), (cvm::sqrt(2.0 * PI) / 2.0)),
+  for (i = 0; i < num_variables(); i++) {
-             parse_silent);
+    // For backward compatibility, read the widths instead of the sigmas
+    h_sigmas[i] /= 2.0;
+  }
 
   std::string h_replica = "";
   if (comm != single_replica) {
-    get_keyval(data, "replicaID", h_replica, replica_id, parse_silent);
+    get_keyval(data, "replicaID", h_replica, replica_id, parse_restart);
     if (h_replica != replica_id)
       cvm::fatal_error("Error: trying to read a hill created by replica \""+h_replica+
                        "\" for replica \""+replica_id+
@@ -1531,7 +1558,7 @@ std::istream & colvarbias_meta::read_hill(std::istream &is)
   }
 
   hill_iter const hills_end = hills.end();
-  hills.push_back(hill(h_it, h_weight, h_centers, h_widths, h_replica));
+  hills.push_back(hill(h_it, h_weight, h_centers, h_sigmas, h_replica));
   if (new_hills_begin == hills_end) {
     // if new_hills_begin is unset, set it for the first time
     new_hills_begin = hills.end();
@@ -1753,8 +1780,8 @@ void colvarbias_meta::write_pmf()
     if (ebmeta) {
       int nt_points=pmf->number_of_points();
       for (int i = 0; i < nt_points; i++) {
-         cvm:: real pmf_val=0.0;
+         cvm::real pmf_val=0.0;
-         cvm:: real target_val=target_dist->value(i);
+         cvm::real target_val=target_dist->value(i);
          if (target_val>0) {
            pmf_val=pmf->value(i);
            pmf_val=pmf_val+cvm::temperature() * cvm::boltzmann() * cvm::logn(target_val);
@@ -1794,8 +1821,8 @@ void colvarbias_meta::write_pmf()
     if (ebmeta) {
       int nt_points=pmf->number_of_points();
       for (int i = 0; i < nt_points; i++) {
-         cvm:: real pmf_val=0.0;
+         cvm::real pmf_val=0.0;
-         cvm:: real target_val=target_dist->value(i);
+         cvm::real target_val=target_dist->value(i);
          if (target_val>0) {
            pmf_val=pmf->value(i);
            pmf_val=pmf_val+cvm::temperature() * cvm::boltzmann() * cvm::logn(target_val);
@@ -1889,10 +1916,10 @@ std::string colvarbias_meta::hill::output_traj()
   }
 
   os << "  ";
-  for (i = 0; i < widths.size(); i++) {
+  for (i = 0; i < sigmas.size(); i++) {
     os << " ";
     os << std::setprecision(cvm::cv_prec)
-       << std::setw(cvm::cv_width) << widths[i];
+       << std::setw(cvm::cv_width) << sigmas[i];
   }
 
   os << "  ";
@@ -1903,6 +1930,47 @@ std::string colvarbias_meta::hill::output_traj()
 }
 
 
+colvarbias_meta::hill::hill(cvm::step_number it_in,
+                            cvm::real W_in,
+                            std::vector<colvarvalue> const &cv_values,
+                            std::vector<cvm::real> const &cv_sigmas,
+                            std::string const &replica_in)
+  : it(it_in),
+    sW(1.0),
+    W(W_in),
+    centers(cv_values.size()),
+    sigmas(cv_values.size()),
+    replica(replica_in)
+{
+  for (size_t i = 0; i < cv_values.size(); i++) {
+    centers[i].type(cv_values[i]);
+    centers[i] = cv_values[i];
+    sigmas[i] = cv_sigmas[i];
+  }
+  if (cvm::debug()) {
+    cvm::log("New hill, applied to "+cvm::to_str(cv_values.size())+
+             " collective variables, with centers "+
+             cvm::to_str(centers)+", sigmas "+
+             cvm::to_str(sigmas)+" and weight "+
+             cvm::to_str(W)+".\n");
+  }
+}
+
+
+colvarbias_meta::hill::hill(colvarbias_meta::hill const &h)
+  : sW(1.0),
+    W(h.W),
+    centers(h.centers),
+    sigmas(h.sigmas),
+    it(h.it),
+    replica(h.replica)
+{}
+
+
+colvarbias_meta::hill::~hill()
+{}
+
+
 std::ostream & operator << (std::ostream &os, colvarbias_meta::hill const &h)
 {
   os.setf(std::ios::scientific, std::ios::floatfield);
@@ -1927,12 +1995,13 @@ std::ostream & operator << (std::ostream &os, colvarbias_meta::hill const &h)
   }
   os << "\n";
 
+  // For backward compatibility, write the widths instead of the sigmas
   os << "  widths  ";
-  for (i = 0; i < (h.widths).size(); i++) {
+  for (i = 0; i < (h.sigmas).size(); i++) {
     os << " "
        << std::setprecision(cvm::cv_prec)
        << std::setw(cvm::cv_width)
-       << h.widths[i];
+       << 2.0 * h.sigmas[i];
   }
   os << "\n";
 

lib/colvars/colvarbias_meta.h

@@ -71,11 +71,14 @@ public:
 
 protected:
 
-  /// \brief width of a hill
+  /// Width of a hill in number of grid points
   ///
   /// The local width of each collective variable, multiplied by this
   /// number, provides the hill width along that direction
-  cvm::real  hill_width;
+  cvm::real hill_width;
+
+  /// The sigma parameters of the Gaussian hills
+  std::vector<cvm::real> colvar_sigmas;
 
   /// \brief Number of simulation steps between two hills
   size_t     new_hill_freq;
@@ -113,11 +116,11 @@ protected:
   /// \brief Add a new hill; if a .hills trajectory is written,
   /// write it there; if there is more than one replica, communicate
   /// it to the others
-  virtual std::list<hill>::const_iterator create_hill(hill const &h);
+  std::list<hill>::const_iterator add_hill(hill const &h);
 
   /// \brief Remove a previously saved hill (returns an iterator for
   /// the next hill in the list)
-  virtual std::list<hill>::const_iterator delete_hill(hill_iter &h);
+  std::list<hill>::const_iterator delete_hill(hill_iter &h);
 
   /// \brief Calculate the values of the hills, incrementing
   /// bias_energy
@@ -266,6 +269,9 @@ class colvarbias_meta::hill {
 
 protected:
 
+  /// Time step at which this hill was added
+  cvm::step_number it;
+
   /// Value of the hill function (ranges between 0 and 1)
   cvm::real hill_value;
 
@@ -275,83 +281,35 @@ protected:
   /// Maximum height in energy of the hill
   cvm::real W;
 
-  /// Center of the hill in the collective variable space
+  /// Centers of the hill in the collective variable space
-  std::vector<colvarvalue>  centers;
+  std::vector<colvarvalue> centers;
 
-  /// Widths of the hill in the collective variable space
+  /// Half-widths of the hill in the collective variable space
-  std::vector<cvm::real>    widths;
+  std::vector<cvm::real> sigmas;
+
+  /// Identity of the replica who added this hill
+  std::string replica;
 
 public:
 
   friend class colvarbias_meta;
 
-  /// Time step at which this hill was added
+  /// Constructor of a hill object
-  cvm::step_number it;
+  /// \param it Step number at which the hill was added
-
+  /// \param W Weight of the hill (energy units)
-  /// Identity of the replica who added this hill (only in multiple replica simulations)
+  /// \param cv_values Array of collective variable values
-  std::string replica;
+  /// \param cv_sigmas Array of collective variable values
-
+  /// \param replica ID of the replica that creates the hill (optional)
-  /// \brief Runtime constructor: data are read directly from
+  hill(cvm::step_number it, cvm::real W,
-  /// collective variables \param weight Weight of the hill \param
+       std::vector<colvarvalue> const &cv_values,
-  /// cv Pointer to the array of collective variables involved \param
+       std::vector<cvm::real> const &cv_sigmas,
-  /// replica (optional) Identity of the replica which creates the
+       std::string const &replica = "");
-  /// hill
-  inline hill(cvm::real             const &W_in,
-              std::vector<colvar *>       &cv,
-              cvm::real             const &hill_width,
-              std::string           const &replica_in = "")
-    : sW(1.0),
-      W(W_in),
-      centers(cv.size()),
-      widths(cv.size()),
-      it(cvm::step_absolute()),
-      replica(replica_in)
-  {
-    for (size_t i = 0; i < cv.size(); i++) {
-      centers[i].type(cv[i]->value());
-      centers[i] = cv[i]->value();
-      widths[i] = cv[i]->width * hill_width;
-    }
-    if (cvm::debug())
-      cvm::log("New hill, applied to "+cvm::to_str(cv.size())+
-                " collective variables, with centers "+
-                cvm::to_str(centers)+", widths "+
-                cvm::to_str(widths)+" and weight "+
-                cvm::to_str(W)+".\n");
-  }
-
-  /// \brief General constructor: all data are explicitly passed as
-  /// arguments (used for instance when reading hills saved on a
-  /// file) \param it Time step of creation of the hill \param
-  /// weight Weight of the hill \param centers Center of the hill
-  /// \param widths Width of the hill around centers \param replica
-  /// (optional) Identity of the replica which creates the hill
-  inline hill(cvm::step_number          const &it_in,
-              cvm::real                 const &W_in,
-              std::vector<colvarvalue>  const &centers_in,
-              std::vector<cvm::real>    const &widths_in,
-              std::string               const &replica_in = "")
-    : sW(1.0),
-      W(W_in),
-      centers(centers_in),
-      widths(widths_in),
-      it(it_in),
-      replica(replica_in)
-  {}
 
   /// Copy constructor
-  inline hill(colvarbias_meta::hill const &h)
+  hill(colvarbias_meta::hill const &h);
-    : sW(1.0),
-      W(h.W),
-      centers(h.centers),
-      widths(h.widths),
-      it(h.it),
-      replica(h.replica)
-  {}
 
   /// Destructor
-  inline ~hill()
+  ~hill();
-  {}
 
   /// Get the energy
   inline cvm::real energy()
@@ -440,8 +398,7 @@ public:
   std::string output_traj();
 
   /// Write the hill to an output stream
-  inline friend std::ostream & operator << (std::ostream &os,
+  friend std::ostream & operator << (std::ostream &os, hill const &h);
-                                            hill const &h);
 
 };
 

lib/colvars/colvarbias_restraint.cpp

@@ -17,6 +17,7 @@
 colvarbias_restraint::colvarbias_restraint(char const *key)
   : colvarbias(key), colvarbias_ti(key)
 {
+  state_keyword = "restraint";
 }
 
 
@@ -705,67 +706,6 @@ std::ostream & colvarbias_restraint_k_moving::write_traj(std::ostream &os)
 
 
 
-// redefined due to legacy state file keyword "harmonic"
-std::istream & colvarbias_restraint::read_state(std::istream &is)
-{
-  size_t const start_pos = is.tellg();
-
-  std::string key, brace, conf;
-  if ( !(is >> key)   || !(key == "restraint" || key == "harmonic") ||
-       !(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(static_cast<size_t>(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(static_cast<size_t>(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::log("brace = "+brace+"\n");
-    cvm::error("Error: corrupt restart information for \""+bias_type+"\" bias \""+
-               this->name+"\": no matching brace at position "+
-               cvm::to_str(static_cast<size_t>(is.tellg()))+" in stream.\n");
-    is.setstate(std::ios::failbit);
-  }
-
-  return is;
-}
-
-
-std::ostream & colvarbias_restraint::write_state(std::ostream &os)
-{
-  os.setf(std::ios::scientific, std::ios::floatfield);
-  os << "restraint {\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;
-}
-
-
-
 colvarbias_restraint_harmonic::colvarbias_restraint_harmonic(char const *key)
   : colvarbias(key),
     colvarbias_ti(key),
@@ -928,6 +868,9 @@ colvarbias_restraint_harmonic_walls::colvarbias_restraint_harmonic_walls(char co
 {
   lower_wall_k = -1.0;
   upper_wall_k = -1.0;
+  // This bias implements the bias_actual_colvars feature (most others do not)
+  provide(f_cvb_bypass_ext_lagrangian);
+  set_enabled(f_cvb_bypass_ext_lagrangian); // Defaults to enabled
 }
 
 
@@ -1075,23 +1018,13 @@ int colvarbias_restraint_harmonic_walls::update()
 }
 
 
-void colvarbias_restraint_harmonic_walls::communicate_forces()
-{
-  for (size_t i = 0; i < num_variables(); i++) {
-    if (cvm::debug()) {
-      cvm::log("Communicating a force to colvar \""+
-               variables(i)->name+"\".\n");
-    }
-    // Impulse-style multiple timestep
-    variables(i)->add_bias_force_actual_value(cvm::real(time_step_factor) * colvar_forces[i]);
-  }
-}
-
-
 cvm::real colvarbias_restraint_harmonic_walls::colvar_distance(size_t i) const
 {
   colvar *cv = variables(i);
-  colvarvalue const &cvv = variables(i)->actual_value();
+
+  colvarvalue const &cvv = is_enabled(f_cvb_bypass_ext_lagrangian) ?
+    variables(i)->actual_value() :
+    variables(i)->value();
 
   // For a periodic colvar, both walls may be applicable at the same time
   // in which case we pick the closer one

lib/colvars/colvarbias_restraint.h

@@ -35,8 +35,6 @@ public:
   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);
 
   virtual std::ostream & write_traj_label(std::ostream &os);
   virtual std::ostream & write_traj(std::ostream &os);
@@ -258,7 +256,6 @@ 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_state_data(std::ostream &os);

lib/colvars/colvardeps.cpp

@@ -428,11 +428,11 @@ void colvardeps::require_feature_alt(int f, int g, int h, int i, int j) {
 
 void colvardeps::print_state() {
   size_t i;
-  cvm::log("Enabled features of \"" + description + "\" (with reference count)");
+  cvm::log("Features of \"" + description + "\" ON/OFF (refcount)");
   for (i = 0; i < feature_states.size(); i++) {
-    if (is_enabled(i))
+    std::string onoff = is_enabled(i) ? "ON" : "OFF";
-      cvm::log("- " + features()[i]->description + " ("
+    cvm::log("- " + features()[i]->description + "   " + onoff + " ("
-        + cvm::to_str(feature_states[i].ref_count) + ")");
+      + cvm::to_str(feature_states[i].ref_count) + ")");
   }
   cvm::increase_depth();
   for (i=0; i<children.size(); i++) {

lib/colvars/colvardeps.h

@@ -229,6 +229,8 @@ public:
     f_cvb_awake,
     /// \brief will apply forces
     f_cvb_apply_force,
+    /// \brief force this bias to act on actual value for extended-Lagrangian coordinates
+    f_cvb_bypass_ext_lagrangian,
     /// \brief requires total forces
     f_cvb_get_total_force,
     /// \brief whether this bias should record the accumulated work

lib/colvars/colvargrid.h

@@ -932,7 +932,7 @@ public:
     std::string key, conf;
     if ((is >> key) && (key == std::string("grid_parameters"))) {
       is.seekg(start_pos, std::ios::beg);
-      is >> colvarparse::read_block("grid_parameters", conf);
+      is >> colvarparse::read_block("grid_parameters", &conf);
       parse_params(conf, colvarparse::parse_silent);
     } else {
       cvm::log("Grid parameters are missing in the restart file, using those from the configuration.\n");

lib/colvars/colvarmodule.cpp

@@ -1284,7 +1284,7 @@ std::istream & colvarmodule::read_restart(std::istream &is)
   {
     // read global restart information
     std::string restart_conf;
-    if (is >> colvarparse::read_block("configuration", restart_conf)) {
+    if (is >> colvarparse::read_block("configuration", &restart_conf)) {
 
       parse->get_keyval(restart_conf, "step",
                         it_restart, static_cast<step_number>(0),
@@ -1329,30 +1329,84 @@ std::istream & colvarmodule::read_restart(std::istream &is)
     parse->clear_keyword_registry();
   }
 
-  // colvars restart
+  print_total_forces_errning(warn_total_forces);
-  cvm::increase_depth();
+
-  for (std::vector<colvar *>::iterator cvi = colvars.begin();
+  read_objects_state(is);
-       cvi != colvars.end();
+
-       cvi++) {
+  return is;
-    if ( !((*cvi)->read_restart(is)) ) {
+}
-      cvm::error("Error: in reading restart configuration for collective variable \""+
+
-                 (*cvi)->name+"\".\n",
+
-                 INPUT_ERROR);
+
+std::istream & colvarmodule::read_objects_state(std::istream &is)
+{
+  size_t pos = 0;
+  std::string word;
+
+  while (is.good()) {
+    pos = is.tellg();
+    word.clear();
+    is >> word;
+
+    if (word.size()) {
+
+      is.seekg(pos, std::ios::beg);
+
+      if (word == "colvar") {
+
+        cvm::increase_depth();
+        for (std::vector<colvar *>::iterator cvi = colvars.begin();
+             cvi != colvars.end();
+             cvi++) {
+          if ( !((*cvi)->read_state(is)) ) {
+            // Here an error signals that the variable is a match, but the
+            // state is corrupt; otherwise, the variable rewinds is silently
+            cvm::error("Error: in reading restart configuration for "
+                       "collective variable \""+(*cvi)->name+"\".\n",
+                       INPUT_ERROR);
+          }
+          if (static_cast<size_t>(is.tellg()) > pos) break; // found it
+        }
+        cvm::decrease_depth();
+
+      } else {
+
+        cvm::increase_depth();
+        for (std::vector<colvarbias *>::iterator bi = biases.begin();
+             bi != biases.end();
+             bi++) {
+          if (((*bi)->state_keyword != word) && (*bi)->bias_type != word) {
+            // Skip biases with different type; state_keyword is used to
+            // support different versions of the state file format
+            continue;
+          }
+          if (!((*bi)->read_state(is))) {
+            // Same as above, an error means a match but the state is incorrect
+            cvm::error("Error: in reading restart configuration for bias \""+
+                       (*bi)->name+"\".\n",
+                       INPUT_ERROR);
+          }
+          if (static_cast<size_t>(is.tellg()) > pos) break; // found it
+        }
+        cvm::decrease_depth();
+      }
     }
+
+    if (static_cast<size_t>(is.tellg()) == pos) {
+      // This block has not been read by any object: discard it and move on
+      // to the next one
+      is >> colvarparse::read_block(word, NULL);
+    }
+
+    if (!is) break;
   }
 
-  // biases restart
+  return is;
-  for (std::vector<colvarbias *>::iterator bi = biases.begin();
+}
-       bi != biases.end();
-       bi++) {
-    if (!((*bi)->read_state(is))) {
-      cvm::error("Error: in reading restart configuration for bias \""+
-                 (*bi)->name+"\".\n",
-                 INPUT_ERROR);
-    }
-  }
-  cvm::decrease_depth();
 
+
+int colvarmodule::print_total_forces_errning(bool warn_total_forces)
+{
   if (warn_total_forces) {
     cvm::log(cvm::line_marker);
     cvm::log("WARNING: The definition of system forces has changed.  Please see:\n");
@@ -1367,10 +1421,11 @@ to:\n\
 and load it to continue this simulation.\n");
     output_prefix() = output_prefix()+".tmp";
     write_restart_file(output_prefix()+".colvars.state");
-    cvm::error("Exiting with error until issue is addressed.\n", INPUT_ERROR);
+    return cvm::error("Exiting with error until issue is addressed.\n",
+                      INPUT_ERROR);
   }
 
-  return is;
+  return COLVARS_OK;
 }
 
 
@@ -1499,7 +1554,7 @@ std::ostream & colvarmodule::write_restart(std::ostream &os)
   for (std::vector<colvar *>::iterator cvi = colvars.begin();
        cvi != colvars.end();
        cvi++) {
-    (*cvi)->write_restart(os);
+    (*cvi)->write_state(os);
   }
 
   for (std::vector<colvarbias *>::iterator bi = biases.begin();

lib/colvars/colvarmodule.h

@@ -429,8 +429,15 @@ public:
   /// (Re)initialize the output trajectory and state file (does not write it yet)
   int setup_output();
 
-  /// Read the input restart file
+  /// Read a restart file
   std::istream & read_restart(std::istream &is);
+
+  /// Read the states of individual objects; allows for changes
+  std::istream & read_objects_state(std::istream &is);
+
+  /// If needed (old restart file), print the warning that cannot be ignored
+  int print_total_forces_errning(bool warn_total_forces);
+
   /// Write the output restart file
   std::ostream & write_restart(std::ostream &os);
 

lib/colvars/colvarparse.cpp

@@ -840,6 +840,17 @@ bool colvarparse::key_lookup(std::string const &conf,
 }
 
 
+colvarparse::read_block::read_block(std::string const &key_in,
+                                    std::string *data_in)
+  : key(key_in), data(data_in)
+{
+}
+
+
+colvarparse::read_block::~read_block()
+{}
+
+
 std::istream & operator>> (std::istream &is, colvarparse::read_block const &rb)
 {
   size_t start_pos = is.tellg();
@@ -856,7 +867,9 @@ std::istream & operator>> (std::istream &is, colvarparse::read_block const &rb)
   }
 
   if (next != "{") {
-    (*rb.data) = next;
+    if (rb.data) {
+      *(rb.data) = next;
+    }
     return is;
   }
 
@@ -870,9 +883,15 @@ std::istream & operator>> (std::istream &is, colvarparse::read_block const &rb)
       br_old = br;
       br++;
     }
-    if (brace_count) (*rb.data).append(line + "\n");
+    if (brace_count) {
+      if (rb.data) {
+        (rb.data)->append(line + "\n");
+      }
+    }
     else {
-      (*rb.data).append(line, 0, br_old);
+      if (rb.data) {
+        (rb.data)->append(line, 0, br_old);
+      }
       break;
     }
   }

lib/colvars/colvarparse.h

@@ -272,14 +272,18 @@ public:
   /// skipping other blocks
   class read_block {
 
-    std::string const   key;
+    /// The keyword that identifies the block
+    std::string const key;
+
+    /// Where to keep the data (may be NULL)
     std::string * const data;
 
   public:
-    inline read_block(std::string const &key_in, std::string &data_in)
+
-      : key(key_in), data(&data_in)
+    read_block(std::string const &key_in, std::string *data_in = NULL);
-    {}
+
-    inline ~read_block() {}
+    ~read_block();
+
     friend std::istream & operator >> (std::istream &is, read_block const &rb);
   };
 

lib/colvars/colvars_version.h

@@ -1,3 +1,3 @@
 #ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2020-01-27"
+#define COLVARS_VERSION "2020-02-25"
 #endif

src/USER-COLVARS/colvarproxy_lammps.cpp

@@ -160,14 +160,20 @@ void colvarproxy_lammps::init(const char *conf_file)
   }
 }
 
-void colvarproxy_lammps::add_config_file(const char *conf_file)
+int colvarproxy_lammps::add_config_file(const char *conf_file)
 {
-  colvars->read_config_file(conf_file);
+  return colvars->read_config_file(conf_file);
 }
 
-void colvarproxy_lammps::add_config_string(const std::string &conf)
+int colvarproxy_lammps::add_config_string(const std::string &conf)
 {
-  colvars->read_config_string(conf);
+  return colvars->read_config_string(conf);
+}
+
+int colvarproxy_lammps::read_state_file(char const *state_filename)
+{
+  input_prefix() = std::string(state_filename);
+  return colvars->setup_input();
 }
 
 colvarproxy_lammps::~colvarproxy_lammps()

src/USER-COLVARS/colvarproxy_lammps.h

@@ -86,10 +86,13 @@ class colvarproxy_lammps : public colvarproxy {
   void write_output_files();
 
   // read additional config from file
-  void add_config_file(char const *config_filename);
+  int add_config_file(char const *config_filename);
 
   // read additional config from string
-  void add_config_string(const std::string &config);
+  int add_config_string(const std::string &config);
+
+  // load a state file
+  int read_state_file(char const *state_filename);
 
   // Request to set the units used internally by Colvars
   int set_unit_system(std::string const &units_in, bool check_only);

src/USER-COLVARS/fix_colvars.cpp

@@ -502,7 +502,7 @@ int FixColvars::modify_param(int narg, char **arg)
     if (me == 0) {
       if (! proxy)
         error->one(FLERR,"Cannot use fix_modify before initialization");
-      proxy->add_config_file(arg[1]);
+      return proxy->add_config_file(arg[1]) == COLVARS_OK ? 0 : 2;
     }
     return 2;
   } else if (strcmp(arg[0],"config") == 0) {
@@ -510,8 +510,16 @@ int FixColvars::modify_param(int narg, char **arg)
     if (me == 0) {
       if (! proxy)
         error->one(FLERR,"Cannot use fix_modify before initialization");
-      std::string conf(arg[1]);
+      std::string const conf(arg[1]);
-      proxy->add_config_string(conf);
+      return proxy->add_config_string(conf) == COLVARS_OK ? 0 : 2;
+    }
+    return 2;
+  } else if (strcmp(arg[0],"load") == 0) {
+    if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
+    if (me == 0) {
+      if (! proxy)
+        error->one(FLERR,"Cannot use fix_modify before initialization");
+      return proxy->read_state_file(arg[1]) == COLVARS_OK ? 0 : 2;
     }
     return 2;
   }