From 5a0a8f1335ff9459875287b8db6f5574ed218f0e Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Mon, 9 Mar 2015 14:50:53 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@13207
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 lib/colvars/colvar.cpp               | 31 ++++++++++----------
 lib/colvars/colvarbias.cpp           |  2 +-
 lib/colvars/colvarbias_abf.cpp       | 16 ++++++-----
 lib/colvars/colvarbias_restraint.cpp |  3 +-
 lib/colvars/colvarmodule.cpp         | 18 ++++++++----
 lib/colvars/colvarmodule.h           |  8 +++---
 lib/colvars/colvarparse.h            | 13 +++++++++
 lib/colvars/colvarscript.cpp         | 42 ++++++++++++++++++++++------
 lib/colvars/colvartypes.cpp          | 18 ++++++++----
 lib/colvars/colvartypes.h            |  2 ++
 lib/colvars/colvarvalue.h            | 26 +++++++++--------
 11 files changed, 118 insertions(+), 61 deletions(-)

diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp
index 8f0850c0ff..80c04f9f8b 100644
--- a/lib/colvars/colvar.cpp
+++ b/lib/colvars/colvar.cpp
@@ -11,6 +11,7 @@
 
 
 colvar::colvar(std::string const &conf)
+  : colvarparse(conf)
 {
   size_t i, j;
   cvm::log("Initializing a new collective variable.\n");
@@ -151,6 +152,7 @@ colvar::colvar(std::string const &conf)
     cvm::error("Error: no valid components were provided "
                       "for this collective variable.\n",
               INPUT_ERROR);
+    return;
   }
 
   cvm::log("All components initialized.\n");
@@ -232,7 +234,6 @@ colvar::colvar(std::string const &conf)
     x.type(cvc_value);
     x_reported.type(cvc_value);
   }
-
   // If using scripted biases, any colvar may receive bias forces
   // and will need its gradient
   if (cvm::scripted_forces()) {
@@ -319,6 +320,10 @@ colvar::colvar(std::string const &conf)
     }
   }
 
+  // at this point, the colvar's type is defined
+  f.type(value());
+  fb.type(value());
+
   get_keyval(conf, "width", width, 1.0);
   if (width <= 0.0) {
     cvm::error("Error: \"width\" must be positive.\n", INPUT_ERROR);
@@ -374,10 +379,10 @@ colvar::colvar(std::string const &conf)
   if (tasks[task_lower_wall] && tasks[task_upper_wall]) {
     if (lower_wall >= upper_wall) {
       cvm::error("Error: the upper wall, "+
-                        cvm::to_str(upper_wall)+
-                        ", is not higher than the lower wall, "+
-                        cvm::to_str(lower_wall)+".\n",
-                INPUT_ERROR);
+                 cvm::to_str(upper_wall)+
+                 ", is not higher than the lower wall, "+
+                 cvm::to_str(lower_wall)+".\n",
+                 INPUT_ERROR);
     }
 
     if (dist2(lower_wall, upper_wall) < 1.0E-12) {
@@ -391,15 +396,15 @@ colvar::colvar(std::string const &conf)
   get_keyval(conf, "expandBoundaries", expand_boundaries, false);
   if (expand_boundaries && periodic_boundaries()) {
     cvm::error("Error: trying to expand boundaries that already "
-                      "cover a whole period of a periodic colvar.\n",
-              INPUT_ERROR);
+               "cover a whole period of a periodic colvar.\n",
+               INPUT_ERROR);
   }
   if (expand_boundaries && hard_lower_boundary && hard_upper_boundary) {
     cvm::error("Error: inconsistent configuration "
-                      "(trying to expand boundaries with both "
-                      "hardLowerBoundary and hardUpperBoundary enabled).\n",
-              INPUT_ERROR);
-}
+               "(trying to expand boundaries with both "
+               "hardLowerBoundary and hardUpperBoundary enabled).\n",
+               INPUT_ERROR);
+  }
 
   {
     bool b_extended_lagrangian;
@@ -722,11 +727,7 @@ int colvar::enable(colvar::task const &t)
   case task_langevin:
   case task_output_energy:
   case task_scripted:
-    break;
-
   case task_gradients:
-    f.type(value());
-    fb.type(value());
     break;
 
   case task_collect_gradients:
diff --git a/lib/colvars/colvarbias.cpp b/lib/colvars/colvarbias.cpp
index 8f1d85b242..5853e22127 100644
--- a/lib/colvars/colvarbias.cpp
+++ b/lib/colvars/colvarbias.cpp
@@ -6,7 +6,7 @@
 
 
 colvarbias::colvarbias(std::string const &conf, char const *key)
-  : colvarparse(), has_data(false)
+  : colvarparse(conf), bias_energy(0.), has_data(false)
 {
   cvm::log("Initializing a new \""+std::string(key)+"\" instance.\n");
 
diff --git a/lib/colvars/colvarbias_abf.cpp b/lib/colvars/colvarbias_abf.cpp
index d9a56794fd..10490152d7 100644
--- a/lib/colvars/colvarbias_abf.cpp
+++ b/lib/colvars/colvarbias_abf.cpp
@@ -190,13 +190,6 @@ cvm::real colvarbias_abf::update()
     // initialization stuff (file operations relying on n_abf_biases
     // compute current value of colvars
 
-    if ( cvm::n_abf_biases == 1 && cvm::n_meta_biases == 0 ) {
-      // This is the only ABF bias
-      output_prefix = cvm::output_prefix;
-    } else {
-      output_prefix = cvm::output_prefix + "." + this->name;
-    }
-
     for (size_t i=0; i<colvars.size(); i++) {
       bin[i] = samples->current_bin_scalar(i);
     }
@@ -260,10 +253,19 @@ cvm::real colvarbias_abf::update()
     }
   }
 
