Update Colvars to version 2018-12-14 63/4863/2
authorGiacomo Fiorin <giacomo.fiorin@gmail.com>
Fri, 14 Dec 2018 20:22:02 +0000 (15:22 -0500)
committerDavid Hardy <dhardy@ks.uiuc.edu>
Fri, 14 Dec 2018 21:49:24 +0000 (15:49 -0600)
Bugfixes concerning corrFunc/runAve:
https://github.com/Colvars/colvars/pull/194

Fix for the dependency handling in orientationAngle, orientationProj, tilt,
abnd spinAngle (an error condition was raised when it shouldn't have):
https://github.com/Colvars/colvars/pull/199

The "cv save" command now triggers writes of output files from the ABF bias:
https://github.com/Colvars/colvars/pull/188

Relevant commit messages:
394145d 2018-12-14 Split off init() function for orientation derivatives [Giacomo Fiorin]
df94a7c 2018-12-12 Silence compiler warnings [Giacomo Fiorin]
026ef05 2018-11-16 Update version strings [Giacomo Fiorin]
4d8244f 2018-11-16 Make internal numbering of default colvar names consistent [Giacomo Fiorin]
331a205 2018-11-16 Add end_of_step() function to variables and biases [Giacomo Fiorin]
e71730c 2018-11-15 Fix declaration of local variable [Giacomo Fiorin]
07b9ff4 2018-11-15 Check that second variable is defined before accessing it [Giacomo Fiorin]
563f081 2018-11-15 Improve output format for runAve and corrFunc, cosmetic changes [Giacomo Fiorin]
475f493 2018-11-15 Use consistent return types for colvar-class analysis functions [Giacomo Fiorin]
acf8cc2 2018-11-15 Do not recompute finite-difference velocity; make it unavailable only during simulation [Giacomo Fiorin]
24b2db4 2018-11-15 Use properly name of second variable for correlation function [Giacomo Fiorin]
e5397e6 2018-11-15 Fix default file names for corrFunc and runAve output files in the doc [Giacomo Fiorin]
9fe9ff9 2018-11-15 Initialize corrFunc, runAve output files when cvm::output_prefix() is available [Giacomo Fiorin]
7e8bbaf 2018-10-30 Implement cv save in ABF [Jérôme Hénin]
d388035 2018-10-29 Add missing hyperlink [Giacomo Fiorin]

Change-Id: I22b3b0ac000bee69d4dc7fe3d41d657f426a3550

12 files changed:
colvars/src/colvar.cpp
colvars/src/colvar.h
colvars/src/colvarbias.cpp
colvars/src/colvarbias.h
colvars/src/colvarbias_abf.cpp
colvars/src/colvarbias_abf.h
colvars/src/colvardeps.cpp
colvars/src/colvarmodule.cpp
colvars/src/colvarmodule.h
colvars/src/colvars_version.h
colvars/src/colvarscript.cpp
ug/ug_colvars.tex

index 22befdf..8b28eaa 100644 (file)
 
 
 colvar::colvar()
-  : prev_timestep(-1)
 {
-  // Initialize static array once and for all
   runave_os = NULL;
+
+  prev_timestep = -1;
+  after_restart = false;
+  kinetic_energy = 0.0;
+  potential_energy = 0.0;
+
   init_cv_requires();
 }
 
@@ -38,6 +42,7 @@ namespace {
   }
 }
 
+
 int colvar::init(std::string const &conf)
 {
   cvm::log("Initializing a new collective variable.\n");
@@ -48,7 +53,7 @@ int colvar::init(std::string const &conf)
   colvarmodule *cv = cvm::main();
 
   get_keyval(conf, "name", this->name,
-             (std::string("colvar")+cvm::to_str(cv->variables()->size()+1)));
+             (std::string("colvar")+cvm::to_str(cv->variables()->size())));
 
   if ((cvm::colvar_by_name(this->name) != NULL) &&
       (cvm::colvar_by_name(this->name) != this)) {
@@ -63,9 +68,6 @@ int colvar::init(std::string const &conf)
 
   this->description = "colvar " + this->name;
 
-  kinetic_energy = 0.0;
-  potential_energy = 0.0;
-
   error_code |= init_components(conf);
   if (error_code != COLVARS_OK) {
     return cvm::get_error();
@@ -260,7 +262,6 @@ int colvar::init(std::string const &conf)
   f_old.reset();
 
   x_restart.type(value());
-  after_restart = false;
 
   reset_bias_force();
 
@@ -282,8 +283,7 @@ int colvar::init(std::string const &conf)
   // Now that the children are defined we can solve dependencies
   enable(f_cv_active);
 
-  if (cvm::b_analysis)
-    parse_analysis(conf);
+  error_code |= parse_analysis(conf);
 
   if (cvm::debug())
     cvm::log("Done initializing collective variable \""+this->name+"\".\n");
@@ -881,19 +881,7 @@ int colvar::parse_analysis(std::string const &conf)
       cvm::error("Error: runAveStride must be commensurate with the restart frequency.\n", INPUT_ERROR);
     }
 
-    get_keyval(conf, "runAveOutputFile", runave_outfile,
-                std::string(cvm::output_prefix()+"."+
-                             this->name+".runave.traj"));
-
-    size_t const this_cv_width = x.output_width(cvm::cv_width);
-    cvm::proxy->backup_file(runave_outfile);
-    runave_os = cvm::proxy->output_stream(runave_outfile);
-    *runave_os << "# " << cvm::wrap_string("step", cvm::it_width-2)
-               << "  "
-               << cvm::wrap_string("running average", this_cv_width)
-               << " "
-               << cvm::wrap_string("running stddev", this_cv_width)
-               << "\n";
+    get_keyval(conf, "runAveOutputFile", runave_outfile, runave_outfile);
   }
 
   acf_length = 0;
@@ -902,7 +890,6 @@ int colvar::parse_analysis(std::string const &conf)
 
     enable(f_cv_corrfunc);
 
