git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@15315 f3b2605a-c512-4ea7-a41b-209d697bcdaa

lib/colvars/colvarbias.cpp

@@ -5,60 +5,95 @@
 #include "colvarbias.h"
 
 
-colvarbias::colvarbias(std::string const &conf, char const *key)
+colvarbias::colvarbias(char const *key)
-  : colvarparse(conf), bias_energy(0.), has_data(false)
+  : bias_type(to_lower_cppstr(key))
 {
-  cvm::log("Initializing a new \""+std::string(key)+"\" instance.\n");
-
   init_cvb_requires();
 
-  size_t rank = 1;
+  rank = 1;
-  std::string const key_str(key);
 
-  if (to_lower_cppstr(key_str) == std::string("abf")) {
+  if (bias_type == std::string("abf")) {
     rank = cvm::n_abf_biases+1;
   }
-  if (to_lower_cppstr(key_str) == std::string("harmonic") ||
+  if (bias_type == std::string("harmonic") ||
-      to_lower_cppstr(key_str) == std::string("linear")) {
+      bias_type == std::string("linear")) {
     rank = cvm::n_rest_biases+1;
   }
-  if (to_lower_cppstr(key_str) == std::string("histogram")) {
+  if (bias_type == std::string("histogram")) {
     rank = cvm::n_histo_biases+1;
   }
-  if (to_lower_cppstr(key_str) == std::string("metadynamics")) {
+  if (bias_type == std::string("metadynamics")) {
     rank = cvm::n_meta_biases+1;
   }
 
-  get_keyval(conf, "name", name, key_str+cvm::to_str(rank));
+  has_data = false;
+  b_output_energy = false;
+  reset();
 
-  if (cvm::bias_by_name(this->name) != NULL) {
-    cvm::error("Error: this bias cannot have the same name, \""+this->name+
-                "\", as another bias.\n", INPUT_ERROR);
-    return;
-  }
-
-  description = "bias " + name;
-
-  // lookup the associated colvars
-  std::vector<std::string> colvars_str;
-  if (get_keyval(conf, "colvars", colvars_str)) {
-    for (size_t i = 0; i < colvars_str.size(); i++) {
-      add_colvar(colvars_str[i]);
-    }
-  }
-  if (!colvars.size()) {
-    cvm::error("Error: no collective variables specified.\n");
-    return;
-  }
-  for (size_t i=0; i<colvars.size(); i++) {
-    // All biases need at least the value of colvars
-    // although possibly not at all timesteps
-    add_child(colvars[i]);
-  }
   // Start in active state by default
   enable(f_cvb_active);
+}
 
-  get_keyval(conf, "outputEnergy", b_output_energy, false);
+
+int colvarbias::init(std::string const &conf)
+{
+  colvarparse::init(conf);
+
+  if (name.size() == 0) {
+    cvm::log("Initializing a new \""+bias_type+"\" instance.\n");
+    get_keyval(conf, "name", name, bias_type+cvm::to_str(rank));
+
+    {
+      colvarbias *bias_with_name = cvm::bias_by_name(this->name);
+      if (bias_with_name != NULL) {
+        if ((bias_with_name->rank != this->rank) ||
+            (bias_with_name->bias_type != this->bias_type)) {
+          cvm::error("Error: this bias cannot have the same name, \""+this->name+
+                     "\", as another bias.\n", INPUT_ERROR);
+          return INPUT_ERROR;
+        }
+      }
+    }
+
+    description = "bias " + name;
+
+    {
+      // lookup the associated colvars
+      std::vector<std::string> colvar_names;
+      if (get_keyval(conf, "colvars", colvar_names)) {
+        if (colvars.size()) {
+          cvm::error("Error: cannot redefine the colvars that a bias was already defined on.\n",
+                     INPUT_ERROR);
+          return INPUT_ERROR;
+        }
+        for (size_t i = 0; i < colvar_names.size(); i++) {
+          add_colvar(colvar_names[i]);
+        }
+      }
+    }
+
+    if (!colvars.size()) {
+      cvm::error("Error: no collective variables specified.\n", INPUT_ERROR);
+      return INPUT_ERROR;
+    }
+
+  } else {
+    cvm::log("Reinitializing bias \""+name+"\".\n");
+  }
+
+  get_keyval(conf, "outputEnergy", b_output_energy, b_output_energy);
+
+  return COLVARS_OK;
+}
+
+
+int colvarbias::reset()
+{
+  bias_energy = 0.0;
+  for (size_t i = 0; i < colvars.size(); i++) {
+    colvar_forces[i].reset();
+  }
+  return COLVARS_OK;
 }
 
 
@@ -66,7 +101,14 @@ colvarbias::colvarbias()
   : colvarparse(), has_data(false)
 {}
 
+
 colvarbias::~colvarbias()
+{
+  colvarbias::clear();
+}
+
+
+int colvarbias::clear()
 {
   // Remove references to this bias from colvars
   for (std::vector<colvar *>::iterator cvi = colvars.begin();
@@ -81,6 +123,7 @@ colvarbias::~colvarbias()
       }
     }
   }
+
   // ...and from the colvars module
   for (std::vector<colvarbias *>::iterator bi = cvm::biases.begin();
        bi != cvm::biases.end();
@@ -90,31 +133,46 @@ colvarbias::~colvarbias()
       break;
     }
   }
+
+  return COLVARS_OK;
 }
 
-void colvarbias::add_colvar(std::string const &cv_name)
+
+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())
+    if (cvm::debug()) {
       cvm::log("Applying this bias to collective variable \""+
-                cv->name+"\".\n");
+               cv->name+"\".\n");
+    }
     colvars.push_back(cv);
