From 6e40300d26bf0456a3fd86d94e94166e11b58f68 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Wed, 22 Jul 2015 14:36:59 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@13670
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 lib/colvars/colvar.cpp               | 46 +++++++++++---
 lib/colvars/colvar.h                 |  2 +
 lib/colvars/colvarbias.cpp           |  7 +++
 lib/colvars/colvarbias.h             |  2 +-
 lib/colvars/colvarbias_histogram.cpp | 37 +++++++++--
 lib/colvars/colvarbias_histogram.h   |  8 +--
 lib/colvars/colvarcomp.cpp           |  8 ++-
 lib/colvars/colvarcomp.h             |  9 ++-
 lib/colvars/colvargrid.cpp           | 47 ++++++++++++++
 lib/colvars/colvargrid.h             | 91 +++++++++++++++++++---------
 lib/colvars/colvarmodule.cpp         |  1 +
 lib/colvars/colvarmodule.h           |  2 +-
 lib/colvars/colvarparse.cpp          | 14 ++++-
 lib/colvars/colvarparse.h            |  5 +-
 lib/colvars/colvarscript.cpp         | 23 +++++++
 15 files changed, 246 insertions(+), 56 deletions(-)

diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp
index b60c685395..5e90636e1d 100644
--- a/lib/colvars/colvar.cpp
+++ b/lib/colvars/colvar.cpp
@@ -845,6 +845,7 @@ void colvar::calc()
   // prepare atom groups for calculation
   if (cvm::debug())
     cvm::log("Collecting data from atom groups.\n");
+
   for (i = 0; i < cvcs.size(); i++) {
     for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
       cvm::atom_group &atoms = *(cvcs[i]->atom_groups[ig]);
@@ -856,6 +857,7 @@ void colvar::calc()
       // each atom group will take care of its own ref_pos_group, if defined
     }
   }
+
 ////  Don't try to get atom velocities, as no back-end currently implements it
 //   if (tasks[task_output_velocity] && !tasks[task_fdiff_velocity]) {
 //     for (i = 0; i < cvcs.size(); i++) {
@@ -864,6 +866,7 @@ void colvar::calc()
 //       }
 //     }
 //   }
