Update Colvars to version 2021-08-03

lib/colvars/Makefile.common

@@ -30,6 +30,7 @@ COLVARS_SRCS = \
         colvarbias_histogram.cpp \
         colvarbias_meta.cpp \
         colvarbias_restraint.cpp \
+        colvarcomp_alchlambda.cpp \
         colvarcomp_angles.cpp \
         colvarcomp_apath.cpp \
         colvarcomp_coordnums.cpp \

lib/colvars/Makefile.deps

@@ -86,6 +86,20 @@ $(COLVARS_OBJ_DIR)colvarbias_restraint.o: colvarbias_restraint.cpp \
  lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \
  lepton/include/lepton/Exception.h \
  lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h
+$(COLVARS_OBJ_DIR)colvarcomp_alchlambda.o: colvarcomp_alchlambda.cpp \
+ colvarmodule.h colvars_version.h colvar.h colvarvalue.h colvartypes.h \
+ colvarparse.h colvarparams.h colvardeps.h lepton/include/Lepton.h \
+ lepton/include/lepton/CompiledExpression.h \
+ lepton/include/lepton/ExpressionTreeNode.h \
+ lepton/include/lepton/windowsIncludes.h \
+ lepton/include/lepton/CustomFunction.h \
+ lepton/include/lepton/ExpressionProgram.h \
+ lepton/include/lepton/ExpressionTreeNode.h \
+ lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \
+ lepton/include/lepton/Exception.h \
+ lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \
+ colvarcomp.h colvaratoms.h colvarproxy.h colvarproxy_tcl.h \
+ colvarproxy_volmaps.h colvar_arithmeticpath.h colvar_geometricpath.h
 $(COLVARS_OBJ_DIR)colvarcomp_angles.o: colvarcomp_angles.cpp \
  colvarmodule.h colvars_version.h colvar.h colvarvalue.h colvartypes.h \
  colvarparse.h colvarparams.h colvardeps.h lepton/include/Lepton.h \
@@ -285,7 +299,7 @@ $(COLVARS_OBJ_DIR)colvarproxy.o: colvarproxy.cpp colvarmodule.h \
  lepton/include/lepton/Exception.h \
  lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \
  colvarscript_commands.h colvarscript_commands_colvar.h \
- colvarscript_commands_bias.h colvaratoms.h
+ colvarscript_commands_bias.h colvaratoms.h colvarmodule_utils.h
 $(COLVARS_OBJ_DIR)colvarproxy_replicas.o: colvarproxy_replicas.cpp \
  colvarmodule.h colvars_version.h colvarproxy.h colvartypes.h \
  colvarvalue.h colvarproxy_tcl.h colvarproxy_volmaps.h
@@ -294,7 +308,8 @@ $(COLVARS_OBJ_DIR)colvarproxy_tcl.o: colvarproxy_tcl.cpp colvarmodule.h \
  colvarproxy_tcl.h colvarproxy_volmaps.h colvaratoms.h colvarparse.h \
  colvarparams.h colvardeps.h
 $(COLVARS_OBJ_DIR)colvarproxy_volmaps.o: colvarproxy_volmaps.cpp \
- colvarmodule.h colvars_version.h colvarproxy_volmaps.h
+ colvarmodule.h colvars_version.h colvarproxy_volmaps.h \
+ colvarmodule_utils.h
 $(COLVARS_OBJ_DIR)colvarscript.o: colvarscript.cpp colvarproxy.h \
  colvarmodule.h colvars_version.h colvartypes.h colvarvalue.h \
  colvarproxy_tcl.h colvarproxy_volmaps.h colvardeps.h colvarparse.h \

lib/colvars/colvar.cpp

