update colvars library to the state of 16 Oct 2014

lib/colvars/colvar.cpp

@@ -53,11 +53,17 @@ colvar::colvar (std::string const &conf)
       if (cvcp != NULL) {                                               \
         cvcs.push_back (cvcp);                                          \
         cvcp->check_keywords (def_conf, def_config_key);                \
+        if (cvm::get_error()) {                                         \
+          cvm::error("Error: in setting up component \""                \
+                      def_config_key"\".\n");                           \
+          return;                                                       \
+        }                                                               \
         cvm::decrease_depth();                                          \
       } else {                                                          \
-        cvm::error ("Error: in allocating component \""                 \
+        cvm::error("Error: in allocating component \""                  \
                           def_config_key"\".\n",                        \
                           MEMORY_ERROR);                                \
+        return;                                                         \
       }                                                                 \
       if ( (cvcp->period != 0.0) || (cvcp->wrap_center != 0.0) ) {      \
         if ( (cvcp->function_type != std::string ("distance_z")) &&     \
@@ -70,6 +76,7 @@ colvar::colvar (std::string const &conf)
                             "Period: "+cvm::to_str(cvcp->period) +      \
                         " wrapAround: "+cvm::to_str(cvcp->wrap_center), \
                         INPUT_ERROR);                                   \
+          return;                                                       \
         }                                                               \
       }                                                                 \
       if ( ! cvcs.back()->name.size())                                  \
@@ -150,7 +157,7 @@ colvar::colvar (std::string const &conf)
     "", colvarparse::parse_silent)) {
 
     enable(task_scripted);
-    cvm::log("This colvar is a scripted function.");
+    cvm::log("This colvar uses scripted function \"" + scripted_function + "\".");
 
     std::string type_str;
     get_keyval (conf, "scriptedFunctionType", type_str, "scalar");
@@ -187,18 +194,20 @@ colvar::colvar (std::string const &conf)
     }
   }
 