+  // update the output prefix; TODO: move later to setup_output() function
+  if ( cvm::n_abf_biases == 1 && cvm::n_meta_biases == 0 ) {
+    // This is the only ABF bias
+    output_prefix = cvm::output_prefix;
+  } else {
+    output_prefix = cvm::output_prefix + "." + this->name;
+  }
+
   if (output_freq && (cvm::step_absolute() % output_freq) == 0) {
     if (cvm::debug()) cvm::log("ABF bias trying to write gradients and samples to disk");
     write_gradients_samples(output_prefix);
   }
+
   if (b_history_files && (cvm::step_absolute() % history_freq) == 0) {
     cvm::log("ABFHISTORYFILE "+cvm::to_str(cvm::step_absolute()));
     // file already exists iff cvm::step_relative() > 0
diff --git a/lib/colvars/colvarbias_restraint.cpp b/lib/colvars/colvarbias_restraint.cpp
index 76586291e2..5ac7b1318f 100644
--- a/lib/colvars/colvarbias_restraint.cpp
+++ b/lib/colvars/colvarbias_restraint.cpp
@@ -153,8 +153,7 @@ cvm::real colvarbias_restraint::energy_difference(std::string const &conf)
   }
   if (get_keyval(conf, "centers", alt_colvar_centers, colvar_centers)) {
     for (i = 0; i < colvars.size(); i++) {
-      colvar_centers[i].type(colvars[i]->value());
-      colvar_centers[i].apply_constraints();
+      alt_colvar_centers[i].apply_constraints();
     }
   }
 
diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp
index f737a14fc7..372cc3b402 100644
--- a/lib/colvars/colvarmodule.cpp
+++ b/lib/colvars/colvarmodule.cpp
@@ -41,7 +41,10 @@ colvarmodule::colvarmodule(colvarproxy *proxy_in)
   colvarmodule::use_scripted_forces = false;
 
   colvarmodule::b_analysis = false;
+
   colvarmodule::debug_gradients_step_size = 1.0e-07;
+
+  colvarmodule::rotation::monitor_crossings = false;
   colvarmodule::rotation::crossing_threshold = 1.0e-02;
 
   colvarmodule::cv_traj_freq = 100;
