Update Colvars to version 2018-12-14
[namd.git] / colvars / src / colvar.cpp
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