-  // this is set false if any of the components has an exponent
+  if (!tasks[task_scripted]) {
-  // different from 1 in the polynomial
+    // this is set false if any of the components has an exponent
-  b_linear = true;
+    // different from 1 in the polynomial
-
+    b_linear = true;
-  // these will be set to false if any of the cvcs has them false
+    // these will be set to false if any of the cvcs has them false
-  b_inverse_gradients = true;
+    b_inverse_gradients = true;
-  b_Jacobian_force    = true;
+    b_Jacobian_force    = true;
+  }
 
   // Decide whether the colvar is periodic
   // Used to wrap extended DOF if extendedLagrangian is on
   if (cvcs.size() == 1 && (cvcs[0])->b_periodic && (cvcs[0])->sup_np == 1
-                                                && (cvcs[0])->sup_coeff == 1.0 ) {
+                                                && (cvcs[0])->sup_coeff == 1.0
+                        && !tasks[task_scripted]) {
     this->b_periodic = true;
     this->period = (cvcs[0])->period;
     // TODO write explicit wrap() function for colvars to allow for

lib/colvars/colvaratoms.cpp

@@ -138,6 +138,7 @@ int cvm::atom_group::parse (std::string const &conf,
         for (size_t i = 0; i < atom_indexes.size(); i++) {
           this->push_back (cvm::atom (atom_indexes[i]));
         }
+        if (cvm::get_error()) return COLVARS_ERROR;
       } else {
         cvm::error ("Error: no numbers provided for \""
                     "atomNumbers\".\n", INPUT_ERROR);
@@ -166,6 +167,7 @@ int cvm::atom_group::parse (std::string const &conf,
       for (size_t i = 0; i < index_groups_i->size(); i++) {
         this->push_back (cvm::atom ((*index_groups_i)[i]));
       }
+      if (cvm::get_error()) return COLVARS_ERROR;
     }
   }
 
@@ -185,6 +187,7 @@ int cvm::atom_group::parse (std::string const &conf,
           for (int anum = initial; anum <= final; anum++) {
             this->push_back (cvm::atom (anum));
           }
+          if (cvm::get_error()) return COLVARS_ERROR;
           range_conf = "";
           continue;
         }
@@ -235,6 +238,7 @@ int cvm::atom_group::parse (std::string const &conf,
           for (int resid = initial; resid <= final; resid++) {
             this->push_back (cvm::atom (resid, atom_name, psf_segid));
           }
+          if (cvm::get_error()) return COLVARS_ERROR;
           range_conf = "";
         } else {
           cvm::error ("Error: cannot parse definition for \""
@@ -272,6 +276,9 @@ int cvm::atom_group::parse (std::string const &conf,
     }
   }
 
+  // Catch any errors from all the initialization steps above
+  if (cvm::get_error()) return COLVARS_ERROR;
+
   for (std::vector<cvm::atom>::iterator a1 = this->begin();
        a1 != this->end(); ++a1) {
     std::vector<cvm::atom>::iterator a2 = a1;
@@ -433,6 +440,10 @@ int cvm::atom_group::parse (std::string const &conf,
               std::string (key)+"\".\n");
 
   this->check_keywords (group_conf, key);
+  if (cvm::get_error()) {
+    cvm::error("Error setting up atom group \""+std::string (key)+"\".");
+    return COLVARS_ERROR;
+  }
 
   cvm::log ("Atom group \""+std::string (key)+"\" defined, "+
             cvm::to_str (this->size())+" atoms initialized: total mass = "+

lib/colvars/colvarcomp.cpp

@@ -47,11 +47,16 @@ void colvar::cvc::parse_group (std::string const &conf,
                                bool optional)
 {
   if (key_lookup (conf, group_key)) {
-    group.parse (conf, group_key);
+    if (group.parse (conf, group_key) != COLVARS_OK) {
+      cvm::error ("Error parsing definition for atom group \""+
+                         std::string (group_key)+"\".\n");
+      return;
+    }
   } else {
     if (! optional) {
-      cvm::fatal_error ("Error: definition for atom group \""+
+      cvm::error ("Error: definition for atom group \""+
                       std::string (group_key)+"\" not found.\n");
+      return;
     }
   }
 }
@@ -137,9 +142,7 @@ void colvar::cvc::debug_gradients (cvm::atom_group &group)
                              21, 14)+"\n");
       cvm::log ("|dx(actual) - dx(interp)|/|dx(actual)| = "+
                 cvm::to_str (std::fabs (x_1 - x_0 - dx_pred) /
-                             std::fabs (x_1 - x_0),
+                             std::fabs (x_1 - x_0), 12, 5)+"\n");
-                             12, 5)+
-                ".\n");
     }
   }
 

lib/colvars/colvarcomp_distances.cpp

@@ -725,7 +725,7 @@ colvar::rmsd::rmsd (std::string const &conf)
     return;
   }
 
-  if (atoms.ref_pos_group != NULL) {
+  if (atoms.ref_pos_group != NULL && b_Jacobian_derivative) {
     cvm::log ("The option \"refPositionsGroup\" (alternative group for fitting) was enabled: "
               "Jacobian derivatives of the RMSD will not be calculated.\n");
     b_Jacobian_derivative = false;
@@ -788,10 +788,21 @@ colvar::rmsd::rmsd (std::string const &conf)
     atoms.ref_pos = ref_pos;
     atoms.center_ref_pos();
 
+    cvm::log ("This is a standard minimum RMSD, derivatives of the optimal rotation "
+              "will not be computed as they cancel out in the gradients.");
+    atoms.b_fit_gradients = false;
+  }
+
+  if (atoms.b_rotate) {
+    // TODO: finer-grained control of this would require exposing a
+    // "request_Jacobian_derivative()" method to the colvar, and the same
+    // from the colvar to biases
+    // TODO: this should not be enabled here anyway, as it is not specific of the
+    // component - instead it should be decided in a generic way by the atom group
+
     // request the calculation of the derivatives of the rotation defined by the atom group
     atoms.rot.request_group1_gradients (atoms.size());
     // request derivatives of optimal rotation wrt reference coordinates for Jacobian:
-    // this is only required for ABF, but we do both groups here for better caching
     atoms.rot.request_group2_gradients (atoms.size());
   }
 }

lib/colvars/colvarmodule.cpp

@@ -99,7 +99,7 @@ int colvarmodule::config (std::string &conf)
   // parse global options
   error_code |= parse_global_params (conf);
 