-    std::string acf_colvar_name;
     get_keyval(conf, "corrFuncWithColvar", acf_colvar_name, this->name);
     if (acf_colvar_name == this->name) {
       cvm::log("Calculating auto-correlation function.\n");
@@ -918,8 +905,12 @@ int colvar::parse_analysis(std::string const &conf)
     } else if (acf_type_str == to_lower_cppstr(std::string("velocity"))) {
       acf_type = acf_vel;
       enable(f_cv_fdiff_velocity);
-      if (acf_colvar_name.size())
-        (cvm::colvar_by_name(acf_colvar_name))->enable(f_cv_fdiff_velocity);
+      colvar *cv2 = cvm::colvar_by_name(acf_colvar_name);
+      if (cv2 == NULL) {
+        return cvm::error("Error: collective variable \""+acf_colvar_name+
+                          "\" is not defined at this time.\n", INPUT_ERROR);
+      }
+      cv2->enable(f_cv_fdiff_velocity);
     } else if (acf_type_str == to_lower_cppstr(std::string("coordinate_p2"))) {
       acf_type = acf_p2coor;
     } else {
@@ -937,9 +928,7 @@ int colvar::parse_analysis(std::string const &conf)
     }
 
     get_keyval(conf, "corrFuncNormalize", acf_normalize, true);
-    get_keyval(conf, "corrFuncOutputFile", acf_outfile,
-                std::string(cvm::output_prefix()+"."+this->name+
-                             ".corrfunc.dat"));
+    get_keyval(conf, "corrFuncOutputFile", acf_outfile, acf_outfile);
   }
   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
 }
@@ -1389,9 +1378,9 @@ int colvar::calc_colvar_properties()
 {
   if (is_enabled(f_cv_fdiff_velocity)) {
     // calculate the velocity by finite differences
-    if (cvm::step_relative() == 0)
+    if (cvm::step_relative() == 0) {
       x_old = x;
-    else {
+    else {
       v_fdiff = fdiff_velocity(x_old, x);
       v_reported = v_fdiff;
     }
@@ -1486,7 +1475,6 @@ cvm::real colvar::update_forces_energy()
         return 0.;
       }
     }
-    prev_timestep = cvm::step_relative();
 
     // Integrate with slow timestep (if time_step_factor != 1)
     cvm::real dt = cvm::dt() * cvm::real(time_step_factor);
@@ -1547,8 +1535,18 @@ cvm::real colvar::update_forces_energy()
   // bypass the extended Lagrangian mass)
   f += fb_actual;
 
+  if (cvm::debug())
+    cvm::log("Done updating colvar \""+this->name+"\".\n");
+  return (potential_energy + kinetic_energy);
+}
+
+
+int colvar::end_of_step()
+{
+  if (cvm::debug())
+    cvm::log("End of step for colvar \""+this->name+"\".\n");
+
   if (is_enabled(f_cv_fdiff_velocity)) {
-    // set it for the next step
     x_old = x;
   }
 
@@ -1556,9 +1554,9 @@ cvm::real colvar::update_forces_energy()
     f_old = f;
   }
 
-  if (cvm::debug())
-    cvm::log("Done updating colvar \""+this->name+"\".\n");
-  return (potential_energy + kinetic_energy);
+  prev_timestep = cvm::step_relative();
+
+  return COLVARS_OK;
 }
 
 
@@ -1966,6 +1964,10 @@ std::ostream & colvar::write_restart(std::ostream &os) {
 
   os << "}\n\n";
 
+  if (runave_os) {
+    cvm::main()->proxy->flush_output_stream(runave_os);
+  }
+
   return os;
 }
 
@@ -2075,55 +2077,61 @@ std::ostream & colvar::write_traj(std::ostream &os)
   return os;
 }
 
+
 int colvar::write_output_files()
 {
-  if (cvm::b_analysis) {
+  int error_code = COLVARS_OK;
 
+  if (is_enabled(f_cv_corrfunc)) {
     if (acf.size()) {
-      cvm::log("Writing acf to file \""+acf_outfile+"\".\n");
-
+      if (acf_outfile.size() == 0) {
+        acf_outfile = std::string(cvm::output_prefix()+"."+this->name+
+                                  ".corrfunc.dat");
+      }
+      cvm::log("Writing correlation function to file \""+acf_outfile+"\".\n");
       cvm::backup_file(acf_outfile.c_str());
       std::ostream *acf_os = cvm::proxy->output_stream(acf_outfile);
       if (!acf_os) return cvm::get_error();
-      write_acf(*acf_os);
+      error_code |= write_acf(*acf_os);
       cvm::proxy->close_output_stream(acf_outfile);
     }
-
-    if (runave_os) {
-      cvm::proxy->close_output_stream(runave_outfile);
-      runave_os = NULL;
-    }
   }
-  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
+
+  return error_code;
 }
 
 
 
 // ******************** ANALYSIS FUNCTIONS ********************
 
-void colvar::analyze()
+int colvar::analyze()
 {
+  int error_code = COLVARS_OK;
+
   if (is_enabled(f_cv_runave)) {
-    calc_runave();
+    error_code |= calc_runave();
   }
 
   if (is_enabled(f_cv_corrfunc)) {
-    calc_acf();
+    error_code |= calc_acf();
   }
+
+  return error_code;
 }
 
 
 inline void history_add_value(size_t const           &history_length,
-                               std::list<colvarvalue> &history,
-                               colvarvalue const      &new_value)
+                              std::list<colvarvalue> &history,
+                              colvarvalue const      &new_value)
 {
   history.push_front(new_value);
   if (history.size() > history_length)
     history.pop_back();
 }
 