+
     colvar_forces.push_back(colvarvalue());
-    colvar_forces.back().type(cv->value()); // make sure each forces is initialized to zero
+    colvar_forces.back().type(cv->value()); // make sure each force is initialized to zero
     colvar_forces.back().reset();
+
     cv->biases.push_back(this); // add back-reference to this bias to colvar
+
+    // Add dependency link.
+    // All biases need at least the value of each colvar
+    // although possibly not at all timesteps
+    add_child(cv);
+
   } else {
     cvm::error("Error: cannot find a colvar named \""+
-               cv_name+"\".\n");
+               cv_name+"\".\n", INPUT_ERROR);
+    return INPUT_ERROR;
   }
+  return COLVARS_OK;
 }
 
 
 int colvarbias::update()
 {
   // Note: if anything is added here, it should be added also in the SMP block of calc_biases()
+  // TODO move here debug msg of bias update
   has_data = true;
   return COLVARS_OK;
 }

lib/colvars/colvarbias.h

@@ -13,10 +13,16 @@ class colvarbias : public colvarparse, public cvm::deps {
 public:
 
   /// Name of this bias
-  std::string    name;
+  std::string name;
+
+  /// Type of this bias
+  std::string bias_type;
+
+  /// If there is more than one bias of this type, record its rank
+  int rank;
 
   /// Add a new collective variable to this bias
-  void add_colvar(std::string const &cv_name);
+  int add_colvar(std::string const &cv_name);
 
   /// Retrieve colvar values and calculate their biasing forces
   /// Return bias energy
@@ -46,11 +52,29 @@ public:
   void communicate_forces();
 
   /// \brief Constructor
-  colvarbias(std::string const &conf, char const *key);
+  colvarbias(char const *key);
+
+  /// \brief Parse config string and (re)initialize
+  virtual int init(std::string const &conf);
+
+  /// \brief Set to zero all mutable data
+  virtual int reset();
+
+protected:
 
   /// Default constructor
   colvarbias();
 
+private:
+
+  /// Copy constructor
+  colvarbias(colvarbias &);
+
+public:
+
+  /// \brief Delete everything
+  virtual int clear();
+
   /// Destructor
   virtual ~colvarbias();
 
@@ -103,7 +127,7 @@ protected:
   bool                     b_output_energy;
 
   /// \brief Whether this bias has already accumulated information
-  /// (when relevant)
+  /// (for history-dependent biases)
   bool                     has_data;
 
 };

lib/colvars/colvarbias_abf.cpp

@@ -4,16 +4,24 @@
 #include "colvar.h"
 #include "colvarbias_abf.h"
 
-/// ABF bias constructor; parses the config file
 
-colvarbias_abf::colvarbias_abf(std::string const &conf, char const *key)
+colvarbias_abf::colvarbias_abf(char const *key)
-  : colvarbias(conf, key),
+  : colvarbias(key),
     force(NULL),
     gradients(NULL),
     samples(NULL),
     last_gradients(NULL),
     last_samples(NULL)
 {
+}
+
+
+int colvarbias_abf::init(std::string const &conf)
+{
+  colvarbias::init(conf);
+
+  provide(f_cvb_history_dependent);
+
   // TODO relax this in case of VMD plugin
   if (cvm::temperature() == 0.0)
     cvm::log("WARNING: ABF should not be run without a thermostat or at 0 Kelvin!\n");
@@ -21,10 +29,18 @@ colvarbias_abf::colvarbias_abf(std::string const &conf, char const *key)
   // ************* parsing general ABF options ***********************
 
   get_keyval(conf, "applyBias",  apply_bias, true);
-  if (!apply_bias) cvm::log("WARNING: ABF biases will *not* be applied!\n");
+  if (apply_bias) {
+    enable(f_cvb_apply_force);
+  } else {
+    cvm::log("WARNING: ABF biases will *not* be applied!\n");
+  }
 
   get_keyval(conf, "updateBias",  update_bias, true);
-  if (!update_bias) cvm::log("WARNING: ABF biases will *not* be updated!\n");
+  if (update_bias) {
+    enable(f_cvb_history_dependent);
+  } else {
+    cvm::log("WARNING: ABF biases will *not* be updated!\n");
+  }
 
   get_keyval(conf, "hideJacobian", hide_Jacobian, false);
   if (hide_Jacobian) {
@@ -38,7 +54,7 @@ colvarbias_abf::colvarbias_abf(std::string const &conf, char const *key)
   min_samples = full_samples / 2;
   // full_samples - min_samples >= 1 is guaranteed
 
-  get_keyval(conf, "inputPrefix",  input_prefix, std::vector<std::string> ());
+  get_keyval(conf, "inputPrefix",  input_prefix, std::vector<std::string>());
   get_keyval(conf, "outputFreq", output_freq, cvm::restart_out_freq);
   get_keyval(conf, "historyFreq", history_freq, 0);
   b_history_files = (history_freq > 0);
@@ -63,10 +79,10 @@ colvarbias_abf::colvarbias_abf(std::string const &conf, char const *key)
 
   if (update_bias) {
   // Request calculation of system force (which also checks for availability)
-    if(enable(f_cvb_get_system_force)) return;
+    if(enable(f_cvb_get_system_force)) return cvm::get_error();
   }
   if (apply_bias) {
-    if(enable(f_cvb_apply_force)) return;
+    if(enable(f_cvb_apply_force)) return cvm::get_error();
   }
 
   for (size_t i = 0; i < colvars.size(); i++) {
@@ -126,6 +142,8 @@ colvarbias_abf::colvarbias_abf(std::string const &conf, char const *key)
   }
 
   cvm::log("Finished ABF setup.\n");
+
+  return COLVARS_OK;
 }
 
 /// Destructor
@@ -277,6 +295,7 @@ int colvarbias_abf::update()
   return COLVARS_OK;
 }
 