-  if (error_code != COLVARS_OK) {
+  if (error_code != COLVARS_OK || cvm::get_error()) {
     set_error_bits(INPUT_ERROR);
     return COLVARS_ERROR;
   }
@@ -107,13 +107,23 @@ int colvarmodule::config (std::string &conf)
   // parse the options for collective variables
   error_code |= parse_colvars (conf);
 
+  if (error_code != COLVARS_OK || cvm::get_error()) {
+    set_error_bits(INPUT_ERROR);
+    return COLVARS_ERROR;
+  }
+
   // parse the options for biases
   error_code |= parse_biases (conf);
 
+  if (error_code != COLVARS_OK || cvm::get_error()) {
+    set_error_bits(INPUT_ERROR);
+    return COLVARS_ERROR;
+  }
+
   // done parsing known keywords, check that all keywords found were valid ones
   error_code |= parse->check_keywords (conf, "colvarmodule");
 
-  if (error_code != COLVARS_OK) {
+  if (error_code != COLVARS_OK || cvm::get_error()) {
     set_error_bits(INPUT_ERROR);
     return COLVARS_ERROR;
   }
@@ -179,12 +189,11 @@ int colvarmodule::parse_colvars (std::string const &conf)
       cvm::log (cvm::line_marker);
       cvm::increase_depth();
       colvars.push_back (new colvar (colvar_conf));
-      if (cvm::get_error()) {
+      if (cvm::get_error() ||
-        delete colvars.back();
+          ((colvars.back())->check_keywords (colvar_conf, "colvar") != COLVARS_OK)) {
-        colvars.pop_back();
+            cvm::log("Error while constructing colvar number " +
-        return COLVARS_ERROR;
+        cvm::to_str(colvars.size()) + " : deleting.");
-      }
+        delete colvars.back();  // the colvar destructor updates the colvars array
-      if ((colvars.back())->check_keywords (colvar_conf, "colvar") != COLVARS_OK) {
         return COLVARS_ERROR;
       }
       cvm::decrease_depth();
@@ -192,6 +201,7 @@ int colvarmodule::parse_colvars (std::string const &conf)
       cvm::error("Error: \"colvar\" keyword found without any configuration.\n", INPUT_ERROR);
       return COLVARS_ERROR;
     }
+    cvm::decrease_depth();
     colvar_conf = "";
   }
 
@@ -208,6 +218,17 @@ int colvarmodule::parse_colvars (std::string const &conf)
   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
 }
 