+
 inline void history_incr(std::list< std::list<colvarvalue> >           &history,
-                          std::list< std::list<colvarvalue> >::iterator &history_p)
+                         std::list< std::list<colvarvalue> >::iterator &history_p)
 {
   if ((++history_p) == history.end())
     history_p = history.begin();
@@ -2133,18 +2141,21 @@ inline void history_incr(std::list< std::list<colvarvalue> >           &history,
 int colvar::calc_acf()
 {
   // using here an acf_stride-long list of vectors for either
-  // coordinates(acf_x_history) or velocities (acf_v_history); each vector can
+  // coordinates (acf_x_history) or velocities (acf_v_history); each vector can
   // contain up to acf_length values, which are contiguous in memory
   // representation but separated by acf_stride in the time series;
   // the pointer to each vector is changed at every step
 
+  colvar const *cfcv = cvm::colvar_by_name(acf_colvar_name);
+  if (cfcv == NULL) {
+    return cvm::error("Error: collective variable \""+acf_colvar_name+
+                      "\" is not defined at this time.\n", INPUT_ERROR);
+  }
+
   if (acf_x_history.empty() && acf_v_history.empty()) {
 
     // first-step operations
 
-    colvar *cfcv = (acf_colvar_name.size() ?
-                    cvm::colvar_by_name(acf_colvar_name) :
-                    this);
     if (colvarvalue::check_types(cfcv->value(), value())) {
       cvm::error("Error: correlation function between \""+cfcv->name+
                  "\" and \""+this->name+"\" cannot be calculated, "
@@ -2153,7 +2164,8 @@ int colvar::calc_acf()
     }
     acf_nframes = 0;
 
-    cvm::log("Colvar \""+this->name+"\": initializing ACF calculation.\n");
+    cvm::log("Colvar \""+this->name+"\": initializing correlation function "
+             "calculation.\n");
 
     if (acf.size() < acf_length+1)
       acf.resize(acf_length+1, 0.0);
@@ -2182,41 +2194,31 @@ int colvar::calc_acf()
       break;
     }
 
-  } else {
-
-    colvar *cfcv = (acf_colvar_name.size() ?
-                    cvm::colvar_by_name(acf_colvar_name) :
-                    this);
+  } else if (cvm::step_relative() > prev_timestep) {
 
     switch (acf_type) {
 
     case acf_vel:
 
-      if (is_enabled(f_cv_fdiff_velocity)) {
-        // calc() should do this already, but this only happens in a
-        // simulation; better do it again in case a trajectory is
-        // being read
-        v_reported = v_fdiff = fdiff_velocity(x_old, cfcv->value());
-      }
-
       calc_vel_acf((*acf_v_history_p), cfcv->velocity());
-      // store this value in the history
-      history_add_value(acf_length+acf_offset, *acf_v_history_p, cfcv->velocity());
-      // if stride is larger than one, cycle among different histories
+      history_add_value(acf_length+acf_offset, *acf_v_history_p,
+                        cfcv->velocity());
       history_incr(acf_v_history, acf_v_history_p);
       break;
 
     case acf_coor:
 
       calc_coor_acf((*acf_x_history_p), cfcv->value());
-      history_add_value(acf_length+acf_offset, *acf_x_history_p, cfcv->value());
+      history_add_value(acf_length+acf_offset, *acf_x_history_p,
+                        cfcv->value());
       history_incr(acf_x_history, acf_x_history_p);
       break;
 
     case acf_p2coor:
 
       calc_p2coor_acf((*acf_x_history_p), cfcv->value());
-      history_add_value(acf_length+acf_offset, *acf_x_history_p, cfcv->value());
+      history_add_value(acf_length+acf_offset, *acf_x_history_p,
+                        cfcv->value());
       history_incr(acf_x_history, acf_x_history_p);
       break;
 
@@ -2225,18 +2227,14 @@ int colvar::calc_acf()
     }
   }
 
-  if (is_enabled(f_cv_fdiff_velocity)) {
-    // set it for the next step
-    x_old = x;
-  }
-  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
+  return COLVARS_OK;
 }
 
 
-int colvar::calc_vel_acf(std::list<colvarvalue> &v_list,
-                           colvarvalue const      &v)
+void colvar::calc_vel_acf(std::list<colvarvalue> &v_list,
+                          colvarvalue const      &v)
 {
-  // loop over stored velocities and add to the ACF, but only the
+  // loop over stored velocities and add to the ACF, but only if the
   // length is sufficient to hold an entire row of ACF values
   if (v_list.size() >= acf_length+acf_offset) {
     std::list<colvarvalue>::iterator  vs_i = v_list.begin();
@@ -2255,7 +2253,6 @@ int colvar::calc_vel_acf(std::list<colvarvalue> &v_list,
 
     acf_nframes++;
   }
-  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
 }
 
 
@@ -2280,7 +2277,7 @@ void colvar::calc_coor_acf(std::list<colvarvalue> &x_list,
 
 
 void colvar::calc_p2coor_acf(std::list<colvarvalue> &x_list,
-                              colvarvalue const      &x)
+                             colvarvalue const      &x)
 {
   // same as above but with second order Legendre polynomial instead
   // of just the scalar product
@@ -2301,20 +2298,46 @@ void colvar::calc_p2coor_acf(std::list<colvarvalue> &x_list,
 }
 
 
-void colvar::write_acf(std::ostream &os)
+int colvar::write_acf(std::ostream &os)
 {
-  if (!acf_nframes)
-    cvm::log("Warning: ACF was not calculated (insufficient frames).\n");
+  if (!acf_nframes) {
+    return COLVARS_OK;
+  }
+
   os.setf(std::ios::scientific, std::ios::floatfield);
-  os << "# Autocorrelation function for collective variable \""
-     << this->name << "\"\n";
-  // one frame is used for normalization, the statistical sample is
-  // hence decreased
-  os << "# nframes = " << (acf_normalize ?
-                           acf_nframes - 1 :
-                           acf_nframes) << "\n";
+  os << "# ";
+  switch (acf_type) {
+  case acf_vel:
+    os << "Velocity";
+    break;
+  case acf_coor:
+    os << "Coordinate";
+    break;
+  case acf_p2coor:
+    os << "Coordinate (2nd Legendre poly)";
+    break;
+  }
+
+  if (acf_colvar_name == name) {
+    os << " autocorrelation function for variable \""
+       << this->name << "\"\n";
+  } else {
+    os << " correlation function between variables \"" //
+       << this->name << "\" and \"" << acf_colvar_name << "\"\n";
+  }
+
+  os << "# Number of samples = ";
+  if (acf_normalize) {
+    os << (acf_nframes-1) << " (one DoF is used for normalization)\n";
+  } else {
+    os << acf_nframes << "\n";
+  }
+
+  os << "# " << cvm::wrap_string("step", cvm::it_width-2) << " "
+     << cvm::wrap_string("corrfunc(step)", cvm::cv_width) << "\n";
 
   cvm::real const acf_norm = acf.front() / cvm::real(acf_nframes);
+
   std::vector<cvm::real>::iterator acf_i;
   size_t it = acf_offset;
   for (acf_i = acf.begin(); acf_i != acf.end(); ++acf_i) {
@@ -2325,11 +2348,15 @@ void colvar::write_acf(std::ostream &os)
             (*acf_i)/(acf_norm * cvm::real(acf_nframes)) :
             (*acf_i)/(cvm::real(acf_nframes)) ) << "\n";
   }
+
+  return os.good() ? COLVARS_OK : FILE_ERROR;
 }
 
 
-void colvar::calc_runave()
+int colvar::calc_runave()
 {
+  int error_code = COLVARS_OK;
+
   if (x_history.empty()) {
 
     runave.type(value().type());
@@ -2348,10 +2375,29 @@ void colvar::calc_runave()
 
   } else {
 
-    if ( (cvm::step_relative() % runave_stride) == 0) {
+    if ( (cvm::step_relative() % runave_stride) == 0 &&
+         (cvm::step_relative() > prev_timestep) ) {
 
       if ((*x_history_p).size() >= runave_length-1) {
 
+        if (runave_os == NULL) {
+          if (runave_outfile.size() == 0) {
+            runave_outfile = std::string(cvm::output_prefix()+"."+
+                                         this->name+".runave.traj");
+          }
+
+          size_t const this_cv_width = x.output_width(cvm::cv_width);
+          cvm::proxy->backup_file(runave_outfile);
+          runave_os = cvm::proxy->output_stream(runave_outfile);
+          runave_os->setf(std::ios::scientific, std::ios::floatfield);
+          *runave_os << "# " << cvm::wrap_string("step", cvm::it_width-2)
+                     << "   "
+                     << cvm::wrap_string("running average", this_cv_width)
+                     << " "
+                     << cvm::wrap_string("running stddev", this_cv_width)
+                     << "\n";
+        }
+
         runave = x;
         std::list<colvarvalue>::iterator xs_i;
         for (xs_i = (*x_history_p).begin();
@@ -2370,7 +2416,7 @@ void colvar::calc_runave()
         runave_variance *= 1.0 / cvm::real(runave_length-1);
 
         *runave_os << std::setw(cvm::it_width) << cvm::step_relative()
-                   << "  "
+                   << "   "
                    << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
                    << runave << " "
                    << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
@@ -2381,6 +2427,7 @@ void colvar::calc_runave()
     }
   }
 
+  return error_code;
 }
 
 // Static members
index 989d551..a67749d 100644 (file)
@@ -291,6 +291,9 @@ public:
   /// \brief Calculate the colvar's value and related quantities
   int calc();
 
+  /// Carry out operations needed before next step is run
+  int end_of_step();
+
   /// \brief Calculate a subset of the colvar components (CVCs) currently active
   /// (default: all active CVCs)
   /// Note: both arguments refer to the sect of *active* CVCs, not all CVCs
@@ -410,8 +413,9 @@ public:
 
   /// Read the analysis tasks
   int parse_analysis(std::string const &conf);
+
   /// Perform analysis tasks
-  void analyze();
+  int analyze();
 
 
   /// Read the value from a collective variable trajectory file
@@ -489,23 +493,23 @@ protected:
   acf_type_e             acf_type;
 
   /// \brief Velocity ACF, scalar product between v(0) and v(t)
-  int calc_vel_acf(std::list<colvarvalue> &v_history,
-                     colvarvalue const      &v);
+  void calc_vel_acf(std::list<colvarvalue> &v_history,
+                    colvarvalue const      &v);
 
   /// \brief Coordinate ACF, scalar product between x(0) and x(t)
   /// (does not work with scalar numbers)
   void calc_coor_acf(std::list<colvarvalue> &x_history,
-                      colvarvalue const      &x);
+                     colvarvalue const      &x);
 
   /// \brief Coordinate ACF, second order Legendre polynomial between
   /// x(0) and x(t) (does not work with scalar numbers)
   void calc_p2coor_acf(std::list<colvarvalue> &x_history,
-                        colvarvalue const      &x);
+                       colvarvalue const      &x);
 
   /// Calculate the auto-correlation function (ACF)
   int calc_acf();
   /// Save the ACF to a file
-  void write_acf(std::ostream &os);
+  int write_acf(std::ostream &os);
 
   /// Length of running average series
   size_t         runave_length;
@@ -521,7 +525,7 @@ protected:
   cvm::real      runave_variance;
 
   /// Calculate the running average and its standard deviation
-  void calc_runave();
+  int calc_runave();
 
   /// If extended Lagrangian active: colvar energies (kinetic and harmonic potential)
   cvm::real kinetic_energy;
index 29620fb..9363fcd 100644 (file)
@@ -236,6 +236,12 @@ void colvarbias::communicate_forces()
 }
 
 
+int colvarbias::end_of_step()
+{
+  return COLVARS_OK;
+}
+
+
 int colvarbias::change_configuration(std::string const &conf)
 {
   cvm::error("Error: change_configuration() not implemented.\n",
index 083b9d7..391826e 100644 (file)
@@ -66,6 +66,9 @@ public:
   /// Send forces to the collective variables
   virtual void communicate_forces();
 
+  /// Carry out operations needed before next step is run
+  virtual int end_of_step();
+
   /// Load new configuration - force constant and/or centers only
   virtual int change_configuration(std::string const &conf);
 
index 771955a..e5edd92 100644 (file)
@@ -564,6 +564,8 @@ int colvarbias_abf::replica_share() {
   return COLVARS_OK;
 }
 
+
+
 void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool append)
 {
   std::string  samples_out_name = prefix + ".count";
@@ -572,10 +574,7 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
 
   std::ostream *samples_os =
     cvm::proxy->output_stream(samples_out_name, mode);
-  if (!samples_os) {
-    cvm::error("Error opening ABF samples file " + samples_out_name + " for writing");
-    return;
-  }
+  if (!samples_os) return;
   samples->write_multicol(*samples_os);
   cvm::proxy->close_output_stream(samples_out_name);
 
@@ -583,10 +582,7 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
   if (num_variables() > 2) {
     std::string  samples_dx_out_name = prefix + ".count.dx";
     std::ostream *samples_dx_os = cvm::proxy->output_stream(samples_dx_out_name, mode);
-    if (!samples_os) {
-      cvm::error("Error opening samples file " + samples_dx_out_name + " for writing");
-      return;
-    }
+    if (!samples_os) return;
     samples->write_opendx(*samples_dx_os);
     *samples_dx_os << std::endl;
     cvm::proxy->close_output_stream(samples_dx_out_name);
@@ -594,10 +590,7 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
 
   std::ostream *gradients_os =
     cvm::proxy->output_stream(gradients_out_name, mode);
-  if (!gradients_os) {
-    cvm::error("Error opening ABF gradient file " + gradients_out_name + " for writing");
-    return;
-  }
+  if (!gradients_os) return;
   gradients->write_multicol(*gradients_os);
   cvm::proxy->close_output_stream(gradients_out_name);
 
@@ -609,20 +602,14 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
 
     std::string  pmf_out_name = prefix + ".pmf";
     std::ostream *pmf_os = cvm::proxy->output_stream(pmf_out_name, mode);
-    if (!pmf_os) {
-      cvm::error("Error opening pmf file " + pmf_out_name + " for writing");
-      return;
-    }
+    if (!pmf_os) return;
     pmf->write_multicol(*pmf_os);
 
     // In dimension higher than 2, dx is easier to handle and visualize
     if (num_variables() > 2) {
       std::string  pmf_dx_out_name = prefix + ".pmf.dx";
       std::ostream *pmf_dx_os = cvm::proxy->output_stream(pmf_dx_out_name, mode);
-      if (!pmf_dx_os) {
-        cvm::error("Error opening pmf file " + pmf_dx_out_name + " for writing");
-        return;
-      }
+      if (!pmf_dx_os) return;
       pmf->write_opendx(*pmf_dx_os);
       *pmf_dx_os << std::endl;
       cvm::proxy->close_output_stream(pmf_dx_out_name);
@@ -639,10 +626,7 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
 
     std::ostream *z_samples_os =
       cvm::proxy->output_stream(z_samples_out_name, mode);
-    if (!z_samples_os) {
-      cvm::error("Error opening eABF z-histogram file " + z_samples_out_name + " for writing");
-      return;
-    }
+    if (!z_samples_os) return;
     z_samples->write_multicol(*z_samples_os);
     cvm::proxy->close_output_stream(z_samples_out_name);
 
@@ -651,10 +635,7 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
 
       std::ostream *z_gradients_os =
         cvm::proxy->output_stream(z_gradients_out_name, mode);
-      if (!z_gradients_os) {
-        cvm::error("Error opening eABF z-gradient file " + z_gradients_out_name + " for writing");
-        return;
-      }
+      if (!z_gradients_os) return;
       z_gradients->write_multicol(*z_gradients_os);
       cvm::proxy->close_output_stream(z_gradients_out_name);
     }
@@ -672,10 +653,7 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
 
     std::ostream *czar_gradients_os =
       cvm::proxy->output_stream(czar_gradients_out_name, mode);
-    if (!czar_gradients_os) {
-      cvm::error("Error opening CZAR gradient file " + czar_gradients_out_name + " for writing");
-      return;
-    }
+    if (!czar_gradients_os) return;
     czar_gradients->write_multicol(*czar_gradients_os);
     cvm::proxy->close_output_stream(czar_gradients_out_name);
 
@@ -688,20 +666,14 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
 
       std::string  czar_pmf_out_name = prefix + ".czar.pmf";
       std::ostream *czar_pmf_os = cvm::proxy->output_stream(czar_pmf_out_name, mode);
-      if (!czar_pmf_os) {
-        cvm::error("Error opening CZAR pmf file " + czar_pmf_out_name + " for writing");
-        return;
-      }
+      if (!czar_pmf_os) return;
       czar_pmf->write_multicol(*czar_pmf_os);
 
       // In dimension higher than 2, dx is easier to handle and visualize
       if (num_variables() > 2) {
         std::string  czar_pmf_dx_out_name = prefix + ".czar.pmf.dx";
         std::ostream *czar_pmf_dx_os = cvm::proxy->output_stream(czar_pmf_dx_out_name, mode);
-        if (!czar_pmf_dx_os) {
-          cvm::error("Error opening CZAR pmf file " + czar_pmf_dx_out_name + " for writing");
-          return;
-        }
+        if (!czar_pmf_dx_os) return;
         czar_pmf->write_opendx(*czar_pmf_dx_os);
         *czar_pmf_dx_os << std::endl;
         cvm::proxy->close_output_stream(czar_pmf_dx_out_name);
@@ -854,3 +826,9 @@ std::istream & colvarbias_abf::read_state_data(std::istream& is)
 
   return is;
 }
+
+int colvarbias_abf::write_output_files()
+{
+  write_gradients_samples(output_prefix);
+  return COLVARS_OK;
+}
index 0260401..52bf2df 100644 (file)
@@ -136,13 +136,13 @@ private:
 
   /// Write human-readable FE gradients and sample count, and DX file in dim > 2
   void write_gradients_samples(const std::string &prefix, bool append = false);
-  void write_last_gradients_samples(const std::string &prefix, bool append = false);
 
   /// Read human-readable FE gradients and sample count (if not using restart)
   void read_gradients_samples();
 
-  std::istream& read_state_data(std::istream&);
-  std::ostream& write_state_data(std::ostream&);
+  virtual std::istream& read_state_data(std::istream&);
+  virtual std::ostream& write_state_data(std::ostream&);
+  virtual int write_output_files();
 };
 
 #endif
index 80bd667..d20ee6e 100644 (file)
@@ -567,6 +567,9 @@ void colvardeps::init_cv_requires() {
     // Most features are available, so we set them so
     // and list exceptions below
    }
+
+  feature_states[f_cv_fdiff_velocity].available =
+    cvm::main()->proxy->simulation_running();
 }
 
 
index 87b08b1..d88a97a 100644 (file)
@@ -62,8 +62,6 @@ colvarmodule::colvarmodule(colvarproxy *proxy_in)
   use_scripted_forces = false;
   scripting_after_biases = false;
 
-  b_analysis = false;
-
   colvarmodule::debug_gradients_step_size = 1.0e-07;
 
   colvarmodule::rotation::monitor_crossings = false;
@@ -274,7 +272,12 @@ int colvarmodule::parse_global_params(std::string const &conf)
     }
   }
 
-  parse->get_keyval(conf, "analysis", b_analysis, b_analysis);
+  bool b_analysis = true;
+  if (parse->get_keyval(conf, "analysis", b_analysis, true,
+                        colvarparse::parse_silent)) {
+    cvm::log("Warning: keyword \"analysis\" is deprecated: it is now set "
+             "to true; individual analyses are performed only if requested.");
+  }
 
   parse->get_keyval(conf, "debugGradientsStepSize", debug_gradients_step_size,
                     debug_gradients_step_size,
@@ -715,9 +718,7 @@ int colvarmodule::calc()
   error_code |= calc_biases();
   error_code |= update_colvar_forces();
 
-  if (cvm::b_analysis) {
-    error_code |= analyze();
-  }
+  error_code |= analyze();
 
   // write trajectory files, if needed
   if (cv_traj_freq && cv_traj_name.size()) {
@@ -736,6 +737,8 @@ int colvarmodule::calc()
     write_output_files();
   }
 
+  error_code |= end_of_step();
+
   return error_code;
 }
 
@@ -1056,6 +1059,33 @@ int colvarmodule::analyze()
 }
 
 
+int colvarmodule::end_of_step()
+{
+  if (cvm::debug()) {
+    cvm::log("colvarmodule::end_of_step(), step = "+cvm::to_str(it)+".\n");
+  }
+
+  for (std::vector<colvar *>::iterator cvi = variables_active()->begin();
+       cvi != variables_active()->end();
+       cvi++) {
+    cvm::increase_depth();
+    (*cvi)->end_of_step();
+    cvm::decrease_depth();
+  }
+
+  // perform bias-specific analysis
+  for (std::vector<colvarbias *>::iterator bi = biases.begin();
+       bi != biases.end();
+       bi++) {
+    cvm::increase_depth();
+    (*bi)->end_of_step();
+    cvm::decrease_depth();
+  }
+
+  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
+}
+
+
 int colvarmodule::setup()
 {
   if (this->size() == 0) return cvm::get_error();
@@ -1895,7 +1925,6 @@ long      colvarmodule::it = 0;
 long      colvarmodule::it_restart = 0;
 size_t    colvarmodule::restart_out_freq = 0;
 size_t    colvarmodule::cv_traj_freq = 0;
-bool      colvarmodule::b_analysis = false;
 bool      colvarmodule::use_scripted_forces = false;
 bool      colvarmodule::scripting_after_biases = true;
 
index 3d93798..99b7976 100644 (file)
@@ -422,6 +422,9 @@ public:
   /// Perform analysis
   int analyze();
 
+  /// Carry out operations needed before next step is run
+  int end_of_step();
+
   /// \brief Read a collective variable trajectory (post-processing
   /// only, not called at runtime)
   int read_traj(char const *traj_filename,
@@ -546,9 +549,6 @@ public:
   /// Frequency for collective variables trajectory output
   static size_t cv_traj_freq;
 
-  /// \brief True if only analysis is performed and not a run
-  static bool   b_analysis;
-
   /// Frequency for saving output restarts
   static size_t restart_out_freq;
   /// Output restart file name
index a37d264..bd84d07 100644 (file)
@@ -1,5 +1,5 @@
 #ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2018-10-16"
+#define COLVARS_VERSION "2018-11-16"
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
index c9fe049..3ecd660 100644 (file)
@@ -137,6 +137,10 @@ int colvarscript::run(int objc, unsigned char *const objv[])
 
   if (cmd == "update") {
     error_code |= proxy->update_input();
+    if (error_code) {
+      result += "Error updating the Colvars module.\n";
+      return error_code;
+    }
     error_code |= colvars->calc();
     error_code |= proxy->update_output();
     if (error_code) {
index db8c83a..ffff828 100644 (file)
@@ -71,8 +71,6 @@ Detailed explanations of the design of the Colvars module are provided in refere
 
 
 \cvsec{A crash course}{sec:colvars_crash_course}
-% FIXME: the example uses NAMD syntax to select Calphas, which may confuse
-% LAMMPS users
 
 Suppose that we want to run a steered MD experiment where a small molecule is pulled away from a protein binding site.
 In Colvars terms, this is done by applying a moving restraint to the distance between the two objects.
@@ -81,30 +79,35 @@ The configuration will contain two blocks, one defining the distance variable (s
 \bigskip
 % verbatim can't appear within commands
 {\noindent\ttfamily
-colvar \{\\
+\cvlammpsonly{indexFile index.ndx\\}colvar \{\\
 \-~~name dist\\
 \-~~distance \{\\
 \-~~~~group1 \{ atomNumbersRange 42-55 \}\\
-\-~~~~group2 \{\\
+\cvnamebasedonly{\-~~~~group2 \{\\
 \-~~~~~~psfSegID  PR\\
-\-~~~~~~atomNameResidueRange  CA 15-30 \}\\
-\-~~~~\}\\
+\-~~~~~~atomNameResidueRange  CA 15-30\\
+\-~~~~\}\\}
+\cvlammpsonly{\-~~~~group2 \{ indexGroup C-alpha\_15-30 \}\\}
 \-~~\}\\
 \}\\
 \\
 harmonic \{\\
 \-~~colvars dist\\
 \-~~forceConstant 20.0\\
-\-~~centers 4.~~~~~~~~~~\# initial distance\\
-\-~~targetCenters 15.~~~\# final distance\\
+\-~~centers 4.0~~~~~~~~~\# initial distance\\
+\-~~targetCenters 15.0~~\# final distance\\
 \-~~targetNumSteps 500000\\
 \}\\
 }
 
-Reading this input in plain English: the variable here named \emph{dist} consists in a distance function between the centers of two groups: the ligand (atoms 42 to 55) and the alpha carbon atoms (CA) of residues 15 to 30 in the protein (segment name PR).
-The atom selection syntax is detailed in section \ref{sec:colvar_atom_groups}.
+Reading this input in plain English: the variable here named \emph{dist} consists in a distance function between the centers of two groups: the ligand (atoms 42 to 55) and the $\alpha$-carbon atoms of residues 15 to 30 in the protein \cvnamebasedonly{(segment name PR)}\cvlammpsonly{(selected from the index group ``\texttt{C-alpha\_15-30}'')}.
+To the ``\emph{dist}'' variable, we apply a \texttt{harmonic} potential of force constant 20~\cvnamdonly{kcal/mol}\cvvmdonly{kcal/mol}\cvlammpsonly{(energy unit)}/\cvnamdonly{\AA}\cvvmdonly{\AA}\cvlammpsonly{(length unit)}$^2$, initially centered around a value of 4~\cvnamdonly{\AA}\cvvmdonly{\AA}\cvlammpsonly{(length unit)}, which will increase to 15~\cvnamdonly{\AA}\cvvmdonly{\AA}\cvlammpsonly{(length unit)} over 500,000 simulation steps.
+
+The atom selection keywords is detailed in section \ref{sec:colvar_atom_groups}.
+\cvlammpsonly{If the selection is too complex to implement only via internal keywords, an external index file may be created following the NDX format used in GROMACS:\\
+\url{http://manual.gromacs.org/documentation/current/user-guide/file-formats.html}\\
+or by using the \texttt{group2ndx} LAMMPS command.}
 
-To the ``\emph{dist}'' variable, we apply a \texttt{harmonic} potential of force constant 20~\cvnamdonly{kcal/mol}\cvvmdonly{kcal/mol}\cvlammpsonly{energy units}/\AA$^2$, initially centered around a value of 4~\AA, which will increase to 15~\AA{} over 500,000 simulation steps.
 
 
 \cvsec{General parameters}{sec:colvarmodule}
@@ -430,18 +433,6 @@ The following keywords are available in the global context of the Colvars config
     boolean}{%
     \texttt{on}}{%
     If this flag is enabled (default), SMP parallelism over threads will be used to compute variables and biases, provided that this is supported by the \MDENGINE{} build in use.}
-
-\item %
-  \keydef
-    {analysis}{%
-    global}{%
-    Turn on run-time statistical analysis}{%
-    boolean}{%
-    \texttt{off}}{%
-    If this flag is enabled, each colvar is instructed to perform
-    whatever run-time statistical analysis it is configured to, such as
-    correlation functions, or running averages and standard deviations.
-    See section~\ref{sec:colvar_acf} for details.}
     
 \end{itemize}
 
@@ -454,7 +445,7 @@ The corresponding configuration is given below. The options within the \texttt{c
   \centering
 \cvnamdugonly{\includegraphics[width=12cm]{figures/colvars_diagram}}
 \cvvmdugonly{\includegraphics[width=12cm]{pictures/colvars_diagram}}
-\cvrefmanonly{\includegraphics[width=12cm]{colvars_diagram}}
+\cvrefmanonly{\ifdefined\HCode{\HCode{<img class="diagram" src="colvars_diagram.png" alt="Colvars diagram">}}\else\includegraphics[width=12cm]{colvars_diagram}\fi}
   \caption[Graphical representation of a Colvars configuration.]{Graphical representation of a Colvars configuration\cvlammpsonly{ \textbf{(note:} \emph{currently, the $\alpha$-helical content colvar is unavailable in LAMMPS)}}.
     The colvar called ``$d$'' is defined as the difference between two distances: the first distance ($d_{1}$) is taken between the center of mass of atoms 1 and 2 and that of atoms 3 to 5, the second ($d_{2}$) between atom 7 and the center of mass of atoms 8 to 10.
 The difference $d = d_{1} - d_{2}$ is obtained by multiplying the two by a coefficient $C = +1$ or $C = -1$, respectively.
@@ -2585,15 +2576,14 @@ to perform eABF simulations (\ref{sec:eABF}).
 
 \cvsubsec{Statistical analysis}{sec:colvar_acf}
 
-When the global keyword \texttt{analysis} is defined in the
-configuration file, run-time calculations of statistical properties for
-individual colvars can be performed.  At the moment, several types of
-time correlation functions, running averages and running standard
-deviations are available.
+Run-time calculations of statistical properties that depend explicitly on time can be performed for individual collective variables.
+Currently, several types of time correlation functions, running averages and running standard deviations are implemented.
+For run-time computation of histograms, please see the histogram bias (\ref{sec:colvarbias_histogram}).
 
 \begin{itemize}
 
 \item %
+  \labelkey{colvar|corrFunc}
   \keydef
     {corrFunc}{%
     \texttt{colvar}}{%
@@ -2679,10 +2669,11 @@ deviations are available.
     \texttt{colvar}}{%
     Output file for the time correlation function}{%
     UNIX filename}{%
-    \texttt{$<$name$>$.corrfunc.dat}}{%
+    \outputName\texttt{.$<$name$>$.corrfunc.dat}}{%
     The time correlation function is saved in this file.}
 
 \item %
+  \labelkey{colvar|runAve}
   \keydef
     {runAve}{%
     \texttt{colvar}}{%
@@ -2716,7 +2707,7 @@ deviations are available.
     \texttt{colvar}}{%
     Output file for the running average and standard deviation}{%
     UNIX filename}{%
-    \texttt{$<$name$>$.runave.dat}}{%
+    \outputName\texttt{.$<$name$>$.runave.traj}}{%
     The running average and standard deviation are saved in this file.}
 
 \end{itemize}
@@ -3170,7 +3161,7 @@ In addition, restraint biases (\ref{sec:colvarbias_harmonic}, \ref{sec:colvarbia
 \begin{itemize}
 
 \item %
-  \labelkey{writeTIPMF}{colvarbias|writeTIPMF}
+  \labelkey{colvarbias|writeTIPMF}
   \keydef
     {writeTIPMF}{%
     colvar bias}{%
@@ -3184,7 +3175,7 @@ In addition, restraint biases (\ref{sec:colvarbias_harmonic}, \ref{sec:colvarbia
 
 
 \item %
-  \labelkey{writeTISamples}{colvarbias|writeTISamples}
+  \labelkey{colvarbias|writeTISamples}
   \keydef
     {writeTISamples}{%
     colvar bias}{%
@@ -3232,7 +3223,7 @@ $\bm{F}_\xi$ exerted on $\bm{\xi}$, taken over an iso-$\bm{\xi}$ surface:
 
 \begin{equation}
   \label{eq:gradient}
-  \bm{\nabla}_\xi A(\bm{\xi}) = \left\langle -\bm{F}_\xi \right\rangle_\bm{\xi}
+  \bm{\nabla}_\xi A(\bm{\xi}) = \left\langle -\bm{F}_\xi \right\rangle_{\bm{\xi}}
 \end{equation}
 
 Several formulae that take the form of~(\ref{eq:gradient}) have been
@@ -3242,7 +3233,7 @@ originating in a work by Ruiz-Montero et al.~\cite{Ruiz-Montero1997},
 generalized by den Otter~\cite{denOtter2000} and extended to multiple
 variables by Ciccotti et al.~\cite{Ciccotti2005}.  Consider a system
 subject to constraints of the form $\sigma_{k}(\vx) = 0$.  Let
-($\bm{v}_{i})_{i\in[1,n]}$ be arbitrarily chosen vector fields
+$(\bm{v}_{i})_{i\in[1,n]}$ be arbitrarily chosen vector fields
 ($\mathbb{R}^{3N}\rightarrow\mathbb{R}^{3N}$) verifying, for all $i$,
 $j$, and $k$:
 
@@ -3258,7 +3249,7 @@ then the following holds~\cite{Ciccotti2005}:
 \begin{equation}
 \label{eq:gradient_vector}
 \frac{\partial A}{\partial \xi_{i}} = \left\langle \bm{v}_{i} \cdot \gradx V
-- k_B T \gradx \cdot \bm{v}_{i} \right\rangle_\bm{\xi}
+- k_B T \gradx \cdot \bm{v}_{i} \right\rangle_{\bm{\xi}}
 \end{equation}
 
 where $V$ is the potential energy function.
@@ -3602,7 +3593,7 @@ $(k/2) (z - \lambda)^2$.
 Under eABF dynamics, the adaptive bias on $\lambda$ is
 the running estimate of the average spring force:
 \begin{equation}
-F^\mathrm{bias}(\lambda^*) = \left\langle k(\lambda - z) \right\rangle_{\lambda^*}
+F^{\mathrm{bias}}(\lambda^{*}) = \left\langle k(\lambda{} - z) \right\rangle_{\lambda^{*}}
 \end{equation}
 where the angle brackets indicate a canonical average conditioned by $\lambda=\lambda^*$.
 At long simulation times, eABF produces a flat histogram of the extended variable $\lambda$,
@@ -4870,6 +4861,9 @@ please be aware that \emph{several of these tutorials are not actively maintaine
 \item \textbf{Colvars version 2017-01-09 or later \cvnamdonly{(NAMD version 2.13b1 or later)}\cvvmdonly{(VMD version 1.9.4 or later)}\cvlammpsonly{(LAMMPS version 10Mar2017 or later)}.}\\
   A new type of restraint, \texttt{harmonicWalls} (see~\ref{sec:colvarbias_harmonic_walls}), replaces and improves upon the legacy keywords \texttt{lowerWall} and \texttt{upperWall}: these are still supported as short-hands.
 
+\item \textbf{Colvars version 2018-11-15 or later \cvnamdonly{(NAMD version 2.14b1 or later)}\cvvmdonly{(VMD version 1.9.4 or later)}\cvlammpsonly{(LAMMPS version 23Nov2018 or later)}.}\\
+  The global \texttt{analysis} keyword has been discontinued: specific analysis tasks are controlled directly by the keywords \refkey{corrFunc}{colvar|corrFunc} and \refkey{runAve}{colvar|runAve}, which continue to remain \texttt{off} by default.
+
 \item \textbf{Deprecation warning for calculations including wall potentials.}\\
   The legacy keywords \texttt{lowerWall} and \texttt{upperWall} will stop having default values and will need to be set explicitly (preferably as part of the \texttt{harmonicWalls} restraint).
   When using an ABF bias, it is recommended to set the two walls equal to \refkey{lowerBoundary}{colvar|lowerBoundary} and \refkey{upperBoundary}{colvar|upperBoundary}, respectively.