@@ -52,7 +55,7 @@ colvarmodule::colvarmodule(colvarproxy *proxy_in)
 }
 
 
-int colvarmodule::config_file(char const  *config_filename)
+int colvarmodule::read_config_file(char const  *config_filename)
 {
   cvm::log(cvm::line_marker);
   cvm::log("Reading new configuration from file \""+
@@ -75,11 +78,11 @@ int colvarmodule::config_file(char const  *config_filename)
   }
   config_s.close();
 
-  return config(conf);
+  return parse_config(conf);
 }
 
 
-int colvarmodule::config_string(std::string const &config_str)
+int colvarmodule::read_config_string(std::string const &config_str)
 {
   cvm::log(cvm::line_marker);
   cvm::log("Reading new configuration:\n");
@@ -91,10 +94,10 @@ int colvarmodule::config_string(std::string const &config_str)
   while (colvarparse::getline_nocomments(config_s, line)) {
     conf.append(line+"\n");
   }
-  return config(conf);
+  return parse_config(conf);
 }
 
-int colvarmodule::config(std::string &conf)
+int colvarmodule::parse_config(std::string &conf)
 {
   int error_code = 0;
 
@@ -156,6 +159,10 @@ int colvarmodule::parse_global_params(std::string const &conf)
                     debug_gradients_step_size,
                     colvarparse::parse_silent);
 