+
   if (tasks[task_system_force]) {
     for (i = 0; i < cvcs.size(); i++) {
       for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
@@ -880,6 +883,7 @@ void colvar::calc()
 
   // First, update component values
   for (i = 0; i < cvcs.size(); i++) {
+    if (!cvcs[i]->b_enabled) continue;
     cvm::increase_depth();
     (cvcs[i])->calc_value();
     cvm::decrease_depth();
@@ -905,6 +909,7 @@ void colvar::calc()
   } else if (x.type() == colvarvalue::type_scalar) {
     // polynomial combination allowed
     for (i = 0; i < cvcs.size(); i++) {
+      if (!cvcs[i]->b_enabled) continue;
       x += (cvcs[i])->sup_coeff *
       ( ((cvcs[i])->sup_np != 1) ?
         std::pow((cvcs[i])->value().real_value, (cvcs[i])->sup_np) :
@@ -913,6 +918,7 @@ void colvar::calc()
   } else {
     // only linear combination allowed
     for (i = 0; i < cvcs.size(); i++) {
+      if (!cvcs[i]->b_enabled) continue;
       x += (cvcs[i])->sup_coeff * (cvcs[i])->value();
     }
   }
@@ -928,17 +934,15 @@ void colvar::calc()
 
     for (i = 0; i < cvcs.size(); i++) {
       // calculate the gradients of each component
+      if (!cvcs[i]->b_enabled) continue;
       cvm::increase_depth();
-
       (cvcs[i])->calc_gradients();
-
       // if requested, propagate (via chain rule) the gradients above
       // to the atoms used to define the roto-translation
       for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
         if (cvcs[i]->atom_groups[ig]->b_fit_gradients)
           cvcs[i]->atom_groups[ig]->calc_fit_gradients();
       }
-
       cvm::decrease_depth();
     }
 
@@ -958,6 +962,7 @@ void colvar::calc()
         atomic_gradients[a].reset();
       }
       for (i = 0; i < cvcs.size(); i++) {
+        if (!cvcs[i]->b_enabled) continue;
         // Coefficient: d(a * x^n) = a * n * x^(n-1) * dx
         cvm::real coeff = (cvcs[i])->sup_coeff * cvm::real((cvcs[i])->sup_np) *
           std::pow((cvcs[i])->value().real_value, (cvcs[i])->sup_np-1);
@@ -1123,19 +1128,23 @@ cvm::real colvar::update()
 
   if (tasks[task_Jacobian_force]) {
     size_t i;
-    cvm::increase_depth();
     for (i = 0; i < cvcs.size(); i++) {
+      if (!cvcs[i]->b_enabled) continue;
+      cvm::increase_depth();
       (cvcs[i])->calc_Jacobian_derivative();
+      cvm::decrease_depth();
     }
-    cvm::decrease_depth();
 
+    size_t ncvc = 0;
     fj.reset();
     for (i = 0; i < cvcs.size(); i++) {
+      if (!cvcs[i]->b_enabled) continue;
       // linear combination is assumed
-      fj += 1.0 / ( cvm::real(cvcs.size()) *  cvm::real((cvcs[i])->sup_coeff) ) *
+      fj += 1.0 / cvm::real((cvcs[i])->sup_coeff) *
         (cvcs[i])->Jacobian_derivative();
+      ncvc++;
     }
-    fj *= cvm::boltzmann() * cvm::temperature();
+    fj *= (1.0/cvm::real(ncvc)) * cvm::boltzmann() * cvm::temperature();
 
     // the instantaneous Jacobian force was not included in the reported system force;
     // instead, it is subtracted from the applied force (silent Jacobian correction)
@@ -1196,6 +1205,7 @@ void colvar::communicate_forces()
   if (tasks[task_scripted]) {
     std::vector<cvm::matrix2d<cvm::real> > func_grads;
     for (i = 0; i < cvcs.size(); i++) {
+      if (!cvcs[i]->b_enabled) continue;
       func_grads.push_back(cvm::matrix2d<cvm::real> (x.size(),
                                                      cvcs[i]->value().size()));
     }
@@ -1211,6 +1221,7 @@ void colvar::communicate_forces()
     }
 
     for (i = 0; i < cvcs.size(); i++) {
+      if (!cvcs[i]->b_enabled) continue;
       cvm::increase_depth();
       // cvc force is colvar force times colvar/cvc Jacobian
       // (vector-matrix product)
@@ -1221,6 +1232,7 @@ void colvar::communicate_forces()
   } else if (x.type() == colvarvalue::type_scalar) {
 
     for (i = 0; i < cvcs.size(); i++) {
+      if (!cvcs[i]->b_enabled) continue;
       cvm::increase_depth();
       (cvcs[i])->apply_force(f * (cvcs[i])->sup_coeff *
                               cvm::real((cvcs[i])->sup_np) *
@@ -1232,6 +1244,7 @@ void colvar::communicate_forces()
   } else {
 
     for (i = 0; i < cvcs.size(); i++) {
+      if (!cvcs[i]->b_enabled) continue;
       cvm::increase_depth();
       (cvcs[i])->apply_force(f * (cvcs[i])->sup_coeff);
       cvm::decrease_depth();
@@ -1243,6 +1256,25 @@ void colvar::communicate_forces()
 }
 
 
+int colvar::set_cvc_flags(std::vector<bool> const &flags) {
+
+  size_t i;
+  if (flags.size() != cvcs.size()) {
+    cvm::error("ERROR: Wrong number of CVC flags provided.");
+    return COLVARS_ERROR;
+  }
+  bool e = false;
+  for (i = 0; i < cvcs.size(); i++) {
+    cvcs[i]->b_enabled = flags[i];
+    e = e || flags[i];
+  }
+  if (!e) {
+    cvm::error("ERROR: All CVCs are disabled for this colvar.");
+    return COLVARS_ERROR;
+  }
+  return COLVARS_OK;
+}
+
 
 // ******************** METRIC FUNCTIONS ********************
 // Use the metrics defined by \link cvc \endlink objects
diff --git a/lib/colvars/colvar.h b/lib/colvars/colvar.h
index beadc52f51..79b00c026b 100644
--- a/lib/colvars/colvar.h
+++ b/lib/colvars/colvar.h
@@ -329,6 +329,8 @@ public:
   /// colvar::update()) to the external degrees of freedom
   void communicate_forces();
 
+  /// \brief Enables and disables individual CVCs based on flags
+  int set_cvc_flags(std::vector<bool> const &flags);
 
   /// \brief Use the internal metrics (as from \link cvc
   /// \endlink objects) to calculate square distances and gradients
diff --git a/lib/colvars/colvarbias.cpp b/lib/colvars/colvarbias.cpp
index 5853e22127..cb2037f0e5 100644
--- a/lib/colvars/colvarbias.cpp
+++ b/lib/colvars/colvarbias.cpp
@@ -100,6 +100,13 @@ void colvarbias::add_colvar(std::string const &cv_name)
 }
 
 
+cvm::real colvarbias::update()
+{
+  has_data = true;
+  return 0.0;
+}
+
+
 void colvarbias::communicate_forces()
 {
   for (size_t i = 0; i < colvars.size(); i++) {
diff --git a/lib/colvars/colvarbias.h b/lib/colvars/colvarbias.h
index e4cc89f463..3fe2b97833 100644
--- a/lib/colvars/colvarbias.h
+++ b/lib/colvars/colvarbias.h
@@ -19,7 +19,7 @@ public:
 
   /// Retrieve colvar values and calculate their biasing forces
   /// Return bias energy
-  virtual cvm::real update() = 0;
+  virtual cvm::real update();
 
   // TODO: move update_bias here (share with metadynamics)
 
diff --git a/lib/colvars/colvarbias_histogram.cpp b/lib/colvars/colvarbias_histogram.cpp
index 29f3c9bd9f..a131456663 100644
--- a/lib/colvars/colvarbias_histogram.cpp
+++ b/lib/colvars/colvarbias_histogram.cpp
@@ -11,7 +11,9 @@ colvarbias_histogram::colvarbias_histogram(std::string const &conf, char const *
     grid(NULL), out_name("")
 {
   get_keyval(conf, "outputFile", out_name, std::string(""));
+  get_keyval(conf, "outputFileDX", out_name_dx, std::string(""));
   get_keyval(conf, "outputFreq", output_freq, cvm::restart_out_freq);
+
   /// with VMD, this may not be an error
   // if ( output_freq == 0 ) {
   //   cvm::error("User required histogram with zero output frequency");
@@ -74,9 +76,6 @@ colvarbias_histogram::colvarbias_histogram(std::string const &conf, char const *
 /// Destructor
 colvarbias_histogram::~colvarbias_histogram()
 {
-  if (grid_os.is_open())
-    grid_os.close();
-
   if (grid) {
     delete grid;
     grid = NULL;
@@ -89,6 +88,9 @@ colvarbias_histogram::~colvarbias_histogram()
 /// Update the grid
 cvm::real colvarbias_histogram::update()
 {
+  // update base class
+  colvarbias::update();
+
   if (cvm::debug()) {
     cvm::log("Updating histogram bias " + this->name);
   }
@@ -105,6 +107,13 @@ cvm::real colvarbias_histogram::update()
     }
   }
 
+  if (out_name_dx.size() == 0) {
+    if (cvm::step_relative() == 0) {
+      out_name_dx = cvm::output_prefix + "." + this->name + ".dx";
+      cvm::log("Histogram " + this->name + " will be written to file \"" + out_name_dx + "\"");
+    }
+  }
+
   if (colvar_array_size == 0) {
     // update indices for scalar values
     size_t i;
@@ -132,23 +141,39 @@ cvm::real colvarbias_histogram::update()
   if (output_freq && (cvm::step_absolute() % output_freq) == 0) {
     write_output_files();
   }
+
   return 0.0; // no bias energy for histogram
 }
 
 
 int colvarbias_histogram::write_output_files()
 {
+  if (!has_data) {
+    // nothing to write
+    return COLVARS_OK;
+  }
+
   if (out_name.size()) {
     cvm::log("Writing the histogram file \""+out_name+"\".\n");
-
-    grid_os.open(out_name.c_str());
+    cvm::ofstream grid_os(out_name.c_str());
     if (!grid_os.is_open()) {
-      cvm::error("Error opening histogram file " + out_name + " for writing");
+      cvm::error("Error opening histogram file " + out_name + " for writing.\n", FILE_ERROR);
     }
     // TODO add return code here
     grid->write_multicol(grid_os);
     grid_os.close();
   }
+
+  if (out_name_dx.size()) {
+    cvm::log("Writing the histogram file \""+out_name_dx+"\".\n");
+    cvm::ofstream grid_os(out_name_dx.c_str());
+    if (!grid_os.is_open()) {
+      cvm::error("Error opening histogram file " + out_name_dx + " for writing.\n", FILE_ERROR);
+    }
+    // TODO add return code here
+    grid->write_opendx(grid_os);
+    grid_os.close();
+  }
   return COLVARS_OK;
 }
 
diff --git a/lib/colvars/colvarbias_histogram.h b/lib/colvars/colvarbias_histogram.h
index b9d07151d3..90fb5441c8 100644
--- a/lib/colvars/colvarbias_histogram.h
+++ b/lib/colvars/colvarbias_histogram.h
@@ -27,9 +27,8 @@ protected:
 
   /// n-dim histogram
   colvar_grid_scalar *grid;
-  std::vector<int>  bin;
-  std::string	  out_name;
-
+  std::vector<int> bin;
+  std::string out_name, out_name_dx;
   size_t output_freq;
 
   /// If one or more of the variables are \link type_vector \endlink, treat them as arrays of this length
@@ -37,9 +36,6 @@ protected:
   /// If colvar_array_size is larger than 1, weigh each one by this number before accumulating the histogram
   std::vector<cvm::real> weights;
 
-  void		  write_grid();
-  cvm::ofstream	  grid_os;  /// Stream for writing grid to disk
-
   std::istream& read_restart(std::istream&);
   std::ostream& write_restart(std::ostream&);
 };
diff --git a/lib/colvars/colvarcomp.cpp b/lib/colvars/colvarcomp.cpp
index fc31c44ab5..c05ab3f98b 100644
--- a/lib/colvars/colvarcomp.cpp
+++ b/lib/colvars/colvarcomp.cpp
@@ -8,7 +8,9 @@
 
 
 colvar::cvc::cvc()
-  : sup_coeff(1.0), sup_np(1),
+  : sup_coeff(1.0),
+    sup_np(1),
+    b_enabled(true),
     b_periodic(false),
     b_inverse_gradients(false),
     b_Jacobian_derivative(false),
@@ -17,7 +19,9 @@ colvar::cvc::cvc()
 
 
 colvar::cvc::cvc(std::string const &conf)
-  : sup_coeff(1.0), sup_np(1),
+  : sup_coeff(1.0),
+    sup_np(1),
+    b_enabled(true),
     b_periodic(false),
     b_inverse_gradients(false),
     b_Jacobian_derivative(false),
diff --git a/lib/colvars/colvarcomp.h b/lib/colvars/colvarcomp.h
index 62bd8a274e..0a66253b1d 100644
--- a/lib/colvars/colvarcomp.h
+++ b/lib/colvars/colvarcomp.h
@@ -84,14 +84,19 @@ public:
   /// \brief Exponent in the polynomial combination (default: 1)
   int       sup_np;
 
+  /// \brief This defaults to true; setting it to false disables
+  /// update of this cvc to save compute time (useful with scriptedFunction)
+  bool b_enabled;
+
+  /// \brief Is this a periodic component?
+  bool b_periodic;
+
   /// \brief Period of this cvc value, (default: 0.0, non periodic)
   cvm::real period;
 
   /// \brief If the component is periodic, wrap around this value (default: 0.0)
   cvm::real wrap_center;
 
-  bool b_periodic;
-
   /// \brief Constructor
   ///
   /// At least one constructor which reads a string should be defined
diff --git a/lib/colvars/colvargrid.cpp b/lib/colvars/colvargrid.cpp
index 76f7688b7c..6cc667b51c 100644
--- a/lib/colvars/colvargrid.cpp
+++ b/lib/colvars/colvargrid.cpp
@@ -103,6 +103,53 @@ std::ostream & colvar_grid_scalar::write_restart(std::ostream &os)
 }
 
 
+cvm::real colvar_grid_scalar::maximum_value() const
+{
+  cvm::real max = data[0];
+  for (size_t i = 0; i < nt; i++) {
+    if (data[i] > max) max = data[i];
+  }
+  return max;
+}
+
+
+cvm::real colvar_grid_scalar::minimum_value() const
+{
+  cvm::real min = data[0];
+  for (size_t i = 0; i < nt; i++) {
+    if (data[i] < min) min = data[i];
+  }
+  return min;
+}
+
+
+cvm::real colvar_grid_scalar::integral() const
+{
+  cvm::real sum = 0.0;
+  for (size_t i = 0; i < nt; i++) {
+    sum += data[i];
+  }
+  cvm::real bin_volume = 1.0;
+  for (size_t id = 0; id < widths.size(); id++) {
+    bin_volume *= widths[id];
+  }
+  return bin_volume * sum;
+}
+
+
+cvm::real colvar_grid_scalar::entropy() const
+{
+  cvm::real sum = 0.0;
+  for (size_t i = 0; i < nt; i++) {
+    sum += -1.0 * data[i] * std::log(data[i]);
+  }
+  cvm::real bin_volume = 1.0;
+  for (size_t id = 0; id < widths.size(); id++) {
+    bin_volume *= widths[id];
+  }
+  return bin_volume * sum;
+}
+
 
 colvar_grid_gradient::colvar_grid_gradient()
   : colvar_grid<cvm::real>(), samples(NULL)
diff --git a/lib/colvars/colvargrid.h b/lib/colvars/colvargrid.h
index be2d7cfc8f..29d21feb35 100644
--- a/lib/colvars/colvargrid.h
+++ b/lib/colvars/colvargrid.h
@@ -1025,6 +1025,49 @@ public:
     has_data = true;
     return is;
   }
+
+  /// \brief Write the grid data without labels, as they are
+  /// represented in memory
+  /// \param buf_size Number of values per line
+  std::ostream & write_opendx(std::ostream &os)
+  {
+    // write the header
+    os << "object 1 class gridpositions counts";
+    int icv;
+    for (icv = 0; icv < number_of_colvars(); icv++) {
+      os << " " << number_of_points(icv);
+    }
+    os << "\n";
+
+    os << "origin";
+    for (icv = 0; icv < number_of_colvars(); icv++) {
+      os << " " << (lower_boundaries[icv].real_value + 0.5 * widths[icv]);
+    }
+    os << "\n";
+
+    for (icv = 0; icv < number_of_colvars(); icv++) {
+      os << "delta";
+      for (size_t icv2 = 0; icv2 < number_of_colvars(); icv2++) {
+        if (icv == icv2) os << " " << widths[icv];
+        else os << " " << 0.0;
+      }
+      os << "\n";
+    }
+
+    os << "object 2 class gridconnections counts";
+    for (icv = 0; icv < number_of_colvars(); icv++) {
+      os << " " << number_of_points(icv);
+    }
+    os << "\n";
+
+    os << "object 3 class array type double rank 0 items "
+       << number_of_points() << " data follows\n";
+
+    write_raw(os);
+
+    os << "object \"collective variables scalar field\" class field\n";
+    return os;
+  }
 };
 
 
@@ -1114,12 +1157,12 @@ public:
 
   /// Constructor from a vector of colvars
   colvar_grid_scalar(std::vector<colvar *> &colvars,
-                      bool margin = 0);
+                     bool margin = 0);
 
   /// Accumulate the value
   inline void acc_value(std::vector<int> const &ix,
-                         cvm::real const &new_value,
-                         size_t const &imult = 0)
+                        cvm::real const &new_value,
+                        size_t const &imult = 0)
   {
     // only legal value of imult here is 0
     data[address(ix)] += new_value;
@@ -1154,32 +1197,33 @@ public:
   /// \brief Return the value of the function at ix divided by its
   /// number of samples (if the count grid is defined)
   virtual cvm::real value_output(std::vector<int> const &ix,
-                                  size_t const &imult = 0)
+                                 size_t const &imult = 0)
   {
     if (imult > 0) {
       cvm::error("Error: trying to access a component "
-                  "larger than 1 in a scalar data grid.\n");
+                 "larger than 1 in a scalar data grid.\n");
       return 0.;
     }
-    if (samples)
+    if (samples) {
       return (samples->value(ix) > 0) ?
         (data[address(ix)] / cvm::real(samples->value(ix))) :
         0.0;
-    else
+    } else {
       return data[address(ix)];
+    }
   }
 
   /// \brief Get the value from a formatted output and transform it
   /// into the internal representation (it may have been rescaled or
   /// manipulated)
   virtual void value_input(std::vector<int> const &ix,
-                            cvm::real const &new_value,
-                            size_t const &imult = 0,
-                            bool add = false)
+                           cvm::real const &new_value,
+                           size_t const &imult = 0,
+                           bool add = false)
   {
     if (imult > 0) {
       cvm::error("Error: trying to access a component "
-                  "larger than 1 in a scalar data grid.\n");
+                 "larger than 1 in a scalar data grid.\n");
       return;
     }
     if (add) {
@@ -1203,24 +1247,17 @@ public:
   std::ostream & write_restart(std::ostream &os);
 
   /// \brief Return the highest value
-  inline cvm::real maximum_value()
-  {
-    cvm::real max = data[0];
-    for (size_t i = 0; i < nt; i++) {
-      if (data[i] > max) max = data[i];
-    }
-    return max;
-  }
+  cvm::real maximum_value() const;
 
   /// \brief Return the lowest value
-  inline cvm::real minimum_value()
-  {
-    cvm::real min = data[0];
-    for (size_t i = 0; i < nt; i++) {
-      if (data[i] < min) min = data[i];
-    }
-    return min;
-  }
+  cvm::real minimum_value() const;
+
+  /// \brief Calculates the integral of the map (uses widths if they are defined)
+  cvm::real integral() const;
+
+  /// \brief Assuming that the map is a normalized probability density,
+  ///        calculates the entropy (uses widths if they are defined)
+  cvm::real entropy() const;
 
 private:
   // gradient
diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp
index 99a32ca8ff..2d62cace56 100644
--- a/lib/colvars/colvarmodule.cpp
+++ b/lib/colvars/colvarmodule.cpp
@@ -782,6 +782,7 @@ std::istream & colvarmodule::read_restart(std::istream &is)
       }
     }
     is.clear();
+    parse->clear_keyword_registry();
   }
 
   // colvars restart
diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h
index 25d2ce817b..e1d25b1d97 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-04-22"
+#define COLVARS_VERSION "2015-07-21"
 #endif
 
 #ifndef COLVARS_DEBUG
diff --git a/lib/colvars/colvarparse.cpp b/lib/colvars/colvarparse.cpp
index b564ec331b..a97d95e36e 100644
--- a/lib/colvars/colvarparse.cpp
+++ b/lib/colvars/colvarparse.cpp
@@ -355,6 +355,14 @@ void colvarparse::strip_values(std::string &conf)
 }
 
 
+void colvarparse::clear_keyword_registry()
+{
+  allowed_keywords.clear();
+  data_begin_pos.clear();
+  data_end_pos.clear();
+}
+
+
 int colvarparse::check_keywords(std::string &conf, char const *key)
 {
   if (cvm::debug())
@@ -396,9 +404,9 @@ int colvarparse::check_keywords(std::string &conf, char const *key)
       return COLVARS_ERROR;
     }
   }
