From b69be175f03afb2109f69d77e165b45c43c1ec83 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Thu, 14 Aug 2014 23:29:42 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@12326
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 lib/colvars/colvarbias_alb.cpp | 426 +++++++++++++++++++++++++++++++++
 lib/colvars/colvarbias_alb.h   |  85 +++++++
 2 files changed, 511 insertions(+)
 create mode 100644 lib/colvars/colvarbias_alb.cpp
 create mode 100644 lib/colvars/colvarbias_alb.h

diff --git a/lib/colvars/colvarbias_alb.cpp b/lib/colvars/colvarbias_alb.cpp
new file mode 100644
index 0000000000..41e01d9695
--- /dev/null
+++ b/lib/colvars/colvarbias_alb.cpp
@@ -0,0 +1,426 @@
+#include "colvarmodule.h"
+#include "colvarbias_alb.h"
+#include "colvarbias.h"
+#include <stdio.h>
+#include <stdlib.h>
+
+/* Note about nomenclature. Force constant is called a coupling
+ * constant here to emphasize its changing in the code. Outwards,
+ * everything is called a force constant to keep it consistent with
+ * the rest of colvars.
+ *
+ */
+
+colvarbias_alb::colvarbias_alb(std::string const &conf, char const *key) :
+  colvarbias(conf, key), update_calls(0), b_equilibration(true) {
+
+  // get the initial restraint centers
+  colvar_centers.resize (colvars.size());
+
+  means.resize(colvars.size());
+  ssd.resize(colvars.size()); //sum of squares of differences from mean
+
+  //setup force vectors
+  max_coupling_range.resize(colvars.size());
+  max_coupling_rate.resize(colvars.size());
+  coupling_accum.resize(colvars.size());
+  set_coupling.resize(colvars.size());
+  current_coupling.resize(colvars.size());
+  coupling_rate.resize(colvars.size());
+
+
+  for (size_t i = 0; i < colvars.size(); i++) {
+    colvar_centers[i].type (colvars[i]->type());
+    //zero moments
+    means[i] = ssd[i] = 0;
+
+    //zero force some of the force vectors that aren't initialized
+    coupling_accum[i] = current_coupling[i] = 0;
+
+  }
+  if (get_keyval (conf, "centers", colvar_centers, colvar_centers)) {
+    for (size_t i = 0; i < colvars.size(); i++) {
+      colvar_centers[i].apply_constraints();
+    }
+  } else {
+    colvar_centers.clear();
+    cvm::fatal_error ("Error: must define the initial centers of adaptive linear bias .\n");
+  }
+
+  if (colvar_centers.size() != colvars.size())
+    cvm::fatal_error ("Error: number of centers does not match "
+                      "that of collective variables.\n");
+
+  if(!get_keyval (conf, "UpdateFrequency", update_freq, 0))
+    cvm::fatal_error("Error: must set updateFrequency for apadtive linear bias.\n");
+
+  //we split the time between updating and equilibrating
+  update_freq /= 2;
+
+  if(update_freq == 0)
+    cvm::fatal_error("Error: must set updateFrequency to greater than 2.\n");
+
+  get_keyval (conf, "outputCenters", b_output_centers, false);
+  get_keyval (conf, "outputGradient", b_output_grad, false);
+  get_keyval (conf, "outputCoupling", b_output_coupling, true);
+  get_keyval (conf, "hardForceRange", b_hard_coupling_range, true);
+
+  //initial guess
+  if(!get_keyval (conf, "forceConstant", set_coupling, set_coupling))
+    for(size_t i =0 ; i < colvars.size(); i++)
+      set_coupling[i] = 0.;
+
+  //how we're going to increase to that point
+  for(size_t i = 0; i < colvars.size(); i++)
+    coupling_rate[i] = (set_coupling[i] - current_coupling[i]) / update_freq;
+
+
+  if(!get_keyval (conf, "forceRange", max_coupling_range, max_coupling_range)) {
+    //set to default
+    for(size_t i = 0; i < colvars.size(); i++) {
+      if(cvm::temperature() > 0)
+	max_coupling_range[i] =   3 * cvm::temperature() * cvm::boltzmann();
+      else
+	max_coupling_range[i] =   3 * cvm::boltzmann();
+    }
+  }
+
+
+
+  if(!get_keyval (conf, "rateMax", max_coupling_rate, max_coupling_rate)) {
+    //set to default
+    for(size_t i = 0; i < colvars.size(); i++) {
+      max_coupling_rate[i] =   max_coupling_range[i] / (10 * update_freq);
+    }
+  }
+
+
+  if (cvm::debug())
+    cvm::log (" bias.\n");
+
+}
+
+colvarbias_alb::~colvarbias_alb() {
+
+  if (cvm::n_rest_biases > 0)
+    cvm::n_rest_biases -= 1;
+
+}
+
+cvm::real colvarbias_alb::update() {
+
+  bias_energy = 0.0;
+  update_calls++;
+
+  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;
+  cvm::real delta;
+  for (size_t i = 0; i < colvars.size(); i++) {
+    colvar_forces[i] = -restraint_force(restraint_convert_k(current_coupling[i], colvars[i]->width),
+					colvars[i],
+					colvar_centers[i]);
+    bias_energy += restraint_potential(restraint_convert_k(current_coupling[i], colvars[i]->width),
+				       colvars[i],
+				       colvar_centers[i]);
+
+    if(!b_equilibration) {
+      //Welford, West, and Hanso online variance method
+
+      delta = static_cast<cvm::real>(colvars[i]->value())  - means[i];
+      means[i] += delta / update_calls;
+      ssd[i] += delta * (static_cast<cvm::real>(colvars[i]->value())  - means[i]);
+
+    } else {
+      //check if we've reached the setpoint
+      if(coupling_rate[i] == 0 || pow(current_coupling[i] - set_coupling[i],2)  < pow(coupling_rate[i],2)) {
+	finished_equil_flag &= 1; //we continue equilibrating as long as we haven't reached all the set points
+      }
+      else {
+	current_coupling[i] += coupling_rate[i];
+	finished_equil_flag = 0;
+      }
+
+
+      //update max_coupling_range
+      if(!b_hard_coupling_range && fabs(current_coupling[i]) > max_coupling_range[i]) {
+	std::ostringstream logStream;
+	logStream << "Coupling constant for "
+		  << colvars[i]->name
+		  << " has exceeded coupling range of "
+		  << max_coupling_range[i]
+		  << ".\n";
+
+	max_coupling_range[i] *= 1.25;
+	logStream << "Expanding coupling range to "  << max_coupling_range[i] << ".\n";
+	cvm::log(logStream.str());
+      }
+
+
+    }
+  }
+
+  if(b_equilibration && finished_equil_flag) {
+    b_equilibration = false;
+    update_calls = 0;
+  }
+
+
+  //now we update coupling constant, if necessary
+  if(!b_equilibration && update_calls == update_freq) {
+
+    //use estimated variance to take a step
+    cvm::real step_size = 0;
+    cvm::real temp;
+
+    //reset means and sum of squares of differences
+    for(size_t i = 0; i < colvars.size(); i++) {
+
+      temp = 2. * (means[i] / (static_cast<cvm::real> (colvar_centers[i])) - 1) * ssd[i] / (update_calls - 1);
+
+      if(cvm::temperature() > 0)
+	step_size = temp / (cvm::temperature()  * cvm::boltzmann());
+      else
+	step_size = temp / cvm::boltzmann();
+
+      means[i] = 0;
+      ssd[i] = 0;
+
+      //stochastic if we do that update or not
+      if(colvars.size() == 1 || rand() < RAND_MAX / colvars.size()) {
+	coupling_accum[i] += step_size * step_size;
+	current_coupling[i] = set_coupling[i];
+	set_coupling[i] += max_coupling_range[i] / sqrt(coupling_accum[i]) * step_size;
+	coupling_rate[i] = (set_coupling[i] - current_coupling[i]) / update_freq;
+	//set to the minimum rate and then put the sign back on it
+	coupling_rate[i] = copysign(fmin(fabs(coupling_rate[i]), max_coupling_rate[i]), coupling_rate[i]);
+      } else {
+	coupling_rate[i] = 0;
+      }
+
+    }
+
+    update_calls = 0;
+    b_equilibration = true;
+
+  }
+
+  return bias_energy;
+
+}
+
+
+std::istream & colvarbias_alb::read_restart (std::istream &is)
+{
+  size_t const start_pos = is.tellg();
+
+  cvm::log ("Restarting adaptive linear bias \""+
+            this->name+"\".\n");
+
+  std::string key, brace, conf;
+  if ( !(is >> key)   || !(key == "ALB") ||
+       !(is >> brace) || !(brace == "{") ||
+       !(is >> colvarparse::read_block ("configuration", conf)) ) {
+
+    cvm::log ("Error: in reading restart configuration for restraint bias \""+
+              this->name+"\" at position "+
+              cvm::to_str (is.tellg())+" in stream.\n");
+    is.clear();
+    is.seekg (start_pos, std::ios::beg);
+    is.setstate (std::ios::failbit);
+    return is;
+  }
+
+  std::string name = "";
+  if ( (colvarparse::get_keyval (conf, "name", name, std::string (""), colvarparse::parse_silent)) &&
+       (name != this->name) )
+    cvm::fatal_error ("Error: in the restart file, the "
+                      "\"ALB\" block has a wrong name\n");
+  if (name.size() == 0) {
+    cvm::fatal_error ("Error: \"ALB\" block in the restart file "
+                      "has no identifiers.\n");
+  }
+
+  if (!get_keyval (conf, "setCoupling", set_coupling))
+    cvm::fatal_error ("Error: current setCoupling  is missing from the restart.\n");
+
+  if (!get_keyval (conf, "currentCoupling", current_coupling))
+    cvm::fatal_error ("Error: current setCoupling  is missing from the restart.\n");
+
+  if (!get_keyval (conf, "maxCouplingRange", max_coupling_range))
+    cvm::fatal_error ("Error: maxCouplingRange  is missing from the restart.\n");
+
+
+  if (!get_keyval (conf, "couplingRate", coupling_rate))
+    cvm::fatal_error ("Error: current setCoupling  is missing from the restart.\n");
+
+  if (!get_keyval (conf, "couplingAccum", coupling_accum))
+    cvm::fatal_error ("Error: couplingAccum is missing from the restart.\n");
+
+
+  if (!get_keyval (conf, "mean", means))
+    cvm::fatal_error ("Error: current mean is missing from the restart.\n");
+
+  if (!get_keyval (conf, "ssd", ssd))
+    cvm::fatal_error ("Error: current ssd is missing from the restart.\n");
+
+  if (!get_keyval (conf, "updateCalls", update_calls))
+    cvm::fatal_error ("Error: current updateCalls is missing from the restart.\n");
+
+  if (!get_keyval (conf, "b_equilibration", b_equilibration))
+    cvm::fatal_error ("Error: current updateCalls is missing from the restart.\n");
+
+  is >> brace;
+  if (brace != "}") {
+    cvm::fatal_error ("Error: corrupt restart information for adaptive linear bias \""+
+                      this->name+"\": no matching brace at position "+
+                      cvm::to_str (is.tellg())+" in the restart file.\n");
+    is.setstate (std::ios::failbit);
+  }
+
+  return is;
+}
+
+
+std::ostream & colvarbias_alb::write_restart (std::ostream &os)
+{
+  os << "ALB {\n"
+     << "  configuration {\n"
+     << "    name " << this->name << "\n";
+  os << "    setCoupling ";
+  for(size_t i = 0; i < colvars.size(); i++) {
+    os << std::setprecision (cvm::en_prec)
+       << std::setw (cvm::en_width) << set_coupling[i] << "\n";
+  }
+  os << "    currentCoupling ";
+  for(size_t i = 0; i < colvars.size(); i++) {
+    os << std::setprecision (cvm::en_prec)
+       << std::setw (cvm::en_width) << current_coupling[i] << "\n";
+  }
+  os << "    maxCouplingRange ";
+  for(size_t i = 0; i < colvars.size(); i++) {
+    os << std::setprecision (cvm::en_prec)
+       << std::setw (cvm::en_width) << max_coupling_range[i] << "\n";
+  }
+  os << "    couplingRate ";
+  for(size_t i = 0; i < colvars.size(); i++) {
+    os << std::setprecision (cvm::en_prec)
+       << std::setw (cvm::en_width) << coupling_rate[i] << "\n";
+  }
+  os << "    couplingAccum ";
+  for(size_t i = 0; i < colvars.size(); i++) {
+    os << std::setprecision (cvm::en_prec)
+       << std::setw (cvm::en_width) << coupling_accum[i] << "\n";
+  }
+  os << "    mean ";
+  for(size_t i = 0; i < colvars.size(); i++) {
+    os << std::setprecision (cvm::en_prec)
+       << std::setw (cvm::en_width) << means[i] << "\n";
+  }
+  os << "    ssd ";
+  for(size_t i = 0; i < colvars.size(); i++) {
+    os << std::setprecision (cvm::en_prec)
+       << std::setw (cvm::en_width) << ssd[i] << "\n";
+  }
+  os << "    updateCalls " << update_calls << "\n";
+  if(b_equilibration)
+    os << "    b_equilibration yes\n";
+  else
+    os << "    b_equilibration no\n";
+
+  os << "  }\n"
+     << "}\n\n";
+
+  return os;
+}
+
+
+std::ostream & colvarbias_alb::write_traj_label (std::ostream &os)
+{
+  os << " ";
+
+  if (b_output_energy)
+    os << " E_"
+       << cvm::wrap_string (this->name, cvm::en_width-2);
+
+  if (b_output_coupling)
+    for(size_t i = 0; i < current_coupling.size(); i++) {
+      os << " ForceConst_" << i
+	 <<std::setw(cvm::en_width - 6 - (i / 10 + 1))
+	 << "";
+    }
+
+  if(b_output_grad)
+    for(size_t i = 0; i < means.size(); i++) {
+      os << "Grad_"
+	 << cvm::wrap_string(colvars[i]->name, cvm::cv_width - 4);
+    }
+
+  if (b_output_centers)
+    for (size_t i = 0; i < colvars.size(); i++) {
+      size_t const this_cv_width = (colvars[i]->value()).output_width (cvm::cv_width);
+      os << " x0_"
+         << cvm::wrap_string (colvars[i]->name, this_cv_width-3);
+    }
+
+  return os;
+}
+
+
+std::ostream & colvarbias_alb::write_traj (std::ostream &os)
+{
+  os << " ";
+
+  if (b_output_energy)
+    os << " "
+       << std::setprecision (cvm::en_prec) << std::setw (cvm::en_width)
+       << bias_energy;
+
+  if(b_output_coupling)
+    for(size_t i = 0; i < current_coupling.size(); i++) {
+      os << " "
+	 << std::setprecision (cvm::en_prec) << std::setw (cvm::en_width)
+	 << current_coupling[i];
+    }
+
+
+  if (b_output_centers)
+    for (size_t i = 0; i < colvars.size(); i++) {
+      os << " "
+         << std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width)
+         << colvar_centers[i];
+    }
+
+  if(b_output_grad)
+    for(size_t i = 0; i < means.size(); i++) {
+      os << " "
+	 << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
+	 << -2. * (means[i] / (static_cast<cvm::real> (colvar_centers[i])) - 1) * ssd[i] / (fmax(update_calls,2) - 1);
+
+    }
+
+  return os;
+}
+
+
+cvm::real colvarbias_alb::restraint_potential(cvm::real k,  const colvar* x,  const colvarvalue &xcenter) const
+{
+  return k * (x->value() - xcenter);
+}
+
+colvarvalue colvarbias_alb::restraint_force(cvm::real k,  const colvar* x,  const colvarvalue &xcenter) const
+{
+  return k;
+}
+
+cvm::real colvarbias_alb::restraint_convert_k(cvm::real k, cvm::real dist_measure) const
+{
+  return k / dist_measure;
+}
+
diff --git a/lib/colvars/colvarbias_alb.h b/lib/colvars/colvarbias_alb.h
new file mode 100644
index 0000000000..20ad20cf75
--- /dev/null
+++ b/lib/colvars/colvarbias_alb.h
@@ -0,0 +1,85 @@
+#ifndef COLVARBIAS_ALB_H
+#define COLVARBIAS_ALB_H
+
+#include "colvar.h"
+#include "colvarbias.h"
+
+class colvarbias_alb : public colvarbias {
+
+public:
+  colvarbias_alb(std::string const &conf, char const *key);
+  virtual ~colvarbias_alb();
+
+
+  virtual cvm::real update();
+
+  /// Read the bias configuration from a restart file
+  virtual std::istream & read_restart (std::istream &is);
+
+  /// Write the bias configuration to a restart file
+  virtual std::ostream & write_restart (std::ostream &os);
+
+  /// 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:
+
+  /// \brief Restraint centers
+  std::vector<colvarvalue> colvar_centers;
+
+  /// \brief colvar parameters, used for calculating the gradient/variance
+  std::vector<cvm::real> means;
+  std::vector<cvm::real> ssd; // SSD = sum of squares of differences from mean
+  int update_calls;
+
+  ///\brief how often to update coupling constant
+  int update_freq;
+
+  ///\brief Estimated range of coupling constant values in kT
+  std::vector<cvm::real> max_coupling_range;
+
+  //\brief Estimated max for how quickly the rate can change in kT / time
+  std::vector<cvm::real> max_coupling_rate;
+
+  /// \brief accumated couping force; used in stochastic online gradient descent algorithm
+  std::vector<cvm::real> coupling_accum;
+
+  /// \brief coupling constant
+  std::vector<cvm::real> set_coupling;
+
+  /// \brief current coupling constant, which is ramped up during equilibration to coupling
+  std::vector<cvm::real> current_coupling;
+
+  /// \brief how quickly to change the coupling constant
+  std::vector<cvm::real> coupling_rate;
+
+  // \brief if we're equilibrating our estimates or collecting data
+  bool b_equilibration;
+
+  // \brief If the coupling range should be increased
+  bool b_hard_coupling_range;
+
+
+  /// \brief flag for outputting colvar centers
+  bool b_output_centers;
+
+  /// \brief flag for outputting current gradient
+  bool b_output_grad;
+
+  /// \brief flag for outputting coupling constant
+  bool b_output_coupling;
+
+  cvm::real restraint_potential(cvm::real k,  const colvar*  x, const colvarvalue& xcenter) const;
+
+  /// \brief Force function
+  colvarvalue restraint_force(cvm::real k,  const colvar* x,  const colvarvalue& xcenter) const;
+
+  ///\brief Unit scaling
+  cvm::real restraint_convert_k(cvm::real k, cvm::real dist_measure) const;
+
+};
+
+#endif