@@ -96,7 +96,7 @@ int colvar::init(std::string const &conf)
     "", colvarparse::parse_silent)) {
 
     enable(f_cv_scripted);
-    cvm::log("This colvar uses scripted function \"" + scripted_function + "\".");
+    cvm::log("This colvar uses scripted function \"" + scripted_function + "\".\n");
 
     std::string type_str;
     get_keyval(conf, "scriptedFunctionType", type_str, "scalar");
@@ -134,7 +134,7 @@ int colvar::init(std::string const &conf)
     std::sort(cvcs.begin(), cvcs.end(), compare);
 
     if(cvcs.size() > 1) {
-      cvm::log("Sorted list of components for this scripted colvar:");
+      cvm::log("Sorted list of components for this scripted colvar:\n");
       for (i = 0; i < cvcs.size(); i++) {
         cvm::log(cvm::to_str(i+1) + " " + cvcs[i]->name);
       }
@@ -298,6 +298,11 @@ int colvar::init(std::string const &conf)
   error_code |= init_extended_Lagrangian(conf);
   error_code |= init_output_flags(conf);
 
+  // Detect if we have one component that is an alchemical lambda
+  if (is_enabled(f_cv_single_cvc) && cvcs[0]->function_type == "alch_lambda") {
+    enable(f_cv_external);
+  }
+
   // Now that the children are defined we can solve dependencies
   enable(f_cv_active);
 
@@ -636,7 +641,6 @@ int colvar::init_extended_Lagrangian(std::string const &conf)
     x_ext.type(colvarvalue::type_notset);
     v_ext.type(value());
     fr.type(value());
-
     const bool found = get_keyval(conf, "extendedTemp", temp, cvm::temperature());
     if (temp <= 0.0) {
       if (found)
@@ -654,7 +658,7 @@ int colvar::init_extended_Lagrangian(std::string const &conf)
       return INPUT_ERROR;
     }
     ext_force_k = cvm::boltzmann() * temp / (tolerance * tolerance);
-    cvm::log("Computed extended system force constant: " + cvm::to_str(ext_force_k) + " [E]/U^2");
+    cvm::log("Computed extended system force constant: " + cvm::to_str(ext_force_k) + " [E]/U^2\n");
 
     get_keyval(conf, "extendedTimeConstant", extended_period, 200.0);
     if (extended_period <= 0.0) {
@@ -662,7 +666,7 @@ int colvar::init_extended_Lagrangian(std::string const &conf)
     }
     ext_mass = (cvm::boltzmann() * temp * extended_period * extended_period)
       / (4.0 * PI * PI * tolerance * tolerance);
-    cvm::log("Computed fictitious mass: " + cvm::to_str(ext_mass) + " [E]/(U/fs)^2   (U: colvar unit)");
+    cvm::log("Computed fictitious mass: " + cvm::to_str(ext_mass) + " [E]/(U/fs)^2   (U: colvar unit)\n");
 
     {
       bool b_output_energy;
@@ -770,7 +774,9 @@ template<typename def_class_name> int colvar::init_components_type(std::string c
       if ( (cvcp->function_type != std::string("distance_z")) &&
            (cvcp->function_type != std::string("dihedral")) &&
            (cvcp->function_type != std::string("polar_phi")) &&
-           (cvcp->function_type != std::string("spin_angle")) ) {
+           (cvcp->function_type != std::string("spin_angle")) &&
+           (cvcp->function_type != std::string("euler_phi")) &&
+           (cvcp->function_type != std::string("euler_psi"))) {
         cvm::error("Error: invalid use of period and/or "
                    "wrapAround in a \""+
                    std::string(def_config_key)+
@@ -860,6 +866,7 @@ int colvar::init_components(std::string const &conf)
     "inertia", "inertia");
   error_code |= init_components_type<inertia_z>(conf, "moment of inertia around an axis", "inertiaZ");
   error_code |= init_components_type<eigenvector>(conf, "eigenvector", "eigenvector");
+  error_code |= init_components_type<alch_lambda>(conf, "alchemical coupling parameter", "alchLambda");
   error_code |= init_components_type<gspath>(conf, "geometrical path collective variables (s)", "gspath");
   error_code |= init_components_type<gzpath>(conf, "geometrical path collective variables (z)", "gzpath");
   error_code |= init_components_type<linearCombination>(conf, "linear combination of other collective variables", "linearCombination");
@@ -867,6 +874,9 @@ int colvar::init_components(std::string const &conf)
   error_code |= init_components_type<gzpathCV>(conf, "geometrical path collective variables (z) for other CVs", "gzpathCV");
   error_code |= init_components_type<aspathCV>(conf, "arithmetic path collective variables (s) for other CVs", "aspathCV");
   error_code |= init_components_type<azpathCV>(conf, "arithmetic path collective variables (s) for other CVs", "azpathCV");
+  error_code |= init_components_type<euler_phi>(conf, "euler phi angle of the optimal orientation", "eulerPhi");
+  error_code |= init_components_type<euler_psi>(conf, "euler psi angle of the optimal orientation", "eulerPsi");
+  error_code |= init_components_type<euler_theta>(conf, "euler theta angle of the optimal orientation", "eulerTheta");
 
   error_code |= init_components_type<map_total>(conf, "total value of atomic map", "mapTotal");
 
@@ -1074,6 +1084,7 @@ int colvar::init_dependencies() {
 
     init_feature(f_cv_hide_Jacobian, "hide_Jacobian_force", f_type_user);
     require_feature_self(f_cv_hide_Jacobian, f_cv_Jacobian); // can only hide if calculated
+    exclude_feature_self(f_cv_hide_Jacobian, f_cv_extended_Lagrangian);
 
     init_feature(f_cv_extended_Lagrangian, "extended_Lagrangian", f_type_user);
     require_feature_self(f_cv_extended_Lagrangian, f_cv_scalar);
@@ -1082,6 +1093,9 @@ int colvar::init_dependencies() {
     init_feature(f_cv_Langevin, "Langevin_dynamics", f_type_user);
     require_feature_self(f_cv_Langevin, f_cv_extended_Lagrangian);
 
+    init_feature(f_cv_external, "external", f_type_user);
+    require_feature_self(f_cv_external, f_cv_single_cvc);
+
     init_feature(f_cv_single_cvc, "single_component", f_type_static);
 
     init_feature(f_cv_linear, "linear", f_type_static);
@@ -1195,6 +1209,21 @@ std::vector<std::vector<int> > colvar::get_atom_lists()
 }
 
 
+std::vector<int> const &colvar::get_volmap_ids()
+{
+  volmap_ids_.resize(cvcs.size());
+  for (size_t i = 0; i < cvcs.size(); i++) {
+    if (cvcs[i]->param_exists("mapID") == COLVARS_OK) {
+      volmap_ids_[i] =
+        *(reinterpret_cast<int const *>(cvcs[i]->get_param_ptr("mapID")));
+    } else {
+      volmap_ids_[i] = -1;
+    }
+  }
+  return volmap_ids_;
+}
+
+
 colvar::~colvar()
 {
   // There is no need to call free_children_deps() here
@@ -1552,9 +1581,11 @@ int colvar::collect_cvc_total_forces()
       }
     }
 
-    if (!is_enabled(f_cv_hide_Jacobian)) {
+    if (!(is_enabled(f_cv_hide_Jacobian) && is_enabled(f_cv_subtract_applied_force))) {
       // add the Jacobian force to the total force, and don't apply any silent
       // correction internally: biases such as colvarbias_abf will handle it
+      // If f_cv_hide_Jacobian is enabled, a force of -fj is present in ft due to the
+      // Jacobian-compensating force
       ft += fj;
     }
   }
@@ -1687,14 +1718,16 @@ cvm::real colvar::update_forces_energy()
 
   // add the biases' force, which at this point should already have
   // been summed over each bias using this colvar
+  // fb is already multiplied by the relevant time step factor for each bias
   f += fb;
 
   if (is_enabled(f_cv_Jacobian)) {
     // the instantaneous Jacobian force was not included in the reported total force;
     // instead, it is subtracted from the applied force (silent Jacobian correction)
     // This requires the Jacobian term for the *current* timestep
+    // Need to scale it for impulse MTS
     if (is_enabled(f_cv_hide_Jacobian))
-      f -= fj;
+      f -= fj * cvm::real(time_step_factor);
   }
 
   // At this point f is the force f from external biases that will be applied to the
@@ -1727,6 +1760,25 @@ cvm::real colvar::update_forces_energy()
       colvarvalue f_ext(fr.type()); // force acting on the extended variable
       f_ext.reset();
 
+      if (is_enabled(f_cv_external)) {
+        // There are no forces on the "actual colvar" bc there is no gradient wrt atomic coordinates
+        // So we apply this to the extended DOF
+        f += fb_actual;
+      }
+
+      fr    = f;
+      // External force has been scaled for a 1-timestep impulse, scale it back because we will
+      // integrate it with the colvar's own timestep factor
+      f_ext = f / cvm::real(time_step_factor);
+
+      colvarvalue f_system(fr.type()); // force exterted by the system on the extended DOF
+
+      if (is_enabled(f_cv_external)) {
+        // Add "alchemical" force from external variable
+        f_system = cvcs[0]->total_force();
+        // f is now irrelevant because we are not applying atomic forces in the simulation
+        // just driving the external variable lambda
+      } else {
         // the total force is applied to the fictitious mass, while the
         // atoms only feel the harmonic force + wall force
         // fr: bias force on extended variable (without harmonic spring), for output in trajectory
@@ -1734,19 +1786,18 @@ cvm::real colvar::update_forces_energy()
         // f: - initially, external biasing force
         //    - after this code block, colvar force to be applied to atomic coordinates
         //      ie. spring force (fb_actual will be added just below)
-      fr    = f;
+        f_system = (-0.5 * ext_force_k) * this->dist2_lgrad(x_ext, x);
-      // External force has been scaled for a 1-timestep impulse, scale it back because we will
+        f        = -1.0 * f_system;
-      // integrate it with the colvar's own timestep factor
-      f_ext = f / cvm::real(time_step_factor);
-      f_ext += (-0.5 * ext_force_k) * this->dist2_lgrad(x_ext, x);
-      f      = (-0.5 * ext_force_k) * this->dist2_rgrad(x_ext, x);
         // Coupling force is a slow force, to be applied to atomic coords impulse-style
+        // over a single MD timestep
         f *= cvm::real(time_step_factor);
+      }
+      f_ext += f_system;
 
       if (is_enabled(f_cv_subtract_applied_force)) {
         // Report a "system" force without the biases on this colvar
-        // that is, just the spring force
+        // that is, just the spring force (or alchemical force)
-        ft_reported = (-0.5 * ext_force_k) * this->dist2_lgrad(x_ext, x);
+        ft_reported = f_system;
       } else {
         // The total force acting on the extended variable is f_ext
         // This will be used in the next timestep
@@ -1778,10 +1829,19 @@ cvm::real colvar::update_forces_energy()
           (is_enabled(f_cv_reflecting_upper_boundary) && (delta = x_ext - upper_boundary) > 0)) {
         x_ext -= 2.0 * delta;
         v_ext *= -1.0;
+        if ((is_enabled(f_cv_reflecting_lower_boundary) && (delta = x_ext - lower_boundary) < 0) ||
+            (is_enabled(f_cv_reflecting_upper_boundary) && (delta = x_ext - upper_boundary) > 0)) {
+          cvm::error("Error: extended coordinate value " + cvm::to_str(x_ext) + " is still outside boundaries after reflection.\n");
+        }
       }
 
       x_ext.apply_constraints();
       this->wrap(x_ext);
+      if (is_enabled(f_cv_external)) {
+        // Colvar value is constrained to the extended value
+        x = x_ext;
+        cvcs[0]->set_value(x_ext);
+      }
     } else {
       // If this is a postprocessing run (eg. in VMD), the extended DOF
       // is equal to the actual coordinate
@@ -1789,9 +1849,11 @@ cvm::real colvar::update_forces_energy()
     }
   }
 
+  if (!is_enabled(f_cv_external)) {
     // Now adding the force on the actual colvar (for those biases that
     // bypass the extended Lagrangian mass)
     f += fb_actual;
+  }
 
   if (cvm::debug())
     cvm::log("Done updating colvar \""+this->name+"\".\n");
@@ -1954,7 +2016,7 @@ int colvar::update_cvc_flags()
 
 int colvar::update_cvc_config(std::vector<std::string> const &confs)
 {
-  cvm::log("Updating configuration for colvar \""+name+"\"");
+  cvm::log("Updating configuration for colvar \""+name+"\"\n");
 
   if (confs.size() != cvcs.size()) {
     return cvm::error("Error: Wrong number of CVC config strings.  "
@@ -2076,6 +2138,15 @@ bool colvar::periodic_boundaries() const
 cvm::real colvar::dist2(colvarvalue const &x1,
                          colvarvalue const &x2) const
 {
+  if ( is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function) ) {
+    if (is_enabled(f_cv_periodic) && is_enabled(f_cv_scalar)) {
+      cvm::real diff = x1.real_value - x2.real_value;
+      const cvm::real period_lower_boundary = wrap_center - period / 2.0;
+      const cvm::real period_upper_boundary = wrap_center + period / 2.0;
+      diff = (diff < period_lower_boundary ? diff + period : ( diff > period_upper_boundary ? diff - period : diff));
+      return diff * diff;
+    }
+  }
   if (is_enabled(f_cv_homogeneous)) {
     return (cvcs[0])->dist2(x1, x2);
   } else {
@@ -2086,6 +2157,15 @@ cvm::real colvar::dist2(colvarvalue const &x1,
 colvarvalue colvar::dist2_lgrad(colvarvalue const &x1,
                                  colvarvalue const &x2) const
 {
+  if ( is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function) ) {
+    if (is_enabled(f_cv_periodic) && is_enabled(f_cv_scalar)) {
+      cvm::real diff = x1.real_value - x2.real_value;
+      const cvm::real period_lower_boundary = wrap_center - period / 2.0;
+      const cvm::real period_upper_boundary = wrap_center + period / 2.0;
+      diff = (diff < period_lower_boundary ? diff + period : ( diff > period_upper_boundary ? diff - period : diff));
+      return 2.0 * diff;
+    }
+  }
   if (is_enabled(f_cv_homogeneous)) {
     return (cvcs[0])->dist2_lgrad(x1, x2);
   } else {
@@ -2096,6 +2176,15 @@ colvarvalue colvar::dist2_lgrad(colvarvalue const &x1,
 colvarvalue colvar::dist2_rgrad(colvarvalue const &x1,
                                  colvarvalue const &x2) const
 {
+  if ( is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function) ) {
+    if (is_enabled(f_cv_periodic) && is_enabled(f_cv_scalar)) {
+      cvm::real diff = x1.real_value - x2.real_value;
+      const cvm::real period_lower_boundary = wrap_center - period / 2.0;
+      const cvm::real period_upper_boundary = wrap_center + period / 2.0;
+      diff = (diff < period_lower_boundary ? diff + period : ( diff > period_upper_boundary ? diff - period : diff));
+      return (-2.0) * diff;
+    }
+  }
   if (is_enabled(f_cv_homogeneous)) {
     return (cvcs[0])->dist2_rgrad(x1, x2);
   } else {
@@ -2302,7 +2391,7 @@ std::ostream & colvar::write_traj_label(std::ostream & os)
     os << " "
        << cvm::wrap_string(this->name, this_cv_width);
 
-    if (is_enabled(f_cv_extended_Lagrangian)) {
+    if (is_enabled(f_cv_extended_Lagrangian) && !is_enabled(f_cv_external)) {
       // extended DOF
       os << " r_"
          << cvm::wrap_string(this->name, this_cv_width-2);
@@ -2314,7 +2403,7 @@ std::ostream & colvar::write_traj_label(std::ostream & os)
     os << " v_"
        << cvm::wrap_string(this->name, this_cv_width-2);
 
-    if (is_enabled(f_cv_extended_Lagrangian)) {
+    if (is_enabled(f_cv_extended_Lagrangian) && !is_enabled(f_cv_external)) {
       // extended DOF
       os << " vr_"
          << cvm::wrap_string(this->name, this_cv_width-3);
@@ -2347,7 +2436,7 @@ std::ostream & colvar::write_traj(std::ostream &os)
   os << " ";
   if (is_enabled(f_cv_output_value)) {
 
-    if (is_enabled(f_cv_extended_Lagrangian)) {
+    if (is_enabled(f_cv_extended_Lagrangian) && !is_enabled(f_cv_external)) {
       os << " "
          << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
          << x;
@@ -2360,7 +2449,7 @@ std::ostream & colvar::write_traj(std::ostream &os)
 
   if (is_enabled(f_cv_output_velocity)) {
 
-    if (is_enabled(f_cv_extended_Lagrangian)) {
+    if (is_enabled(f_cv_extended_Lagrangian) && !is_enabled(f_cv_external)) {
       os << " "
          << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
          << v_fdiff;

lib/colvars/colvar.h

@@ -192,7 +192,7 @@ protected:
   /// Amplitude of Gaussian white noise for Langevin extended dynamics
   cvm::real ext_sigma;
 
-  /// \brief Harmonic restraint force
+  /// \brief Applied force on extended DOF, for output (unscaled if using MTS)
   colvarvalue fr;
 
   /// \brief Jacobian force, when Jacobian_force is enabled
@@ -591,6 +591,7 @@ public:
   class alpha_dihedrals;
   class alpha_angles;
   class dihedPC;
+  class alch_lambda;
   class componentDisabled;
   class CartesianBasedPath;
   class gspath;
@@ -601,6 +602,9 @@ public:
   class gzpathCV;
   class aspathCV;
   class azpathCV;
+  class euler_phi;
+  class euler_psi;
+  class euler_theta;
 
   // non-scalar components
   class distance_vec;
@@ -658,6 +662,9 @@ protected:
   static std::map<std::string, std::function<colvar::cvc* (const std::string& subcv_conf)>> global_cvc_map;
 #endif
 
+  /// Volmap numeric IDs, one for each CVC (-1 if not available)
+  std::vector<int> volmap_ids_;
+
 public:
 
   /// \brief Sorted array of (zero-based) IDs for all atoms involved
@@ -670,6 +677,10 @@ public:
 
   /// \brief Get vector of vectors of atom IDs for all atom groups
   virtual std::vector<std::vector<int> > get_atom_lists();
+
+  /// Volmap numeric IDs, one for each CVC (-1 if not available)
+  std::vector<int> const &get_volmap_ids();
+
 };
 
 

lib/colvars/colvar_arithmeticpath.h

@@ -7,6 +7,7 @@
 #include <cmath>
 #include <iostream>
 #include <limits>
+#include <string>
 
 namespace ArithmeticPathCV {
 
@@ -24,6 +25,7 @@ public:
     virtual void computeValue();
     virtual void computeDerivatives();
     virtual void compute();
+    virtual void reComputeLambda(const vector<scalar_type>& rmsd_between_refs);
 protected:
     scalar_type lambda;
     vector<scalar_type> weights;
@@ -124,6 +126,16 @@ void ArithmeticPathBase<element_type, scalar_type, path_type>::computeDerivative
     }
 }
 
+template <typename element_type, typename scalar_type, path_sz path_type>
+void ArithmeticPathBase<element_type, scalar_type, path_type>::reComputeLambda(const vector<scalar_type>& rmsd_between_refs) {
+    scalar_type mean_square_displacements = 0.0;
+    for (size_t i_frame = 1; i_frame < total_frames; ++i_frame) {
+        cvm::log(std::string("Distance between frame ") + cvm::to_str(i_frame) + " and " + cvm::to_str(i_frame + 1) + " is " + cvm::to_str(rmsd_between_refs[i_frame - 1]) + std::string("\n"));
+        mean_square_displacements += rmsd_between_refs[i_frame - 1] * rmsd_between_refs[i_frame - 1];
+    }
+    mean_square_displacements /= scalar_type(total_frames - 1);
+    lambda = 1.0 / mean_square_displacements;
+}
 }
 
 #endif // ARITHMETICPATHCV_H

lib/colvars/colvaratoms.cpp

@@ -37,6 +37,7 @@ cvm::atom::atom(int atom_number)
   }
   id = p->get_atom_id(index);
   update_mass();
+  update_charge();
   reset_data();
 }
 
@@ -53,6 +54,7 @@ cvm::atom::atom(cvm::residue_id const &residue,
   }
   id = p->get_atom_id(index);
   update_mass();
+  update_charge();
   reset_data();
 }
 
@@ -62,6 +64,7 @@ cvm::atom::atom(atom const &a)
 {
   id = (cvm::proxy)->get_atom_id(index);
   update_mass();
+  update_charge();
   reset_data();
 }
 
@@ -214,8 +217,6 @@ int cvm::atom_group::init()
   index = -1;
 
   b_dummy = false;
-  b_center = false;
-  b_rotate = false;
   b_user_defined_fit = false;
   fitting_group = NULL;
 
@@ -240,8 +241,9 @@ int cvm::atom_group::init_dependencies() {
     }
 
     init_feature(f_ag_active, "active", f_type_dynamic);
-    init_feature(f_ag_center, "translational_fit", f_type_static);
+    init_feature(f_ag_center, "center_to_reference", f_type_user);
-    init_feature(f_ag_rotate, "rotational_fit", f_type_static);
+    init_feature(f_ag_center_origin, "center_to_origin", f_type_user);
+    init_feature(f_ag_rotate, "rotate_to_origin", f_type_user);
     init_feature(f_ag_fitting_group, "fitting_group", f_type_static);
     init_feature(f_ag_explicit_gradient, "explicit_atom_gradient", f_type_dynamic);
     init_feature(f_ag_fit_gradients, "fit_gradients", f_type_user);
@@ -272,6 +274,10 @@ int cvm::atom_group::init_dependencies() {
 
   // Features that are implemented (or not) by all atom groups
   feature_states[f_ag_active].available = true;
+  feature_states[f_ag_center].available = true;
+  feature_states[f_ag_center_origin].available = true;
+  feature_states[f_ag_rotate].available = true;
+
   // f_ag_scalable_com is provided by the CVC iff it is COM-based
   feature_states[f_ag_scalable_com].available = false;
   // TODO make f_ag_scalable depend on f_ag_scalable_com (or something else)
@@ -317,6 +323,9 @@ void cvm::atom_group::update_total_mass()
       total_mass += ai->mass;
     }
   }
+  if (total_mass < 1e-15) {
+    cvm::error("ERROR: " + description + " has zero total mass.\n");
+  }
 }
 
 
@@ -381,12 +390,22 @@ int cvm::atom_group::parse(std::string const &group_conf)
 
   // We need to know about fitting to decide whether the group is scalable
   // and we need to know about scalability before adding atoms
-  bool b_defined_center = get_keyval(group_conf, "centerReference", b_center, false);
+  bool b_defined_center = get_keyval_feature(this, group_conf, "centerToOrigin", f_ag_center_origin, false);
-  bool b_defined_rotate = get_keyval(group_conf, "rotateReference", b_rotate, false);
+  // Legacy alias
+  b_defined_center |= get_keyval_feature(this, group_conf, "centerReference", f_ag_center, is_enabled(f_ag_center_origin), parse_deprecated);
+  b_defined_center |= get_keyval_feature(this, group_conf, "centerToReference", f_ag_center, is_enabled(f_ag_center));
+
+  if (is_enabled(f_ag_center_origin) && ! is_enabled(f_ag_center)) {
+    return cvm::error("centerToReference may not be disabled if centerToOrigin is enabled.\n");
+  }
+  // Legacy alias
+  bool b_defined_rotate = get_keyval_feature(this, group_conf, "rotateReference", f_ag_rotate, false, parse_deprecated);
+  b_defined_rotate |= get_keyval_feature(this, group_conf, "rotateToReference", f_ag_rotate, is_enabled(f_ag_rotate));
+
   // is the user setting explicit options?
   b_user_defined_fit = b_defined_center || b_defined_rotate;
 
-  if (is_available(f_ag_scalable_com) && !b_rotate && !b_center) {
+  if (is_available(f_ag_scalable_com) && !is_enabled(f_ag_rotate) && !is_enabled(f_ag_center)) {
     enable(f_ag_scalable_com);
     enable(f_ag_scalable);
   }
@@ -507,12 +526,8 @@ int cvm::atom_group::parse(std::string const &group_conf)
 
     // whether these atoms will ever receive forces or not
     bool enable_forces = true;
-    // disableForces is deprecated
+    get_keyval(group_conf, "enableForces", enable_forces, true, colvarparse::parse_silent);
-    if (get_keyval(group_conf, "enableForces", enable_forces, true)) {
     noforce = !enable_forces;
-    } else {
-      get_keyval(group_conf, "disableForces", noforce, false, colvarparse::parse_silent);
-    }
   }
 
   // Now that atoms are defined we can parse the detailed fitting options
@@ -760,10 +775,10 @@ std::string const cvm::atom_group::print_atom_ids() const
 
 int cvm::atom_group::parse_fitting_options(std::string const &group_conf)
 {
-  if (b_center || b_rotate) {
+  if (is_enabled(f_ag_center) || is_enabled(f_ag_rotate)) {
 
     if (b_dummy)
-      cvm::error("Error: centerReference or rotateReference "
+      cvm::error("Error: centerToReference or rotateToReference "
                  "cannot be defined for a dummy atom.\n");
 
     bool b_ref_pos_group = false;
@@ -827,7 +842,7 @@ int cvm::atom_group::parse_fitting_options(std::string const &group_conf)
 
     if (ref_pos.size()) {
 
-      if (b_rotate) {
+      if (is_enabled(f_ag_rotate)) {
         if (ref_pos.size() != group_for_fit->size())
           cvm::error("Error: the number of reference positions provided("+
                      cvm::to_str(ref_pos.size())+
@@ -846,7 +861,7 @@ int cvm::atom_group::parse_fitting_options(std::string const &group_conf)
       return COLVARS_ERROR;
     }
 
-    if (b_rotate && !noforce) {
+    if (is_enabled(f_ag_rotate) && !noforce) {
       cvm::log("Warning: atom group \""+key+
                "\" will be aligned to a fixed orientation given by the reference positions provided.  "
                "If the internal structure of the group changes too much (i.e. its RMSD is comparable "
@@ -865,7 +880,7 @@ int cvm::atom_group::parse_fitting_options(std::string const &group_conf)
     bool b_fit_gradients;
     get_keyval(group_conf, "enableFitGradients", b_fit_gradients, true);
 
-    if (b_fit_gradients && (b_center || b_rotate)) {
+    if (b_fit_gradients && (is_enabled(f_ag_center) || is_enabled(f_ag_rotate))) {
       enable(f_ag_fit_gradients);
     }
   }
@@ -879,7 +894,7 @@ void cvm::atom_group::do_feature_side_effects(int id)
   // If enabled features are changed upstream, the features below should be refreshed
   switch (id) {
     case f_ag_fit_gradients:
-      if (b_center || b_rotate) {
+      if (is_enabled(f_ag_center) || is_enabled(f_ag_rotate)) {
         atom_group *group_for_fit = fitting_group ? fitting_group : this;
         group_for_fit->fit_gradients.assign(group_for_fit->size(), cvm::atom_pos(0.0, 0.0, 0.0));
         rot.request_group1_gradients(group_for_fit->size());
@@ -971,7 +986,7 @@ int cvm::atom_group::calc_required_properties()
   calc_center_of_geometry();
 
   if (!is_enabled(f_ag_scalable)) {
-    if (b_center || b_rotate) {
+    if (is_enabled(f_ag_center) || is_enabled(f_ag_rotate)) {
       if (fitting_group) {
         fitting_group->calc_center_of_geometry();
       }
@@ -1001,7 +1016,7 @@ void cvm::atom_group::calc_apply_roto_translation()
     fitting_group->cog_orig = fitting_group->center_of_geometry();
   }
 
-  if (b_center) {
+  if (is_enabled(f_ag_center)) {
     // center on the origin first
     cvm::atom_pos const rpg_cog = fitting_group ?
       fitting_group->center_of_geometry() : this->center_of_geometry();
@@ -1011,9 +1026,9 @@ void cvm::atom_group::calc_apply_roto_translation()
     }
   }
 
-  if (b_rotate) {
+  if (is_enabled(f_ag_rotate)) {
-    // rotate the group (around the center of geometry if b_center is
+    // rotate the group (around the center of geometry if f_ag_center is
-    // true, around the origin otherwise)
+    // enabled, around the origin otherwise)
     rot.calc_optimal_rotation(fitting_group ?
                               fitting_group->positions() :
                               this->positions(),
@@ -1030,7 +1045,7 @@ void cvm::atom_group::calc_apply_roto_translation()
     }
   }
 
-  if (b_center) {
+  if (is_enabled(f_ag_center) && !is_enabled(f_ag_center_origin)) {
     // align with the center of geometry of ref_pos
     apply_translation(ref_pos_cog);
     if (fitting_group) {
@@ -1063,7 +1078,7 @@ void cvm::atom_group::read_velocities()
 {
   if (b_dummy) return;
 
-  if (b_rotate) {
+  if (is_enabled(f_ag_rotate)) {
 
     for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
       ai->read_velocity();
@@ -1084,7 +1099,7 @@ void cvm::atom_group::read_total_forces()
 {
   if (b_dummy) return;
 
-  if (b_rotate) {
+  if (is_enabled(f_ag_rotate)) {
 
     for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
       ai->read_total_force();
@@ -1172,14 +1187,14 @@ void cvm::atom_group::calc_fit_gradients()
 
   cvm::atom_group *group_for_fit = fitting_group ? fitting_group : this;
 
-  if (b_center) {
+  if (is_enabled(f_ag_center)) {
     // add the center of geometry contribution to the gradients
     cvm::rvector atom_grad;
 
     for (size_t i = 0; i < this->size(); i++) {
       atom_grad += atoms[i].grad;
     }
-    if (b_rotate) atom_grad = (rot.inverse()).rotate(atom_grad);
+    if (is_enabled(f_ag_rotate)) atom_grad = (rot.inverse()).rotate(atom_grad);
     atom_grad *= (-1.0)/(cvm::real(group_for_fit->size()));
 
     for (size_t j = 0; j < group_for_fit->size(); j++) {
@@ -1187,7 +1202,7 @@ void cvm::atom_group::calc_fit_gradients()
     }
   }
 
-  if (b_rotate) {
+  if (is_enabled(f_ag_rotate)) {
 
     // add the rotation matrix contribution to the gradients
     cvm::rotation const rot_inv = rot.inverse();
@@ -1196,7 +1211,7 @@ void cvm::atom_group::calc_fit_gradients()
 
       // compute centered, unrotated position
       cvm::atom_pos const pos_orig =
-        rot_inv.rotate((b_center ? (atoms[i].pos - ref_pos_cog) : (atoms[i].pos)));
+        rot_inv.rotate((is_enabled(f_ag_center) ? (atoms[i].pos - ref_pos_cog) : (atoms[i].pos)));
 
       // calculate \partial(R(q) \vec{x}_i)/\partial q) \cdot \partial\xi/\partial\vec{x}_i
       cvm::quaternion const dxdq =
@@ -1340,7 +1355,7 @@ void cvm::atom_group::apply_colvar_force(cvm::real const &force)
     return;
   }
 
-  if (b_rotate) {
+  if (is_enabled(f_ag_rotate)) {
 
     // rotate forces back to the original frame
     cvm::rotation const rot_inv = rot.inverse();
@@ -1355,7 +1370,7 @@ void cvm::atom_group::apply_colvar_force(cvm::real const &force)
     }
   }
 
-  if ((b_center || b_rotate) && is_enabled(f_ag_fit_gradients)) {
+  if ((is_enabled(f_ag_center) || is_enabled(f_ag_rotate)) && is_enabled(f_ag_fit_gradients)) {
 
     atom_group *group_for_fit = fitting_group ? fitting_group : this;
 
@@ -1386,7 +1401,7 @@ void cvm::atom_group::apply_force(cvm::rvector const &force)
     return;
   }
 
-  if (b_rotate) {
+  if (is_enabled(f_ag_rotate)) {
 
     cvm::rotation const rot_inv = rot.inverse();
     for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {

lib/colvars/colvaratoms.h

@@ -101,19 +101,15 @@ public:
   inline void update_mass()
   {
     colvarproxy *p = cvm::proxy;
-    if (p->updated_masses()) {
     mass = p->get_atom_mass(index);
   }
-  }
 
   /// Get the latest value of the charge
   inline void update_charge()
   {
     colvarproxy *p = cvm::proxy;
-    if (p->updated_charges()) {
     charge = p->get_atom_charge(index);
   }
-  }
 
   /// Get the current position
   inline void read_position()
@@ -328,35 +324,23 @@ public:
   /// If yes, returns 1-based number of a common atom; else, returns 0
   static int overlap(const atom_group &g1, const atom_group &g2);
 
-  /// \brief When updating atomic coordinates, translate them to align with the
+  /// The rotation calculated automatically if f_ag_rotate is defined
-  /// center of mass of the reference coordinates
-  bool b_center;
-
-  /// \brief When updating atom coordinates (and after
-  /// centering them if b_center is set), rotate the group to
-  /// align with the reference coordinates.
-  ///
-  /// Note: gradients will be calculated in the rotated frame: when
-  /// forces will be applied, they will rotated back to the original
-  /// frame
-  bool b_rotate;
-  /// The rotation calculated automatically if b_rotate is defined
   cvm::rotation rot;
 
-  /// \brief Indicates that the user has explicitly set centerReference or
+  /// \brief Indicates that the user has explicitly set centerToReference or
   /// rotateReference, and the corresponding reference:
   /// cvc's (eg rmsd, eigenvector) will not override the user's choice
   bool b_user_defined_fit;
 
-  /// \brief use reference coordinates for b_center or b_rotate
+  /// \brief use reference coordinates for f_ag_center or f_ag_rotate
   std::vector<cvm::atom_pos> ref_pos;
 
   /// \brief Center of geometry of the reference coordinates; regardless
-  /// of whether b_center is true, ref_pos is centered to zero at
+  /// of whether f_ag_center is true, ref_pos is centered to zero at
   /// initialization, and ref_pos_cog serves to center the positions
   cvm::atom_pos              ref_pos_cog;
 
-  /// \brief If b_center or b_rotate is true, use this group to
+  /// \brief If f_ag_center or f_ag_rotate is true, use this group to
   /// define the transformation (default: this group itself)
   atom_group                *fitting_group;
 
@@ -395,12 +379,12 @@ public:
   void apply_translation(cvm::rvector const &t);
 
   /// \brief Get the current velocities; this must be called always
-  /// *after* read_positions(); if b_rotate is defined, the same
+  /// *after* read_positions(); if f_ag_rotate is defined, the same
   /// rotation applied to the coordinates will be used
   void read_velocities();
 
   /// \brief Get the current total_forces; this must be called always
-  /// *after* read_positions(); if b_rotate is defined, the same
+  /// *after* read_positions(); if f_ag_rotate is defined, the same
   /// rotation applied to the coordinates will be used
   void read_total_forces();
 

lib/colvars/colvarbias_abf.cpp

@@ -136,10 +136,11 @@ int colvarbias_abf::init(std::string const &conf)
       colvars[i]->enable(f_cv_hide_Jacobian);
     }
 
-    // If any colvar is extended-system, we need to collect the extended
+    // If any colvar is extended-system (restrained style, not external with constraint), we are running eABF
-    // system gradient
+    if (colvars[i]->is_enabled(f_cv_extended_Lagrangian)
-    if (colvars[i]->is_enabled(f_cv_extended_Lagrangian))
+        && !colvars[i]->is_enabled(f_cv_external)) {
       b_extended = true;
+    }
 
     // Cannot mix and match coarse time steps with ABF because it gives
     // wrong total force averages - total force needs to be averaged over
@@ -475,7 +476,7 @@ int colvarbias_abf::update()
     last_gradients->copy_grid(*gradients);
     last_samples->copy_grid(*samples);
     shared_last_step = cvm::step_absolute();
-    cvm::log("Prepared sample and gradient buffers at step "+cvm::to_str(cvm::step_absolute())+".");
+    cvm::log("Prepared sample and gradient buffers at step "+cvm::to_str(cvm::step_absolute())+".\n");
   }
 
   // update UI estimator every step
@@ -812,8 +813,10 @@ int colvarbias_abf::write_output_files()
     cvm::log("ABF bias trying to write gradients and samples to disk");
   }
 
-  if (shared_on && cvm::main()->proxy->replica_index() > 0) {
+  if (shared_on && cvm::main()->proxy->replica_index() > 0
+    && ! (b_CZAR_estimator || b_UI_estimator) ) {
     // No need to report the same data as replica 0, let it do the I/O job
+    // except if using an eABF FE estimator
     return COLVARS_OK;
   }
 

lib/colvars/colvarbias_abf.h

@@ -121,6 +121,8 @@ private:
     } else {
       system_force[i] = colvars[i]->total_force().real_value
         - colvar_forces[i].real_value;
+        // If hideJacobian is active then total_force has an extra term of -fj
+        // which is the Jacobian-compensating force at the colvar level
     }
     if (cvm::debug())
       cvm::log("ABF System force calc: cv " + cvm::to_str(i) +

lib/colvars/colvarbias_histogram.cpp

@@ -129,14 +129,14 @@ int colvarbias_histogram::update()
     // output_prefix is unset during the constructor
     if (cvm::step_relative() == 0) {
       out_name = cvm::output_prefix() + "." + this->name + ".dat";
-      cvm::log("Histogram " + this->name + " will be written to file \"" + out_name + "\"");
+      cvm::log("Histogram " + this->name + " will be written to file \"" + out_name + "\"\n");
     }
   }
 
   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 + "\"");
+      cvm::log("Histogram " + this->name + " will be written to file \"" + out_name_dx + "\"\n");
     }
   }
 

lib/colvars/colvarbias_histogram.h

@@ -35,7 +35,6 @@ protected:
   colvar_grid_scalar *grid;
   std::vector<int> bin;
   std::string out_name, out_name_dx;
-  size_t output_freq;
 
   /// If one or more of the variables are \link colvarvalue::type_vector \endlink, treat them as arrays of this length
   size_t colvar_array_size;

lib/colvars/colvarbias_meta.cpp

@@ -50,6 +50,7 @@ colvarbias_meta::colvarbias_meta(char const *key)
 
   dump_fes = true;
   keep_hills = false;
+  restart_keep_hills = false;
   dump_fes_save = false;
   dump_replica_fes = false;
 
@@ -902,7 +903,7 @@ void colvarbias_meta::project_hills(colvarbias_meta::hill_iter  h_first,
           count++) {
       size_t i;
       for (i = 0; i < num_variables(); i++) {
-        new_colvar_values[i] = hills_energy->bin_to_value_scalar(he_ix[i], i);
+        new_colvar_values[i] = he->bin_to_value_scalar(he_ix[i], i);
       }
 
       // loop over the hills and increment the energy grid locally
@@ -934,40 +935,12 @@ void colvarbias_meta::project_hills(colvarbias_meta::hill_iter  h_first,
     }
 
   } else {
-
+    cvm::error("No grid object provided in metadynamics::project_hills()\n",
-    // TODO delete this (never used)
+               BUG_ERROR);
-
-    // simpler version, with just the energy
-
-    for ( ; (he->index_ok(he_ix)); ) {
-
-      for (size_t i = 0; i < num_variables(); i++) {
-        new_colvar_values[i] = hills_energy->bin_to_value_scalar(he_ix[i], i);
-      }
-
-      hills_energy_here = 0.0;
-      calc_hills(h_first, h_last, hills_energy_here, &new_colvar_values);
-      he->acc_value(he_ix, hills_energy_here);
-
-      he->incr(he_ix);
-
-      count++;
-      if ((count % print_frequency) == 0) {
-        if (print_progress) {
-          cvm::real const progress = cvm::real(count) / cvm::real(he->number_of_points());
-          std::ostringstream os;
-          os.setf(std::ios::fixed, std::ios::floatfield);
-          os << std::setw(6) << std::setprecision(2)
-             << 100.0 * progress
-             << "% done.";
-          cvm::log(os.str());
-        }
-      }
-    }
   }
 
   if (print_progress) {
-    cvm::log("100.00% done.");
+    cvm::log("100.00% done.\n");
   }
 
   if (! keep_hills) {
@@ -1281,9 +1254,25 @@ int colvarbias_meta::set_state_params(std::string const &state_conf)
     return error_code;
   }
 
+  colvarparse::get_keyval(state_conf, "keepHills", restart_keep_hills, false,
+                          colvarparse::parse_restart);
+
+  if ((!restart_keep_hills) && (cvm::main()->restart_version_number() < 20210604)) {
+    if (keep_hills) {
+      cvm::log("Warning: could not ensure that keepHills was enabled when "
+               "this state file was written; because it is enabled now, "
+               "it is assumed that it was also then, but please verify.\n");
+      restart_keep_hills = true;
+    }
+  } else {
+    if (restart_keep_hills) {
+      cvm::log("This state file/stream contains explicit hills.\n");
+    }
+  }
+
   std::string check_replica = "";
   if (colvarparse::get_keyval(state_conf, "replicaID", check_replica,
-                              std::string(""), colvarparse::parse_silent) &&
+                              std::string(""), colvarparse::parse_restart) &&
       (check_replica != this->replica_id)) {
     return cvm::error("Error: in the state file , the "
                       "\"metadynamics\" block has a different replicaID ("+
@@ -1297,8 +1286,6 @@ int colvarbias_meta::set_state_params(std::string const &state_conf)
 
 std::istream & colvarbias_meta::read_state_data(std::istream& is)
 {
-  bool grids_from_restart_file = use_grids;
-
   if (use_grids) {
 
     if (expand_grids) {
@@ -1343,7 +1330,6 @@ std::istream & colvarbias_meta::read_state_data(std::istream& is)
                !(hills_energy->read_restart(is))) {
       is.clear();
       is.seekg(hills_energy_pos, std::ios::beg);
-      grids_from_restart_file = false;
       if (!rebin_grids) {
         if (hills_energy_backup == NULL)
           cvm::fatal_error("Error: couldn't read the free energy grid for metadynamics bias \""+
@@ -1381,7 +1367,6 @@ std::istream & colvarbias_meta::read_state_data(std::istream& is)
                !(hills_energy_gradients->read_restart(is))) {
       is.clear();
       is.seekg(hills_energy_gradients_pos, std::ios::beg);
-      grids_from_restart_file = false;
       if (!rebin_grids) {
         if (hills_energy_backup == NULL)
           cvm::fatal_error("Error: couldn't read the free energy gradients grid for metadynamics bias \""+
@@ -1419,28 +1404,29 @@ std::istream & colvarbias_meta::read_state_data(std::istream& is)
     }
   }
 
+  // Save references to the end of the list of existing hills, so that it can
+  // be cleared if hills are read successfully state
   bool const existing_hills = !hills.empty();
   size_t const old_hills_size = hills.size();
   hill_iter old_hills_end = hills.end();
   hill_iter old_hills_off_grid_end = hills_off_grid.end();
 
-  // read the hills explicitly written (if there are any)
+  // Read any hills following the grid data (if any)
   while (read_hill(is)) {
-    if (cvm::debug())
+    if (cvm::debug()) {
       cvm::log("Read a previously saved hill under the "
                "metadynamics bias \""+
                this->name+"\", created at step "+
                cvm::to_str((hills.back()).it)+".\n");
     }
+  }
   is.clear();
   new_hills_begin = hills.end();
-  if (grids_from_restart_file) {
+  cvm::log("Read "+cvm::to_str(hills.size() - old_hills_size)+" hills.\n");
-    if (hills.size() > old_hills_size)
+
-      cvm::log("Read "+cvm::to_str(hills.size())+
+  if (existing_hills) {
-               " hills in addition to the grids.\n");
+    hills.erase(hills.begin(), old_hills_end);
-  } else {
+    hills_off_grid.erase(hills_off_grid.begin(), old_hills_off_grid_end);
-    if (!hills.empty())
-      cvm::log("Read "+cvm::to_str(hills.size())+" hills.\n");
   }
 
   if (rebin_grids) {
@@ -1454,7 +1440,16 @@ std::istream & colvarbias_meta::read_state_data(std::istream& is)
     colvar_grid_gradient *new_hills_energy_gradients =
       new colvar_grid_gradient(colvars);
 
-    if (!grids_from_restart_file || (keep_hills && !hills.empty())) {
+    if (cvm::debug()) {
+      std::ostringstream tmp_os;
+      tmp_os << "hills_energy parameters:\n";
+      hills_energy->write_params(tmp_os);
+      tmp_os << "new_hills_energy parameters:\n";
+      new_hills_energy->write_params(tmp_os);
+      cvm::log(tmp_os.str());
+    }
+
+    if (restart_keep_hills && !hills.empty()) {
       // if there are hills, recompute the new grids from them
       cvm::log("Rebinning the energy and forces grids from "+
                cvm::to_str(hills.size())+" hills (this may take a while)...\n");
@@ -1494,11 +1489,6 @@ std::istream & colvarbias_meta::read_state_data(std::istream& is)
   if (cvm::debug())
     cvm::log("colvarbias_meta::read_restart() done\n");
 
-  if (existing_hills) {
-    hills.erase(hills.begin(), old_hills_end);
-    hills_off_grid.erase(hills_off_grid.begin(), old_hills_off_grid_end);
-  }
-
   has_data = true;
 
   if (comm != single_replica) {
@@ -1509,6 +1499,15 @@ std::istream & colvarbias_meta::read_state_data(std::istream& is)
 }
 
 
+inline std::istream & reset_istream(std::istream &is, size_t start_pos)
+{
+  is.clear();
+  is.seekg(start_pos, std::ios::beg);
+  is.setstate(std::ios::failbit);
+  return is;
+}
+
+
 std::istream & colvarbias_meta::read_hill(std::istream &is)
 {
   if (!is) return is; // do nothing if failbit is set
@@ -1518,45 +1517,72 @@ std::istream & colvarbias_meta::read_hill(std::istream &is)
 
   std::string data;
   if ( !(is >> read_block("hill", &data)) ) {
-    is.clear();
+    return reset_istream(is, start_pos);
-    is.seekg(start_pos, std::ios::beg);
-    is.setstate(std::ios::failbit);
-    return is;
   }
 
+  std::istringstream data_is(data);
+
   cvm::step_number h_it = 0L;
-  get_keyval(data, "step", h_it, h_it, parse_restart);
+  cvm::real h_weight;
-  if (h_it <= state_file_step) {
+  std::vector<colvarvalue> h_centers(num_variables());
+  for (i = 0; i < num_variables(); i++) {
+    h_centers[i].type(variables(i)->value());
+  }
+  std::vector<cvm::real> h_sigmas(num_variables());
+  std::string h_replica;
+
+  std::string keyword;
+  while (data_is >> keyword) {
+
+    if (keyword == "step") {
+      if ( !(data_is >> h_it)) {
+        return reset_istream(is, start_pos);
+      }
+      if ((h_it <= state_file_step) && !restart_keep_hills) {
         if (cvm::debug())
           cvm::log("Skipping a hill older than the state file for metadynamics bias \""+
                    this->name+"\""+
                    ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+"\n");
         return is;
       }
-
-  cvm::real h_weight;
-  get_keyval(data, "weight", h_weight, hill_weight, parse_restart);
-
-  std::vector<colvarvalue> h_centers(num_variables());
-  for (i = 0; i < num_variables(); i++) {
-    h_centers[i].type(variables(i)->value());
     }
-  get_keyval(data, "centers", h_centers, h_centers, parse_restart);
 
-  std::vector<cvm::real> h_sigmas(num_variables());
+    if (keyword == "weight") {
-  get_keyval(data, "widths", h_sigmas, h_sigmas, parse_restart);
+      if ( !(data_is >> h_weight)) {
+        return reset_istream(is, start_pos);
+      }
+    }
+
+    if (keyword == "centers") {
       for (i = 0; i < num_variables(); i++) {
+        if ( !(data_is >> h_centers[i])) {
+          return reset_istream(is, start_pos);
+        }
+      }
+    }
+
+    if (keyword == "widths") {
+      for (i = 0; i < num_variables(); i++) {
+        if ( !(data_is >> h_sigmas[i])) {
+          return reset_istream(is, start_pos);
+        }
         // For backward compatibility, read the widths instead of the sigmas
         h_sigmas[i] /= 2.0;
       }
+    }
 
-  std::string h_replica = "";
     if (comm != single_replica) {
-    get_keyval(data, "replicaID", h_replica, replica_id, parse_restart);
+      if (keyword == "replicaID") {
-    if (h_replica != replica_id)
+        if ( !(data_is >> h_replica)) {
-      cvm::fatal_error("Error: trying to read a hill created by replica \""+h_replica+
+          return reset_istream(is, start_pos);
-                       "\" for replica \""+replica_id+
+        }
-                       "\"; did you swap output files?\n");
+        if (h_replica != replica_id) {
+          cvm::error("Error: trying to read a hill created by replica \""+
+                     h_replica+"\" for replica \""+replica_id+
+                     "\"; did you swap output files?\n", INPUT_ERROR);
+        }
+      }
+    }
   }
 
   hill_iter const hills_end = hills.end();
@@ -1699,8 +1725,12 @@ std::string const colvarbias_meta::hills_traj_file_name() const
 std::string const colvarbias_meta::get_state_params() const
 {
   std::ostringstream os;
-  if (this->comm != single_replica)
+  if (keep_hills) {
+    os << "keepHills on" << "\n";
+  }
+  if (this->comm != single_replica) {
     os << "replicaID " << this->replica_id << "\n";
+  }
   return (colvarbias::get_state_params() + os.str());
 }
 
@@ -1944,6 +1974,7 @@ colvarbias_meta::hill::hill(cvm::step_number it_in,
     sigmas(cv_values.size()),
     replica(replica_in)
 {
+  hill_value = 0.0;
   for (size_t i = 0; i < cv_values.size(); i++) {
     centers[i].type(cv_values[i]);
     centers[i] = cv_values[i];
@@ -1967,7 +1998,9 @@ colvarbias_meta::hill::hill(colvarbias_meta::hill const &h)
     centers(h.centers),
     sigmas(h.sigmas),
     replica(h.replica)
-{}
+{
+  hill_value = 0.0;
+}
 
 
 colvarbias_meta::hill::~hill()

lib/colvars/colvarbias_meta.h

@@ -155,10 +155,12 @@ protected:
   /// \brief How often the hills should be projected onto the grids
   size_t     grids_freq;
 
-  /// \brief Whether to keep the hills in the restart file (e.g. to do
+  /// Keep hills in the restart file (e.g. to accurately rebin later)
-  /// meaningful accurate rebinning afterwards)
   bool       keep_hills;
 
+  /// value of keepHills saved in the most recent restart file
+  bool restart_keep_hills;
+
   /// \brief Dump the free energy surface (.pmf file) every restartFrequency
   bool       dump_fes;
 

lib/colvars/colvarcomp.cpp

@@ -51,7 +51,7 @@ int colvar::cvc::init(std::string const &conf)
   std::string const old_name(name);
 
   if (name.size() > 0) {
-    cvm::log("Updating configuration for component \""+name+"\"");
+    cvm::log("Updating configuration for component \""+name+"\"\n");
   }
 
   if (get_keyval(conf, "name", name, name)) {
@@ -112,7 +112,7 @@ int colvar::cvc::init_total_force_params(std::string const &conf)
   }
   if (get_keyval_feature(this, conf, "oneSiteTotalForce",
                          f_cvc_one_site_total_force, is_enabled(f_cvc_one_site_total_force))) {
-    cvm::log("Computing total force on group 1 only");
+    cvm::log("Computing total force on group 1 only\n");
   }
 
   if (! is_enabled(f_cvc_one_site_total_force)) {
@@ -426,7 +426,7 @@ void colvar::cvc::collect_gradients(std::vector<int> const &atom_ids, std::vecto
 
     // If necessary, apply inverse rotation to get atomic
     // gradient in the laboratory frame
-    if (ag.b_rotate) {
+    if (ag.is_enabled(f_ag_rotate)) {
       cvm::rotation const rot_inv = ag.rot.inverse();
 
       for (size_t k = 0; k < ag.size(); k++) {
@@ -505,7 +505,7 @@ void colvar::cvc::debug_gradients()
     cvm::atom_pos fit_gradient_sum, gradient_sum;
 
     // print the values of the fit gradients
-    if (group->b_rotate || group->b_center) {
+    if (group->is_enabled(f_ag_center) || group->is_enabled(f_ag_rotate)) {
       if (group->is_enabled(f_ag_fit_gradients)) {
         size_t j;
 
@@ -514,7 +514,7 @@ void colvar::cvc::debug_gradients()
         for (j = 0; j < group_for_fit->fit_gradients.size(); j++) {
           cvm::log((group->fitting_group ? std::string("refPosGroup") : group->key) +
                   "[" + cvm::to_str(j) + "] = " +
-                  (group->b_rotate ?
+                  (group->is_enabled(f_ag_rotate) ?
                     cvm::to_str(rot_0.rotate(group_for_fit->fit_gradients[j])) :
                     cvm::to_str(group_for_fit->fit_gradients[j])));
         }
@@ -525,7 +525,7 @@ void colvar::cvc::debug_gradients()
     for (size_t ia = 0; ia < group->size(); ia++) {
 
       // tests are best conducted in the unrotated (simulation) frame
-      cvm::rvector const atom_grad = (group->b_rotate ?
+      cvm::rvector const atom_grad = (group->is_enabled(f_ag_rotate) ?
                                       rot_inv.rotate((*group)[ia].grad) :
                                       (*group)[ia].grad);
       gradient_sum += atom_grad;

lib/colvars/colvarcomp.h

@@ -271,6 +271,12 @@ public:
   /// \brief Whether or not this CVC will be computed in parallel whenever possible
   bool b_try_scalable;
 
+  /// Forcibly set value of CVC - useful for driving an external coordinate,
+  /// eg. lambda dynamics
+  inline void set_value(colvarvalue const &new_value) {
+    x = new_value;
+  }
+
 protected:
 
   /// \brief Cached value
@@ -1344,6 +1350,71 @@ public:
 };
 
 
+class colvar::euler_phi
+  : public colvar::orientation
+{
+public:
+  euler_phi(std::string const &conf);
+  euler_phi();
+  virtual int init(std::string const &conf);
+  virtual ~euler_phi() {}
+  virtual void calc_value();
+  virtual void calc_gradients();
+  virtual void apply_force(colvarvalue const &force);
+  virtual cvm::real dist2(colvarvalue const &x1,
+                          colvarvalue const &x2) const;
+  virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
+                                  colvarvalue const &x2) const;
+  virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
+                                  colvarvalue const &x2) const;
+  /// Redefined to handle the 2*PI periodicity
+  virtual void wrap(colvarvalue &x_unwrapped) const;
+};
+
+
+class colvar::euler_psi
+  : public colvar::orientation
+{
+public:
+  euler_psi(std::string const &conf);
+  euler_psi();
+  virtual int init(std::string const &conf);
+  virtual ~euler_psi() {}
+  virtual void calc_value();
+  virtual void calc_gradients();
+  virtual void apply_force(colvarvalue const &force);
+  virtual cvm::real dist2(colvarvalue const &x1,
+                          colvarvalue const &x2) const;
+  virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
+                                  colvarvalue const &x2) const;
+  virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
+                                  colvarvalue const &x2) const;
+  /// Redefined to handle the 2*PI periodicity
+  virtual void wrap(colvarvalue &x_unwrapped) const;
+};
+
+
+class colvar::euler_theta
+  : public colvar::orientation
+{
+public:
+  euler_theta(std::string const &conf);
+  euler_theta();
+  virtual int init(std::string const &conf);
+  virtual ~euler_theta() {}
+  virtual void calc_value();
+  virtual void calc_gradients();
+  virtual void apply_force(colvarvalue const &force);
+  // theta angle is a scalar variable and not periodic
+  // we need to override the virtual functions from orientation
+  virtual cvm::real dist2(colvarvalue const &x1,
+                          colvarvalue const &x2) const;
+  virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
+                                  colvarvalue const &x2) const;
+  virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
+                                  colvarvalue const &x2) const;
+};
+
 
 /// \brief Colvar component: root mean square deviation (RMSD) of a
 /// group with respect to a set of reference coordinates; uses \link
@@ -1406,6 +1477,28 @@ public:
 };
 
 
+// \brief Colvar component: alch_lambda
+// To communicate value with back-end in lambda-dynamics
+class colvar::alch_lambda
+  : public colvar::cvc
+{
+protected:
+  // No atom groups needed
+public:
+  alch_lambda(std::string const &conf);
+  alch_lambda();
+  virtual ~alch_lambda() {}
+  virtual void calc_value();
+  virtual void calc_gradients();
+  virtual void apply_force(colvarvalue const &force);
+  virtual cvm::real dist2(colvarvalue const &x1,
+                          colvarvalue const &x2) const;
+  virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
+                                  colvarvalue const &x2) const;
+  virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
+                                  colvarvalue const &x2) const;
+};
+
 
 class colvar::componentDisabled
   : public colvar::cvc
@@ -1681,11 +1774,20 @@ public:
 
 protected:
 
-  /// Identifier of the map object (as used by the simulation engine)
+  /// String identifier of the map object (as used by the simulation engine)
-  std::string map_name;
+  std::string volmap_name;
+
+  /// Numeric identifier of the map object (as used by the simulation engine)
+  int volmap_id;
 
   /// Index of the map objet in the proxy arrays
   int volmap_index;
+
+  /// Group of atoms selected internally (optional)
+  cvm::atom_group *atoms;
+
+  /// Weights assigned to each atom (default: uniform weights)
+  std::vector<cvm::real> atom_weights;
 };
 
 

lib/colvars/colvarcomp_alchlambda.cpp

@@ -0,0 +1,56 @@
+// -*- c++ -*-
+
+// 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
+
+#include <algorithm>
+
+#include "colvarmodule.h"
+#include "colvarvalue.h"
+#include "colvarparse.h"
+#include "colvar.h"
+#include "colvarcomp.h"
+
+
+colvar::alch_lambda::alch_lambda(std::string const &conf)
+  : cvc(conf)
+{
+  function_type = "alch_lambda";
+
+  disable(f_cvc_explicit_gradient);
+  disable(f_cvc_gradient);
+
+  x.type(colvarvalue::type_scalar);
+  // Query initial value from back-end
+  cvm::proxy->get_alch_lambda(&x.real_value);
+}
+
+
+void colvar::alch_lambda::calc_value()
+{
+  // Special workflow:
+  // at the beginning of the timestep we get a force instead of calculating the value
+
+  cvm::proxy->get_dE_dLambda(&ft.real_value);
+  ft.real_value *= -1.0; // Energy derivative to force
+}
+
+
+void colvar::alch_lambda::calc_gradients()
+{
+}
+
+
+void colvar::alch_lambda::apply_force(colvarvalue const &force)
+{
+  // Special workflow:
+  // at the end of the time step we send a new value
+  cvm::proxy->set_alch_lambda(&x.real_value);
+}
+
+simple_scalar_dist_functions(alch_lambda)

lib/colvars/colvarcomp_apath.cpp

@@ -19,17 +19,8 @@ colvar::aspathCV::aspathCV(std::string const &conf): CVBasedPath(conf) {
     get_keyval(conf, "weights", p_weights, std::vector<cvm::real>(cv.size(), 1.0));
     x.type(colvarvalue::type_scalar);
     use_explicit_gradients = true;
-    std::vector<cvm::real> rmsd_between_refs(total_reference_frames - 1, 0.0);
-    computeDistanceBetweenReferenceFrames(rmsd_between_refs);
-    cvm::real mean_square_displacements = 0.0;
-    for (size_t i_frame = 1; i_frame < total_reference_frames; ++i_frame) {
-        cvm::log(std::string("Distance between frame ") + cvm::to_str(i_frame) + " and " + cvm::to_str(i_frame + 1) + " is " + cvm::to_str(rmsd_between_refs[i_frame - 1]) + std::string("\n"));
-        mean_square_displacements += rmsd_between_refs[i_frame - 1] * rmsd_between_refs[i_frame - 1];
-    }
-    mean_square_displacements /= cvm::real(total_reference_frames - 1);
-    cvm::real suggested_lambda = 1.0 / mean_square_displacements;
     cvm::real p_lambda;
-    get_keyval(conf, "lambda", p_lambda, suggested_lambda);
+    get_keyval(conf, "lambda", p_lambda, -1.0);
     ArithmeticPathCV::ArithmeticPathBase<colvarvalue, cvm::real, ArithmeticPathCV::path_sz::S>::initialize(cv.size(), total_reference_frames, p_lambda, ref_cv[0], p_weights);
     cvm::log(std::string("Lambda is ") + cvm::to_str(lambda) + std::string("\n"));
     for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
@@ -58,6 +49,16 @@ void colvar::aspathCV::updateDistanceToReferenceFrames() {
 }
 
 void colvar::aspathCV::calc_value() {
+    if (lambda < 0) {
+        // this implies that the user may not set a valid lambda value
+        // so recompute it by the suggested value in Parrinello's paper
+        cvm::log("A non-positive value of lambda is detected, which implies that it may not set in the configuration.\n");
+        cvm::log("This component (aspathCV) will recompute a value for lambda following the suggestion in the origin paper.\n");
+        std::vector<cvm::real> rmsd_between_refs(total_reference_frames - 1, 0.0);
+        computeDistanceBetweenReferenceFrames(rmsd_between_refs);
+        reComputeLambda(rmsd_between_refs);
+        cvm::log("Ok, the value of lambda is updated to " + cvm::to_str(lambda));
+    }
     computeValue();
     x = s;
 }
@@ -107,17 +108,8 @@ colvar::azpathCV::azpathCV(std::string const &conf): CVBasedPath(conf) {
     get_keyval(conf, "weights", p_weights, std::vector<cvm::real>(cv.size(), 1.0));
     x.type(colvarvalue::type_scalar);
     use_explicit_gradients = true;
-    std::vector<cvm::real> rmsd_between_refs(total_reference_frames - 1, 0.0);
-    computeDistanceBetweenReferenceFrames(rmsd_between_refs);
-    cvm::real mean_square_displacements = 0.0;
-    for (size_t i_frame = 1; i_frame < total_reference_frames; ++i_frame) {
-        cvm::log(std::string("Distance between frame ") + cvm::to_str(i_frame) + " and " + cvm::to_str(i_frame + 1) + " is " + cvm::to_str(rmsd_between_refs[i_frame - 1]) + std::string("\n"));
-        mean_square_displacements += rmsd_between_refs[i_frame - 1] * rmsd_between_refs[i_frame - 1];
-    }
-    mean_square_displacements /= cvm::real(total_reference_frames - 1);
-    cvm::real suggested_lambda = 1.0 / mean_square_displacements;
     cvm::real p_lambda;
-    get_keyval(conf, "lambda", p_lambda, suggested_lambda);
+    get_keyval(conf, "lambda", p_lambda, -1.0);
     ArithmeticPathCV::ArithmeticPathBase<colvarvalue, cvm::real, ArithmeticPathCV::path_sz::Z>::initialize(cv.size(), total_reference_frames, p_lambda, ref_cv[0], p_weights);
     cvm::log(std::string("Lambda is ") + cvm::to_str(lambda) + std::string("\n"));
     for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
@@ -146,6 +138,16 @@ void colvar::azpathCV::updateDistanceToReferenceFrames() {
 }
 
 void colvar::azpathCV::calc_value() {
+    if (lambda < 0) {
+        // this implies that the user may not set a valid lambda value
+        // so recompute it by the suggested value in Parrinello's paper
+        cvm::log("A non-positive value of lambda is detected, which implies that it may not set in the configuration.\n");
+        cvm::log("This component (azpathCV) will recompute a value for lambda following the suggestion in the origin paper.\n");
+        std::vector<cvm::real> rmsd_between_refs(total_reference_frames - 1, 0.0);
+        computeDistanceBetweenReferenceFrames(rmsd_between_refs);
+        reComputeLambda(rmsd_between_refs);
+        cvm::log("Ok, the value of lambda is updated to " + cvm::to_str(lambda));
+    }
     computeValue();
     x = z;
 }

lib/colvars/colvarcomp_distances.cpp

@@ -194,10 +194,10 @@ colvar::distance_z::distance_z(std::string const &conf)
   ref2 = parse_group(conf, "ref2", true);
 
   if ( ref2 ) {
-    cvm::log("Using axis joining the centers of mass of groups \"ref\" and \"ref2\"");
+    cvm::log("Using axis joining the centers of mass of groups \"ref\" and \"ref2\"\n");
     fixed_axis = false;
     if (key_lookup(conf, "axis"))
-      cvm::log("Warning: explicit axis definition will be ignored!");
+      cvm::log("Warning: explicit axis definition will be ignored!\n");
   } else {
     if (get_keyval(conf, "axis", axis, cvm::rvector(0.0, 0.0, 1.0))) {
       if (axis.norm2() == 0.0) {
@@ -808,9 +808,9 @@ colvar::gyration::gyration(std::string const &conf)
   atoms = parse_group(conf, "atoms");
 
   if (atoms->b_user_defined_fit) {
-    cvm::log("WARNING: explicit fitting parameters were provided for atom group \"atoms\".");
+    cvm::log("WARNING: explicit fitting parameters were provided for atom group \"atoms\".\n");
   } else {
-    atoms->b_center = true;
+    atoms->enable(f_ag_center);
     atoms->ref_pos.assign(1, cvm::atom_pos(0.0, 0.0, 0.0));
     atoms->fit_gradients.assign(atoms->size(), cvm::rvector(0.0, 0.0, 0.0));
   }
@@ -1025,13 +1025,13 @@ colvar::rmsd::rmsd(std::string const &conf)
   }
 
   if (atoms->b_user_defined_fit) {
-    cvm::log("WARNING: explicit fitting parameters were provided for atom group \"atoms\".");
+    cvm::log("WARNING: explicit fitting parameters were provided for atom group \"atoms\".\n");
   } else {
     // Default: fit everything
-    cvm::log("Enabling \"centerReference\" and \"rotateReference\", to minimize RMSD before calculating it as a variable: "
+    cvm::log("Enabling \"centerToReference\" and \"rotateToReference\", to minimize RMSD before calculating it as a variable: "
               "if this is not the desired behavior, disable them explicitly within the \"atoms\" block.\n");
-    atoms->b_center = true;
+    atoms->enable(f_ag_center);
-    atoms->b_rotate = true;
+    atoms->enable(f_ag_rotate);
     // default case: reference positions for calculating the rmsd are also those used
     // for fitting
     atoms->ref_pos = ref_pos;
@@ -1156,7 +1156,7 @@ void colvar::rmsd::calc_Jacobian_derivative()
   cvm::real rotation_term = 0.0;
 
   // The rotation term only applies is coordinates are rotated
-  if (atoms->b_rotate) {
+  if (atoms->is_enabled(f_ag_rotate)) {
 
     // gradient of the rotation matrix
     cvm::matrix2d<cvm::rvector> grad_rot_mat(3, 3);
@@ -1202,7 +1202,7 @@ void colvar::rmsd::calc_Jacobian_derivative()
   }
 
   // The translation term only applies is coordinates are centered
-  cvm::real translation_term = atoms->b_center ? 3.0 : 0.0;
+  cvm::real translation_term = atoms->is_enabled(f_ag_center) ? 3.0 : 0.0;
 
   jd.real_value = x.real_value > 0.0 ?
     (3.0 * atoms->size() - 1.0 - translation_term - rotation_term) / x.real_value :
@@ -1284,10 +1284,10 @@ colvar::eigenvector::eigenvector(std::string const &conf)
     cvm::log("WARNING: explicit fitting parameters were provided for atom group \"atoms\".\n");
   } else {
     // default: fit everything
-    cvm::log("Enabling \"centerReference\" and \"rotateReference\", to minimize RMSD before calculating the vector projection: "
+    cvm::log("Enabling \"centerToReference\" and \"rotateToReference\", to minimize RMSD before calculating the vector projection: "
               "if this is not the desired behavior, disable them explicitly within the \"atoms\" block.\n");
-    atoms->b_center = true;
+    atoms->enable(f_ag_center);
-    atoms->b_rotate = true;
+    atoms->enable(f_ag_rotate);
     atoms->ref_pos = ref_pos;
     atoms->center_ref_pos();
     atoms->disable(f_ag_fit_gradients); // cancel out if group is fitted on itself
@@ -1355,14 +1355,14 @@ colvar::eigenvector::eigenvector(std::string const &conf)
 
   if (b_difference_vector) {
 
-    if (atoms->b_center) {
+    if (atoms->is_enabled(f_ag_center)) {
       // both sets should be centered on the origin for fitting
       for (size_t i = 0; i < atoms->size(); i++) {
         eigenvec[i] -= eig_center;
         ref_pos[i]  -= ref_pos_center;
       }
     }
-    if (atoms->b_rotate) {
+    if (atoms->is_enabled(f_ag_rotate)) {
       atoms->rot.calc_optimal_rotation(eigenvec, ref_pos);
       for (size_t i = 0; i < atoms->size(); i++) {
         eigenvec[i] = atoms->rot.rotate(eigenvec[i]);
@@ -1372,7 +1372,7 @@ colvar::eigenvector::eigenvector(std::string const &conf)
     for (size_t i = 0; i < atoms->size(); i++) {
       eigenvec[i] -= ref_pos[i];
     }
-    if (atoms->b_center) {
+    if (atoms->is_enabled(f_ag_center)) {
       // bring back the ref positions to where they were
       for (size_t i = 0; i < atoms->size(); i++) {
         ref_pos[i] += ref_pos_center;
@@ -1521,7 +1521,8 @@ colvar::cartesian::cartesian(std::string const &conf)
 
   x.type(colvarvalue::type_vector);
   disable(f_cvc_explicit_gradient);
-  x.vector1d_value.resize(atoms->size() * axes.size());
+  // Don't try to access atoms if creation of the atom group failed
+  if (atoms != NULL) x.vector1d_value.resize(atoms->size() * axes.size());
 }
 
 

lib/colvars/colvarcomp_gpath.cpp

@@ -65,8 +65,8 @@ colvar::CartesianBasedPath::CartesianBasedPath(std::string const &conf): cvc(con
         cvm::atom_group* tmp_atoms = parse_group(conf, "atoms");
         if (!has_user_defined_fitting) {
             // Swipe from the rmsd class
-            tmp_atoms->b_center = true;
+            tmp_atoms->enable(f_ag_center);
-            tmp_atoms->b_rotate = true;
+            tmp_atoms->enable(f_ag_rotate);
             tmp_atoms->ref_pos = reference_frames[i_frame];
             tmp_atoms->center_ref_pos();
             tmp_atoms->enable(f_ag_fit_gradients);
@@ -87,8 +87,8 @@ colvar::CartesianBasedPath::CartesianBasedPath(std::string const &conf): cvc(con
             std::vector<cvm::atom_pos> reference_fitting_position(tmp_fitting_atoms->size());
             cvm::load_coords(reference_position_filename.c_str(), &reference_fitting_position, tmp_fitting_atoms, reference_column, reference_column_value);
             // setup the atom group for calculating
-            tmp_atoms->b_center = true;
+            tmp_atoms->enable(f_ag_center);
-            tmp_atoms->b_rotate = true;
+            tmp_atoms->enable(f_ag_rotate);
             tmp_atoms->b_user_defined_fit = true;
             tmp_atoms->disable(f_ag_scalable);
             tmp_atoms->disable(f_ag_scalable_com);

lib/colvars/colvarcomp_rotations.cpp

@@ -442,3 +442,335 @@ void colvar::spin_angle::wrap(colvarvalue &x_unwrapped) const
 
   return;
 }
+
+
+colvar::euler_phi::euler_phi(std::string const &conf)
+  : orientation()
+{
+  function_type = "euler_phi";
+  period = 360.0;
+  enable(f_cvc_periodic);
+  enable(f_cvc_explicit_gradient);
+  x.type(colvarvalue::type_scalar);
+  init(conf);
+}
+
+
+colvar::euler_phi::euler_phi()
+  : orientation()
+{
+  function_type = "euler_phi";
+  period = 360.0;
+  enable(f_cvc_periodic);
+  enable(f_cvc_explicit_gradient);
+  x.type(colvarvalue::type_scalar);
+}
+
+
+int colvar::euler_phi::init(std::string const &conf)
+{
+  int error_code = COLVARS_OK;
+  error_code |= orientation::init(conf);
+  return error_code;
+}
+
+
+void colvar::euler_phi::calc_value()
+{
+  atoms_cog = atoms->center_of_geometry();
+
+  rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog));
+
+  const cvm::real& q0 = rot.q.q0;
+  const cvm::real& q1 = rot.q.q1;
+  const cvm::real& q2 = rot.q.q2;
+  const cvm::real& q3 = rot.q.q3;
+  const cvm::real tmp_y = 2 * (q0 * q1 + q2 * q3);
+  const cvm::real tmp_x = 1 - 2 * (q1 * q1 + q2 * q2);
+  x.real_value = cvm::atan2(tmp_y, tmp_x) * (180.0/PI);
+}
+
+
+void colvar::euler_phi::calc_gradients()
+{
+  const cvm::real& q0 = rot.q.q0;
+  const cvm::real& q1 = rot.q.q1;
+  const cvm::real& q2 = rot.q.q2;
+  const cvm::real& q3 = rot.q.q3;
+  const cvm::real denominator = (2 * q0 * q1 + 2 * q2 * q3) * (2 * q0 * q1 + 2 * q2 * q3) + (-2 * q1 * q1 - 2 * q2 * q2 + 1) * (-2 * q1 * q1 - 2 * q2 * q2 + 1);
+  const cvm::real dxdq0 = (180.0/PI) * 2 * q1 * (-2 * q1 * q1 - 2 * q2 * q2 + 1) / denominator;
+  const cvm::real dxdq1 = (180.0/PI) * (2 * q0 * (-2 * q1 * q1 - 2 * q2 * q2 + 1) - 4 * q1 * (-2 * q0 * q1 - 2 * q2 * q3)) / denominator;
+  const cvm::real dxdq2 = (180.0/PI) * (-4 * q2 * (-2 * q0 * q1 - 2 * q2 * q3) + 2 * q3 * (-2 * q1 * q1 - 2 * q2 * q2 + 1)) / denominator;
+  const cvm::real dxdq3 = (180.0/PI) * 2 * q2 * (-2 * q1 * q1 - 2 * q2 * q2 + 1) / denominator;
+  for (size_t ia = 0; ia < atoms->size(); ia++) {
+    (*atoms)[ia].grad = (dxdq0 * (rot.dQ0_2[ia])[0]) +
+                        (dxdq1 * (rot.dQ0_2[ia])[1]) +
+                        (dxdq2 * (rot.dQ0_2[ia])[2]) +
+                        (dxdq3 * (rot.dQ0_2[ia])[3]);
+  }
+}
+
+
+void colvar::euler_phi::apply_force(colvarvalue const &force)
+{
+  cvm::real const &fw = force.real_value;
+  if (!atoms->noforce) {
+    atoms->apply_colvar_force(fw);
+  }
+}
+
+
+cvm::real colvar::euler_phi::dist2(colvarvalue const &x1,
+                                  colvarvalue const &x2) const
+{
+  cvm::real diff = x1.real_value - x2.real_value;
+  diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
+  return diff * diff;
+}
+
+
+colvarvalue colvar::euler_phi::dist2_lgrad(colvarvalue const &x1,
+                                          colvarvalue const &x2) const
+{
+  cvm::real diff = x1.real_value - x2.real_value;
+  diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
+  return 2.0 * diff;
+}
+
+
+colvarvalue colvar::euler_phi::dist2_rgrad(colvarvalue const &x1,
+                                          colvarvalue const &x2) const
+{
+  cvm::real diff = x1.real_value - x2.real_value;
+  diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
+  return (-2.0) * diff;
+}
+
+
+void colvar::euler_phi::wrap(colvarvalue &x_unwrapped) const
+{
+  if ((x_unwrapped.real_value - wrap_center) >= 180.0) {
+    x_unwrapped.real_value -= 360.0;
+    return;
+  }
+
+  if ((x_unwrapped.real_value - wrap_center) < -180.0) {
+    x_unwrapped.real_value += 360.0;
+    return;
+  }
+
+  return;
+}
+
+
+colvar::euler_psi::euler_psi(std::string const &conf)
+  : orientation()
+{
+  function_type = "euler_psi";
+  period = 360.0;
+  enable(f_cvc_periodic);
+  enable(f_cvc_explicit_gradient);
+  x.type(colvarvalue::type_scalar);
+  init(conf);
+}
+
+
+colvar::euler_psi::euler_psi()
+  : orientation()
+{
+  function_type = "euler_psi";
+  period = 360.0;
+  enable(f_cvc_periodic);
+  enable(f_cvc_explicit_gradient);
+  x.type(colvarvalue::type_scalar);
+}
+
+
+int colvar::euler_psi::init(std::string const &conf)
+{
+  int error_code = COLVARS_OK;
+  error_code |= orientation::init(conf);
+  return error_code;
+}
+
+
+void colvar::euler_psi::calc_value()
+{
+  atoms_cog = atoms->center_of_geometry();
+
+  rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog));
+
+  const cvm::real& q0 = rot.q.q0;
+  const cvm::real& q1 = rot.q.q1;
+  const cvm::real& q2 = rot.q.q2;
+  const cvm::real& q3 = rot.q.q3;
+  const cvm::real tmp_y = 2 * (q0 * q3 + q1 * q2);
+  const cvm::real tmp_x = 1 - 2 * (q2 * q2 + q3 * q3);
+  x.real_value = cvm::atan2(tmp_y, tmp_x) * (180.0/PI);
+}
+
+
+void colvar::euler_psi::calc_gradients()
+{
+  const cvm::real& q0 = rot.q.q0;
+  const cvm::real& q1 = rot.q.q1;
+  const cvm::real& q2 = rot.q.q2;
+  const cvm::real& q3 = rot.q.q3;
+  const cvm::real denominator = (2 * q0 * q3 + 2 * q1 * q2) * (2 * q0 * q3 + 2 * q1 * q2) + (-2 * q2 * q2 - 2 * q3 * q3 + 1) * (-2 * q2 * q2 - 2 * q3 * q3 + 1);
+  const cvm::real dxdq0 = (180.0/PI) * 2 * q3 * (-2 * q2 * q2 - 2 * q3 * q3 + 1) / denominator;
+  const cvm::real dxdq1 = (180.0/PI) * 2 * q2 * (-2 * q2 * q2 - 2 * q3 * q3 + 1) / denominator;
+  const cvm::real dxdq2 = (180.0/PI) * (2 * q1 * (-2 * q2 * q2 - 2 * q3 * q3 + 1) - 4 * q2 * (-2 * q0 * q3 - 2 * q1 * q2)) / denominator;
+  const cvm::real dxdq3 = (180.0/PI) * (2 * q0 * (-2 * q2 * q2 - 2 * q3 * q3 + 1) - 4 * q3 * (-2 * q0 * q3 - 2 * q1 * q2)) / denominator;
+  for (size_t ia = 0; ia < atoms->size(); ia++) {
+    (*atoms)[ia].grad = (dxdq0 * (rot.dQ0_2[ia])[0]) +
+                        (dxdq1 * (rot.dQ0_2[ia])[1]) +
+                        (dxdq2 * (rot.dQ0_2[ia])[2]) +
+                        (dxdq3 * (rot.dQ0_2[ia])[3]);
+  }
+}
+
+
+void colvar::euler_psi::apply_force(colvarvalue const &force)
+{
+  cvm::real const &fw = force.real_value;
+  if (!atoms->noforce) {
+    atoms->apply_colvar_force(fw);
+  }
+}
+
+
+cvm::real colvar::euler_psi::dist2(colvarvalue const &x1,
+                                  colvarvalue const &x2) const
+{
+  cvm::real diff = x1.real_value - x2.real_value;
+  diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
+  return diff * diff;
+}
+
+
+colvarvalue colvar::euler_psi::dist2_lgrad(colvarvalue const &x1,
+                                          colvarvalue const &x2) const
+{
+  cvm::real diff = x1.real_value - x2.real_value;
+  diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
+  return 2.0 * diff;
+}
+
+
+colvarvalue colvar::euler_psi::dist2_rgrad(colvarvalue const &x1,
+                                          colvarvalue const &x2) const
+{
+  cvm::real diff = x1.real_value - x2.real_value;
+  diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
+  return (-2.0) * diff;
+}
+
+
+void colvar::euler_psi::wrap(colvarvalue &x_unwrapped) const
+{
+  if ((x_unwrapped.real_value - wrap_center) >= 180.0) {
+    x_unwrapped.real_value -= 360.0;
+    return;
+  }
+
+  if ((x_unwrapped.real_value - wrap_center) < -180.0) {
+    x_unwrapped.real_value += 360.0;
+    return;
+  }
+
+  return;
+}
+
+
+colvar::euler_theta::euler_theta(std::string const &conf)
+  : orientation()
+{
+  function_type = "euler_theta";
+  enable(f_cvc_explicit_gradient);
+  x.type(colvarvalue::type_scalar);
+  init(conf);
+}
+
+
+colvar::euler_theta::euler_theta()
+  : orientation()
+{
+  function_type = "euler_theta";
+  enable(f_cvc_explicit_gradient);
+  x.type(colvarvalue::type_scalar);
+}
+
+
+int colvar::euler_theta::init(std::string const &conf)
+{
+  int error_code = COLVARS_OK;
+  error_code |= orientation::init(conf);
+  return error_code;
+}
+
+
+void colvar::euler_theta::calc_value()
+{
+  atoms_cog = atoms->center_of_geometry();
+
+  rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog));
+
+  const cvm::real& q0 = rot.q.q0;
+  const cvm::real& q1 = rot.q.q1;
+  const cvm::real& q2 = rot.q.q2;
+  const cvm::real& q3 = rot.q.q3;
+  x.real_value = cvm::asin(2 * (q0 * q2 - q3 * q1)) * (180.0/PI);
+}
+
+
+void colvar::euler_theta::calc_gradients()
+{
+  const cvm::real& q0 = rot.q.q0;
+  const cvm::real& q1 = rot.q.q1;
+  const cvm::real& q2 = rot.q.q2;
+  const cvm::real& q3 = rot.q.q3;
+  const cvm::real denominator = cvm::sqrt(1 - (2 * q0 * q2 - 2 * q1 * q3) * (2 * q0 * q2 - 2 * q1 * q3));
+  const cvm::real dxdq0 = (180.0/PI) * 2 * q2 / denominator;
+  const cvm::real dxdq1 = (180.0/PI) * -2 * q3 / denominator;
+  const cvm::real dxdq2 = (180.0/PI) * 2 * q0 / denominator;
+  const cvm::real dxdq3 = (180.0/PI) * -2 * q1 / denominator;
+  for (size_t ia = 0; ia < atoms->size(); ia++) {
+    (*atoms)[ia].grad = (dxdq0 * (rot.dQ0_2[ia])[0]) +
+                        (dxdq1 * (rot.dQ0_2[ia])[1]) +
+                        (dxdq2 * (rot.dQ0_2[ia])[2]) +
+                        (dxdq3 * (rot.dQ0_2[ia])[3]);
+  }
+}
+
+
+void colvar::euler_theta::apply_force(colvarvalue const &force)
+{
+  cvm::real const &fw = force.real_value;
+  if (!atoms->noforce) {
+    atoms->apply_colvar_force(fw);
+  }
+}
+
+
+cvm::real colvar::euler_theta::dist2(colvarvalue const &x1,
+                                     colvarvalue const &x2) const
+{
+  // theta angle is not periodic
+  return cvc::dist2(x1, x2);
+}
+
+
+colvarvalue colvar::euler_theta::dist2_lgrad(colvarvalue const &x1,
+                                             colvarvalue const &x2) const
+{
+  // theta angle is not periodic
+  return cvc::dist2_lgrad(x1, x2);
+}
+
+
+colvarvalue colvar::euler_theta::dist2_rgrad(colvarvalue const &x1,
+                                             colvarvalue const &x2) const
+{
+  // theta angle is not periodic
+  return cvc::dist2_rgrad(x1, x2);
+}

lib/colvars/colvarcomp_volmaps.cpp

@@ -19,6 +19,9 @@ colvar::map_total::map_total()
   : cvc(), volmap_index(-1)
 {
   function_type = "map_total";
+  volmap_id = -1;
+  volmap_index = -1;
+  atoms = NULL;
   x.type(colvarvalue::type_scalar);
 }
 
@@ -27,6 +30,9 @@ colvar::map_total::map_total(std::string const &conf)
   : cvc(), volmap_index(-1)
 {
   function_type = "map_total";
+  volmap_id = -1;
+  volmap_index = -1;
+  atoms = NULL;
   x.type(colvarvalue::type_scalar);
   map_total::init(conf);
 }
@@ -35,26 +41,101 @@ colvar::map_total::map_total(std::string const &conf)
 int colvar::map_total::init(std::string const &conf)
 {
   int error_code = cvc::init(conf);
-  get_keyval(conf, "mapName", map_name, map_name);
+  colvarproxy *proxy = cvm::main()->proxy;
-  volmap_index = (cvm::proxy)->init_volmap(map_name);
+  get_keyval(conf, "mapName", volmap_name, volmap_name);
+  get_keyval(conf, "mapID", volmap_id, volmap_id);
+  register_param("mapID", reinterpret_cast<void *>(&volmap_id));
+
+  if ((volmap_name.size() > 0) && (volmap_id >= 0)) {
+    error_code |=
+      cvm::error("Error: mapName and mapID are mutually exclusive.\n");
+  }
+
+  // Parse optional group
+  atoms = parse_group(conf, "atoms", true);
+  if (atoms != NULL) {
+
+    // Using internal selection
+    if (volmap_name.size()) {
+      error_code |= proxy->check_volmap_by_name(volmap_name);
+    }
+    if (volmap_id >= 0) {
+      error_code |= proxy->check_volmap_by_id(volmap_id);
+    }
+
+  } else {
+
+    // Using selection from the MD engine
+    if (volmap_name.size()) {
+      volmap_index = proxy->init_volmap_by_name(volmap_name);
+    }
+    if (volmap_id >= 0) {
+      volmap_index = proxy->init_volmap_by_id(volmap_id);
+    }
     error_code |= volmap_index > 0 ? COLVARS_OK : INPUT_ERROR;
+  }
+
+  if (get_keyval(conf, "atomWeights", atom_weights, atom_weights)) {
+    if (atoms == NULL) {
+      error_code |= cvm::error("Error: weights can only be assigned when atoms "
+                               "are selected explicitly in Colvars.\n",
+                               INPUT_ERROR);
+    } else {
+      if (atoms->size() != atom_weights.size()) {
+        error_code |= cvm::error("Error: if defined, the number of weights ("+
+                                 cvm::to_str(atom_weights.size())+
+                                 ") must equal the number of atoms ("+
+                                 cvm::to_str(atoms->size())+
+                                 ").\n", INPUT_ERROR);
+      }
+    }
+  }
+
+  if (volmap_name.size() > 0) {
+    volmap_id = proxy->get_volmap_id_from_name(volmap_name.c_str());
+  }
+
   return error_code;
 }
 
 
 void colvar::map_total::calc_value()
 {
-  x.real_value = (cvm::proxy)->get_volmap_value(volmap_index);
+  colvarproxy *proxy = cvm::main()->proxy;
+  int flags = is_enabled(f_cvc_gradient) ? colvarproxy::volmap_flag_gradients :
+    colvarproxy::volmap_flag_gradients;
+
+  if (atoms != NULL) {
+    // Compute the map inside Colvars
+    x.real_value = 0.0;
+
+    cvm::real *w = NULL;
+    if (atom_weights.size() > 0) {
+      flags |= colvarproxy::volmap_flag_use_atom_field;
+      w = &(atom_weights[0]);
+    }
+    proxy->compute_volmap(flags, volmap_id, atoms->begin(), atoms->end(),
+                          &(x.real_value), w);
+  } else {
+    // Get the externally computed value
+    x.real_value = proxy->get_volmap_value(volmap_index);
+  }
 }
 
 
 void colvar::map_total::calc_gradients()
 {
-  // Atomic coordinates are not available here
+  // Computed in calc_value() or by the MD engine
 }
 
 
 void colvar::map_total::apply_force(colvarvalue const &force)
 {
-  (cvm::proxy)->apply_volmap_force(volmap_index, force.real_value);
+  colvarproxy *proxy = cvm::main()->proxy;
+  if (atoms) {
+    if (!atoms->noforce)
+      atoms->apply_colvar_force(force.real_value);
+  } else {
+    proxy->apply_volmap_force(volmap_index, force.real_value);
+  }
 }

lib/colvars/colvardeps.cpp

@@ -49,7 +49,7 @@ void colvardeps::free_children_deps() {
   // Cannot be in the base class destructor because it needs the derived class features()
   size_t i,j,fid;
 
-  if (cvm::debug()) cvm::log("DEPS: freeing children deps for " + description);
+  if (cvm::debug()) cvm::log("DEPS: freeing children deps for " + description + "\n");
 
   cvm::increase_depth();
   for (fid = 0; fid < feature_states.size(); fid++) {
@@ -58,7 +58,7 @@ void colvardeps::free_children_deps() {
         int g = features()[fid]->requires_children[i];
         for (j=0; j<children.size(); j++) {
           if (cvm::debug()) cvm::log("DEPS: dereferencing children's "
-            + children[j]->features()[g]->description);
+            + children[j]->features()[g]->description + "\n");
           children[j]->decr_ref_count(g);
         }
       }
@@ -80,7 +80,7 @@ void colvardeps::restore_children_deps() {
         int g = features()[fid]->requires_children[i];
         for (j=0; j<children.size(); j++) {
           if (cvm::debug()) cvm::log("DEPS: re-enabling children's "
-            + children[j]->features()[g]->description);
+            + children[j]->features()[g]->description + "\n");
           children[j]->enable(g, false, false);
         }
       }
@@ -135,7 +135,7 @@ int colvardeps::enable(int feature_id,
   if (cvm::debug()) {
     cvm::log("DEPS: " + description +
       (dry_run ? " testing " : " enabling ") +
-      "\"" + f->description +"\"");
+      "\"" + f->description +"\"\n");
   }
 
   if (fs->enabled) {
@@ -144,7 +144,7 @@ int colvardeps::enable(int feature_id,
       // as requirement is enabled
       fs->ref_count++;
       if (cvm::debug())
-        cvm::log("DEPS: bumping ref_count to " + cvm::to_str(fs->ref_count));
+        cvm::log("DEPS: bumping ref_count to " + cvm::to_str(fs->ref_count) + "\n");
     }
     // Do not try to further resolve deps
     return COLVARS_OK;
@@ -243,7 +243,7 @@ int colvardeps::enable(int feature_id,
           enable(g, false, false); // Just for printing error output
         }
         cvm::decrease_depth();
-        cvm::log("-----------------------------------------");
+        cvm::log("-----------------------------------------\n");
         if (toplevel) {
           cvm::error("Error: Failed dependency in " + description + ".");
         }
@@ -285,7 +285,7 @@ int colvardeps::enable(int feature_id,
     do_feature_side_effects(feature_id);
     if (cvm::debug())
       cvm::log("DEPS: feature \"" + f->description + "\" in "
-        + description + " enabled, ref_count = 1.");
+        + description + " enabled, ref_count = 1." + "\n");
   }
   return COLVARS_OK;
 }
@@ -297,7 +297,7 @@ int colvardeps::disable(int feature_id) {
   feature_state *fs = &feature_states[feature_id];
 
   if (cvm::debug()) cvm::log("DEPS: disabling feature \""
-      + f->description + "\" in " + description);
+      + f->description + "\" in " + description + "\n");
 
   if (fs->enabled == false) {
     return COLVARS_OK;
@@ -313,14 +313,14 @@ int colvardeps::disable(int feature_id) {
   // internal deps (self)
   for (i=0; i<f->requires_self.size(); i++) {
     if (cvm::debug()) cvm::log("DEPS: dereferencing self "
-      + features()[f->requires_self[i]]->description);
+      + features()[f->requires_self[i]]->description + "\n");
     decr_ref_count(f->requires_self[i]);
   }
 
   // alternates
   for (i=0; i<fs->alternate_refs.size(); i++) {
     if (cvm::debug()) cvm::log("DEPS: dereferencing alt "
-      + features()[fs->alternate_refs[i]]->description);
+      + features()[fs->alternate_refs[i]]->description + "\n");
     decr_ref_count(fs->alternate_refs[i]);
   }
   // Forget these, now that they are dereferenced
@@ -337,7 +337,7 @@ int colvardeps::disable(int feature_id) {
       int g = f->requires_children[i];
       for (j=0; j<children.size(); j++) {
         if (cvm::debug()) cvm::log("DEPS: dereferencing children's "
-          + children[j]->features()[g]->description);
+          + children[j]->features()[g]->description + "\n");
         children[j]->decr_ref_count(g);
       }
     }
@@ -430,11 +430,13 @@ void colvardeps::require_feature_alt(int f, int g, int h, int i, int j) {
 
 void colvardeps::print_state() {
   size_t i;
-  cvm::log("Features of \"" + description + "\" ON/OFF (refcount)");
+  cvm::log("Features of \"" + description + "\" (refcount)\n");
   for (i = 0; i < feature_states.size(); i++) {
-    std::string onoff = is_enabled(i) ? "ON" : "OFF";
+    std::string onoff = is_enabled(i) ? "ON " : "   ";
-    cvm::log("- " + features()[i]->description + "   " + onoff + " ("
+    // Only display refcount if non-zero for less clutter
-      + cvm::to_str(feature_states[i].ref_count) + ")");
+    std::string refcount = feature_states[i].ref_count != 0 ?
+      " (" + cvm::to_str(feature_states[i].ref_count) + ") " : "";
+    cvm::log("- " + onoff + features()[i]->description + refcount + "\n");
   }
   cvm::increase_depth();
   for (i=0; i<children.size(); i++) {
@@ -460,7 +462,7 @@ void colvardeps::add_child(colvardeps *child) {
       for (i=0; i<features()[fid]->requires_children.size(); i++) {
         int g = features()[fid]->requires_children[i];
         if (cvm::debug()) cvm::log("DEPS: re-enabling children's "
-          + child->features()[g]->description);
+          + child->features()[g]->description + "\n");
         child->enable(g, false, false);
       }
     }

lib/colvars/colvardeps.h

@@ -283,6 +283,10 @@ public:
     /// center with fictitious mass; bias forces will be applied to
     /// the center
     f_cv_extended_Lagrangian,
+    /// \brief An extended variable that sets an external variable in the
+    /// back-end (eg. an alchemical coupling parameter for lambda-dynamics)
+    /// Can have a single component
+    f_cv_external,
     /// \brief The extended system coordinate undergoes Langevin dynamics
     f_cv_Langevin,
     /// \brief Output the potential and kinetic energies
@@ -375,6 +379,7 @@ public:
   enum features_atomgroup {
     f_ag_active,
     f_ag_center,
+    f_ag_center_origin,
     f_ag_rotate,
     f_ag_fitting_group,
     /// Perform a standard minimum msd fit for given atoms

lib/colvars/colvargrid.cpp

@@ -122,17 +122,83 @@ cvm::real colvar_grid_scalar::entropy() const
 
 
 colvar_grid_gradient::colvar_grid_gradient()
-  : colvar_grid<cvm::real>(), samples(NULL)
+  : colvar_grid<cvm::real>(),
+    samples(NULL),
+    weights(NULL)
 {}
 
 colvar_grid_gradient::colvar_grid_gradient(std::vector<int> const &nx_i)
-  : colvar_grid<cvm::real>(nx_i, 0.0, nx_i.size()), samples(NULL)
+  : colvar_grid<cvm::real>(nx_i, 0.0, nx_i.size()),
+    samples(NULL),
+    weights(NULL)
 {}
 
 colvar_grid_gradient::colvar_grid_gradient(std::vector<colvar *> &colvars)
-  : colvar_grid<cvm::real>(colvars, 0.0, colvars.size()), samples(NULL)
+  : colvar_grid<cvm::real>(colvars, 0.0, colvars.size()),
+    samples(NULL),
+    weights(NULL)
 {}
 
+
+colvar_grid_gradient::colvar_grid_gradient(std::string &filename)
+  : colvar_grid<cvm::real>(),
+    samples(NULL),
+    weights(NULL)
+{
+  std::ifstream is;
+  is.open(filename.c_str());
+  if (!is.is_open()) {
+    cvm::error("Error opening multicol gradient file " + filename + " for reading.\n");
+    return;
+  }
+
+  // Data in the header: nColvars, then for each
+  // xiMin, dXi, nPoints, periodic flag
+
+  std::string  hash;
+  size_t i;
+
+  if ( !(is >> hash) || (hash != "#") ) {
+    cvm::error("Error reading grid at position "+
+                cvm::to_str(static_cast<size_t>(is.tellg()))+
+                " in stream(read \"" + hash + "\")\n");
+    return;
+  }
+
+  is >> nd;
+  mult = nd;
+  std::vector<cvm::real> lower_in(nd), widths_in(nd);
+  std::vector<int>       nx_in(nd);
+  std::vector<int>       periodic_in(nd);
+
+  for (i = 0; i < nd; i++ ) {
+    if ( !(is >> hash) || (hash != "#") ) {
+      cvm::error("Error reading grid at position "+
+                  cvm::to_str(static_cast<size_t>(is.tellg()))+
+                  " in stream(read \"" + hash + "\")\n");
+      return;
+    }
+
+    is >> lower_in[i] >> widths_in[i] >> nx_in[i] >> periodic_in[i];
+  }
+
+  this->setup(nx_in, 0., mult);
+
+  widths = widths_in;
+
+  for (i = 0; i < nd; i++ ) {
+    lower_boundaries.push_back(colvarvalue(lower_in[i]));
+    periodic.push_back(static_cast<bool>(periodic_in[i]));
+  }
+
+  // Reset the istream for read_multicol, which expects the whole file
+  is.clear();
+  is.seekg(0);
+  read_multicol(is);
+  is.close();
+}
+
+
 void colvar_grid_gradient::write_1D_integral(std::ostream &os)
 {
   cvm::real bin, min, integral;
@@ -202,7 +268,7 @@ integrate_potential::integrate_potential(std::vector<colvar *> &colvars, colvar_
     // Compute inverse of Laplacian diagonal for Jacobi preconditioning
     // For now all code related to preconditioning is commented out
     // until a method better than Jacobi is implemented
-//     cvm::log("Preparing inverse diagonal for preconditioning...");
+//     cvm::log("Preparing inverse diagonal for preconditioning...\n");
 //     inv_lap_diag.resize(nt);
 //     std::vector<cvm::real> id(nt), lap_col(nt);
 //     for (int i = 0; i < nt; i++) {
@@ -213,7 +279,30 @@ integrate_potential::integrate_potential(std::vector<colvar *> &colvars, colvar_
 //       id[i] = 0.;
 //       inv_lap_diag[i] = 1. / lap_col[i];
 //     }
-//     cvm::log("Done.");
+//     cvm::log("Done.\n");
+  }
+}
+
+
+integrate_potential::integrate_potential(colvar_grid_gradient * gradients)
+  : gradients(gradients)
+{
+  nd = gradients->num_variables();
+  nx = gradients->number_of_points_vec();
+  widths = gradients->widths;
+  periodic = gradients->periodic;
+
+  // Expand grid by 1 bin in non-periodic dimensions
+  for (size_t i = 0; i < nd; i++ ) {
+    if (!periodic[i]) nx[i]++;
+    // Shift the grid by half the bin width (values at edges instead of center of bins)
+    lower_boundaries.push_back(gradients->lower_boundaries[i].real_value - 0.5 * widths[i]);
+  }
+
+  setup(nx);
+
+  if (nd > 1) {
+    divergence.resize(nt);
   }
 }
 
@@ -246,7 +335,7 @@ int integrate_potential::integrate(const int itmax, const cvm::real &tol, cvm::r
   } else if (nd <= 3) {
 
     nr_linbcg_sym(divergence, data, tol, itmax, iter, err);
-    cvm::log("Integrated in " + cvm::to_str(iter) + " steps, error: " + cvm::to_str(err));
+    cvm::log("Integrated in " + cvm::to_str(iter) + " steps, error: " + cvm::to_str(err) + "\n");
 
   } else {
     cvm::error("Cannot integrate PMF in dimension > 3\n");

lib/colvars/colvargrid.h

@@ -102,6 +102,12 @@ public:
     return nd;
   }
 
+  /// Return the numbers of points in all dimensions
+  inline std::vector<int> const &number_of_points_vec() const
+  {
+    return nx;
+  }
+
   /// Return the number of points in the i-th direction, if provided, or
   /// the total number
   inline size_t number_of_points(int const icv = -1) const
@@ -199,6 +205,7 @@ public:
   {
     nd = nt = 0;
     mult = 1;
+    has_parent_data = false;
     this->setup();
   }
 
@@ -222,9 +229,9 @@ public:
                                          hard_lower_boundaries(g.hard_lower_boundaries),
                                          hard_upper_boundaries(g.hard_upper_boundaries),
                                          widths(g.widths),
+                                         has_parent_data(false),
                                          has_data(false)
-  {
+  {}
-  }
 
   /// \brief Constructor from explicit grid sizes \param nx_i Number
   /// of grid points along each dimension \param t Initial value for
@@ -233,7 +240,7 @@ public:
   colvar_grid(std::vector<int> const &nx_i,
               T const &t = T(),
               size_t mult_i = 1)
-    : has_data(false)
+    : has_parent_data(false), has_data(false)
   {
     this->setup(nx_i, t, mult_i);
   }
@@ -245,7 +252,7 @@ public:
               T const &t = T(),
               size_t mult_i = 1,
               bool add_extra_bin = false)
-    : has_data(false)
+    : has_parent_data(false), has_data(false)
   {
     this->init_from_colvars(colvars, t, mult_i, add_extra_bin);
   }
@@ -1066,8 +1073,8 @@ public:
     std::vector<int>      nx_read;
     std::vector<int>      bin;
 
-    if ( cv.size() != nd ) {
+    if ( cv.size() > 0 && cv.size() != nd ) {
-      cvm::error("Cannot read grid file: missing reference to colvars.");
+      cvm::error("Cannot read grid file: number of variables in file differs from number referenced by grid.\n");
       return is;
     }
 
@@ -1525,6 +1532,9 @@ public:
   /// Constructor from a vector of colvars
   colvar_grid_gradient(std::vector<colvar *>  &colvars);
 
+  /// Constructor from a multicol file
+  colvar_grid_gradient(std::string &filename);
+
   /// \brief Get a vector with the binned value(s) indexed by ix, normalized if applicable
   inline void vector_value(std::vector<int> const &ix, std::vector<cvm::real> &v) const
   {
@@ -1658,10 +1668,13 @@ class integrate_potential : public colvar_grid_scalar
   {}
 
   /// Constructor from a vector of colvars + gradient grid
-  integrate_potential (std::vector<colvar *> &colvars, colvar_grid_gradient * gradients);
+  integrate_potential(std::vector<colvar *> &colvars, colvar_grid_gradient * gradients);
+
+  /// Constructor from a gradient grid (for processing grid files without a Colvars config)
+  integrate_potential(colvar_grid_gradient * gradients);
 
   /// \brief Calculate potential from divergence (in 2D); return number of steps
-  int integrate (const int itmax, const cvm::real & tol, cvm::real & err);
+  int integrate(const int itmax, const cvm::real & tol, cvm::real & err);
 
   /// \brief Update matrix containing divergence and boundary conditions
   /// based on new gradient point value, in neighboring bins

lib/colvars/colvarmodule.cpp

@@ -34,6 +34,9 @@ colvarmodule::colvarmodule(colvarproxy *proxy_in)
 
   xyz_reader_use_count = 0;
 
+  restart_version_str.clear();
+  restart_version_int = 0;
+
   if (proxy == NULL) {
     proxy = proxy_in; // Pointer to the proxy object
     parse = new colvarparse(); // Parsing object for global options
@@ -48,7 +51,7 @@ colvarmodule::colvarmodule(colvarproxy *proxy_in)
 
   cvm::log(cvm::line_marker);
   cvm::log("Initializing the collective variables module, version "+
-           cvm::to_str(COLVARS_VERSION)+".\n");
+           version()+".\n");
   cvm::log("Please cite Fiorin et al, Mol Phys 2013:\n "
            "https://dx.doi.org/10.1080/00268976.2013.813594\n"
            "in any publication based on this calculation.\n");
@@ -58,7 +61,7 @@ colvarmodule::colvarmodule(colvarproxy *proxy_in)
   }
 
 #if (__cplusplus >= 201103L)
-  cvm::log("This version was built with the C++11 standard or higher.");
+  cvm::log("This version was built with the C++11 standard or higher.\n");
 #else
   cvm::log("This version was built without the C++11 standard: some features are disabled.\n"
     "Please see the following link for details:\n"
@@ -186,6 +189,7 @@ std::istream & colvarmodule::getline(std::istream &is, std::string &line)
     size_t const sz = l.size();
     if (sz > 0) {
       if (l[sz-1] == '\r' ) {
+        // Replace Windows newlines with Unix newlines
         line = l.substr(0, sz-1);
       } else {
         line = l;
@@ -200,6 +204,7 @@ std::istream & colvarmodule::getline(std::istream &is, std::string &line)
 
 int colvarmodule::parse_config(std::string &conf)
 {
+  // Auto-generated additional configuration
   extra_conf.clear();
 
   // Check that the input has matching braces
@@ -208,6 +213,9 @@ int colvarmodule::parse_config(std::string &conf)
                       INPUT_ERROR);
   }
 
+  // Check that the input has only ASCII characters, and warn otherwise
+  colvarparse::check_ascii(conf);
+
   // Parse global options
   if (catch_input_errors(parse_global_params(conf))) {
     return get_error();
@@ -472,7 +480,7 @@ int colvarmodule::parse_biases(std::string const &conf)
   if (use_scripted_forces) {
     cvm::log(cvm::line_marker);
     cvm::increase_depth();
-    cvm::log("User forces script will be run at each bias update.");
+    cvm::log("User forces script will be run at each bias update.\n");
     cvm::decrease_depth();
   }
 
@@ -754,6 +762,9 @@ int colvarmodule::calc()
 
   error_code |= end_of_step();
 
+  // TODO move this to a base-class proxy method that calls this function
+  error_code |= proxy->end_of_step();
+
   return error_code;
 }
 
@@ -1311,21 +1322,23 @@ std::istream & colvarmodule::read_restart(std::istream &is)
                         colvarparse::parse_restart);
       it = it_restart;
 
-      std::string restart_version;
+      restart_version_str.clear();
-      int restart_version_int = 0;
+      restart_version_int = 0;
       parse->get_keyval(restart_conf, "version",
-                        restart_version, std::string(""),
+                        restart_version_str, std::string(""),
                         colvarparse::parse_restart);
-      if (restart_version.size()) {
+      if (restart_version_str.size()) {
-        if (restart_version != std::string(COLVARS_VERSION)) {
+        // Initialize integer version number of this restart file
-          cvm::log("This state file was generated with version "+
-                   restart_version+"\n");
-        }
         restart_version_int =
-          proxy->get_version_from_string(restart_version.c_str());
+          proxy->get_version_from_string(restart_version_str.c_str());
       }
 
-      if (restart_version_int < 20160810) {
+      if (restart_version() != version()) {
+        cvm::log("This state file was generated with version "+
+                 restart_version()+"\n");
+      }
+
+      if (restart_version_number() < 20160810) {
         // check for total force change
         if (proxy->total_forces_enabled()) {
           warn_total_forces = true;
@@ -1769,6 +1782,8 @@ int cvm::read_index_file(char const *filename)
     cvm::error("Error: in opening index file \""+
                std::string(filename)+"\".\n",
                FILE_ERROR);
+  } else {
+    index_file_names.push_back(std::string(filename));
   }
 
   while (is.good()) {
@@ -1861,6 +1876,7 @@ int colvarmodule::reset_index_groups()
   }
   index_group_names.clear();
   index_groups.clear();
+  index_file_names.clear();
   return COLVARS_OK;
 }
 
@@ -1924,48 +1940,75 @@ int cvm::load_coords_xyz(char const *filename,
   std::string line;
   cvm::real x = 0.0, y = 0.0, z = 0.0;
 
+  std::string const error_msg("Error: cannot parse XYZ file \""+
+                              std::string(filename)+"\".\n");
+
   if ( ! (xyz_is >> natoms) ) {
-    cvm::error("Error: cannot parse XYZ file "
+    return cvm::error(error_msg, INPUT_ERROR);
-               + std::string(filename) + ".\n", INPUT_ERROR);
   }
 
   ++xyz_reader_use_count;
   if (xyz_reader_use_count < 2) {
-    cvm::log("Warning: beginning from 2019-11-26 the XYZ file reader assumes Angstrom units.");
+    cvm::log("Warning: beginning from 2019-11-26 the XYZ file reader assumes Angstrom units.\n");
   }
 
+  if (xyz_is.good()) {
     // skip comment line
     cvm::getline(xyz_is, line);
     cvm::getline(xyz_is, line);
     xyz_is.width(255);
-  std::vector<atom_pos>::iterator pos_i = pos->begin();
+  } else {
+    return cvm::error(error_msg, INPUT_ERROR);
+  }
 
+  std::vector<atom_pos>::iterator pos_i = pos->begin();
+  size_t xyz_natoms = 0;
   if (pos->size() != natoms) { // Use specified indices
     int next = 0; // indices are zero-based
     std::vector<int>::const_iterator index = atoms->sorted_ids().begin();
-    for ( ; pos_i != pos->end() ; pos_i++, index++) {
 
+    for ( ; pos_i != pos->end() ; pos_i++, index++) {
       while ( next < *index ) {
         cvm::getline(xyz_is, line);
         next++;
       }
+      if (xyz_is.good()) {
         xyz_is >> symbol;
         xyz_is >> x >> y >> z;
         // XYZ files are assumed to be in Angstrom (as eg. VMD will)
         (*pos_i)[0] = proxy->angstrom_to_internal(x);
         (*pos_i)[1] = proxy->angstrom_to_internal(y);
         (*pos_i)[2] = proxy->angstrom_to_internal(z);
+        xyz_natoms++;
+      } else {
+        return cvm::error(error_msg, INPUT_ERROR);
       }
+    }
+
   } else {          // Use all positions
+
     for ( ; pos_i != pos->end() ; pos_i++) {
+      if (xyz_is.good()) {
         xyz_is >> symbol;
         xyz_is >> x >> y >> z;
         (*pos_i)[0] = proxy->angstrom_to_internal(x);
         (*pos_i)[1] = proxy->angstrom_to_internal(y);
         (*pos_i)[2] = proxy->angstrom_to_internal(z);
+        xyz_natoms++;
+      } else {
+        return cvm::error(error_msg, INPUT_ERROR);
       }
     }
-  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
+  }
+
+  if (xyz_natoms != pos->size()) {
+    return cvm::error("Error: The number of positions read from file \""+
+                      std::string(filename)+"\" does not match the number of "+
+                      "positions required: "+cvm::to_str(xyz_natoms)+" vs. "+
+                      cvm::to_str(pos->size())+".\n", INPUT_ERROR);
+  }
+
+  return COLVARS_OK;
 }
 
 

lib/colvars/colvarmodule.h

@@ -81,6 +81,12 @@ private:
 
 public:
 
+  /// Get the version string (YYYY-MM-DD format)
+  std::string version() const
+  {
+    return std::string(COLVARS_VERSION);
+  }
+
   /// Get the version number (higher = more recent)
   int version_number() const
   {
@@ -150,6 +156,12 @@ public:
     return ::cos(static_cast<double>(x));
   }
 
+  /// Reimplemented to work around MS compiler issues
+  static inline real asin(real const &x)
+  {
+    return ::asin(static_cast<double>(x));
+  }
+
   /// Reimplemented to work around MS compiler issues
   static inline real acos(real const &x)
   {
@@ -685,6 +697,9 @@ public:
   static rvector position_distance(atom_pos const &pos1,
                                    atom_pos const &pos2);
 
+  /// \brief Names of .ndx files that have been loaded
+  std::vector<std::string> index_file_names;
+
   /// \brief Names of groups from one or more Gromacs .ndx files
   std::vector<std::string> index_group_names;
 
@@ -758,7 +773,11 @@ protected:
   /// Write labels at the next iteration
   bool cv_traj_write_labels;
 
-private:
+  /// Version of the most recent state file read
+  std::string restart_version_str;
+
+  /// Integer version of the most recent state file read
+  int restart_version_int;
 
   /// Counter for the current depth in the object hierarchy (useg e.g. in output)
   size_t depth_s;
@@ -771,6 +790,18 @@ private:
 
 public:
 
+  /// Version of the most recent state file read
+  inline std::string restart_version() const
+  {
+    return restart_version_str;
+  }
+
+  /// Integer version of the most recent state file read
+  inline int restart_version_number() const
+  {
+    return restart_version_int;
+  }
+
   /// Get the current object depth in the hierarchy
   static size_t & depth();
 

lib/colvars/colvarmodule_utils.h

@@ -0,0 +1,80 @@
+// -*- c++ -*-
+
+// 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
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
+
+#ifndef COLVARMODULE_UTILS_H
+#define COLVARMODULE_UTILS_H
+
+
+#include "colvarmodule.h"
+
+
+template <typename T>
+cvm::real get_force_norm2(T const &x)
+{
+  return x.norm2();
+}
+
+
+template <>
+inline cvm::real get_force_norm2(cvm::real const &x)
+{
+  return x*x;
+}
+
+
+template <typename T, int flag, bool get_index>
+cvm::real compute_norm2_stats(std::vector<T> const &v,
+                              int *minmax_index = NULL)
+{
+  cvm::real result = 0.0;
+  if (flag == -1) {
+    // Initialize for minimum search, using approx. largest float32 value
+    result = 1.0e38;
+  }
+
+  typename std::vector<T>::const_iterator xi = v.begin();
+  size_t i = 0;
+
+  if (get_index) *minmax_index = -1; // Let's not assume minmax_index is initialized to -1
+
+  for ( ; xi != v.end(); xi++, i++) {
+    cvm::real const norm2 = get_force_norm2<T>(*xi);
+    if (flag == 0) {
+      result += norm2;
+    }
+    if (flag == 1) {
+      if (norm2 > result) {
+        result = norm2;
+        if (get_index) *minmax_index = i;
+      }
+    }
+    if (flag == -1) {
+      if (norm2 < result) {
+        result = norm2;
+        if (get_index) *minmax_index = i;
+      }
+    }
+  }
+
+  size_t const n = v.size();
+
+  if (flag == 0) {
+    if (n > 0) {
+      result /= cvm::real(n);
+    }
+  }
+
+  result = cvm::sqrt(result);
+
+  return result;
+}
+
+
+#endif

lib/colvars/colvarparse.cpp

@@ -125,6 +125,10 @@ void colvarparse::mark_key_set_user(std::string const &key_str,
     cvm::log("# "+key_str+" = "+cvm::to_str(value)+"\n",
              cvm::log_user_params());
   }
+  if (parse_mode & parse_deprecation_warning) {
+    cvm::log("Warning: keyword "+key_str+
+      " is deprecated. Check the documentation for the current equivalent.\n");
+  }
 }
 
 
@@ -919,6 +923,26 @@ int colvarparse::check_braces(std::string const &conf,
   return (brace_count != 0) ? INPUT_ERROR : COLVARS_OK;
 }
 
+
+int colvarparse::check_ascii(std::string const &conf)
+{
+  // Check for non-ASCII characters
+  std::string line;
+  std::istringstream is(conf);
+  while (cvm::getline(is, line)) {
+    unsigned char const * const uchars =
+      reinterpret_cast<unsigned char const *>(line.c_str());
+    for (size_t i = 0; i < line.size(); i++) {
+      if (uchars[i] & 0x80U) {
+        cvm::log("Warning: non-ASCII character detected in this line: \""+
+                 line+"\".\n");
+      }
+    }
+  }
+  return COLVARS_OK;
+}
+
+
 void colvarparse::split_string(const std::string& data, const std::string& delim, std::vector<std::string>& dest) {
     size_t index = 0, new_index = 0;
     std::string tmpstr;

lib/colvars/colvarparse.h

@@ -56,6 +56,8 @@ public:
     parse_echo = (1<<1),
     /// Print the default value of a keyword, if it is NOT given
     parse_echo_default = (1<<2),
+    /// Print a deprecation warning if the keyword is given
+    parse_deprecation_warning = (1<<3),
     /// Do not print the keyword
     parse_silent = 0,
     /// Raise error if the keyword is not provided
@@ -66,7 +68,9 @@ public:
     /// The call is being executed from a read_restart() function
     parse_restart = (1<<18),
     /// Alias for old default behavior (should be phased out)
-    parse_normal = (1<<2) | (1<<1) | (1<<17)
+    parse_normal = (1<<2) | (1<<1) | (1<<17),
+    /// Settings for a deprecated keyword
+    parse_deprecated = (1<<1) | (1<<3) | (1<<17)
   };
 
   /// \brief Check that all the keywords within "conf" are in the list
@@ -317,6 +321,10 @@ public:
   /// from this position
   static int check_braces(std::string const &conf, size_t const start_pos);
 
+  /// \brief Check that a config string contains non-ASCII characters
+  /// \param conf The configuration string
+  static int check_ascii(std::string const &conf);
+
   /// \brief Split a string with a specified delimiter into a vector
   /// \param data The string to be splitted
   /// \param delim A delimiter

lib/colvars/colvarproxy.cpp

@@ -24,12 +24,15 @@
 #include "colvarproxy.h"
 #include "colvarscript.h"
 #include "colvaratoms.h"
+#include "colvarmodule_utils.h"
 
 
 
 colvarproxy_system::colvarproxy_system()
 {
   angstrom_value = 0.0;
+  kcal_mol_value = 0.0;
+  boundaries_type = boundaries_unsupported;
   total_force_requested = false;
   reset_pbc_lattice();
 }
@@ -38,6 +41,46 @@ colvarproxy_system::colvarproxy_system()
 colvarproxy_system::~colvarproxy_system() {}
 
 
+int colvarproxy_system::set_unit_system(std::string const & /* units */,
+                                        bool /* check_only */)
+{
+  return COLVARS_NOT_IMPLEMENTED;
+}
+
+
+cvm::real colvarproxy_system::backend_angstrom_value()
+{
+  return 1.0;
+}
+
+
+cvm::real colvarproxy_system::boltzmann()
+{
+  return 0.001987191;
+}
+
+
+cvm::real colvarproxy_system::temperature()
+{
+  // TODO define, document and implement a user method to set the value of this
+  return 300.0;
+}
+
+
+cvm::real colvarproxy_system::dt()
+{
+  // TODO define, document and implement a user method to set the value of this
+  return 1.0;
+}
+
+
+cvm::real colvarproxy_system::rand_gaussian()
+{
+  // TODO define, document and implement a user method to set the value of this
+  return 0.0;
+}
+
+
 void colvarproxy_system::add_energy(cvm::real /* energy */) {}
 
 
@@ -139,9 +182,31 @@ int colvarproxy_system::get_molid(int &)
 }
 
 
+int colvarproxy_system::get_alch_lambda(cvm::real* lambda)
+{
+  return cvm::error("Error in get_alch_lambda: alchemical lambda dynamics is not supported by this build.",
+    COLVARS_NOT_IMPLEMENTED);
+}
+
+
+int colvarproxy_system::set_alch_lambda(cvm::real* lambda)
+{
+  return cvm::error("Error in set_alch_lambda: alchemical lambda dynamics is not supported by this build.",
+    COLVARS_NOT_IMPLEMENTED);
+}
+
+
+int colvarproxy_system::get_dE_dLambda(cvm::real* force)
+{
+  return cvm::error("Error in get_dE_dLambda: alchemical lambda dynamics is not supported by this build.",
+    COLVARS_NOT_IMPLEMENTED);
+}
+
 
 colvarproxy_atoms::colvarproxy_atoms()
 {
+  atoms_rms_applied_force_ = atoms_max_applied_force_ = 0.0;
+  atoms_max_applied_force_id_ = -1;
   updated_masses_ = updated_charges_ = false;
 }
 
@@ -178,6 +243,18 @@ int colvarproxy_atoms::add_atom_slot(int atom_id)
 }
 
 
+int colvarproxy_atoms::init_atom(int /* atom_number */)
+{
+  return COLVARS_NOT_IMPLEMENTED;
+}
+
+
+int colvarproxy_atoms::check_atom_id(int /* atom_number */)
+{
+  return COLVARS_NOT_IMPLEMENTED;
+}
+
+
 int colvarproxy_atoms::init_atom(cvm::residue_id const & /* residue */,
                                  std::string const     & /* atom_name */,
                                  std::string const     & /* segment_id */)
@@ -232,8 +309,39 @@ int colvarproxy_atoms::load_coords(char const * /* filename */,
 }
 
 
+void colvarproxy_atoms::compute_rms_atoms_applied_force()
+{
+  atoms_rms_applied_force_ =
+    compute_norm2_stats<cvm::rvector, 0, false>(atoms_new_colvar_forces);
+}
 
-colvarproxy_atom_groups::colvarproxy_atom_groups() {}
+
+void colvarproxy_atoms::compute_max_atoms_applied_force()
+{
+  int minmax_index = -1;
+  size_t const n_atoms_ids = atoms_ids.size();
+  if ((n_atoms_ids > 0) && (n_atoms_ids == atoms_new_colvar_forces.size())) {
+    atoms_max_applied_force_ =
+      compute_norm2_stats<cvm::rvector, 1, true>(atoms_new_colvar_forces,
+                                                 &minmax_index);
+    if (minmax_index >= 0) {
+      atoms_max_applied_force_id_ = atoms_ids[minmax_index];
+    } else {
+      atoms_max_applied_force_id_ = -1;
+    }
+  } else {
+    atoms_max_applied_force_ =
+      compute_norm2_stats<cvm::rvector, 1, false>(atoms_new_colvar_forces);
+    atoms_max_applied_force_id_ = -1;
+  }
+}
+
+
+
+colvarproxy_atom_groups::colvarproxy_atom_groups()
+{
+  atom_groups_rms_applied_force_ = atom_groups_max_applied_force_ = 0.0;
+}
 
 
 colvarproxy_atom_groups::~colvarproxy_atom_groups()
@@ -296,6 +404,20 @@ void colvarproxy_atom_groups::clear_atom_group(int index)
 }
 
 
+void colvarproxy_atom_groups::compute_rms_atom_groups_applied_force()
+{
+  atom_groups_rms_applied_force_ =
+    compute_norm2_stats<cvm::rvector, 0, false>(atom_groups_new_colvar_forces);
+}
+
+
+void colvarproxy_atom_groups::compute_max_atom_groups_applied_force()
+{
+  atom_groups_max_applied_force_ =
+    compute_norm2_stats<cvm::rvector, 1, false>(atom_groups_new_colvar_forces);
+}
+
+
 
 colvarproxy_smp::colvarproxy_smp()
 {
@@ -464,28 +586,14 @@ int colvarproxy_smp::smp_unlock()
 colvarproxy_script::colvarproxy_script()
 {
   script = NULL;
+  force_script_defined = false;
+  have_scripts = false;
 }
 
 
 colvarproxy_script::~colvarproxy_script() {}
 
 
-char const *colvarproxy_script::script_obj_to_str(unsigned char *obj)
-{
-  cvm::error("Error: trying to print a script object without a scripting "
-             "language interface.\n", BUG_ERROR);
-  return reinterpret_cast<char *>(obj);
-}
-
-
-std::vector<std::string> colvarproxy_script::script_obj_to_str_vector(unsigned char * /* obj */)
-{
-  cvm::error("Error: trying to print a script object without a scripting "
-             "language interface.\n", BUG_ERROR);
-  return std::vector<std::string>();
-}
-
-
 int colvarproxy_script::run_force_callback()
 {
   return COLVARS_NOT_IMPLEMENTED;
@@ -512,6 +620,7 @@ int colvarproxy_script::run_colvar_gradient_callback(std::string const & /* name
 colvarproxy_io::colvarproxy_io()
 {
   input_buffer_ = NULL;
+  restart_frequency_engine = 0;
 }
 
 
@@ -660,6 +769,23 @@ int colvarproxy::update_output()
 }
 
 
+int colvarproxy::end_of_step()
+{
+  // Disable flags that Colvars doesn't need any more
+  updated_masses_ = updated_charges_ = false;
+
+  // Compute force statistics
+  compute_rms_atoms_applied_force();
+  compute_max_atoms_applied_force();
+  compute_rms_atom_groups_applied_force();
+  compute_max_atom_groups_applied_force();
+  compute_rms_volmaps_applied_force();
+  compute_max_volmaps_applied_force();
+
+  return COLVARS_OK;
+}
+
+
 int colvarproxy::post_run()
 {
   int error_code = COLVARS_OK;
@@ -672,6 +798,19 @@ int colvarproxy::post_run()
 }
 
 
+void colvarproxy::log(std::string const &message)
+{
+  fprintf(stdout, "colvars: %s", message.c_str());
+}
+
+
+void colvarproxy::error(std::string const &message)
+{
+  // TODO handle errors?
+  colvarproxy::log(message);
+}
+
+
 void colvarproxy::add_error_msg(std::string const &message)
 {
   std::istringstream is(message);

lib/colvars/colvarproxy.h

@@ -59,7 +59,7 @@ public:
   std::string units;
 
   /// \brief Request to set the units used internally by Colvars
-  virtual int set_unit_system(std::string const &units, bool check_only) = 0;
+  virtual int set_unit_system(std::string const &units, bool check_only);
 
   /// \brief Value of 1 Angstrom in the internal (front-end) Colvars unit for atomic coordinates
   /// * defaults to 0. in the base class; derived proxy classes must set it
@@ -68,7 +68,7 @@ public:
   cvm::real angstrom_value;
 
   /// \brief Value of 1 Angstrom in the backend's unit for atomic coordinates
-  virtual cvm::real backend_angstrom_value() = 0;
+  virtual cvm::real backend_angstrom_value();
 
   /// \brief Value of 1 kcal/mol in the internal Colvars unit for energy
   cvm::real kcal_mol_value;
@@ -79,6 +79,12 @@ public:
     return l * angstrom_value;
   }
 
+  /// \brief Convert a length from internal to Angstrom
+  inline cvm::real internal_to_angstrom(cvm::real l) const
+  {
+    return l / angstrom_value;
+  }
+
   // /// \brief Convert a length from back-end unit to internal
   // inline cvm::real back_end_to_internal_unit(cvm::real l) {
   //   if (angstrom_value == 0.) {
@@ -88,19 +94,19 @@ public:
   // }
 
   /// \brief Boltzmann constant in internal Colvars units
-  virtual cvm::real boltzmann() = 0;
+  virtual cvm::real boltzmann();
 
   /// \brief Target temperature of the simulation (K units)
-  virtual cvm::real temperature() = 0;
+  virtual cvm::real temperature();
 
   /// \brief Time step of the simulation (fs)
-  virtual cvm::real dt() = 0;
+  virtual cvm::real dt();
 
   /// \brief Pseudo-random number with Gaussian distribution
-  virtual cvm::real rand_gaussian(void) = 0;
+  virtual cvm::real rand_gaussian(void);
 
   /// Pass restraint energy value for current timestep to MD engine
-  virtual void add_energy(cvm::real energy) = 0;
+  virtual void add_energy(cvm::real energy);
 
   /// \brief Get the PBC-aware distance vector between two positions
   virtual cvm::rvector position_distance(cvm::atom_pos const &pos1,
@@ -126,6 +132,15 @@ public:
   /// \param molid Set this argument equal to the current VMD molid
   virtual int get_molid(int &molid);
 
+  /// Get value of alchemical lambda parameter from back-end (if available)
+  virtual int get_alch_lambda(cvm::real* lambda);
+
+  /// Set value of alchemical lambda parameter in back-end (if available)
+  virtual int set_alch_lambda(cvm::real* lambda);
+
+  /// Get energy derivative with respect to lambda (if available)
+  virtual int get_dE_dLambda(cvm::real* force);
+
 protected:
 
   /// Whether the total forces have been requested
@@ -167,11 +182,11 @@ public:
 
   /// Prepare this atom for collective variables calculation, selecting it by
   /// numeric index (1-based)
-  virtual int init_atom(int atom_number) = 0;
+  virtual int init_atom(int atom_number);
 
   /// Check that this atom number is valid, but do not initialize the
   /// corresponding atom yet
-  virtual int check_atom_id(int atom_number) = 0;
+  virtual int check_atom_id(int atom_number);
 
   /// Select this atom for collective variables calculation, using name and
   /// residue number.  Not all programs support this: leave this function as
@@ -262,11 +277,16 @@ public:
     return cvm::rvector(0.0);
   }
 
-  inline std::vector<int> *modify_atom_ids()
+  inline std::vector<int> const *get_atom_ids() const
   {
     return &atoms_ids;
   }
 
+  inline std::vector<cvm::real> const *get_atom_masses() const
+  {
+    return &atoms_masses;
+  }
+
   inline std::vector<cvm::real> *modify_atom_masses()
   {
     // assume that we are requesting masses to change them
@@ -274,6 +294,11 @@ public:
     return &atoms_masses;
   }
 
+  inline std::vector<cvm::real> const *get_atom_charges()
+  {
+    return &atoms_charges;
+  }
+
   inline std::vector<cvm::real> *modify_atom_charges()
   {
     // assume that we are requesting charges to change them
@@ -281,21 +306,60 @@ public:
     return &atoms_charges;
   }
 
+  inline std::vector<cvm::rvector> const *get_atom_positions() const
+  {
+    return &atoms_positions;
+  }
+
   inline std::vector<cvm::rvector> *modify_atom_positions()
   {
     return &atoms_positions;
   }
 
+  inline std::vector<cvm::rvector> const *get_atom_total_forces() const
+  {
+    return &atoms_total_forces;
+  }
+
   inline std::vector<cvm::rvector> *modify_atom_total_forces()
   {
     return &atoms_total_forces;
   }
 
-  inline std::vector<cvm::rvector> *modify_atom_new_colvar_forces()
+  inline std::vector<cvm::rvector> const *get_atom_applied_forces() const
   {
     return &atoms_new_colvar_forces;
   }
 
+  inline std::vector<cvm::rvector> *modify_atom_applied_forces()
+  {
+    return &atoms_new_colvar_forces;
+  }
+
+  /// Compute the root-mean-square of the applied forces
+  void compute_rms_atoms_applied_force();
+
+  /// Compute the maximum norm among all applied forces
+  void compute_max_atoms_applied_force();
+
+  /// Get the root-mean-square of the applied forces
+  inline cvm::real rms_atoms_applied_force() const
+  {
+    return atoms_rms_applied_force_;
+  }
+
+  /// Get the maximum norm among all applied forces
+  inline cvm::real max_atoms_applied_force() const
+  {
+    return atoms_max_applied_force_;
+  }
+
+  /// Get the atom ID with the largest applied force
+  inline int max_atoms_applied_force_id() const
+  {
+    return atoms_max_applied_force_id_;
+  }
+
   /// Record whether masses have been updated
   inline bool updated_masses() const
   {
@@ -326,6 +390,15 @@ protected:
   /// \brief Forces applied from colvars, to be communicated to the MD integrator
   std::vector<cvm::rvector> atoms_new_colvar_forces;
 
+  /// Root-mean-square of the applied forces
+  cvm::real atoms_rms_applied_force_;
+
+  /// Maximum norm among all applied forces
+  cvm::real atoms_max_applied_force_;
+
+  /// ID of the atom with the maximum norm among all applied forces
+  int atoms_max_applied_force_id_;
+
   /// Whether the masses and charges have been updated from the host code
   bool updated_masses_, updated_charges_;
 
@@ -404,6 +477,56 @@ public:
     return cvm::rvector(0.0);
   }
 
+  inline std::vector<int> const *get_atom_group_ids() const
+  {
+    return &atom_groups_ids;
+  }
+
+  inline std::vector<cvm::real> *modify_atom_group_masses()
+  {
+    // TODO updated_masses
+    return &atom_groups_masses;
+  }
+
+  inline std::vector<cvm::real> *modify_atom_group_charges()
+  {
+    // TODO updated masses
+    return &atom_groups_charges;
+  }
+
+  inline std::vector<cvm::rvector> *modify_atom_group_positions()
+  {
+    return &atom_groups_coms;
+  }
+
+  inline std::vector<cvm::rvector> *modify_atom_group_total_forces()
+  {
+    return &atom_groups_total_forces;
+  }
+
+  inline std::vector<cvm::rvector> *modify_atom_group_applied_forces()
+  {
+    return &atom_groups_new_colvar_forces;
+  }
+
+  /// Compute the root-mean-square of the applied forces
+  void compute_rms_atom_groups_applied_force();
+
+  /// Compute the maximum norm among all applied forces
+  void compute_max_atom_groups_applied_force();
+
+  /// Get the root-mean-square of the applied forces
+  inline cvm::real rms_atom_groups_applied_force() const
+  {
+    return atom_groups_rms_applied_force_;
+  }
+
+  /// Get the maximum norm among all applied forces
+  inline cvm::real max_atom_groups_applied_force() const
+  {
+    return atom_groups_max_applied_force_;
+  }
+
 protected:
 
   /// \brief Array of 0-based integers used to uniquely associate atom groups
@@ -422,6 +545,12 @@ protected:
   /// \brief Forces applied from colvars, to be communicated to the MD integrator
   std::vector<cvm::rvector> atom_groups_new_colvar_forces;
 
+  /// Root-mean-square of the applied group forces
+  cvm::real atom_groups_rms_applied_force_;
+
+  /// Maximum norm among all applied group forces
+  cvm::real atom_groups_max_applied_force_;
+
   /// Used by all init_atom_group() functions: create a slot for an atom group not requested yet
   int add_atom_group_slot(int atom_group_id);
 };
@@ -519,12 +648,6 @@ public:
   /// Destructor
   virtual ~colvarproxy_script();
 
-  /// Convert a script object (Tcl or Python call argument) to a C string
-  virtual char const *script_obj_to_str(unsigned char *obj);
-
-  /// Convert a script object (Tcl or Python call argument) to a vector of strings
-  virtual std::vector<std::string> script_obj_to_str_vector(unsigned char *obj);
-
   /// Pointer to the scripting interface object
   /// (does not need to be allocated in a new interface)
   colvarscript *script;
@@ -706,11 +829,14 @@ public:
   /// \brief Update data based from the results of a module update (e.g. send forces)
   virtual int update_output();
 
+  /// Carry out operations needed before next step is run
+  int end_of_step();
+
   /// Print a message to the main log
-  virtual void log(std::string const &message) = 0;
+  virtual void log(std::string const &message);
 
   /// Print a message to the main log and/or let the host code know about it
-  virtual void error(std::string const &message) = 0;
+  virtual void error(std::string const &message);
 
   /// Record error message (used by VMD to collect them after a script call)
   void add_error_msg(std::string const &message);

lib/colvars/colvarproxy_tcl.cpp

@@ -11,6 +11,8 @@
 
 #if defined(NAMD_TCL) || defined(VMDTCL)
 #define COLVARS_TCL
+#endif
+#ifdef COLVARS_TCL
 #include <tcl.h>
 #endif
 

lib/colvars/colvarproxy_tcl.h

@@ -27,7 +27,7 @@ public:
   /// Is Tcl available? (trigger initialization if needed)
   int tcl_available();
 
-  /// Tcl implementation of script_obj_to_str()
+  /// Get a string representation of the Tcl object pointed to by obj
   char const *tcl_get_str(void *obj);
 
   /// Tcl implementation of run_force_callback()
@@ -51,6 +51,12 @@ public:
     return tcl_interp_;
   }
 
+  /// Set the pointer to the Tcl interpreter
+  inline void set_tcl_interp(void *interp)
+  {
+    tcl_interp_ = interp;
+  }
+
 protected:
 
   /// Pointer to Tcl interpreter object

lib/colvars/colvarproxy_volmaps.cpp

@@ -9,9 +9,13 @@
 
 #include "colvarmodule.h"
 #include "colvarproxy_volmaps.h"
+#include "colvarmodule_utils.h"
 
 
-colvarproxy_volmaps::colvarproxy_volmaps() {}
+colvarproxy_volmaps::colvarproxy_volmaps()
+{
+  volmaps_rms_applied_force_ = volmaps_max_applied_force_ = 0.0;
+}
 
 
 colvarproxy_volmaps::~colvarproxy_volmaps() {}
@@ -46,25 +50,41 @@ int colvarproxy_volmaps::add_volmap_slot(int volmap_id)
 }
 
 
-int colvarproxy_volmaps::init_volmap(int volmap_id)
+int colvarproxy_volmaps::check_volmap_by_id(int /* volmap_id */)
 {
-  return cvm::error("Error: access to volumetric maps is unavailable "
+  return cvm::error("Error: selecting volumetric maps is not available.\n",
-                    "in this build.\n",
                     COLVARS_NOT_IMPLEMENTED);
 }
 
 
-int colvarproxy_volmaps::init_volmap(const char *volmap_name)
+int colvarproxy_volmaps::check_volmap_by_name(const char * /* volmap_name */)
 {
-  return cvm::error("Error: access to volumetric maps is unavailable "
+  return cvm::error("Error: selecting volumetric maps by name is not "
-                    "in this build.\n",
+                    "available.\n", COLVARS_NOT_IMPLEMENTED);
-                    COLVARS_NOT_IMPLEMENTED);
 }
 
 
-int colvarproxy_volmaps::init_volmap(const std::string &volmap_name)
+int colvarproxy_volmaps::init_volmap_by_name(char const *volmap_name)
 {
-  return init_volmap(volmap_name.c_str());
+  return -1;
+}
+
+
+int colvarproxy_volmaps::init_volmap_by_id(int volmap_id)
+{
+  return -1;
+}
+
+
+int colvarproxy_volmaps::init_volmap_by_name(std::string const &volmap_name)
+{
+  return init_volmap_by_name(volmap_name.c_str());
+}
+
+
+int colvarproxy_volmaps::check_volmap_by_name(std::string const &volmap_name)
+{
+  return check_volmap_by_name(volmap_name.c_str());
 }
 
 
@@ -79,3 +99,36 @@ void colvarproxy_volmaps::clear_volmap(int index)
     volmaps_ncopies[index] -= 1;
   }
 }
+
+
+int colvarproxy_volmaps::get_volmap_id_from_name(char const *volmap_name)
+{
+  // Raise error
+  colvarproxy_volmaps::check_volmap_by_name(volmap_name);
+  return -1;
+}
+
+
+int colvarproxy_volmaps::compute_volmap(int /* flags */,
+                                        int /* volmap_id */,
+                                        cvm::atom_iter /* atom_begin */,
+                                        cvm::atom_iter /* atom_end */,
+                                        cvm::real * /* value */,
+                                        cvm::real * /* atom_field */)
+{
+  return COLVARS_NOT_IMPLEMENTED;
+}
+
+
+void colvarproxy_volmaps::compute_rms_volmaps_applied_force()
+{
+  volmaps_rms_applied_force_ =
+    compute_norm2_stats<cvm::real, 0, false>(volmaps_new_colvar_forces);
+}
+
+
+void colvarproxy_volmaps::compute_max_volmaps_applied_force()
+{
+  volmaps_max_applied_force_ =
+    compute_norm2_stats<cvm::real, 1, false>(volmaps_new_colvar_forces);
+}

lib/colvars/colvarproxy_volmaps.h

@@ -25,17 +25,37 @@ public:
   int add_volmap_slot(int volmap_id);
 
   /// Request and prepare this volumetric map for use by Colvars
-  virtual int init_volmap(int volmap_id);
+  /// \param volmap_id Numeric ID used by the MD engine
+  /// \returns Index of the map in the colvarproxy arrays
+  virtual int init_volmap_by_id(int volmap_id);
 
   /// Request and prepare this volumetric map for use by Colvars
-  virtual int init_volmap(char const *volmap_name);
+  /// \param volmap_name Name used by the MD engine
+  /// \returns Index of the map in the colvarproxy arrays
+  virtual int init_volmap_by_name(char const *volmap_name);
+
+  /// Check that the given volmap ID is valid (return COLVARS_OK if it is)
+  /// \param volmap_id Numeric ID used by the MD engine
+  /// \returns Error code
+  virtual int check_volmap_by_id(int volmap_id);
+
+  /// Check that the given volmap name is valid (return COLVARS_OK if it is)
+  /// \param volmap_name Name used by the MD engine
+  /// \returns Error code
+  virtual int check_volmap_by_name(char const *volmap_name);
 
   /// Request and prepare this volumetric map for use by Colvars
-  int init_volmap(std::string const &volmap_name);
+  int init_volmap_by_name(std::string const &volmap_name);
+
+  /// Check that the given volmap name is valid (return COLVARS_OK if it is)
+  int check_volmap_by_name(std::string const &volmap_name);
 
   /// \brief Used by the CVC destructors
   virtual void clear_volmap(int index);
 
+  /// Get the numeric ID of the given volumetric map (for the MD program)
+  virtual int get_volmap_id_from_name(char const *volmap_name);
+
   /// Get the numeric ID of the given volumetric map (for the MD program)
   inline int get_volmap_id(int index) const
   {
@@ -54,6 +74,32 @@ public:
     volmaps_new_colvar_forces[index] += new_force;
   }
 
+  /// Re-weigh an atomic field (e.g. a colvar) by the value of a volumetric map
+  /// \param flags Combination of flags
+  /// \param volmap_id Numeric index of the map (no need to request it)
+  /// \param atom_begin Iterator pointing to first atom
+  /// \param atom_end Iterator pointing past the last atom
+  /// \param value Pointer to location of total to increment
+  /// \param atom_field Array of atomic field values (if NULL, ones are used)
+  virtual int compute_volmap(int flags,
+                             int volmap_id,
+                             cvm::atom_iter atom_begin,
+                             cvm::atom_iter atom_end,
+                             cvm::real *value,
+                             cvm::real *atom_field);
+
+  /// Flags controlling what computation is done on the map
+  enum {
+    volmap_flag_null = 0,
+    volmap_flag_gradients = 1,
+    volmap_flag_use_atom_field = (1<<8)
+  };
+
+  /// Compute the root-mean-square of the applied forces
+  void compute_rms_volmaps_applied_force();
+
+  /// Compute the maximum norm among all applied forces
+  void compute_max_volmaps_applied_force();
 
 protected:
 
@@ -70,6 +116,12 @@ protected:
   /// \brief Forces applied from colvars, to be communicated to the MD
   /// integrator
   std::vector<cvm::real>    volmaps_new_colvar_forces;
+
+  /// Root-mean-square of the the applied forces
+  cvm::real volmaps_rms_applied_force_;
+
+  /// Maximum norm among all applied forces
+  cvm::real volmaps_max_applied_force_;
 };
 
 

lib/colvars/colvars_version.h

@@ -1,3 +1,3 @@
 #ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2020-09-17"
+#define COLVARS_VERSION "2021-08-03"
 #endif

lib/colvars/colvarscript.cpp

@@ -13,6 +13,8 @@
 
 #if defined(NAMD_TCL) || defined(VMDTCL)
 #define COLVARS_TCL
+#endif
+#ifdef COLVARS_TCL
 #include <tcl.h>
 #endif
 
@@ -23,6 +25,19 @@
 
 
 
+#ifdef COLVARS_TCL
+/// Run the script API via Tcl command-line interface
+/// \param clientData Not used
+/// \param my_interp Pointer to Tcl_Interp object (read from Colvars if NULL)
+/// \param objc Number of Tcl command parameters
+/// \param objv Array of command parameters
+/// \return Result of the script command
+extern "C" int tcl_run_colvarscript_command(ClientData clientData,
+                                            Tcl_Interp *interp_in,
+                                            int objc, Tcl_Obj *const objv[]);
+#endif
+
+
 colvarscript::colvarscript(colvarproxy *p)
  : proxy_(p),
    colvars(p->colvars),
@@ -57,9 +72,11 @@ int colvarscript::init_commands()
   }
 
   cmd_help.resize(colvarscript::cv_n_commands);
+  cmd_rethelp.resize(colvarscript::cv_n_commands);
   cmd_n_args_min.resize(colvarscript::cv_n_commands);
   cmd_n_args_max.resize(colvarscript::cv_n_commands);
   cmd_arghelp.resize(colvarscript::cv_n_commands);
+  cmd_full_help.resize(colvarscript::cv_n_commands);
   cmd_fns.resize(colvarscript::cv_n_commands);
 
   if (cmd_names) {
@@ -95,9 +112,24 @@ int colvarscript::init_command(colvarscript::command const &comm,
 {
   cmd_str_map[std::string(name)] = comm;
   cmd_names[comm] = name;
-  cmd_help[comm] = help;
+
+  // Initialize short help string and return-value help string (if present)
+  {
+    std::string const help_str(help);
+    std::istringstream is(help_str);
+    std::string line;
+    std::getline(is, line);
+    cmd_help[comm] = line;
+    cmd_rethelp[comm] = "";
+    while (std::getline(is, line)) {
+      cmd_rethelp[comm] += line + "\n";
+    }
+  }
+
+  // Initialize arguments' help strings
   cmd_n_args_min[comm] = n_args_min;
   cmd_n_args_max[comm] = n_args_max;
+  {
     std::string const arghelp_str(arghelp);
     std::istringstream is(arghelp_str);
     std::string line;
@@ -108,11 +140,32 @@ int colvarscript::init_command(colvarscript::command const &comm,
       }
       cmd_arghelp[comm].push_back(line);
     }
+  }
+
+  cmd_full_help[comm] = cmd_help[comm]+"\n";
+  if (cmd_n_args_min[comm] > 0) {
+    cmd_full_help[comm] += "\nParameters\n";
+    cmd_full_help[comm] += "----------\n\n";
+    size_t i;
+    for (i = 0; i < cmd_n_args_min[comm]; i++) {
+      cmd_full_help[comm] += cmd_arghelp[comm][i]+"\n";
+    }
+    for (i = cmd_n_args_min[comm]; i < cmd_n_args_max[comm]; i++) {
+      cmd_full_help[comm] += cmd_arghelp[comm][i]+" (optional)\n";
+    }
+  }
+  if (cmd_rethelp[comm].size() > 0) {
+    cmd_full_help[comm] += "\nReturns\n";
+    cmd_full_help[comm] += "-------\n\n";
+    cmd_full_help[comm] += cmd_rethelp[comm]+"\n";
+  }
+
   cmd_fns[comm] = fn;
   if (cvm::debug()) {
     cvm::log("Defined command \""+std::string(name)+"\", with help string:\n");
-    cvm::log(get_command_help(name));
+    cvm::log(get_command_full_help(name));
   }
+
   return COLVARS_OK;
 }
 
@@ -133,27 +186,76 @@ std::string colvarscript::get_cmd_prefix(colvarscript::Object_type t)
 }
 
 
-std::string colvarscript::get_command_help(char const *cmd)
+
+char const *colvarscript::get_command_help(char const *cmd)
 {
   if (cmd_str_map.count(cmd) > 0) {
     colvarscript::command const c = cmd_str_map[std::string(cmd)];
-    std::string new_result(cmd_help[c]+"\n");
+    return cmd_help[c].c_str();
-    if (cmd_n_args_max[c] == 0) return new_result;
-    new_result += "\nParameters\n";
-    new_result += "----------\n\n";
-    size_t i;
-    for (i = 0; i < cmd_n_args_min[c]; i++) {
-      new_result += cmd_arghelp[c][i]+"\n";
   }
-    for (i = cmd_n_args_min[c]; i < cmd_n_args_max[c]; i++) {
-      new_result += cmd_arghelp[c][i]+" (optional)\n";
-    }
-    return new_result;
-  }
-
   cvm::error("Error: command "+std::string(cmd)+
              " is not implemented.\n", INPUT_ERROR);
-  return std::string("");
+  return NULL;
+}
+
+
+char const *colvarscript::get_command_rethelp(char const *cmd)
+{
+  if (cmd_str_map.count(cmd) > 0) {
+    colvarscript::command const c = cmd_str_map[std::string(cmd)];
+    return cmd_rethelp[c].c_str();
+  }
+  cvm::error("Error: command "+std::string(cmd)+
+             " is not implemented.\n", INPUT_ERROR);
+  return NULL;
+}
+
+
+char const *colvarscript::get_command_arghelp(char const *cmd, int i)
+{
+  if (cmd_str_map.count(cmd) > 0) {
+    colvarscript::command const c = cmd_str_map[std::string(cmd)];
+    return cmd_arghelp[c][i].c_str();
+  }
+  cvm::error("Error: command "+std::string(cmd)+
+             " is not implemented.\n", INPUT_ERROR);
+  return NULL;
+}
+
+
+int colvarscript::get_command_n_args_min(char const *cmd)
+{
+  if (cmd_str_map.count(cmd) > 0) {
+    colvarscript::command const c = cmd_str_map[std::string(cmd)];
+    return cmd_n_args_min[c];
+  }
+  cvm::error("Error: command "+std::string(cmd)+
+             " is not implemented.\n", INPUT_ERROR);
+  return -1;
+}
+
+
+int colvarscript::get_command_n_args_max(char const *cmd)
+{
+  if (cmd_str_map.count(cmd) > 0) {
+    colvarscript::command const c = cmd_str_map[std::string(cmd)];
+    return cmd_n_args_max[c];
+  }
+  cvm::error("Error: command "+std::string(cmd)+
+             " is not implemented.\n", INPUT_ERROR);
+  return -1;
+}
+
+
+char const *colvarscript::get_command_full_help(char const *cmd)
+{
+  if (cmd_str_map.count(cmd) > 0) {
+    colvarscript::command const c = cmd_str_map[std::string(cmd)];
+    return cmd_full_help[c].c_str();
+  }
+  cvm::error("Error: command "+std::string(cmd)+
+             " is not implemented.\n", INPUT_ERROR);
+  return NULL;
 }
 
 
@@ -234,7 +336,7 @@ std::string colvarscript::get_command_cmdline_help(colvarscript::Object_type t,
   if (cmd_str_map.count(cmdkey) > 0) {
     command const c = cmd_str_map[cmdkey];
     return get_command_cmdline_syntax(t, c)+"\n\n"+
-      get_command_help(cmd_names[c]);
+      get_command_full_help(cmd_names[c]);
   }
   cvm::error("Error: could not find scripting command \""+cmd+"\".",
              INPUT_ERROR);
@@ -244,7 +346,7 @@ std::string colvarscript::get_command_cmdline_help(colvarscript::Object_type t,
 
 int colvarscript::run(int objc, unsigned char *const objv[])
 {
-  result.clear();
+  clear_str_result();
 
   if (cvm::debug()) {
     cvm::log("Called script run with " + cvm::to_str(objc) + " args:");
@@ -346,6 +448,60 @@ int colvarscript::run(int objc, unsigned char *const objv[])
 }
 
 
+char *colvarscript::obj_to_str(unsigned char *obj)
+{
+  char *strobj = reinterpret_cast<char *>(obj);
+  if (cvm::debug()) {
+    cvm::log("Using simple-cast script::obj_to_str(): result = \"" +
+             (strobj ? std::string(strobj) : std::string("(null)")) + "\"");
+  }
+  return strobj;
+}
+
+
+std::vector<std::string> colvarscript::obj_to_str_vector(unsigned char *obj)
+{
+  if (cvm::debug()) {
+    cvm::log("Using simple-cast colvarscript::obj_to_str_vector().\n");
+  }
+
+  std::vector<std::string> new_result;
+  std::string const str(reinterpret_cast<char *>(obj));
+
+  // TODO get rid of this once colvarscript can handle both fix_modify and Tcl?
+  // LAMMPS has a nicer function in the utils class
+
+  for (size_t i = 0; i < str.length(); i++) {
+    char const c = str[i];
+    if (c == '\"') {
+      i++;
+      if (i >= str.length()) {
+        cvm::error("Error: could not split the following string:\n"+
+                   str+"\n", INPUT_ERROR);
+        break;
+      }
+      new_result.push_back(std::string(""));
+      while (str[i] != '\"') {
+        new_result.back().append(1, str[i]);
+        if (i >= str.length()) {
+          cvm::error("Error: could not split the following string:\n"+
+                     str+"\n", INPUT_ERROR);
+          break;
+        } else {
+          i++;
+        }
+      }
+    }
+  }
+
+  if (cvm::debug()) {
+    cvm::log("result = "+cvm::to_str(new_result)+".\n");
+  }
+
+  return new_result;
+}
+
+
 int colvarscript::proc_features(colvardeps *obj,
                                 int objc, unsigned char *const objv[]) {
 
@@ -428,9 +584,9 @@ int colvarscript::set_result_str(std::string const &s)
 {
   if (cvm::get_error() != COLVARS_OK) {
     // Avoid overwriting the error message
-    result += s;
+    modify_str_result() += s;
   } else {
-    result = s;
+    modify_str_result() = s;
   }
   return COLVARS_OK;
 }
@@ -438,17 +594,17 @@ int colvarscript::set_result_str(std::string const &s)
 
 void colvarscript::add_error_msg(std::string const &s)
 {
-  result += s;
+  modify_str_result() += s;
   // Ensure terminating newlines
   if (s[s.size()-1] != '\n') {
-    result += "\n";
+    modify_str_result() += "\n";
   }
 }
 
 
 int colvarscript::clear_str_result()
 {
-  result.clear();
+  modify_str_result().clear();
   return COLVARS_OK;
 }
 
@@ -487,8 +643,24 @@ const char * get_colvarscript_result()
 int tcl_colvars_vmd_init(Tcl_Interp *interp, int molid);
 #endif
 
-extern "C"
+#if !defined(VMDTCL) && !defined(NAMD_TCL)
-int tcl_run_colvarscript_command(ClientData /* clientData */,
+extern "C" {
+  int Colvars_Init(Tcl_Interp *interp) {
+    colvarproxy *proxy = new colvarproxy();
+    colvarmodule *colvars = new colvarmodule(proxy);
+    proxy->set_tcl_interp(reinterpret_cast<void *>(interp));
+    proxy->colvars = colvars;
+    proxy->script = new colvarscript(proxy);
+    Tcl_CreateObjCommand(interp, "cv", tcl_run_colvarscript_command,
+                         (ClientData *) NULL, (Tcl_CmdDeleteProc *) NULL);
+    Tcl_EvalEx(interp, "package provide colvars", -1, 0);
+    return TCL_OK;
+  }
+}
+#endif
+
+
+extern "C" int tcl_run_colvarscript_command(ClientData /* clientData */,
                                             Tcl_Interp *my_interp,
                                             int objc, Tcl_Obj *const objv[])
 {
@@ -496,23 +668,44 @@ int tcl_run_colvarscript_command(ClientData /* clientData */,
 
   if (!colvars) {
 #if defined(VMDTCL)
+
+    if (objc == 2) {
+      if (!strcmp(Tcl_GetString(objv[1]), "molid")) {
+        // return invalid molid
+        Tcl_SetResult(my_interp, (char *) "-1", TCL_STATIC);
+      }
+      if (!strcmp(Tcl_GetString(objv[1]), "delete") ||
+          !strcmp(Tcl_GetString(objv[1]), "reset")) {
+        // nothing to delete or reset
+        Tcl_SetResult(my_interp, NULL, TCL_STATIC);
+      }
+      if (!strcmp(Tcl_GetString(objv[1]), "help")) {
+        // print message
+        Tcl_SetResult(my_interp,
+                      (char *) "First, setup the Colvars module with: "
+                      "cv molid <id>|top", TCL_STATIC);
+      }
+      return TCL_OK;
+    }
+
     if (objc >= 3) {
       // require a molid to create the module
       if (!strcmp(Tcl_GetString(objv[1]), "molid")) {
-        int molid = -1;
+        int molid = -(1<<16); // This value is used to indicate "top"
         if (strcmp(Tcl_GetString(objv[2]), "top")) {
           // If this is not "top", get the integer value
           Tcl_GetIntFromObj(my_interp, objv[2], &molid);
         }
         return tcl_colvars_vmd_init(my_interp, molid);
       } else {
-        // TODO allow calling cv help after this
+        Tcl_SetResult(my_interp, (char *) "Syntax error.  First, setup the Colvars module with cv molid <id>|top", TCL_STATIC);
-        Tcl_SetResult(my_interp, (char *) "Syntax error.", TCL_STATIC);
         return TCL_ERROR;
       }
     }
+
     Tcl_SetResult(my_interp, (char *) "First, setup the Colvars module with: "
-                  "cv molid <molecule id>", TCL_STATIC);
+                  "cv molid <id>|top", TCL_STATIC);
+
 #else
     Tcl_SetResult(my_interp,
                   const_cast<char *>("Error: Colvars module not yet initialized"),
@@ -534,10 +727,19 @@ int tcl_run_colvarscript_command(ClientData /* clientData */,
 
   cvm::clear_error();
 
-  int retval = script->run(objc,
+  unsigned char * arg_pointers_[100];
-                           reinterpret_cast<unsigned char * const *>(objv));
+  if (objc > 100) {
+    std::string const errstr = "Too many positional arguments ("+
+      cvm::to_str(objc)+") passed to the \"cv\" command.\n";
+    Tcl_SetResult(interp, const_cast<char *>(errstr.c_str()), TCL_VOLATILE);
+    return TCL_ERROR;
+  }
+  for (int i = 0; i < objc; i++) {
+    arg_pointers_[i] = reinterpret_cast<unsigned char *>(const_cast<char *>(proxy->tcl_get_str(objv[i])));
+  }
+  int retval = script->run(objc, arg_pointers_);
 
-  std::string result = proxy->get_error_msgs() + script->result;
+  std::string result = proxy->get_error_msgs() + script->str_result();
 
   Tcl_SetResult(interp, const_cast<char *>(result.c_str()),
                 TCL_VOLATILE);
@@ -558,3 +760,162 @@ int tcl_run_colvarscript_command(ClientData /* clientData */,
 }
 
 #endif // #if defined(COLVARS_TCL)
+
+
+
+
+int colvarscript::set_result_text_from_str(std::string const &x_str,
+                                           unsigned char *obj) {
+  if (obj) {
+    strcpy(reinterpret_cast<char *>(obj), x_str.c_str());
+  } else {
+    set_result_str(x_str);
+  }
+  return COLVARS_OK;
+}
+
+// Template to convert everything to string and use the above
+
+template <typename T>
+int colvarscript::set_result_text(T const &x, unsigned char *obj) {
+  std::string const x_str = x.to_simple_string();
+  return set_result_text_from_str(x_str, obj);
+}
+
+
+template <typename T>
+int colvarscript::pack_vector_elements_text(std::vector<T> const &x,
+                                            std::string &x_str) {
+  x_str.clear();
+  for (size_t i = 0; i < x.size(); ++i) {
+    if (i > 0) x_str.append(1, ' ');
+    x_str += cvm::to_str(x[i]);
+  }
+  return COLVARS_OK;
+}
+
+
+// Specializations for plain old data types that don't have a stringifier member
+
+template <>
+int colvarscript::set_result_text(int const &x, unsigned char *obj) {
+  std::string const x_str = cvm::to_str(x);
+  return set_result_text_from_str(x_str, obj);
+}
+
+template <>
+int colvarscript::set_result_text(std::vector<int> const &x,
+                                  unsigned char *obj) {
+  std::string x_str("");
+  pack_vector_elements_text<int>(x, x_str);
+  return set_result_text_from_str(x_str, obj);
+}
+
+
+template <>
+int colvarscript::set_result_text(long int const &x, unsigned char *obj) {
+  std::string const x_str = cvm::to_str(x);
+  return set_result_text_from_str(x_str, obj);
+}
+
+template <>
+int colvarscript::set_result_text(std::vector<long int> const &x,
+                                  unsigned char *obj) {
+  std::string x_str("");
+  pack_vector_elements_text<long int>(x, x_str);
+  return set_result_text_from_str(x_str, obj);
+}
+
+
+template <>
+int colvarscript::set_result_text(cvm::real const &x, unsigned char *obj) {
+  std::string const x_str = cvm::to_str(x);
+  return set_result_text_from_str(x_str, obj);
+}
+
+template <>
+int colvarscript::set_result_text(std::vector<cvm::real> const &x,
+                                  unsigned char *obj) {
+  std::string x_str("");
+  pack_vector_elements_text<cvm::real>(x, x_str);
+  return set_result_text_from_str(x_str, obj);
+}
+
+
+// TODO these can be removed after the Tcl backend is ready (otherwise, the
+// default template syntax may break scripts or the Dashboard)
+
+template <>
+int colvarscript::set_result_text(std::vector<cvm::rvector> const &x,
+                                  unsigned char *obj) {
+  std::string x_str("");
+  for (size_t i = 0; i < x.size(); i++) {
+    if (i > 0) x_str.append(1, ' ');
+    x_str += "{ "+x[i].to_simple_string()+" }";
+  }
+  return set_result_text_from_str(x_str, obj);
+}
+
+template <>
+int colvarscript::set_result_text(std::vector<colvarvalue> const &x,
+                                  unsigned char *obj) {
+  std::string x_str("");
+  for (size_t i = 0; i < x.size(); i++) {
+    if (i > 0) x_str.append(1, ' ');
+    x_str += "{ "+x[i].to_simple_string()+" }";
+  }
+  return set_result_text_from_str(x_str, obj);
+}
+
+
+// Member functions to set script results for each typexc
+
+int colvarscript::set_result_int(int const &x, unsigned char *obj) {
+  return set_result_text<int>(x, obj);
+}
+
+int colvarscript::set_result_int_vec(std::vector<int> const &x,
+                                     unsigned char *obj) {
+  return set_result_text< std::vector<int> >(x, obj);
+}
+
+
+int colvarscript::set_result_long_int(long int const &x, unsigned char *obj) {
+  return set_result_text<long int>(x, obj);
+}
+
+int colvarscript::set_result_long_int_vec(std::vector<long int> const &x,
+                                          unsigned char *obj) {
+  return set_result_text< std::vector<long int> >(x, obj);
+}
+
+
+int colvarscript::set_result_real(cvm::real const &x, unsigned char *obj) {
+  return set_result_text<cvm::real>(x, obj);
+}
+
+int colvarscript::set_result_real_vec(std::vector<cvm::real> const &x,
+                                      unsigned char *obj) {
+  return set_result_text< std::vector<cvm::real> >(x, obj);
+}
+
+
+int colvarscript::set_result_rvector(cvm::rvector const &x, unsigned char *obj) {
+  return set_result_text<cvm::rvector>(x, obj);
+}
+
+int colvarscript::set_result_rvector_vec(std::vector<cvm::rvector> const &x,
+                                         unsigned char *obj) {
+  return set_result_text< std::vector<cvm::rvector> >(x, obj);
+}
+
+
+int colvarscript::set_result_colvarvalue(colvarvalue const &x,
+                                         unsigned char *obj) {
+  return set_result_text<colvarvalue>(x, obj);
+}
+
+int colvarscript::set_result_colvarvalue_vec(std::vector<colvarvalue> const &x,
+                                             unsigned char *obj) {
+  return set_result_text< std::vector<colvarvalue> >(x, obj);
+}

lib/colvars/colvarscript.h

@@ -46,9 +46,8 @@ public:
   /// COLVARSCRIPT_ERROR
   int proxy_error;
 
-  /// If an error is returned by one of the methods, it should set this to the
+  /// String representation of the result of a script call
-  /// error message
+  std::string str_result_;
-  std::string result;
 
   /// Run a script command with space-separated positional arguments (objects)
   int run(int objc, unsigned char *const objv[]);
@@ -56,13 +55,13 @@ public:
   /// Get the string result of the current scripting call
   inline std::string const &str_result() const
   {
-    return result;
+    return str_result_;
   }
 
   /// Modify the string result of the current scripting call
   inline std::string &modify_str_result()
   {
-    return result;
+    return str_result_;
   }
 
   /// Set the return value to the given string
@@ -137,21 +136,38 @@ public:
   template<colvarscript::Object_type T>
   int cmd_arg_shift();
 
-  /// Use scripting language to get the string representation of an object
-  inline char const *obj_to_str(unsigned char *const obj)
-  {
-    return (obj == NULL ? NULL : proxy_->script_obj_to_str(obj));
-  }
-
   /// Get names of all commands
   inline char const **get_command_names() const
   {
     return cmd_names;
   }
 
+  /// Get one-line help summary for a command
+  /// \param cmd Name of the command's function (e.g. "cv_units")
+  char const *get_command_help(char const *cmd);
+
+  /// Get description of the return value of a command
+  /// \param cmd Name of the command's function (e.g. "cv_units")
+  char const *get_command_rethelp(char const *cmd);
+
+  /// Get description of the argument of a command (excluding prefix)
+  /// \param cmd Name of the command's function (e.g. "cv_units")
+  /// \param i Index of the argument; 0 is the first argument after the
+  /// prefix, e.g. "value" has an index of 0 in the array of arguments:
+  /// { "cv", "colvar", "xi", "value" }
+  char const *get_command_arghelp(char const *cmd, int i);
+
+  /// Get number of required arguments (excluding prefix)
+  /// \param cmd Name of the command's function (e.g. "cv_units")
+  int get_command_n_args_min(char const *cmd);
+
+  /// Get number of total arguments (excluding prefix)
+  /// \param cmd Name of the command's function (e.g. "cv_units")
+  int get_command_n_args_max(char const *cmd);
+
   /// Get help string for a command (does not specify how it is launched)
   /// \param cmd Name of the command's function (e.g. "cv_units")
-  std::string get_command_help(char const *cmd);
+  char const *get_command_full_help(char const *cmd);
 
   /// Get summary of command line syntax for all commands of a given context
   /// \param t One of use_module, use_colvar or use_bias
@@ -182,6 +198,53 @@ public:
     return this->proxy_;
   }
 
+  // Input functions - get the string reps of script argument objects
+
+  /// Get the string representation of an object (by default, a simple cast)
+  char *obj_to_str(unsigned char *obj);
+
+  /// Get a list of strings from an object (does not work with a simple cast)
+  std::vector<std::string> obj_to_str_vector(unsigned char *obj);
+
+
+  // Output functions - convert internal objects to representations suitable
+  // for use in the scripting language.  At the moment only conversion to C
+  // strings is supported, and obj is assumed to be a char * pointer.
+
+  /// Copy x into obj if not NULL, or into the script object's result otherwise
+  int set_result_int(int const &x, unsigned char *obj = NULL);
+
+  /// Copy x into obj if not NULL, or into the script object's result otherwise
+  int set_result_int_vec(std::vector<int> const &x, unsigned char *obj = NULL);
+
+  /// Copy x into obj if not NULL, or into the script object's result otherwise
+  int set_result_long_int(long int const &x, unsigned char *obj = NULL);
+
+  /// Copy x into obj if not NULL, or into the script object's result otherwise
+  int set_result_long_int_vec(std::vector<long int> const &x,
+                              unsigned char *obj = NULL);
+
+  /// Copy x into obj if not NULL, or into the script object's result otherwise
+  int set_result_real(cvm::real const &x, unsigned char *obj = NULL);
+
+  /// Copy x into obj if not NULL, or into the script object's result otherwise
+  int set_result_real_vec(std::vector<cvm::real> const &x,
+                          unsigned char *obj = NULL);
+
+  /// Copy x into obj if not NULL, or into the script object's result otherwise
+  int set_result_rvector(cvm::rvector const &x, unsigned char *obj = NULL);
+
+  /// Copy x into obj if not NULL, or into the script object's result otherwise
+  int set_result_rvector_vec(std::vector<cvm::rvector> const &x,
+                             unsigned char *obj = NULL);
+
+  /// Copy x into obj if not NULL, or into the script object's result otherwise
+  int set_result_colvarvalue(colvarvalue const &x, unsigned char *obj = NULL);
+
+  /// Copy x into obj if not NULL, or into the script object's result otherwise
+  int set_result_colvarvalue_vec(std::vector<colvarvalue> const &x,
+                                 unsigned char *obj = NULL);
+
 private:
 
   /// Set up all script API functions
@@ -193,14 +256,6 @@ private:
                    int n_args_min, int n_args_max, char const *arghelp,
                    int (*fn)(void *, int, unsigned char * const *));
 
-  /// Execute a script command
-  inline int exec_command(command c,
-                          void *pobj,
-                          int objc, unsigned char * const *objv)
-  {
-    return (*(cmd_fns[c]))(pobj, objc, objv);
-  }
-
 public: // TODO this function will be removed soon
 
   /// Run subcommands on base colvardeps object (colvar, bias, ...)
@@ -218,6 +273,9 @@ private: // TODO
   /// Help strings for each command
   std::vector<std::string> cmd_help;
 
+  /// Description of the return values of each command (may be empty)
+  std::vector<std::string> cmd_rethelp;
+
   /// Minimum number of arguments for each command
   std::vector<size_t> cmd_n_args_min;
 
@@ -227,6 +285,9 @@ private: // TODO
   /// Help strings for each command argument
   std::vector< std::vector<std::string> > cmd_arghelp;
 
+  /// Full help strings for each command
+  std::vector<std::string> cmd_full_help;
+
   /// Implementations of each command
   std::vector<int (*)(void *, int, unsigned char * const *)> cmd_fns;
 
@@ -241,6 +302,18 @@ private: // TODO
     return NULL;
   }
 
+  /// Set obj equal to x, using its string representation
+  template <typename T>
+  int set_result_text(T const &x, unsigned char *obj);
+
+  /// Code reused by instances of set_result_text()
+  template <typename T>
+  int pack_vector_elements_text(std::vector<T> const &x, std::string &x_str);
+
+  /// Code reused by all instances of set_result_text()
+  int set_result_text_from_str(std::string const &x_str, unsigned char *obj);
+
+
 };
 
 
@@ -305,13 +378,15 @@ int colvarscript::check_cmd_nargs(char const *cmd,
 {
   int const shift = cmd_arg_shift<T>();
   if (objc < shift+n_args_min) {
-    add_error_msg("Missing arguments for script function \""+std::string(cmd)+
+    add_error_msg("Insufficient number of arguments ("+cvm::to_str(objc)+
-                  "\":\n"+get_command_help(cmd));
+                  ") for script function \""+std::string(cmd)+
+                  "\":\n"+get_command_full_help(cmd));
     return COLVARSCRIPT_ERROR;
   }
   if (objc > shift+n_args_max) {
-    add_error_msg("Too many arguments for script function \""+std::string(cmd)+
+    add_error_msg("Too many arguments ("+cvm::to_str(objc)+
-                  "\":\n"+get_command_help(cmd));
+                  ") for script function \""+std::string(cmd)+
+                  "\":\n"+get_command_full_help(cmd));
     return COLVARSCRIPT_ERROR;
   }
   return COLVARSCRIPT_OK;
@@ -364,18 +439,6 @@ int colvarscript::cmd_arg_shift()
 
 extern "C" {
 
-#if defined(COLVARS_TCL)
-  /// Run the script API via Tcl command-line interface
-  /// \param clientData Not used
-  /// \param my_interp Pointer to Tcl_Interp object (read from Colvars if NULL)
-  /// \param objc Number of Tcl command parameters
-  /// \param objv Array of command parameters
-  /// \return Result of the script command
-  int tcl_run_colvarscript_command(ClientData clientData,
-                                   Tcl_Interp *interp_in,
-                                   int objc, Tcl_Obj *const objv[]);
-#endif // #if defined(COLVARS_TCL)
-
   /// Generic wrapper for string-based scripting
   int run_colvarscript_command(int objc, unsigned char *const objv[]);
 

lib/colvars/colvarscript_commands.cpp

@@ -33,6 +33,54 @@ char const **cvscript_command_names()
 }
 
 
+extern "C"
+char const *cvscript_command_help(char const *c)
+{
+  colvarscript *script = colvarscript_obj();
+  return script->get_command_help(c);
+}
+
+
+extern "C"
+char const *cvscript_command_rethelp(char const *c)
+{
+  colvarscript *script = colvarscript_obj();
+  return script->get_command_rethelp(c);
+}
+
+
+extern "C"
+char const *cvscript_command_arghelp(char const *c, int i)
+{
+  colvarscript *script = colvarscript_obj();
+  return script->get_command_arghelp(c, i);
+}
+
+
+extern "C"
+char const *cvscript_command_full_help(char const *c)
+{
+  colvarscript *script = colvarscript_obj();
+  return script->get_command_full_help(c);
+}
+
+
+extern "C"
+int cvscript_command_n_args_min(char const *c)
+{
+  colvarscript *script = colvarscript_obj();
+  return script->get_command_n_args_min(c);
+}
+
+
+extern "C"
+int cvscript_command_n_args_max(char const *c)
+{
+  colvarscript *script = colvarscript_obj();
+  return script->get_command_n_args_max(c);
+}
+
+
 // Instantiate the body of all script commands
 
 #define CVSCRIPT_COMM_FN(COMM,N_ARGS_MIN,N_ARGS_MAX,ARGS,FN_BODY)       \

lib/colvars/colvarscript_commands.h

@@ -25,7 +25,8 @@
 
 // COMM = the id of the command (must be a member of colvarscript::command)
 
-// HELP = a one-line description (C string literal) for the command
+// HELP = short description (C string literal) for the command; the second line
+//        is optional, and documents the return value (if any)
 
 // N_ARGS_MIN = the lowest number of arguments allowed
 
@@ -68,6 +69,33 @@ extern "C" {
   /// Get the names of all commands (array of strings)
   char const ** cvscript_command_names();
 
+  /// Get the help summary of the given command
+  /// \param cmd Name of the command's function (e.g. "cv_units")
+  char const *cvscript_command_help(char const *cmd);
+
+  /// Get description of the return value of a command
+  /// \param cmd Name of the command's function (e.g. "cv_units")
+  char const *cvscript_command_rethelp(char const *cmd);
+
+  /// Get description of the arguments of a command (excluding prefix)
+  /// \param cmd Name of the command's function (e.g. "cv_units")
+  /// \param i Index of the argument; 0 is the first argument after the
+  /// prefix, e.g. "value" has an index of 0 in the array of arguments:
+  /// { "cv", "colvar", "xi", "value" }
+  char const *cvscript_command_arghelp(char const *cmd, int i);
+
+  /// Get the full help string of a command
+  /// \param cmd Name of the command's function (e.g. "cv_units")
+  char const *cvscript_command_full_help(char const *cmd);
+
+  /// Get number of required arguments (excluding prefix)
+  /// \param cmd Name of the command's function (e.g. "cv_units")
+  int cvscript_command_n_args_min(char const *cmd);
+
+  /// Get number of total arguments (excluding prefix)
+  /// \param cmd Name of the command's function (e.g. "cv_units")
+  int cvscript_command_n_args_max(char const *cmd);
+
 }
 
 #endif
@@ -135,7 +163,8 @@ CVSCRIPT(cv_delete,
          )
 
 CVSCRIPT(cv_frame,
-         "Get or set current frame number (VMD only)",
+         "Get or set current frame number (VMD only)\n"
+         "frame : integer - Frame number",
          0, 1,
          "frame : integer - Frame number",
          char const *arg =
@@ -143,7 +172,7 @@ CVSCRIPT(cv_frame,
          if (arg == NULL) {
            long int f = -1;
            if (script->proxy()->get_frame(f) == COLVARS_OK) {
-             script->set_result_str(cvm::to_str(f));
+             script->set_result_long_int(f);
              return COLVARS_OK;
            } else {
              script->add_error_msg("Frame number is not available");
@@ -161,8 +190,90 @@ CVSCRIPT(cv_frame,
          return COLVARS_OK;
          )
 
+CVSCRIPT(cv_getatomappliedforces,
+         "Get the list of forces applied by Colvars to atoms\n"
+         "forces : array of arrays of floats - Atomic forces",
+         0, 0,
+         "",
+         script->set_result_rvector_vec(*(script->proxy()->get_atom_applied_forces()));
+         return COLVARS_OK;
+         )
+
+CVSCRIPT(cv_getatomappliedforcesmax,
+         "Get the maximum norm of forces applied by Colvars to atoms\n"
+         "force : float - Maximum atomic force",
+         0, 0,
+         "",
+         script->set_result_real(script->proxy()->max_atoms_applied_force());
+         return COLVARS_OK;
+         )
+
+CVSCRIPT(cv_getatomappliedforcesmaxid,
+         "Get the atom ID with the largest applied force\n"
+         "id : int - ID of the atom with the maximum atomic force",
+         0, 0,
+         "",
+         script->set_result_int(script->proxy()->max_atoms_applied_force_id());
+         return COLVARS_OK;
+         )
+
+CVSCRIPT(cv_getatomappliedforcesrms,
+         "Get the root-mean-square norm of forces applied by Colvars to atoms\n"
+         "force : float - RMS atomic force",
+         0, 0,
+         "",
+         script->set_result_real(script->proxy()->rms_atoms_applied_force());
+         return COLVARS_OK;
+         )
+
+CVSCRIPT(cv_getatomids,
+         "Get the list of indices of atoms used in Colvars\n"
+         "indices : array of ints - Atomic indices",
+         0, 0,
+         "",
+         script->set_result_int_vec(*(script->proxy()->get_atom_ids()));
+         return COLVARS_OK;
+         )
+
+CVSCRIPT(cv_getatomcharges,
+         "Get the list of charges of atoms used in Colvars\n"
+         "charges : array of floats - Atomic charges",
+         0, 0,
+         "",
+         script->set_result_real_vec(*(script->proxy()->get_atom_charges()));
+         return COLVARS_OK;
+         )
+
+CVSCRIPT(cv_getatommasses,
+         "Get the list of masses of atoms used in Colvars\n"
+         "masses : array of floats - Atomic masses",
+         0, 0,
+         "",
+         script->set_result_real_vec(*(script->proxy()->get_atom_masses()));
+         return COLVARS_OK;
+         )
+
+CVSCRIPT(cv_getatompositions,
+         "Get the list of cached positions of atoms used in Colvars\n"
+         "positions : array of arrays of floats - Atomic positions",
+         0, 0,
+         "",
+         script->set_result_rvector_vec(*(script->proxy()->get_atom_positions()));
+         return COLVARS_OK;
+         )
+
+CVSCRIPT(cv_getatomtotalforces,
+         "Get the list of cached total forces of atoms used in Colvars\n"
+         "forces : array of arrays of floats - Atomic total foces",
+         0, 0,
+         "",
+         script->set_result_rvector_vec(*(script->proxy()->get_atom_total_forces()));
+         return COLVARS_OK;
+         )
+
 CVSCRIPT(cv_getconfig,
-         "Get the module's configuration string read so far",
+         "Get the module's configuration string read so far\n"
+         "conf : string - Current configuration string",
          0, 0,
          "",
          script->set_result_str(cvm::main()->get_config());
@@ -170,15 +281,17 @@ CVSCRIPT(cv_getconfig,
          )
 
 CVSCRIPT(cv_getenergy,
-         "Get the current Colvars energy",
+         "Get the current Colvars energy\n"
+         "E : float - Amount of energy (internal units)",
          0, 0,
          "",
-         script->set_result_str(cvm::to_str(cvm::main()->total_bias_energy));
+         script->set_result_real(cvm::main()->total_bias_energy);
          return COLVARS_OK;
          )
 
 CVSCRIPT(cv_help,
-         "Get the help string of the Colvars scripting interface",
+         "Get the help string of the Colvars scripting interface\n"
+         "help : string - Help string",
          0, 1,
          "command : string - Get the help string of this specific command",
          unsigned char *const cmdobj =
@@ -205,7 +318,8 @@ CVSCRIPT(cv_help,
          )
 
 CVSCRIPT(cv_list,
-         "Return a list of all variables or biases",
+         "Return a list of all variables or biases\n"
+         "list : sequence of strings - List of elements",
          0, 1,
          "param : string - \"colvars\" or \"biases\"; default is \"colvars\"",
          std::string res;
@@ -235,7 +349,8 @@ CVSCRIPT(cv_list,
          )
 
 CVSCRIPT(cv_listcommands,
-         "Get the list of script functions, prefixed with \"cv_\", \"colvar_\" or \"bias_\"",
+         "Get the list of script functions, prefixed with \"cv_\", \"colvar_\" or \"bias_\"\n"
+         "list : sequence of strings - List of commands",
          0, 0,
          "",
          int const n_commands = cvscript_n_commands();
@@ -249,6 +364,20 @@ CVSCRIPT(cv_listcommands,
          return COLVARS_OK;
          )
 
+CVSCRIPT(cv_listindexfiles,
+         "Get a list of the index files loaded in this session",
+         0, 0,
+         "",
+         int const n_files = script->module()->index_file_names.size();
+         std::string result;
+         for (int i = 0; i < n_files; i++) {
+           if (i > 0) result.append(1, ' ');
+           result.append(script->module()->index_file_names[i]);
+         }
+         script->set_result_str(result);
+         return COLVARS_OK;
+         )
+
 CVSCRIPT(cv_load,
          "Load data from a state file into all matching colvars and biases",
          1, 1,
@@ -280,15 +409,16 @@ CVSCRIPT(cv_loadfromstring,
          )
 
 CVSCRIPT(cv_molid,
-         "Get or set the molecule ID on which Colvars is defined (VMD only)",
+         "Get or set the molecule ID on which Colvars is defined (VMD only)\n"
+         "molid : integer - Current molecule ID",
          0, 1,
-         "molid : integer - Molecule ID; -1 means undefined",
+         "molid : integer - New molecule ID; -1 means undefined",
          char const *arg =
            script->obj_to_str(script->get_module_cmd_arg(0, objc, objv));
          if (arg == NULL) {
            int molid = -1;
            script->proxy()->get_molid(molid);
-           script->set_result_str(cvm::to_str(molid));
+           script->set_result_int(molid);
            return COLVARS_OK;
          } else {
            script->add_error_msg("Error: To change the molecule ID in VMD, use cv delete first.");
@@ -297,7 +427,8 @@ CVSCRIPT(cv_molid,
          )
 
 CVSCRIPT(cv_printframe,
-         "Return the values that would be written to colvars.traj",
+         "Return the values that would be written to colvars.traj\n"
+         "values : string - The values\n",
          0, 0,
          "",
          std::ostringstream os;
@@ -307,7 +438,8 @@ CVSCRIPT(cv_printframe,
          )
 
 CVSCRIPT(cv_printframelabels,
-         "Return the labels that would be written to colvars.traj",
+         "Return the labels that would be written to colvars.traj\n"
+         "Labels : string - The labels",
          0, 0,
          "",
          std::ostringstream os;
@@ -348,14 +480,16 @@ CVSCRIPT(cv_save,
          )
 
 CVSCRIPT(cv_savetostring,
-         "Write the Colvars state to a string and return it",
+         "Write the Colvars state to a string and return it\n"
+         "state : string - The saved state",
          0, 0,
          "",
          return script->module()->write_restart_string(script->modify_str_result());
          )
 
 CVSCRIPT(cv_units,
-         "Get or set the current Colvars unit system",
+         "Get or set the current Colvars unit system\n"
+         "units : string - The current unit system",
          0, 1,
          "units : string - The new unit system",
          char const *argstr =
@@ -390,7 +524,8 @@ CVSCRIPT(cv_update,
          )
 
 CVSCRIPT(cv_version,
-         "Get the Colvars Module version number",
+         "Get the Colvars Module version number\n"
+         "version : string - Colvars version",
          0, 0,
          "",
          script->set_result_str(COLVARS_VERSION);

lib/colvars/colvarscript_commands_bias.h

@@ -9,15 +9,17 @@
 
 
 CVSCRIPT(bias_bin,
-         "Get the current grid bin index (1D ABF only for now)",
+         "Get the current grid bin index (1D ABF only for now)\n"
+         "bin : integer - Bin index",
          0, 0,
          "",
-         script->set_result_str(cvm::to_str(this_bias->current_bin()));
+         script->set_result_int(this_bias->current_bin());
          return COLVARS_OK;
          )
 
 CVSCRIPT(bias_bincount,
-         "Get the number of samples at the given grid bin (1D ABF only for now)",
+         "Get the number of samples at the given grid bin (1D ABF only for now)\n"
+         "samples : integer - Number of samples",
          0, 1,
          "index : integer - Grid index; defaults to current bin",
          int index = this_bias->current_bin();
@@ -30,12 +32,13 @@ CVSCRIPT(bias_bincount,
              return COLVARSCRIPT_ERROR;
            }
          }
-         script->set_result_str(cvm::to_str(this_bias->bin_count(index)));
+         script->set_result_int(this_bias->bin_count(index));
          return COLVARS_OK;
          )
 
 CVSCRIPT(bias_binnum,
-         "Get the total number of grid points of this bias (1D ABF only for now)",
+         "Get the total number of grid points of this bias (1D ABF only for now)\n"
+         "Bins : integer - Number of grid points",
          0, 0,
          "",
          int r = this_bias->bin_num();
@@ -44,7 +47,7 @@ CVSCRIPT(bias_binnum,
                                  this_bias->name);
            return COLVARSCRIPT_ERROR;
          }
-         script->set_result_str(cvm::to_str(r));
+         script->set_result_int(r);
          return COLVARS_OK;
          )
 
@@ -57,22 +60,25 @@ CVSCRIPT(bias_delete,
          )
 
 CVSCRIPT(bias_energy,
-         "Get the current energy of this bias",
+         "Get the current energy of this bias\n"
+         "E : float - Energy value",
          0, 0,
          "",
-         script->set_result_str(cvm::to_str(this_bias->get_energy()));
+         script->set_result_real(this_bias->get_energy());
          return COLVARS_OK;
          )
 
 CVSCRIPT(bias_get,
-         "Get the value of the given feature for this bias",
+         "Get the value of the given feature for this bias\n"
+         "state : 1/0 - State of the given feature",
          1, 1,
          "feature : string - Name of the feature",
          return script->proc_features(this_bias, objc, objv);
          )
 
 CVSCRIPT(bias_getconfig,
-         "Return the configuration string of this bias",
+         "Return the configuration string of this bias\n"
+         "conf : string - Current configuration string",
          0, 0,
          "",
          script->set_result_str(this_bias->get_config());
@@ -80,7 +86,8 @@ CVSCRIPT(bias_getconfig,
          )
 
 CVSCRIPT(bias_help,
-         "Get a help summary or the help string of one bias subcommand",
+         "Get a help summary or the help string of one bias subcommand\n"
+         "help : string - Help string",
          0, 1,
          "command : string - Get the help string of this specific command",
          unsigned char *const cmdobj =
@@ -129,7 +136,8 @@ CVSCRIPT(bias_save,
          )
 
 CVSCRIPT(bias_savetostring,
-         "Save data from this bias into a string and return it",
+         "Save data from this bias into a string and return it\n"
+         "state : string - The bias state",
          0, 0,
          "",
          return this_bias->write_state_string(script->modify_str_result());
@@ -156,7 +164,8 @@ CVSCRIPT(bias_share,
          )
 
 CVSCRIPT(bias_state,
-         "Print a string representation of the feature state of this bias",
+         "Print a string representation of the feature state of this bias\n"
+         "state : string - String representation of the bias features",
          0, 0,
          "",
          this_bias->print_state();
@@ -164,10 +173,11 @@ CVSCRIPT(bias_state,
          )
 
 CVSCRIPT(bias_update,
-         "Recompute this bias and return its up-to-date energy",
+         "Recompute this bias and return its up-to-date energy\n"
+         "E : float - Energy value",
          0, 0,
          "",
          this_bias->update();
-         script->set_result_str(cvm::to_str(this_bias->get_energy()));
+         script->set_result_colvarvalue(this_bias->get_energy());
          return COLVARS_OK;
          )

lib/colvars/colvarscript_commands_colvar.h

@@ -9,7 +9,8 @@
 
 
 CVSCRIPT(colvar_addforce,
-         "Apply the given force onto this colvar and return the same",
+         "Apply the given force onto this colvar and return the same\n"
+         "force : float or array - Applied force; matches colvar dimensionality",
          1, 1,
          "force : float or array - Applied force; must match colvar dimensionality",
          std::string const f_str(script->obj_to_str(script->get_colvar_cmd_arg(0, objc, objv)));
@@ -23,7 +24,7 @@ CVSCRIPT(colvar_addforce,
            return COLVARSCRIPT_ERROR;
          }
          this_colvar->add_bias_force(force);
-         script->set_result_str(force.to_simple_string());
+         script->set_result_colvarvalue(force);
          return COLVARS_OK;
          )
 
@@ -56,22 +57,25 @@ CVSCRIPT(colvar_delete,
          )
 
 CVSCRIPT(colvar_get,
-         "Get the value of the given feature for this colvar",
+         "Get the value of the given feature for this colvar\n"
+         "state : 1/0 - State of the given feature",
          1, 1,
          "feature : string - Name of the feature",
          return script->proc_features(this_colvar, objc, objv);
          )
 
 CVSCRIPT(colvar_getappliedforce,
-         "Return the total of the forces applied to this colvar",
+         "Return the total of the forces applied to this colvar\n"
+         "force : float - Applied force; matches the colvar dimensionality",
          0, 0,
          "",
-         script->set_result_str((this_colvar->applied_force()).to_simple_string());
+         script->set_result_colvarvalue(this_colvar->applied_force());
          return COLVARS_OK;
          )
 
 CVSCRIPT(colvar_getatomgroups,
-         "Return the atom indices used by this colvar as a list of lists",
+         "Return the atom indices used by this colvar as a list of lists\n"
+         "groups : array of arrays of ints - Atom indices",
          0, 0,
          "",
          std::string result;
@@ -91,21 +95,17 @@ CVSCRIPT(colvar_getatomgroups,
          )
 
 CVSCRIPT(colvar_getatomids,
-         "Return the list of atom indices used by this colvar",
+         "Return the list of atom indices used by this colvar\n"
+         "indices : array of ints - Atom indices",
          0, 0,
          "",
-         std::string result;
+         script->set_result_int_vec(this_colvar->atom_ids);
-         std::vector<int>::iterator li = this_colvar->atom_ids.begin();
-         for ( ; li != this_colvar->atom_ids.end(); ++li) {
-           result += cvm::to_str(*li);
-           result += " ";
-         }
-         script->set_result_str(result);
          return COLVARS_OK;
          )
 
 CVSCRIPT(colvar_getconfig,
-         "Return the configuration string of this colvar",
+         "Return the configuration string of this colvar\n"
+         "conf : string - Current configuration string",
          0, 0,
          "",
          script->set_result_str(this_colvar->get_config());
@@ -113,35 +113,34 @@ CVSCRIPT(colvar_getconfig,
          )
 
 CVSCRIPT(colvar_getgradients,
-         "Return the atomic gradients of this colvar",
+         "Return the atomic gradients of this colvar\n"
+         "gradients : array of arrays of floats - Atomic gradients",
          0, 0,
          "",
-         std::string result;
+         script->set_result_rvector_vec(this_colvar->atomic_gradients);
-         std::vector<cvm::rvector>::iterator li =
-           this_colvar->atomic_gradients.begin();
-         for ( ; li != this_colvar->atomic_gradients.end(); ++li) {
-           result += "{";
-           int j;
-           for (j = 0; j < 3; ++j) {
-             result += cvm::to_str((*li)[j]);
-             result += " ";
-           }
-           result += "} ";
-         }
-         script->set_result_str(result);
          return COLVARS_OK;
          )
 
 CVSCRIPT(colvar_gettotalforce,
-         "Return the sum of internal and external forces to this colvar",
+         "Return the sum of internal and external forces to this colvar\n"
+         "force : float - Total force; matches the colvar dimensionality",
          0, 0,
          "",
-         script->set_result_str((this_colvar->total_force()).to_simple_string());
+         script->set_result_colvarvalue(this_colvar->total_force());
+         return COLVARS_OK;
+         )
+
+CVSCRIPT(colvar_getvolmapids,
+         "Return the list of volumetric map indices used by this colvar",
+         0, 0,
+         "",
+         script->set_result_int_vec(this_colvar->get_volmap_ids());
          return COLVARS_OK;
          )
 
 CVSCRIPT(colvar_help,
-         "Get a help summary or the help string of one colvar subcommand",
+         "Get a help summary or the help string of one colvar subcommand\n"
+         "help : string - Help string",
          0, 1,
          "command : string - Get the help string of this specific command",
          unsigned char *const cmdobj =
@@ -167,7 +166,7 @@ CVSCRIPT(colvar_modifycvcs,
          "Modify configuration of individual components by passing string arguments",
          1, 1,
          "confs : sequence of strings - New configurations; empty strings are skipped",
-         std::vector<std::string> const confs(script->proxy()->script_obj_to_str_vector(script->get_colvar_cmd_arg(0, objc, objv)));
+         std::vector<std::string> const confs(script->obj_to_str_vector(script->get_colvar_cmd_arg(0, objc, objv)));
          cvm::increase_depth();
          int res = this_colvar->update_cvc_config(confs);
          cvm::decrease_depth();
@@ -180,10 +179,11 @@ CVSCRIPT(colvar_modifycvcs,
          )
 
 CVSCRIPT(colvar_run_ave,
-         "Get the current running average of the value of this colvar",
+         "Get the current running average of the value of this colvar\n"
+         "value : float or array - Averaged value; matches the colvar dimensionality",
          0, 0,
          "",
-         script->set_result_str(this_colvar->run_ave().to_simple_string());
+         script->set_result_colvarvalue(this_colvar->run_ave());
          return COLVARS_OK;
          )
 
@@ -196,7 +196,8 @@ CVSCRIPT(colvar_set,
          )
 
 CVSCRIPT(colvar_state,
-         "Print a string representation of the feature state of this colvar",
+         "Print a string representation of the feature state of this colvar\n"
+         "state : string - The feature state",
          0, 0,
          "",
          this_colvar->print_state();
@@ -204,7 +205,8 @@ CVSCRIPT(colvar_state,
          )
 
 CVSCRIPT(colvar_type,
-         "Get the type description of this colvar",
+         "Get the type description of this colvar\n"
+         "type : string - Type description",
          0, 0,
          "",
          script->set_result_str(this_colvar->value().type_desc(this_colvar->value().value_type));
@@ -212,25 +214,28 @@ CVSCRIPT(colvar_type,
          )
 
 CVSCRIPT(colvar_update,
-         "Recompute this colvar and return its up-to-date value",
+         "Recompute this colvar and return its up-to-date value\n"
+         "value : float or array - Current value; matches the colvar dimensionality",
          0, 0,
          "",
          this_colvar->calc();
          this_colvar->update_forces_energy();
-         script->set_result_str((this_colvar->value()).to_simple_string());
+         script->set_result_colvarvalue(this_colvar->value());
          return COLVARS_OK;
          )
 
 CVSCRIPT(colvar_value,
-         "Get the current value of this colvar",
+         "Get the current value of this colvar\n"
+         "value : float or array - Current value; matches the colvar dimensionality",
          0, 0,
          "",
-         script->set_result_str(this_colvar->value().to_simple_string());
+         script->set_result_colvarvalue(this_colvar->value());
          return COLVARS_OK;
          )
 
 CVSCRIPT(colvar_width,
-         "Get the width of this colvar",
+         "Get the width of this colvar\n"
+         "width : float - Value of the width",
          0, 0,
          "",
          script->set_result_str(cvm::to_str(this_colvar->width, 0,

lib/colvars/colvartypes.cpp

@@ -256,6 +256,7 @@ namespace {
 colvarmodule::rotation::rotation()
 {
   b_debug_gradients = false;
+  lambda = 0.0;
 #ifdef COLVARS_LAMMPS
   jacobi = new_Jacobi_solver(4);
 #else
@@ -268,6 +269,7 @@ colvarmodule::rotation::rotation(cvm::quaternion const &qi)
   : q(qi)
 {
   b_debug_gradients = false;
+  lambda = 0.0;
 #ifdef COLVARS_LAMMPS
   jacobi = new_Jacobi_solver(4);
 #else
@@ -283,6 +285,7 @@ colvarmodule::rotation::rotation(cvm::real angle, cvm::rvector const &axis)
   cvm::real const sina = cvm::sin(angle/2.0);
   q = cvm::quaternion(cvm::cos(angle/2.0),
                       sina * axis_n.x, sina * axis_n.y, sina * axis_n.z);
+  lambda = 0.0;
 #ifdef COLVARS_LAMMPS
   jacobi = new_Jacobi_solver(4);
 #else

lib/colvars/colvarvalue.cpp

@@ -224,7 +224,7 @@ void colvarvalue::is_derivative()
 
 
 colvarvalue::colvarvalue(colvarvalue const &x)
-  : value_type(x.type())
+  : value_type(x.type()), real_value(0.0)
 {
   switch (x.type()) {
   case type_scalar:
@@ -252,6 +252,7 @@ colvarvalue::colvarvalue(colvarvalue const &x)
 
 
 colvarvalue::colvarvalue(cvm::vector1d<cvm::real> const &v, Type vti)
+  : real_value(0.0)
 {
   if ((vti != type_vector) && (v.size() != num_dimensions(vti))) {
     cvm::error("Error: trying to initialize a variable of type \""+type_desc(vti)+
@@ -579,16 +580,7 @@ colvarvalue colvarvalue::dist2_grad(colvarvalue const &x2) const
       cvm::rvector const &v1 = this->rvector_value;
       cvm::rvector const &v2 = x2.rvector_value;
       cvm::real const cos_t = v1 * v2;
-      cvm::real const sin_t = cvm::sqrt(1.0 - cos_t*cos_t);
+      return colvarvalue(2.0 * (cos_t * v1 - v2), colvarvalue::type_unit3vectorderiv);
-      return colvarvalue( 2.0 * sin_t *
-                          cvm::rvector((-1.0) * sin_t * v2.x +
-                                       cos_t/sin_t * (v1.x - cos_t*v2.x),
-                                       (-1.0) * sin_t * v2.y +
-                                       cos_t/sin_t * (v1.y - cos_t*v2.y),
-                                       (-1.0) * sin_t * v2.z +
-                                       cos_t/sin_t * (v1.z - cos_t*v2.z)
-                                       ),
-                          colvarvalue::type_unit3vectorderiv );
     }
   case colvarvalue::type_quaternion:
   case colvarvalue::type_quaternionderiv:

lib/colvars/colvarvalue.h

@@ -124,7 +124,7 @@ public:
 
   /// Constructor from a type specification
   inline colvarvalue(Type const &vti)
-    : value_type(vti)
+    : value_type(vti), real_value(0.0)
   {
     reset();
   }
@@ -138,12 +138,12 @@ public:
   /// by default a type \link type_3vector \endlink , if you want a
   /// \link type_unit3vector \endlink you must set it explicitly)
   inline colvarvalue(cvm::rvector const &v, Type vti = type_3vector)
-    : value_type(vti), rvector_value(v)
+    : value_type(vti), real_value(0.0), rvector_value(v)
   {}
 
   /// \brief Copy constructor from quaternion base type
   inline colvarvalue(cvm::quaternion const &q, Type vti = type_quaternion)
-    : value_type(vti), quaternion_value(q)
+    : value_type(vti), real_value(0.0), quaternion_value(q)
   {}
 
   /// Copy constructor from vector1d base type

src/COLVARS/colvarproxy_lammps_version.h

@@ -1,3 +1,3 @@
 #ifndef COLVARPROXY_VERSION
-#define COLVARPROXY_VERSION "2020-04-07"
+#define COLVARPROXY_VERSION "2021-03-02"
 #endif

src/COLVARS/fix_colvars.cpp

@@ -480,7 +480,7 @@ void FixColvars::one_time_init()
   memory->create(force_buf,3*num_coords,"colvars:force_buf");
 
   if (me == 0) {
-    std::vector<int> &tl = *(proxy->modify_atom_ids());
+    std::vector<int> const &tl = *(proxy->get_atom_ids());
     inthash_t *hashtable=new inthash_t;
     inthash_init(hashtable, num_coords);
     idmap = (void *)hashtable;
@@ -563,7 +563,7 @@ void FixColvars::setup(int vflag)
 
   if (me == 0) {
 
-    std::vector<int>           &id = *(proxy->modify_atom_ids());
+    std::vector<int>     const &id = *(proxy->get_atom_ids());
     std::vector<int>           &tp = *(proxy->modify_atom_types());
     std::vector<cvm::atom_pos> &cd = *(proxy->modify_atom_positions());
     std::vector<cvm::rvector>  &of = *(proxy->modify_atom_total_forces());
@@ -836,7 +836,7 @@ void FixColvars::post_force(int /*vflag*/)
 
   if (me == 0) {
 
-    std::vector<cvm::rvector> &fo = *(proxy->modify_atom_new_colvar_forces());
+    std::vector<cvm::rvector> &fo = *(proxy->modify_atom_applied_forces());
 
     double *fbuf = force_buf;
     for (int j=0; j < num_coords; ++j) {