-  allowed_keywords.clear();
-  data_begin_pos.clear();
-  data_end_pos.clear();
+
+  clear_keyword_registry();
+
   return COLVARS_OK;
 }
 
diff --git a/lib/colvars/colvarparse.h b/lib/colvars/colvarparse.h
index 5a753780a8..9e3811b660 100644
--- a/lib/colvars/colvarparse.h
+++ b/lib/colvars/colvarparse.h
@@ -37,7 +37,7 @@ protected:
   /// \brief Whether or not to accumulate data_begin_pos and
   /// data_end_pos in key_lookup(); it may be useful to disable
   /// this after the constructor is called, because other files may be
-  /// read (e.g. restart) that would mess up with the registry; in any
+  /// read (e.g. restart) that would mess up the registry; in any
   /// case, nothing serious happens until check_keywords() is invoked
   /// (which should happen only right after construction)
   bool save_delimiters;
@@ -144,6 +144,9 @@ public:
   /// then loop over all words
   int check_keywords(std::string &conf, char const *key);
 
+  /// \brief Use this after parsing a config string (note that check_keywords() calls it already)
+  void clear_keyword_registry();
+
 
   /// \brief Return a lowercased copy of the string
   static inline std::string to_lower_cppstr(std::string const &in)
diff --git a/lib/colvars/colvarscript.cpp b/lib/colvars/colvarscript.cpp
index 380f9d6fda..0c23faf519 100644
--- a/lib/colvars/colvarscript.cpp
+++ b/lib/colvars/colvarscript.cpp
@@ -256,6 +256,29 @@ int colvarscript::proc_colvar(int argc, char const *argv[]) {
     return COLVARSCRIPT_OK;
   }
 
+  if (subcmd == "cvcflags") {
+    if (argc < 4) {
+      result = "cvcflags: missing parameter: vector of flags";
+      return COLVARSCRIPT_ERROR;
+    }
+    std::string flags_str = argv[3];
+    std::istringstream is(flags_str);
+    std::vector<bool> flags;
+
+    int flag;
+    while (is >> flag) {
+      flags.push_back(flag != 0);
+    }
+
+    int res = cv->set_cvc_flags(flags);
+    if (res != COLVARS_OK) {
+      result = "Error setting CVC flags";
+      return COLVARSCRIPT_ERROR;
+    }
+    result = "0";
+    return COLVARSCRIPT_OK;
+  }
+
   result = "Syntax error\n" + help_string();
   return COLVARSCRIPT_ERROR;
 }