+  parse->get_keyval(conf, "monitorEigenvalueCrossing",
+                    colvarmodule::rotation::monitor_crossings,
+                    colvarmodule::rotation::monitor_crossings,
+                    colvarparse::parse_silent);
   parse->get_keyval(conf, "eigenvalueCrossingThreshold",
                     colvarmodule::rotation::crossing_threshold,
                     colvarmodule::rotation::crossing_threshold,
@@ -1196,7 +1203,6 @@ size_t    colvarmodule::restart_out_freq = 0;
 size_t    colvarmodule::cv_traj_freq = 0;
 size_t    colvarmodule::depth = 0;
 bool      colvarmodule::b_analysis = false;
-cvm::real colvarmodule::rotation::crossing_threshold = 1.0E-04;
 std::list<std::string> colvarmodule::index_group_names;
 std::list<std::vector<int> > colvarmodule::index_groups;
 bool     colvarmodule::use_scripted_forces = false;
diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h
index d21b4a2885..d2fcb998de 100644
--- a/lib/colvars/colvarmodule.h
+++ b/lib/colvars/colvarmodule.h
@@ -4,7 +4,7 @@
 #define COLVARMODULE_H
 
 #ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2015-02-04"
+#define COLVARS_VERSION "2015-02-27"
 #endif
 
 #ifndef COLVARS_DEBUG
@@ -200,14 +200,14 @@ public:
   int reset();
 
   /// Open a config file, load its contents, and pass it to config_string()
-  int config_file(char const *config_file_name);
+  int read_config_file(char const *config_file_name);
 
   /// \brief Parse a config string assuming it is a complete configuration
   /// (i.e. calling all parse functions)
-  int config_string(std::string const &conf);
+  int read_config_string(std::string const &conf);
 
   /// \brief Parse a "clean" config string (no comments)
-  int config(std::string &conf);
+  int parse_config(std::string &conf);
 
 
   // Parse functions (setup internal data based on a string)
diff --git a/lib/colvars/colvarparse.h b/lib/colvars/colvarparse.h
index 2c7b8e000e..5f7ba26ba3 100644
--- a/lib/colvars/colvarparse.h
+++ b/lib/colvars/colvarparse.h
@@ -48,12 +48,25 @@ protected:
   /// \brief Remove all the values from the config string
   void strip_values(std::string &conf);
 
+  /// \brief Configuration string of the object
+  std::string config_string;
+
 public:
 
   inline colvarparse()
     : save_delimiters(true)
   {}
 
+  /// Constructor that stores the object's config string
+  inline colvarparse(const std::string& conf)
+    : save_delimiters(true), config_string(conf)
+  {}
+
+  inline const std::string& get_config()
+  {
+      return config_string;
+  }
+
   /// How a keyword is parsed in a string
   enum Parse_Mode {
     /// \brief(default) Read the first instance of a keyword (if
diff --git a/lib/colvars/colvarscript.cpp b/lib/colvars/colvarscript.cpp
index 0c3a90efbf..037f284b90 100644
--- a/lib/colvars/colvarscript.cpp
+++ b/lib/colvars/colvarscript.cpp
@@ -32,6 +32,7 @@ Managing the colvars module:\n\
   config <string>             -- read configuration from the given string\n\
   reset                       -- delete all internal configuration\n\
   delete                      -- delete this colvars module instance\n\
+  version                     -- return version of colvars code\n\
   \n\
 Input and output:\n\
   list                        -- return a list of all variables\n\
@@ -49,15 +50,18 @@ Input and output:\n\
 
   result += "\n\
 Accessing collective variables:\n\
-  colvar <name> value         -- return the current value of the colvar <name>\n\
-  colvar <name> update        -- recalculate the colvar <name>\n\
-  colvar <name> delete        -- delete the colvar <name>\n\
-  colvar <name> addforce <F>  -- apply given force on <name>\n\
+  colvar <name> value         -- return the current value of colvar <name>\n\
+  colvar <name> update        -- recalculate colvar <name>\n\
+  colvar <name> type          -- return the type of colvar <name>\n\
+  colvar <name> delete        -- delete colvar <name>\n\
+  colvar <name> addforce <F>  -- apply given force on colvar <name>\n\
+  colvar <name> getconfig     -- return config string of colvar <name>\n\
 \n\
 Accessing biases:\n\
-  bias <name> energy          -- return the current energy of the bias <name>\n\
-  bias <name> update          -- recalculate the bias <name>\n\
-  bias <name> delete          -- delete the bias <name>\n\
+  bias <name> energy          -- return the current energy of bias <name>\n\
+  bias <name> update          -- recalculate bias <name>\n\
+  bias <name> delete          -- delete bias <name>\n\
+  bias <name> getconfig       -- return config string of bias <name>\n\
 \n\
 ";
     return COLVARSCRIPT_OK;
@@ -73,6 +77,11 @@ Accessing biases:\n\
     return proc_bias(argc-1, &(argv[1]));
   }
 
+  if (cmd == "version") {
+    result = COLVARS_VERSION;
+    return COLVARSCRIPT_OK;
+  }
+
   if (cmd == "reset") {
     /// Delete every child object
     colvars->reset();
@@ -119,7 +128,7 @@ Accessing biases:\n\
       result = "Missing arguments";
       return COLVARSCRIPT_ERROR;
     }
-    if (colvars->config_file(argv[2]) == COLVARS_OK) {
+    if (colvars->read_config_file(argv[2]) == COLVARS_OK) {
       return COLVARSCRIPT_OK;
     } else {
       return COLVARSCRIPT_ERROR;
@@ -133,7 +142,7 @@ Accessing biases:\n\
       return COLVARSCRIPT_ERROR;
     }
     std::string conf = argv[2];
-    if (colvars->config_string(conf) == COLVARS_OK) {
+    if (colvars->read_config_string(conf) == COLVARS_OK) {
       return COLVARSCRIPT_OK;
     } else {
       return COLVARSCRIPT_ERROR;
@@ -221,6 +230,11 @@ int colvarscript::proc_colvar(int argc, char const *argv[]) {
     return COLVARSCRIPT_OK;
   }
 
+  if (subcmd == "type") {
+    result = cv->value().type_desc(cv->value().value_type);
+    return COLVARSCRIPT_OK;
+  }
+
   if (subcmd == "update") {
     cv->calc();
     cv->update();
@@ -240,6 +254,11 @@ int colvarscript::proc_colvar(int argc, char const *argv[]) {
     return COLVARSCRIPT_OK;
   }
 
+  if (subcmd == "getconfig") {
+    result = cv->get_config();
+    return COLVARSCRIPT_OK;
+  }
+
   if (subcmd == "addforce") {
     if (argc < 4) {
       result = "addforce: missing parameter: force value";
@@ -290,6 +309,11 @@ int colvarscript::proc_bias(int argc, char const *argv[]) {
     return COLVARSCRIPT_OK;
   }
 
+  if (subcmd == "getconfig") {
+    result = b->get_config();
+    return COLVARSCRIPT_OK;
+  }
+
   // Subcommands for MW ABF
   if (subcmd == "bin") {
     int r = b->current_bin();
diff --git a/lib/colvars/colvartypes.cpp b/lib/colvars/colvartypes.cpp
index 2ec448bd14..8e0aa57dfc 100644
--- a/lib/colvars/colvartypes.cpp
+++ b/lib/colvars/colvartypes.cpp
@@ -8,6 +8,10 @@
 #include "colvarparse.h"
 
 
+bool      colvarmodule::rotation::monitor_crossings = false;
+cvm::real colvarmodule::rotation::crossing_threshold = 1.0E-02;
+
+
 std::string cvm::rvector::to_simple_string() const
 {
   std::ostringstream os;
@@ -329,14 +333,16 @@ void colvarmodule::rotation::calc_optimal_rotation(std::vector<cvm::atom_pos> co
   lambda = L0;
   q = Q0;
 
-  if (q_old.norm2() > 0.0) {
-    q.match(q_old);
-    if (q_old.inner(q) < (1.0 - crossing_threshold)) {
-      cvm::log("Warning: one molecular orientation has changed by more than "+
-               cvm::to_str(crossing_threshold)+": discontinuous rotation ?\n");
+  if (cvm::rotation::monitor_crossings) {
+    if (q_old.norm2() > 0.0) {
+      q.match(q_old);
+      if (q_old.inner(q) < (1.0 - crossing_threshold)) {
+        cvm::log("Warning: one molecular orientation has changed by more than "+
+                 cvm::to_str(crossing_threshold)+": discontinuous rotation ?\n");
+      }
     }
+    q_old = q;
   }
-  q_old = q;
 
   if (cvm::debug()) {
     if (b_debug_gradients) {
diff --git a/lib/colvars/colvartypes.h b/lib/colvars/colvartypes.h
index aeb2a9facc..fcc26c0b98 100644
--- a/lib/colvars/colvartypes.h
+++ b/lib/colvars/colvartypes.h
@@ -1531,6 +1531,8 @@ public:
     }
   }
 
+  /// \brief Whether to test for eigenvalue crossing
+  static bool monitor_crossings;
   /// \brief Threshold for the eigenvalue crossing test
   static cvm::real crossing_threshold;
 
diff --git a/lib/colvars/colvarvalue.h b/lib/colvars/colvarvalue.h
index f91cb840f1..24af824081 100644
--- a/lib/colvars/colvarvalue.h
+++ b/lib/colvars/colvarvalue.h
@@ -171,27 +171,31 @@ public:
   /// Set the type explicitly
   inline void type(Type const &vti)
   {
-    // reset the value based on the previous type
-    reset();
-    if ((value_type == type_vector) && (vti != type_vector)) {
-      vector1d_value.resize(0);
+    if (vti != value_type) {
+      // reset the value based on the previous type
+      reset();
+      if ((value_type == type_vector) && (vti != type_vector)) {
+        vector1d_value.resize(0);
+      }
+      value_type = vti;
     }
-    value_type = vti;
   }
 
   /// Set the type after another \link colvarvalue \endlink
   inline void type(colvarvalue const &x)
   {
-    // reset the value held from based on the previous type
-    reset();
-    if (x.type() == type_vector) {
-      vector1d_value.resize(x.vector1d_value.size());
-    } else {
+    if (x.type() != value_type) {
+      // reset the value based on the previous type
+      reset();
       if (value_type == type_vector) {
         vector1d_value.resize(0);
       }
+      value_type = x.type();
+    }
+
+    if (x.type() == type_vector) {
+      vector1d_value.resize(x.vector1d_value.size());
     }
-    value_type = x.type();
   }
 
   /// Make the type a derivative of the original type