+
 int colvarbias_abf::replica_share() {
   int p;
 

lib/colvars/colvarbias_abf.h

@@ -19,10 +19,10 @@ class colvarbias_abf : public colvarbias {
 
 public:
 
-  colvarbias_abf(std::string const &conf, char const *key);
+  colvarbias_abf(char const *key);
-  ~colvarbias_abf();
+  virtual int init(std::string const &conf);
-
+  virtual ~colvarbias_abf();
-  int update();
+  virtual int update();
 
 private:
 

lib/colvars/colvarbias_alb.cpp

@@ -23,8 +23,18 @@ double fmin(double A, double B) { return ( A < B ? A : B ); }
  *
  */
 
-colvarbias_alb::colvarbias_alb(std::string const &conf, char const *key) :
+colvarbias_alb::colvarbias_alb(char const *key)
-  colvarbias(conf, key), update_calls(0), b_equilibration(true) {
+  : colvarbias(key), update_calls(0), b_equilibration(true)
+{
+}
+
+
+int colvarbias_alb::init(std::string const &conf)
+{
+  colvarbias::init(conf);
+
+  provide(f_cvb_history_dependent);
+
   size_t i;
 
   // get the initial restraint centers
@@ -74,6 +84,8 @@ colvarbias_alb::colvarbias_alb(std::string const &conf, char const *key) :
   if (update_freq == 0)
     cvm::fatal_error("Error: must set updateFrequency to greater than 2.\n");
 
+  enable(f_cvb_history_dependent);
+
   get_keyval(conf, "outputCenters", b_output_centers, false);
   get_keyval(conf, "outputGradient", b_output_grad, false);
   get_keyval(conf, "outputCoupling", b_output_coupling, true);
@@ -99,8 +111,6 @@ colvarbias_alb::colvarbias_alb(std::string const &conf, char const *key) :
     }
   }
 
-
-
   if (!get_keyval(conf, "rateMax", max_coupling_rate, max_coupling_rate)) {
     //set to default
     for (i = 0; i < colvars.size(); i++) {
@@ -112,16 +122,19 @@ colvarbias_alb::colvarbias_alb(std::string const &conf, char const *key) :
   if (cvm::debug())
     cvm::log(" bias.\n");
 
+  return COLVARS_OK;
 }
 
-colvarbias_alb::~colvarbias_alb() {
 
+colvarbias_alb::~colvarbias_alb()
+{
   if (cvm::n_rest_biases > 0)
     cvm::n_rest_biases -= 1;
-
 }
 
-int colvarbias_alb::update() {
+
+int colvarbias_alb::update()
+{
 
   bias_energy = 0.0;
   update_calls++;
@@ -129,9 +142,6 @@ int colvarbias_alb::update() {
   if (cvm::debug())
     cvm::log("Updating the adaptive linear bias \""+this->name+"\".\n");
 
-
-
-
   //log the moments of the CVs
   // Force and energy calculation
   bool finished_equil_flag = 1;
@@ -423,17 +433,24 @@ std::ostream & colvarbias_alb::write_traj(std::ostream &os)
 }
 
 
-cvm::real colvarbias_alb::restraint_potential(cvm::real k,  const colvar* x,  const colvarvalue &xcenter) const
+cvm::real colvarbias_alb::restraint_potential(cvm::real k,
+                                              colvar const *x,
+                                              colvarvalue const &xcenter) const
 {
   return k * (x->value() - xcenter);
 }
 
-colvarvalue colvarbias_alb::restraint_force(cvm::real k,  const colvar* x,  const colvarvalue &xcenter) const
+
+colvarvalue colvarbias_alb::restraint_force(cvm::real k,
+                                            colvar const *x,
+                                            colvarvalue const &xcenter) const
 {
   return k;
 }
 
-cvm::real colvarbias_alb::restraint_convert_k(cvm::real k, cvm::real dist_measure) const
+
+cvm::real colvarbias_alb::restraint_convert_k(cvm::real k,
+                                              cvm::real dist_measure) const
 {
   return k / dist_measure;
 }

lib/colvars/colvarbias_alb.h

@@ -9,22 +9,15 @@
 class colvarbias_alb : public colvarbias {
 
 public:
-  colvarbias_alb(std::string const &conf, char const *key);
+
+  colvarbias_alb(char const *key);
   virtual ~colvarbias_alb();
-
+  virtual int init(std::string const &conf);
-
   virtual int update();
 
-  /// Read the bias configuration from a restart file
   virtual std::istream & read_restart(std::istream &is);
-
-  /// Write the bias configuration to a restart file
   virtual std::ostream & write_restart(std::ostream &os);
-
-  /// 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);
 
 protected:

lib/colvars/colvarbias_histogram.cpp

@@ -6,10 +6,20 @@
 
 /// Histogram "bias" constructor
 
-colvarbias_histogram::colvarbias_histogram(std::string const &conf, char const *key)
+colvarbias_histogram::colvarbias_histogram(char const *key)
-  : colvarbias(conf, key),
+  : colvarbias(key),
     grid(NULL), out_name("")
 {
+}
+
+
+int colvarbias_histogram::init(std::string const &conf)
+{
+  colvarbias::init(conf);
+
+  provide(f_cvb_history_dependent);
+  enable(f_cvb_history_dependent);
+
   size_t i;
 
   get_keyval(conf, "outputFile", out_name, std::string(""));
@@ -30,18 +40,18 @@ colvarbias_histogram::colvarbias_histogram(std::string const &conf, char const *
       for (i = 0; i < colvars.size(); i++) { // should be all vector
         if (colvars[i]->value().type() != colvarvalue::type_vector) {
           cvm::error("Error: used gatherVectorColvars with non-vector colvar.\n", INPUT_ERROR);
-          return;
+          return INPUT_ERROR;
         }
         if (i == 0) {
           colvar_array_size = colvars[i]->value().size();
           if (colvar_array_size < 1) {
             cvm::error("Error: vector variable has dimension less than one.\n", INPUT_ERROR);
-            return;
+            return INPUT_ERROR;
           }
         } else {
           if (colvar_array_size != colvars[i]->value().size()) {
             cvm::error("Error: trying to combine vector colvars of different lengths.\n", INPUT_ERROR);
-            return;
+            return INPUT_ERROR;
           }
         }
       }
@@ -49,7 +59,7 @@ colvarbias_histogram::colvarbias_histogram(std::string const &conf, char const *
       for (i = 0; i < colvars.size(); i++) { // should be all scalar
         if (colvars[i]->value().type() != colvarvalue::type_scalar) {
           cvm::error("Error: only scalar colvars are supported when gatherVectorColvars is off.\n", INPUT_ERROR);
-          return;
+          return INPUT_ERROR;
         }
       }
     }
@@ -75,10 +85,10 @@ colvarbias_histogram::colvarbias_histogram(std::string const &conf, char const *
     }
   }
 
-  cvm::log("Finished histogram setup.\n");
+  return COLVARS_OK;
 }
 
-/// Destructor
+
 colvarbias_histogram::~colvarbias_histogram()
 {
   if (grid) {
@@ -90,7 +100,7 @@ colvarbias_histogram::~colvarbias_histogram()
     cvm::n_histo_biases -= 1;
 }
 
-/// Update the grid
+
 int colvarbias_histogram::update()
 {
   int error_code = COLVARS_OK;
@@ -124,7 +134,7 @@ int colvarbias_histogram::update()
     // update indices for scalar values
     size_t i;
     for (i = 0; i < colvars.size(); i++) {
-      bin[i] = grid->value_to_bin_scalar(colvars[i]->value(), i);
+      bin[i] = grid->current_bin_scalar(i);
     }
 
     if (grid->index_ok(bin)) {
@@ -135,7 +145,7 @@ int colvarbias_histogram::update()
     size_t iv, i;
     for (iv = 0; iv < colvar_array_size; iv++) {
       for (i = 0; i < colvars.size(); i++) {
-        bin[i] = grid->value_to_bin_scalar(colvars[i]->value().vector1d_value[iv], i);
+        bin[i] = grid->current_bin_scalar(i, iv);
       }
 
       if (grid->index_ok(bin)) {

lib/colvars/colvarbias_histogram.h

@@ -16,12 +16,11 @@ class colvarbias_histogram : public colvarbias {
 
 public:
 
-  colvarbias_histogram(std::string const &conf, char const *key);
+  colvarbias_histogram(char const *key);
   ~colvarbias_histogram();
-
+  virtual int init(std::string const &conf);
-  int update();
+  virtual int update();
-
+  virtual int write_output_files();
-  int write_output_files();
 
 protected:
 
@@ -36,8 +35,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;
 
-  std::istream& read_restart(std::istream&);
+  virtual std::istream& read_restart(std::istream&);
-  std::ostream& write_restart(std::ostream&);
+  virtual std::ostream& write_restart(std::ostream&);
 };
 
 #endif

lib/colvars/colvarbias_meta.cpp

@@ -25,30 +25,38 @@
 #include "colvarbias_meta.h"
 
 
-colvarbias_meta::colvarbias_meta()
+colvarbias_meta::colvarbias_meta(char const *key)
-  : colvarbias(),
+  : colvarbias(key),
     new_hills_begin(hills.end()),
     state_file_step(0)
 {
 }
 
 
-colvarbias_meta::colvarbias_meta(std::string const &conf, char const *key)
+int colvarbias_meta::init(std::string const &conf)
-  : colvarbias(conf, key),
-    new_hills_begin(hills.end()),
-    state_file_step(0)
 {
-  if (cvm::n_abf_biases > 0)
+  colvarbias::init(conf);
-    cvm::log("Warning: running ABF and metadynamics together is not recommended unless applyBias is off for ABF.\n");
+
+  provide(f_cvb_history_dependent);
 
   get_keyval(conf, "hillWeight", hill_weight, 0.0);
-  if (hill_weight <= 0.0) {
+  if (hill_weight > 0.0) {
+    enable(f_cvb_apply_force);
+  } else {
     cvm::error("Error: hillWeight must be provided, and a positive number.\n", INPUT_ERROR);
   }
 
   get_keyval(conf, "newHillFrequency", new_hill_freq, 1000);
+  if (new_hill_freq > 0) {
+    enable(f_cvb_history_dependent);
+  }
 
   get_keyval(conf, "hillWidth", hill_width, std::sqrt(2.0 * PI) / 2.0);
+  cvm::log("Half-widths of the Gaussian hills (sigma's):\n");
+  for (size_t i = 0; i < colvars.size(); i++) {
+    cvm::log(colvars[i]->name+std::string(": ")+
+             cvm::to_str(0.5 * colvars[i]->width * hill_width));
+  }
 
   {
     bool b_replicas = false;
@@ -152,6 +160,7 @@ colvarbias_meta::colvarbias_meta(std::string const &conf, char const *key)
              ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n");
 
   save_delimiters = false;
+  return COLVARS_OK;
 }
 
 
@@ -851,7 +860,7 @@ void colvarbias_meta::update_replicas_registry()
         // add this replica to the registry
         cvm::log("Metadynamics bias \""+this->name+"\""+
                  ": accessing replica \""+new_replica+"\".\n");
-        replicas.push_back(new colvarbias_meta());
+        replicas.push_back(new colvarbias_meta("metadynamics"));
         (replicas.back())->replica_id = new_replica;
         (replicas.back())->replica_list_file = new_replica_file;
         (replicas.back())->replica_state_file = "";

lib/colvars/colvarbias_meta.h

@@ -27,23 +27,13 @@ public:
   /// Communication between different replicas
   Communication comm;
 
-  /// Constructor
+  colvarbias_meta(char const *key);
-  colvarbias_meta(std::string const &conf, char const *key);
+  virtual int init(std::string const &conf);
-
-  /// Default constructor
-  colvarbias_meta();
-
-  /// Destructor
   virtual ~colvarbias_meta();
-
   virtual int update();
-
   virtual std::istream & read_restart(std::istream &is);
-
   virtual std::ostream & write_restart(std::ostream &os);
-
   virtual int setup_output();
-
   virtual void write_pmf();
 
   class hill;

lib/colvars/colvarbias_restraint.cpp

@@ -5,12 +5,25 @@
 #include "colvarbias_restraint.h"
 
 
-colvarbias_restraint::colvarbias_restraint(std::string const &conf,
+colvarbias_restraint::colvarbias_restraint(char const *key)
-                                           char const *key)
+  : colvarbias(key)
-  : colvarbias(conf, key),
-    target_nstages(0),
-    target_nsteps(0)
 {
+}
+
+
+int colvarbias_restraint::init(std::string const &conf)
+{
+  colvarbias::init(conf);
+
+  if (cvm::debug())
+    cvm::log("Initializing a new restraint bias.\n");
+
+  // TODO move these initializations to constructor and let get_keyval
+  // only override existing values
+  target_nstages = 0;
+  target_nsteps = 0;
+  force_k = 0.0;
+
   get_keyval(conf, "forceConstant", force_k, 1.0);
 
   {
@@ -124,6 +137,7 @@ colvarbias_restraint::colvarbias_restraint(std::string const &conf,
     cvm::log("Done initializing a new restraint bias.\n");
 }
 
+
 colvarbias_restraint::~colvarbias_restraint()
 {
   if (cvm::n_rest_biases > 0)
@@ -505,55 +519,92 @@ std::ostream & colvarbias_restraint::write_traj(std::ostream &os)
   return os;
 }
 
-colvarbias_restraint_harmonic::colvarbias_restraint_harmonic(std::string const &conf, char const *key) :
+
-  colvarbias_restraint(conf, key) {
+colvarbias_restraint_harmonic::colvarbias_restraint_harmonic(char const *key)
+  : colvarbias_restraint(key)
+{
+}
+
+
+int colvarbias_restraint_harmonic::init(std::string const &conf)
+{
+  colvarbias_restraint::init(conf);
+
   for (size_t i = 0; i < colvars.size(); i++) {
     if (colvars[i]->width != 1.0)
       cvm::log("The force constant for colvar \""+colvars[i]->name+
-                "\" will be rescaled to "+
+               "\" will be rescaled to "+
-                cvm::to_str(restraint_convert_k(force_k, colvars[i]->width))+
+               cvm::to_str(restraint_convert_k(force_k, colvars[i]->width))+
-                " according to the specified width.\n");
+               " according to the specified width.\n");
   }
 }
 
-cvm::real colvarbias_restraint_harmonic::restraint_potential(cvm::real k,  colvar* x,  const colvarvalue &xcenter) const
+
+cvm::real colvarbias_restraint_harmonic::restraint_potential(cvm::real k,
+                                                             colvar const *x,
+                                                             colvarvalue const &xcenter) const
 {
   return 0.5 * k * x->dist2(x->value(), xcenter);
 }
 
-colvarvalue colvarbias_restraint_harmonic::restraint_force(cvm::real k,  colvar* x,  const colvarvalue &xcenter) const
+
+colvarvalue colvarbias_restraint_harmonic::restraint_force(cvm::real k,
+                                                           colvar const *x,
+                                                           colvarvalue const &xcenter) const
 {
   return 0.5 * k * x->dist2_lgrad(x->value(), xcenter);
 }
 
-cvm::real colvarbias_restraint_harmonic::restraint_convert_k(cvm::real k, cvm::real dist_measure) const
+
+cvm::real colvarbias_restraint_harmonic::restraint_convert_k(cvm::real k,
+                                                             cvm::real dist_measure) const
 {
   return k / (dist_measure * dist_measure);
 }
 
 
-colvarbias_restraint_linear::colvarbias_restraint_linear(std::string const &conf, char const *key) :
+
-  colvarbias_restraint(conf, key) {
+colvarbias_restraint_linear::colvarbias_restraint_linear(char const *key)
+  : colvarbias_restraint(key)
+{
+}
+
+
+int colvarbias_restraint_linear::init(std::string const &conf)
+{
+  colvarbias_restraint::init(conf);
+
   for (size_t i = 0; i < colvars.size(); i++) {
     if (colvars[i]->width != 1.0)
       cvm::log("The force constant for colvar \""+colvars[i]->name+
-                "\" will be rescaled to "+
+               "\" will be rescaled to "+
-                cvm::to_str(restraint_convert_k(force_k, colvars[i]->width))+
+               cvm::to_str(restraint_convert_k(force_k, colvars[i]->width))+
-                " according to the specified width.\n");
+               " according to the specified width.\n");
   }
 }
 
-cvm::real colvarbias_restraint_linear::restraint_potential(cvm::real k,  colvar* x,  const colvarvalue &xcenter) const
+
+cvm::real colvarbias_restraint_linear::restraint_potential(cvm::real k,
+                                                           colvar const *x,
+                                                           colvarvalue const &xcenter) const
 {
   return k * (x->value() - xcenter);
 }
 
-colvarvalue colvarbias_restraint_linear::restraint_force(cvm::real k,  colvar* x,  const colvarvalue &xcenter) const
+
+colvarvalue colvarbias_restraint_linear::restraint_force(cvm::real k,
+                                                         colvar const *x,
+                                                         colvarvalue const &xcenter) const
 {
   return k;
 }
 
-cvm::real colvarbias_restraint_linear::restraint_convert_k(cvm::real k, cvm::real dist_measure) const
+
+cvm::real colvarbias_restraint_linear::restraint_convert_k(cvm::real k,
+                                                           cvm::real dist_measure) const
 {
   return k / dist_measure;
 }
+
+
+

lib/colvars/colvarbias_restraint.h

@@ -14,6 +14,7 @@ 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);
 
@@ -33,19 +34,21 @@ public:
   virtual std::ostream & write_traj(std::ostream &os);
 
   /// \brief Constructor
-  colvarbias_restraint(std::string const &conf, char const *key);
+  colvarbias_restraint(char const *key);
 
-  /// Destructor
+  virtual int init(std::string const &conf);
   virtual ~colvarbias_restraint();
 
 
 protected:
 
   /// \brief Potential function
-  virtual cvm::real restraint_potential(cvm::real k, colvar* x, const colvarvalue& xcenter) const = 0;
+  virtual cvm::real restraint_potential(cvm::real k, colvar const *x,
+                                        colvarvalue const &xcenter) const = 0;
 
   /// \brief Force function
-  virtual colvarvalue restraint_force(cvm::real k, colvar* x, const colvarvalue& xcenter) const = 0;
+  virtual colvarvalue restraint_force(cvm::real k, colvar const *x,
+                                      colvarvalue const &xcenter) const = 0;
 
   ///\brief Unit scaling
   virtual cvm::real restraint_convert_k(cvm::real k, cvm::real dist_measure) const = 0;
@@ -113,36 +116,51 @@ protected:
   long target_nsteps;
 };
 
+
 /// \brief Harmonic bias restraint
 /// (implementation of \link colvarbias_restraint \endlink)
 class colvarbias_restraint_harmonic : public colvarbias_restraint {
 
 public:
-  colvarbias_restraint_harmonic(std::string const &conf, char const *key);
 
-protected: /// \brief Potential function
+  colvarbias_restraint_harmonic(char const *key);
-  virtual cvm::real restraint_potential(cvm::real k,  colvar*  x, const colvarvalue& xcenter) const;
+  virtual int init(std::string const &conf);
+  // no additional members, destructor not needed
+
+protected:
+
+  /// \brief Potential function
+  virtual cvm::real restraint_potential(cvm::real k, colvar const *x,
+                                        colvarvalue const &xcenter) const;
 
   /// \brief Force function
-  virtual colvarvalue restraint_force(cvm::real k,  colvar* x,  const colvarvalue& xcenter) const;
+  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;
 
 };
 
+
 /// \brief Linear bias restraint
 /// (implementation of \link colvarbias_restraint \endlink)
 class colvarbias_restraint_linear : public colvarbias_restraint {
 
 public:
-  colvarbias_restraint_linear(std::string const &conf, char const *key);
+  colvarbias_restraint_linear(char const *key);
+  virtual int init(std::string const &conf);
+  // no additional members, destructor not needed
 
-protected: /// \brief Potential function
+protected:
-  virtual cvm::real restraint_potential(cvm::real k,  colvar*  x, const colvarvalue& xcenter) const;
+
+  /// \brief Potential function
+  virtual cvm::real restraint_potential(cvm::real k, colvar const *x,
+                                        colvarvalue const &xcenter) const;
 
   /// \brief Force function
-  virtual colvarvalue restraint_force(cvm::real k,  colvar* x,  const colvarvalue& xcenter) const;
+  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/colvardeps.cpp

@@ -210,6 +210,8 @@ void cvm::deps::init_cvb_requires() {
   f_description(f_cvb_get_system_force, "obtain system force");
   f_req_children(f_cvb_get_system_force, f_cv_system_force);
 
+  f_description(f_cvb_history_dependent, "history-dependent");
+
   // Initialize feature_states for each instance
   feature_states.reserve(f_cvb_ntot);
   for (i = 0; i < f_cvb_ntot; i++) {
@@ -217,6 +219,9 @@ void cvm::deps::init_cvb_requires() {
     // Most features are available, so we set them so
     // and list exceptions below
   }
+
+  // some biases are not history-dependent
+  feature_states[f_cvb_history_dependent]->available = false;
 }
 
 

lib/colvars/colvardeps.h

@@ -147,8 +147,9 @@ public:
   enum features_biases {
     /// \brief Bias is active
     f_cvb_active,
-    f_cvb_apply_force,
+    f_cvb_apply_force, // will apply forces
-    f_cvb_get_system_force,
+    f_cvb_get_system_force, // requires system forces
+    f_cvb_history_dependent, // depends on simulation history
     f_cvb_ntot
   };
 

lib/colvars/colvargrid.h

@@ -368,6 +368,14 @@ 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 item iv in variable i
+  inline int current_bin_scalar(int const i, int const iv) const
+  {
+    return value_to_bin_scalar(actual_value[i] ?
+        cv[i]->actual_value().vector1d_value[iv] :
+        cv[i]->value().vector1d_value[iv], i);
+  }
+
   /// \brief Use the lower boundary and the width to report which bin
   /// the provided value is in
   inline int value_to_bin_scalar(colvarvalue const &value, const int i) const

lib/colvars/colvarmodule.cpp

@@ -32,8 +32,9 @@ colvarmodule::colvarmodule(colvarproxy *proxy_in)
   cvm::log(cvm::line_marker);
   cvm::log("Initializing the collective variables module, version "+
            cvm::to_str(COLVARS_VERSION)+".\n");
-  cvm::log("Please cite Fiorin et al, Mol Phys 2013 in any publication "
+  cvm::log("Please cite Fiorin et al, Mol Phys 2013:\n "
-           "based on this calculation.\n");
+           "http://dx.doi.org/10.1080/00268976.2013.813594\n"
+           "in any publication based on this calculation.\n");
 
   if (proxy->smp_enabled() == COLVARS_OK) {
     cvm::log("SMP parallelism is available.\n");
@@ -250,8 +251,9 @@ int colvarmodule::parse_biases_type(std::string const &conf,
     if (bias_conf.size()) {
       cvm::log(cvm::line_marker);
       cvm::increase_depth();
-      biases.push_back(new bias_type(bias_conf, keyword));
+      biases.push_back(new bias_type(keyword));
-      if (cvm::check_new_bias(bias_conf, keyword)) {
+      biases.back()->init(bias_conf);
+      if (cvm::check_new_bias(bias_conf, keyword) != COLVARS_OK) {
         return COLVARS_ERROR;
       }
       cvm::decrease_depth();
@@ -297,12 +299,27 @@ int colvarmodule::parse_biases(std::string const &conf)
     cvm::decrease_depth();
   }
 
-  for (size_t i = 0; i < biases.size(); i++) {
+  size_t i;
+
+  for (i = 0; i < biases.size(); i++) {
     biases[i]->enable(cvm::deps::f_cvb_active);
     if (cvm::debug())
       biases[i]->print_state();
   }
 
+  size_t n_hist_dep_biases = 0;
+  for (i = 0; i < biases.size(); i++) {
+    if (biases[i]->is_enabled(cvm::deps::f_cvb_apply_force) &&
+        biases[i]->is_enabled(cvm::deps::f_cvb_history_dependent)) {
+      n_hist_dep_biases++;
+    }
+  }
+  if (n_hist_dep_biases) {
+    cvm::log("WARNING: there are "+cvm::to_str(n_hist_dep_biases)+
+             " history-dependent biases with non-zero force parameters; "
+             "please make sure that their forces do not counteract each other.\n");
+  }
+
   if (biases.size() || use_scripted_forces) {
     cvm::log(cvm::line_marker);
     cvm::log("Collective variables biases initialized, "+
@@ -322,7 +339,7 @@ int colvarmodule::catch_input_errors(int result)
   if (result != COLVARS_OK || get_error()) {
     set_error_bit(result);
     set_error_bit(INPUT_ERROR);
-    parse->reset();
+    parse->init();
     return get_error();
   }
   return COLVARS_OK;
@@ -795,6 +812,7 @@ colvarmodule::~colvarmodule()
       (proxy->smp_thread_id() == 0)) {
     reset();
     delete parse;
+    parse = NULL;
     proxy = NULL;
   }
 }
@@ -802,7 +820,7 @@ colvarmodule::~colvarmodule()
 
 int colvarmodule::reset()
 {
-  parse->reset();
+  parse->init();
 
   cvm::log("Resetting the Collective Variables Module.\n");
 

lib/colvars/colvarmodule.h

@@ -4,7 +4,7 @@
 #define COLVARMODULE_H
 
 #ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2016-06-14"
+#define COLVARS_VERSION "2016-07-05"
 #endif
 
 #ifndef COLVARS_DEBUG

lib/colvars/colvarparse.cpp

@@ -524,7 +524,7 @@ int colvarparse::check_keywords(std::string &conf, char const *key)
     }
   }
 
-  clear_keyword_registry();
+  init();
 
   return COLVARS_OK;
 }

lib/colvars/colvarparse.h

@@ -53,14 +53,35 @@ protected:
 
 public:
 
+
   inline colvarparse()
     : save_delimiters(true)
-  {}
+  {
+    init();
+  }
 
   /// Constructor that stores the object's config string
   inline colvarparse(const std::string& conf)
-    : save_delimiters(true), config_string(conf)
+    : save_delimiters(true)
-  {}
+  {
+    init(conf);
+  }
+
+  /// Set the object ready to parse a new configuration string
+  inline void init()
+  {
+    config_string.clear();
+    clear_keyword_registry();
+  }
+
+  /// Set a new config string for this object
+  inline void init(const std::string& conf)
+  {
+    if (! config_string.size()) {
+      init();
+      config_string = conf;
+    }
+  }
 
   inline const std::string& get_config()
   {
@@ -79,6 +100,16 @@ public:
     parse_silent
   };
 
+  /// \brief Check that all the keywords within "conf" are in the list
+  /// of allowed keywords; this will invoke strip_values() first and
+  /// then loop over all words
+  int check_keywords(std::string &conf, char const *key);
+
+  /// \brief Use this after parsing a config string (note that check_keywords() calls it already)
+  void clear_keyword_registry();
+
+public:
+
   /// \fn get_keyval bool const get_keyval (std::string const &conf,
   /// char const *key, _type_ &value, _type_ const &def_value,
   /// Parse_Mode const parse_mode) \brief Helper function to parse
@@ -210,18 +241,9 @@ protected:
                                                    std::vector<TYPE> &values,
                                                    std::vector<TYPE> const &def_values,
                                                    Parse_Mode const parse_mode);
+
 public:
 
-  /// \brief Check that all the keywords within "conf" are in the list
-  /// of allowed keywords; this will invoke strip_values() first and
-  /// then loop over all words
-  int check_keywords(std::string &conf, char const *key);
-
-  /// \brief Use this after parsing a config string (note that check_keywords() calls it already)
-  void clear_keyword_registry();
-
-  inline void reset() { clear_keyword_registry(); }
-
   /// \brief Return a lowercased copy of the string
   static inline std::string to_lower_cppstr(std::string const &in)
   {