+bool colvarmodule::check_new_bias(std::string &conf, char const *key)
+{
+  if (cvm::get_error() ||
+     (biases.back()->check_keywords(conf, key) != COLVARS_OK)) {
+    cvm::log("Error while constructing bias number " +
+        cvm::to_str(biases.size()) + " : deleting.\n");
+    delete biases.back(); // the bias destructor updates the biases array
+    return true;
+  }
+  return false;
+}
 
 int colvarmodule::parse_biases (std::string const &conf)
 {
@@ -223,12 +244,7 @@ int colvarmodule::parse_biases (std::string const &conf)
         cvm::log (cvm::line_marker);
         cvm::increase_depth();
         biases.push_back (new colvarbias_abf (abf_conf, "abf"));
-        if (cvm::get_error()) {
+        if (cvm::check_new_bias(abf_conf, "abf")) {
-          delete biases.back();
-          biases.pop_back();
-          return COLVARS_ERROR;
-        }
-        if ((biases.back())->check_keywords (abf_conf, "abf") != COLVARS_OK) {
           return COLVARS_ERROR;
         }
         cvm::decrease_depth();
@@ -250,12 +266,7 @@ int colvarmodule::parse_biases (std::string const &conf)
         cvm::log (cvm::line_marker);
         cvm::increase_depth();
         biases.push_back (new colvarbias_restraint_harmonic (harm_conf, "harmonic"));
-        if (cvm::get_error()) {
+        if (cvm::check_new_bias(harm_conf, "harmonic")) {
-          delete biases.back();
-          biases.pop_back();
-          return COLVARS_ERROR;
-        }
-        if ((biases.back())->check_keywords (harm_conf, "harmonic") != COLVARS_OK) {
           return COLVARS_ERROR;
         }
         cvm::decrease_depth();
@@ -277,12 +288,7 @@ int colvarmodule::parse_biases (std::string const &conf)
         cvm::log (cvm::line_marker);
         cvm::increase_depth();
         biases.push_back (new colvarbias_restraint_linear (lin_conf, "linear"));
-        if (cvm::get_error()) {
+        if (cvm::check_new_bias(lin_conf, "linear")) {
-          delete biases.back();
-          biases.pop_back();
-          return COLVARS_ERROR;
-        }
-        if ((biases.back())->check_keywords (lin_conf, "linear") != COLVARS_OK) {
           return COLVARS_ERROR;
         }
         cvm::decrease_depth();
@@ -304,12 +310,7 @@ int colvarmodule::parse_biases (std::string const &conf)
         cvm::log (cvm::line_marker);
         cvm::increase_depth();
         biases.push_back (new colvarbias_alb (alb_conf, "ALB"));
-        if (cvm::get_error()) {
+        if (cvm::check_new_bias(alb_conf, "ALB")) {
-          delete biases.back();
-          biases.pop_back();
-          return COLVARS_ERROR;
-        }
-        if ((biases.back())->check_keywords (alb_conf, "ALB") != COLVARS_OK) {
           return COLVARS_ERROR;
         }
         cvm::decrease_depth();
@@ -332,12 +333,7 @@ int colvarmodule::parse_biases (std::string const &conf)
         cvm::log (cvm::line_marker);
         cvm::increase_depth();
         biases.push_back (new colvarbias_histogram (histo_conf, "histogram"));
-        if (cvm::get_error()) {
+        if (cvm::check_new_bias(histo_conf, "histogram")) {
-          delete biases.back();
-          biases.pop_back();
-          return COLVARS_ERROR;
-        }
-        if ((biases.back())->check_keywords (histo_conf, "histogram") != COLVARS_OK) {
           return COLVARS_ERROR;
         }
         cvm::decrease_depth();
@@ -359,12 +355,7 @@ int colvarmodule::parse_biases (std::string const &conf)
         cvm::log (cvm::line_marker);
         cvm::increase_depth();
         biases.push_back (new colvarbias_meta (meta_conf, "metadynamics"));
-        if (cvm::get_error()) {
+        if (cvm::check_new_bias(meta_conf, "metadynamics")) {
-          delete biases.back();
-          biases.pop_back();
-          return COLVARS_ERROR;
-        }
-        if ((biases.back())->check_keywords (meta_conf, "metadynamics") != COLVARS_OK) {
           return COLVARS_ERROR;
         }
         cvm::decrease_depth();
@@ -675,21 +666,20 @@ colvarmodule::~colvarmodule()
 
 int colvarmodule::reset()
 {
-  if (cvm::debug())
+  cvm::log("Resetting the Collective Variables Module.\n");
-    cvm::log ("colvars::reset() called.\n");
+  // Iterate backwards because we are deleting the elements as we go
-  for (std::vector<colvar *>::iterator cvi = colvars.begin();
+  for (std::vector<colvar *>::reverse_iterator cvi = colvars.rbegin();
-       cvi != colvars.end();
+       cvi != colvars.rend();
        cvi++) {
-    delete *cvi;
+    delete *cvi; // the colvar destructor updates the colvars array
-    cvi--;
   }
   colvars.clear();
 
-  for (std::vector<colvarbias *>::iterator bi = biases.begin();
+  // Iterate backwards because we are deleting the elements as we go
-       bi != biases.end();
+  for (std::vector<colvarbias *>::reverse_iterator bi = biases.rbegin();
+       bi != biases.rend();
        bi++) {
-    delete *bi;
+    delete *bi; // the bias destructor updates the biases array
-    bi--;
   }
   biases.clear();
 

lib/colvars/colvarmodule.h

@@ -4,7 +4,7 @@
 #define COLVARMODULE_H
 
 #ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2014-10-07"
+#define COLVARS_VERSION "2014-10-13"
 #endif
 
 #ifndef COLVARS_DEBUG
@@ -219,6 +219,9 @@ public:
   /// Parse and initialize collective variable biases
   int parse_biases (std::string const &conf);
 
+  /// Test error condition and keyword parsing
+  /// on error, delete new bias
+  bool check_new_bias(std::string &conf, char const *key);
 
   // "Setup" functions (change internal data based on related data
   // from the proxy that may change during program execution)

lib/colvars/colvarscript.cpp

@@ -214,10 +214,7 @@ int colvarscript::proc_colvar (int argc, char const *argv[]) {
   }
 
   if (subcmd == "update") {
-    // note: this is not sufficient for a newly created colvar
+    cv->calc();
-    // as atom coordinates may not be properly loaded
-    // a full CVM update is required
-    // or otherwise updating atom positions
     cv->update();
     result = cvm::to_str(cv->value(), 0, cvm::cv_prec);
     return COLVARSCRIPT_OK;

lib/colvars/colvartypes.cpp

@@ -350,22 +350,22 @@ void colvarmodule::rotation::calc_optimal_rotation
 
     // derivative of the S matrix
     ds_1.reset();
-    ds_1[0][0] = cvm::rvector ( a2x,  a2y,  a2z);
+    ds_1[0][0].set( a2x,  a2y,  a2z);
-    ds_1[1][0] = cvm::rvector ( 0.0,  a2z, -a2y);
+    ds_1[1][0].set( 0.0,  a2z, -a2y);
     ds_1[0][1] = ds_1[1][0];
-    ds_1[2][0] = cvm::rvector (-a2z,  0.0,  a2x);
+    ds_1[2][0].set(-a2z,  0.0,  a2x);
     ds_1[0][2] = ds_1[2][0];
-    ds_1[3][0] = cvm::rvector ( a2y, -a2x,  0.0);
+    ds_1[3][0].set( a2y, -a2x,  0.0);
     ds_1[0][3] = ds_1[3][0];
-    ds_1[1][1] = cvm::rvector ( a2x, -a2y, -a2z);
+    ds_1[1][1].set( a2x, -a2y, -a2z);
-    ds_1[2][1] = cvm::rvector ( a2y,  a2x,  0.0);
+    ds_1[2][1].set( a2y,  a2x,  0.0);
     ds_1[1][2] = ds_1[2][1];
-    ds_1[3][1] = cvm::rvector ( a2z,  0.0,  a2x);
+    ds_1[3][1].set( a2z,  0.0,  a2x);
     ds_1[1][3] = ds_1[3][1];
-    ds_1[2][2] = cvm::rvector (-a2x,  a2y, -a2z);
+    ds_1[2][2].set(-a2x,  a2y, -a2z);
-    ds_1[3][2] = cvm::rvector ( 0.0,  a2z,  a2y);
+    ds_1[3][2].set( 0.0,  a2z,  a2y);
     ds_1[2][3] = ds_1[3][2];
-    ds_1[3][3] = cvm::rvector (-a2x, -a2y,  a2z);
+    ds_1[3][3].set(-a2x, -a2y,  a2z);
 
     cvm::rvector              &dl0_1 = dL0_1[ia];
     vector1d<cvm::rvector, 4> &dq0_1 = dQ0_1[ia];
@@ -404,22 +404,22 @@ void colvarmodule::rotation::calc_optimal_rotation
     matrix2d<cvm::rvector, 4, 4> &ds_2 =  dS_2[ia];
 
     ds_2.reset();
-    ds_2[0][0] = cvm::rvector ( a1x,  a1y,  a1z);
+    ds_2[0][0].set( a1x,  a1y,  a1z);
-    ds_2[1][0] = cvm::rvector ( 0.0, -a1z,  a1y);
+    ds_2[1][0].set( 0.0, -a1z,  a1y);
     ds_2[0][1] = ds_2[1][0];
-    ds_2[2][0] = cvm::rvector ( a1z,  0.0, -a1x);
+    ds_2[2][0].set( a1z,  0.0, -a1x);
     ds_2[0][2] = ds_2[2][0];
-    ds_2[3][0] = cvm::rvector (-a1y,  a1x,  0.0);
+    ds_2[3][0].set(-a1y,  a1x,  0.0);
     ds_2[0][3] = ds_2[3][0];
-    ds_2[1][1] = cvm::rvector ( a1x, -a1y, -a1z);
+    ds_2[1][1].set( a1x, -a1y, -a1z);
-    ds_2[2][1] = cvm::rvector ( a1y,  a1x,  0.0);
+    ds_2[2][1].set( a1y,  a1x,  0.0);
     ds_2[1][2] = ds_2[2][1];
-    ds_2[3][1] = cvm::rvector ( a1z,  0.0,  a1x);
+    ds_2[3][1].set( a1z,  0.0,  a1x);
     ds_2[1][3] = ds_2[3][1];
-    ds_2[2][2] = cvm::rvector (-a1x,  a1y, -a1z);
+    ds_2[2][2].set(-a1x,  a1y, -a1z);
-    ds_2[3][2] = cvm::rvector ( 0.0,  a1z,  a1y);
+    ds_2[3][2].set( 0.0,  a1z,  a1y);
     ds_2[2][3] = ds_2[3][2];
-    ds_2[3][3] = cvm::rvector (-a1x, -a1y,  a1z);
+    ds_2[3][3].set(-a1x, -a1y,  a1z);
 
     cvm::rvector              &dl0_2 = dL0_2[ia];
     vector1d<cvm::rvector, 4> &dq0_2 = dQ0_2[ia];

lib/colvars/colvartypes.h

@@ -37,10 +37,19 @@ public:
   {}
 
   /// \brief Set all components to a scalar value
-  inline void set (cvm::real const &value = 0.0) {
+  inline void set (cvm::real const &value) {
     x = y = z = value;
   }
 
+  /// \brief Assign all components
+  inline void set (cvm::real const &x_i,
+                   cvm::real const &y_i,
+                   cvm::real const &z_i) {
+    x = x_i;
+    y = y_i;
+    z = z_i;
+  }
+
   /// \brief Set all components to zero
   inline void reset() {
     x = y = z = 0.0;

src/USER-COLVARS/colvarproxy_lammps.cpp

@@ -90,7 +90,7 @@ colvarproxy_lammps::colvarproxy_lammps(LAMMPS_NS::LAMMPS *lmp,
   // output prefix is always given
   output_prefix_str = std::string(out_name);
   // not so for restarts
-  restart_prefix_str = std::string("rest");
+  restart_output_prefix_str = std::string("rest");
 
   // check if it is possible to save output configuration
   if ((!output_prefix_str.size()) && (!restart_output_prefix_str.size())) {
@@ -101,13 +101,13 @@ colvarproxy_lammps::colvarproxy_lammps(LAMMPS_NS::LAMMPS *lmp,
   // try to extract a restart prefix from a potential restart command.
   LAMMPS_NS::Output *outp = _lmp->output;
   if ((outp->restart_every_single > 0) && (outp->restart1 != 0)) {
-      restart_prefix_str = std::string(outp->restart1);
+      restart_output_prefix_str = std::string(outp->restart1);
   } else if  ((outp->restart_every_double > 0) && (outp->restart2a != 0)) {
-    restart_prefix_str = std::string(outp->restart2a);
+    restart_output_prefix_str = std::string(outp->restart2a);
   }
   // trim off unwanted stuff from the restart prefix
-  if (restart_prefix_str.rfind(".*") != std::string::npos)
+  if (restart_output_prefix_str.rfind(".*") != std::string::npos)
-    restart_prefix_str.erase(restart_prefix_str.rfind(".*"),2);
+    restart_output_prefix_str.erase(restart_output_prefix_str.rfind(".*"),2);
 }
 
 

src/USER-COLVARS/colvarproxy_lammps.h

@@ -44,10 +44,6 @@ class colvarproxy_lammps : public colvarproxy {
   int  restart_every;
   int  previous_step;
 
-  std::string input_prefix_str;
-  std::string output_prefix_str;
-  std::string restart_prefix_str;
-
   bool first_timestep;
   bool system_force_requested;
   bool do_exit;
@@ -102,9 +98,6 @@ class colvarproxy_lammps : public colvarproxy {
   cvm::real temperature() { return t_target; };
   cvm::real dt() { return _lmp->update->dt * _lmp->force->femtosecond; };
 
-  inline std::string input_prefix() { return input_prefix_str; };
-  inline std::string output_prefix() { return output_prefix_str; };
-  inline std::string restart_output_prefix() { return restart_prefix_str; };
   inline size_t restart_frequency() { return restart_every; };
 
   void add_energy (cvm::real energy) { bias_energy = energy; };