Update Colvars library to version 2019-04-26

The following is list of relevant issues fixed and improvements:

Fix forces and missing output of runtime histogram for histogramRestraint
https://github.com/Colvars/colvars/pull/246

Use fix_modify to add configuration to Colvars:
https://github.com/Colvars/colvars/pull/216

Fix componentCoeff and name not working with orientationAngle components:
https://github.com/Colvars/colvars/issues/213

Fix 1-timestep offset with extendedLagrangian:
https://github.com/Colvars/colvars/pull/210

Changes to improve compiler support:
https://github.com/Colvars/colvars/pull/203

Fix ignored anisotropic cutoff3 for groupCoordNum:
https://github.com/Colvars/colvars/pull/202

New dipoleMagnitude variable:
https://github.com/Colvars/colvars/pull/198

Parser improvements:
https://github.com/Colvars/colvars/pull/196

doc/src/fix_colvars.txt

@@ -98,6 +98,16 @@ fix to add the energy change from the biasing force added by the fix
 to the system's potential energy as part of "thermodynamic
 output"_thermo_style.html.
 
+The {fix_modify configfile <config file>} option allows to add settings
+from an additional config file to the colvars module. This option can
+only be used, after the system has been initialized with a "run"_run.html
+command.
+
+The {fix_modify config <quoted string>} option allows to add settings
+from inline strings. Those have to fit on a single line when enclosed
+in a pair of double quotes ("), or can span multiple lines when bracketed
+by a pair of triple double quotes (""", like python embedded documentation).
+
 This fix computes a global scalar which can be accessed by various
 "output commands"_Howto_output.html.  The scalar is the cumulative
 energy change due to this fix.  The scalar value calculated by this

lib/colvars/colvar.cpp

@@ -24,12 +24,13 @@ colvar::colvar()
 {
   runave_os = NULL;
 
-  prev_timestep = -1;
+  prev_timestep = -1L;
   after_restart = false;
   kinetic_energy = 0.0;
   potential_energy = 0.0;
 
-  init_cv_requires();
+  description = "uninitialized colvar";
+  init_dependencies();
 }
 
 
@@ -193,7 +194,7 @@ int colvar::init(std::string const &conf)
   {
     bool homogeneous = is_enabled(f_cv_linear);
     for (i = 0; i < cvcs.size(); i++) {
-      if ((std::fabs(cvcs[i]->sup_coeff) - 1.0) > 1.0e-10) {
+      if ((cvm::fabs(cvcs[i]->sup_coeff) - 1.0) > 1.0e-10) {
         homogeneous = false;
       }
     }
@@ -224,7 +225,7 @@ int colvar::init(std::string const &conf)
   // Allow scripted/custom functions to be defined as periodic
   if ( (is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function)) && is_enabled(f_cv_scalar) ) {
     if (get_keyval(conf, "period", period, 0.)) {
-      set_enabled(f_cv_periodic, true);
+      enable(f_cv_periodic);
       get_keyval(conf, "wrapAround", wrap_center, 0.);
     }
   }
@@ -471,7 +472,7 @@ int colvar::init_grid_parameters(std::string const &conf)
     if (get_keyval(conf, "lowerWallConstant", lower_wall_k, 0.0,
                    parse_silent)) {
       cvm::log("Reading legacy options lowerWall and lowerWallConstant: "
-               "consider using a harmonicWalls restraint.\n");
+               "consider using a harmonicWalls restraint\n(caution: force constant would then be scaled by width^2).\n");
       lower_wall.type(value());
       if (!get_keyval(conf, "lowerWall", lower_wall, lower_boundary)) {
         cvm::log("Warning: lowerWall will need to be "
@@ -485,7 +486,7 @@ int colvar::init_grid_parameters(std::string const &conf)
     if (get_keyval(conf, "upperWallConstant", upper_wall_k, 0.0,
                    parse_silent)) {
       cvm::log("Reading legacy options upperWall and upperWallConstant: "
-               "consider using a harmonicWalls restraint.\n");
+               "consider using a harmonicWalls restraint\n(caution: force constant would then be scaled by width^2).\n");
       upper_wall.type(value());
       if (!get_keyval(conf, "upperWall", upper_wall, upper_boundary)) {
         cvm::log("Warning: upperWall will need to be "
@@ -562,13 +563,13 @@ int colvar::init_extended_Lagrangian(std::string const &conf)
   get_keyval_feature(this, conf, "extendedLagrangian", f_cv_extended_Lagrangian, false);
 
   if (is_enabled(f_cv_extended_Lagrangian)) {
-    cvm::real temp, tolerance, period;
+    cvm::real temp, tolerance, extended_period;
 
     cvm::log("Enabling the extended Lagrangian term for colvar \""+
              this->name+"\".\n");
 
-    xr.type(value());
+    x_ext.type(value());
-    vr.type(value());
+    v_ext.type(value());
     fr.type(value());
 
     const bool found = get_keyval(conf, "extendedTemp", temp, cvm::temperature());
@@ -590,11 +591,11 @@ int colvar::init_extended_Lagrangian(std::string const &conf)
     ext_force_k = cvm::boltzmann() * temp / (tolerance * tolerance);
     cvm::log("Computed extended system force constant: " + cvm::to_str(ext_force_k) + " [E]/U^2");
 
-    get_keyval(conf, "extendedTimeConstant", period, 200.0);
+    get_keyval(conf, "extendedTimeConstant", extended_period, 200.0);
-    if (period <= 0.0) {
+    if (extended_period <= 0.0) {
       cvm::error("Error: \"extendedTimeConstant\" must be positive.\n", INPUT_ERROR);
     }
-    ext_mass = (cvm::boltzmann() * temp * period * period)
+    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)");
 
@@ -615,7 +616,7 @@ int colvar::init_extended_Lagrangian(std::string const &conf)
       enable(f_cv_Langevin);
       ext_gamma *= 1.0e-3; // correct as long as input is required in ps-1 and cvm::dt() is in fs
       // Adjust Langevin sigma for slow time step if time_step_factor != 1
-      ext_sigma = std::sqrt(2.0 * cvm::boltzmann() * temp * ext_gamma * ext_mass / (cvm::dt() * cvm::real(time_step_factor)));
+      ext_sigma = cvm::sqrt(2.0 * cvm::boltzmann() * temp * ext_gamma * ext_mass / (cvm::dt() * cvm::real(time_step_factor)));
     }
   }
 
@@ -761,6 +762,8 @@ int colvar::init_components(std::string const &conf)
     "weighted by inverse power", "distanceInv");
   error_code |= init_components_type<distance_pairs>(conf, "N1xN2-long vector "
     "of pairwise distances", "distancePairs");
+  error_code |= init_components_type<dipole_magnitude>(conf, "dipole magnitude",
+    "dipoleMagnitude");
   error_code |= init_components_type<coordnum>(conf, "coordination "
     "number", "coordNum");
   error_code |= init_components_type<selfcoordnum>(conf, "self-coordination "
@@ -831,22 +834,25 @@ void colvar::build_atom_list(void)
 
   for (size_t i = 0; i < cvcs.size(); i++) {
     for (size_t j = 0; j < cvcs[i]->atom_groups.size(); j++) {
-      cvm::atom_group &ag = *(cvcs[i]->atom_groups[j]);
+      cvm::atom_group const &ag = *(cvcs[i]->atom_groups[j]);
       for (size_t k = 0; k < ag.size(); k++) {
         temp_id_list.push_back(ag[k].id);
       }
+      if (ag.is_enabled(f_ag_fitting_group) && ag.is_enabled(f_ag_fit_gradients)) {
+        cvm::atom_group const &fg = *(ag.fitting_group);
+        for (size_t k = 0; k < fg.size(); k++) {
+          temp_id_list.push_back(fg[k].id);
+        }
+      }
     }
   }
 
   temp_id_list.sort();
   temp_id_list.unique();
 
-  // atom_ids = std::vector<int> (temp_id_list.begin(), temp_id_list.end());
-  unsigned int id_i = 0;
   std::list<int>::iterator li;
   for (li = temp_id_list.begin(); li != temp_id_list.end(); ++li) {
-    atom_ids[id_i] = *li;
+    atom_ids.push_back(*li);
-    id_i++;
   }
 
   temp_id_list.clear();
@@ -934,27 +940,164 @@ int colvar::parse_analysis(std::string const &conf)
 }
 
 
-void colvar::setup() {
+int colvar::init_dependencies() {
-  // loop over all components to reset masses of all groups
+  size_t i;
+  if (features().size() == 0) {
+    for (i = 0; i < f_cv_ntot; i++) {
+      modify_features().push_back(new feature);
+    }
+
+    init_feature(f_cv_active, "active", f_type_dynamic);
+    // Do not require f_cvc_active in children, as some components may be disabled
+    // Colvars must be either a linear combination, or scalar (and polynomial) or scripted/custom
+    require_feature_alt(f_cv_active, f_cv_scalar, f_cv_linear, f_cv_scripted, f_cv_custom_function);
+
+    init_feature(f_cv_awake, "awake", f_type_static);
+    require_feature_self(f_cv_awake, f_cv_active);
+
+    init_feature(f_cv_gradient, "gradient", f_type_dynamic);
+    require_feature_children(f_cv_gradient, f_cvc_gradient);
+
+    init_feature(f_cv_collect_gradient, "collect gradient", f_type_dynamic);
+    require_feature_self(f_cv_collect_gradient, f_cv_gradient);
+    require_feature_self(f_cv_collect_gradient, f_cv_scalar);
+    // The following exlusion could be lifted by implementing the feature
+    exclude_feature_self(f_cv_collect_gradient, f_cv_scripted);
+    require_feature_children(f_cv_collect_gradient, f_cvc_explicit_gradient);
+
+    init_feature(f_cv_fdiff_velocity, "velocity from finite differences", f_type_dynamic);
+
+    // System force: either trivial (spring force); through extended Lagrangian, or calculated explicitly
+    init_feature(f_cv_total_force, "total force", f_type_dynamic);
+    require_feature_alt(f_cv_total_force, f_cv_extended_Lagrangian, f_cv_total_force_calc);
+
+    // Deps for explicit total force calculation
+    init_feature(f_cv_total_force_calc, "total force calculation", f_type_dynamic);
+    require_feature_self(f_cv_total_force_calc, f_cv_scalar);
+    require_feature_self(f_cv_total_force_calc, f_cv_linear);
+    require_feature_children(f_cv_total_force_calc, f_cvc_inv_gradient);
+    require_feature_self(f_cv_total_force_calc, f_cv_Jacobian);
+
+    init_feature(f_cv_Jacobian, "Jacobian derivative", f_type_dynamic);
+    require_feature_self(f_cv_Jacobian, f_cv_scalar);
+    require_feature_self(f_cv_Jacobian, f_cv_linear);
+    require_feature_children(f_cv_Jacobian, f_cvc_Jacobian);
+
+    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
+
+    init_feature(f_cv_extended_Lagrangian, "extended Lagrangian", f_type_user);
+    require_feature_self(f_cv_extended_Lagrangian, f_cv_scalar);
+    require_feature_self(f_cv_extended_Lagrangian, f_cv_gradient);
+
+    init_feature(f_cv_Langevin, "Langevin dynamics", f_type_user);
+    require_feature_self(f_cv_Langevin, f_cv_extended_Lagrangian);
+
+    init_feature(f_cv_linear, "linear", f_type_static);
+
+    init_feature(f_cv_scalar, "scalar", f_type_static);
+
+    init_feature(f_cv_output_energy, "output energy", f_type_user);
+
+    init_feature(f_cv_output_value, "output value", f_type_user);
+
+    init_feature(f_cv_output_velocity, "output velocity", f_type_user);
+    require_feature_self(f_cv_output_velocity, f_cv_fdiff_velocity);
+
+    init_feature(f_cv_output_applied_force, "output applied force", f_type_user);
+
+    init_feature(f_cv_output_total_force, "output total force", f_type_user);
+    require_feature_self(f_cv_output_total_force, f_cv_total_force);
+
+    init_feature(f_cv_subtract_applied_force, "subtract applied force from total force", f_type_user);
+    require_feature_self(f_cv_subtract_applied_force, f_cv_total_force);
+
+    init_feature(f_cv_lower_boundary, "lower boundary", f_type_user);
+    require_feature_self(f_cv_lower_boundary, f_cv_scalar);
+
+    init_feature(f_cv_upper_boundary, "upper boundary", f_type_user);
+    require_feature_self(f_cv_upper_boundary, f_cv_scalar);
+
+    init_feature(f_cv_grid, "grid", f_type_dynamic);
+    require_feature_self(f_cv_grid, f_cv_lower_boundary);
+    require_feature_self(f_cv_grid, f_cv_upper_boundary);
+
+    init_feature(f_cv_runave, "running average", f_type_user);
+
+    init_feature(f_cv_corrfunc, "correlation function", f_type_user);
+
+    init_feature(f_cv_scripted, "scripted", f_type_user);
+
+    init_feature(f_cv_custom_function, "custom function", f_type_user);
+    exclude_feature_self(f_cv_custom_function, f_cv_scripted);
+
+    init_feature(f_cv_periodic, "periodic", f_type_static);
+    require_feature_self(f_cv_periodic, f_cv_scalar);
+    init_feature(f_cv_scalar, "scalar", f_type_static);
+    init_feature(f_cv_linear, "linear", f_type_static);
+    init_feature(f_cv_homogeneous, "homogeneous", f_type_static);
+
+    // because total forces are obtained from the previous time step,
+    // we cannot (currently) have colvar values and total forces for the same timestep
+    init_feature(f_cv_multiple_ts, "multiple timestep colvar", f_type_static);
+    exclude_feature_self(f_cv_multiple_ts, f_cv_total_force_calc);
+
+    // check that everything is initialized
+    for (i = 0; i < colvardeps::f_cv_ntot; i++) {
+      if (is_not_set(i)) {
+        cvm::error("Uninitialized feature " + cvm::to_str(i) + " in " + description);
+      }
+    }
+  }
+
+  // Initialize feature_states for each instance
+  feature_states.reserve(f_cv_ntot);
+  for (i = 0; i < f_cv_ntot; i++) {
+    feature_states.push_back(feature_state(true, false));
+    // Most features are available, so we set them so
+    // and list exceptions below
+   }
+
+  feature_states[f_cv_fdiff_velocity].available =
+    cvm::main()->proxy->simulation_running();
+
+  return COLVARS_OK;
+}
+
+
+void colvar::setup()
+{
+  // loop over all components to update masses and charges of all groups
   for (size_t i = 0; i < cvcs.size(); i++) {
     for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
-      cvm::atom_group &atoms = *(cvcs[i]->atom_groups[ig]);
+      cvm::atom_group *atoms = cvcs[i]->atom_groups[ig];
-      atoms.setup();
+      atoms->setup();
-      atoms.reset_mass(name,i,ig);
+      atoms->print_properties(name, i, ig);
-      atoms.read_positions();
+      atoms->read_positions();
     }
   }
 }
 
 
+std::vector<std::vector<int> > colvar::get_atom_lists()
+{
+  std::vector<std::vector<int> > lists;
+  for (size_t i = 0; i < cvcs.size(); i++) {
+    std::vector<std::vector<int> > li = cvcs[i]->get_atom_lists();
+    lists.insert(lists.end(), li.begin(), li.end());
+  }
+  return lists;
+}
+
+
 colvar::~colvar()
 {
   // There is no need to call free_children_deps() here
   // because the children are cvcs and will be deleted
   // just below
 
-//   Clear references to this colvar's cvcs as children
+  // Clear references to this colvar's cvcs as children
-//   for dependency purposes
+  // for dependency purposes
   remove_all_children();
 
   for (std::vector<cvc *>::reverse_iterator ci = cvcs.rbegin();
@@ -1231,7 +1374,6 @@ int colvar::calc_cvc_gradients(int first_cvc, size_t num_cvcs)
 int colvar::collect_cvc_gradients()
 {
   size_t i;
-
   if (is_enabled(f_cv_collect_gradient)) {
     // Collect the atomic gradients inside colvar object
     for (unsigned int a = 0; a < atomic_gradients.size(); a++) {
@@ -1239,34 +1381,7 @@ int colvar::collect_cvc_gradients()
     }
     for (i = 0; i < cvcs.size(); i++) {
       if (!cvcs[i]->is_enabled()) continue;
-      // Coefficient: d(a * x^n) = a * n * x^(n-1) * dx
+      cvcs[i]->collect_gradients(atom_ids, atomic_gradients);
-      cvm::real coeff = (cvcs[i])->sup_coeff * cvm::real((cvcs[i])->sup_np) *
-        cvm::integer_power((cvcs[i])->value().real_value, (cvcs[i])->sup_np-1);
-
-      for (size_t j = 0; j < cvcs[i]->atom_groups.size(); j++) {
-
-        cvm::atom_group &ag = *(cvcs[i]->atom_groups[j]);
-
-        // If necessary, apply inverse rotation to get atomic
-        // gradient in the laboratory frame
-        if (ag.b_rotate) {
-          cvm::rotation const rot_inv = ag.rot.inverse();
-
-          for (size_t k = 0; k < ag.size(); k++) {
-            size_t a = std::lower_bound(atom_ids.begin(), atom_ids.end(),
-                                        ag[k].id) - atom_ids.begin();
-            atomic_gradients[a] += coeff * rot_inv.rotate(ag[k].grad);
-          }
-
-        } else {
-
-          for (size_t k = 0; k < ag.size(); k++) {
-            size_t a = std::lower_bound(atom_ids.begin(), atom_ids.end(),
-                                        ag[k].id) - atom_ids.begin();
-            atomic_gradients[a] += coeff * ag[k].grad;
-          }
-        }
-      }
     }
   }
   return COLVARS_OK;
@@ -1391,20 +1506,20 @@ int colvar::calc_colvar_properties()
     // initialize the restraint center in the first step to the value
     // just calculated from the cvcs
     if (cvm::step_relative() == 0 && !after_restart) {
-      xr = x;
+      x_ext = x;
-      vr.reset(); // (already 0; added for clarity)
+      v_ext.reset(); // (already 0; added for clarity)
     }
 
     // Special case of a repeated timestep (eg. multiple NAMD "run" statements)
     // revert values of the extended coordinate and velocity prior to latest integration
-    if (cvm::step_relative() == prev_timestep) {
+    if (cvm::proxy->simulation_running() && cvm::step_relative() == prev_timestep) {
-      xr = prev_xr;
+      x_ext = prev_x_ext;
-      vr = prev_vr;
+      v_ext = prev_v_ext;
     }
 
     // report the restraint center as "value"
-    x_reported = xr;
+    x_reported = x_ext;
-    v_reported = vr;
+    v_reported = v_ext;
     // the "total force" with the extended Lagrangian is
     // calculated in update_forces_energy() below
 
@@ -1458,6 +1573,8 @@ cvm::real colvar::update_forces_energy()
   // extended variable if there is one
 
   if (is_enabled(f_cv_extended_Lagrangian)) {
+    if (cvm::proxy->simulation_running()) {
+      // Only integrate the extended equations of motion in running MD simulations
       if (cvm::debug()) {
         cvm::log("Updating extended-Lagrangian degree of freedom.\n");
       }
@@ -1493,15 +1610,15 @@ cvm::real colvar::update_forces_energy()
       // 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);
-    f_ext += (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x);
+      f_ext += (-0.5 * ext_force_k) * this->dist2_lgrad(x_ext, x);
-    f      = (-0.5 * ext_force_k) * this->dist2_rgrad(xr, 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
       f *= cvm::real(time_step_factor);
 
       if (is_enabled(f_cv_subtract_applied_force)) {
         // Report a "system" force without the biases on this colvar
         // that is, just the spring force
-      ft_reported = (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x);
+        ft_reported = (-0.5 * ext_force_k) * this->dist2_lgrad(x_ext, x);
       } else {
         // The total force acting on the extended variable is f_ext
         // This will be used in the next timestep
@@ -1510,25 +1627,32 @@ cvm::real colvar::update_forces_energy()
 
       // backup in case we need to revert this integration timestep
       // if the same MD timestep is re-run
-    prev_xr = xr;
+      prev_x_ext = x_ext;
-    prev_vr = vr;
+      prev_v_ext = v_ext;
 
       // leapfrog: starting from x_i, f_i, v_(i-1/2)
-    vr  += (0.5 * dt) * f_ext / ext_mass;
+      v_ext  += (0.5 * dt) * f_ext / ext_mass;
       // Because of leapfrog, kinetic energy at time i is approximate
-    kinetic_energy = 0.5 * ext_mass * vr * vr;
+      kinetic_energy = 0.5 * ext_mass * v_ext * v_ext;
-    potential_energy = 0.5 * ext_force_k * this->dist2(xr, x);
+      potential_energy = 0.5 * ext_force_k * this->dist2(x_ext, x);
       // leap to v_(i+1/2)
       if (is_enabled(f_cv_Langevin)) {
-      vr -= dt * ext_gamma * vr;
+        v_ext -= dt * ext_gamma * v_ext;
         colvarvalue rnd(x);
         rnd.set_random();
-      vr += dt * ext_sigma * rnd / ext_mass;
+        v_ext += dt * ext_sigma * rnd / ext_mass;
       }
-    vr  += (0.5 * dt) * f_ext / ext_mass;
+      v_ext  += (0.5 * dt) * f_ext / ext_mass;
-    xr  += dt * vr;
+      x_ext  += dt * v_ext;
-    xr.apply_constraints();
+      x_ext.apply_constraints();
-    this->wrap(xr);
+      this->wrap(x_ext);
+    } else {
+      // If this is a postprocessing run (eg. in VMD), the extended DOF
+      // is equal to the actual coordinate
+      x_ext = x;
+    }
+    // Report extended value
+    x_reported = x_ext;
   }
 
   // Now adding the force on the actual colvar (for those biases that
@@ -1730,7 +1854,7 @@ int colvar::update_cvc_config(std::vector<std::string> const &confs)
 
 
 // ******************** METRIC FUNCTIONS ********************
-// Use the metrics defined by \link cvc \endlink objects
+// Use the metrics defined by \link colvar::cvc \endlink objects
 
 
 bool colvar::periodic_boundaries(colvarvalue const &lb, colvarvalue const &ub) const
@@ -1742,7 +1866,7 @@ bool colvar::periodic_boundaries(colvarvalue const &lb, colvarvalue const &ub) c
   }
 
   if (period > 0.0) {
-    if ( ((std::sqrt(this->dist2(lb, ub))) / this->width)
+    if ( ((cvm::sqrt(this->dist2(lb, ub))) / this->width)
          < 1.0E-10 ) {
       return true;
     }
@@ -1792,21 +1916,21 @@ colvarvalue colvar::dist2_rgrad(colvarvalue const &x1,
   }
 }
 
-void colvar::wrap(colvarvalue &x) const
+
+void colvar::wrap(colvarvalue &x_unwrapped) const
 {
-  if ( !is_enabled(f_cv_periodic) ) {
+  if (!is_enabled(f_cv_periodic)) {
     return;
   }
 
   if ( is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function) ) {
     // Scripted functions do their own wrapping, as cvcs might not be periodic
-    cvm::real shift = std::floor((x.real_value - wrap_center) / period + 0.5);
+    cvm::real shift = cvm::floor((x_unwrapped.real_value - wrap_center) /
-    x.real_value -= shift * period;
+                                 period + 0.5);
+    x_unwrapped.real_value -= shift * period;
   } else {
-    cvcs[0]->wrap(x);
+    cvcs[0]->wrap(x_unwrapped);
   }
-
-  return;
 }
 
 
@@ -1852,15 +1976,15 @@ std::istream & colvar::read_restart(std::istream &is)
 
   if (is_enabled(f_cv_extended_Lagrangian)) {
 
-    if ( !(get_keyval(conf, "extended_x", xr,
+    if ( !(get_keyval(conf, "extended_x", x_ext,
                       colvarvalue(x.type()), colvarparse::parse_silent)) &&
-         !(get_keyval(conf, "extended_v", vr,
+         !(get_keyval(conf, "extended_v", v_ext,
                       colvarvalue(x.type()), colvarparse::parse_silent)) ) {
       cvm::log("Error: restart file does not contain "
                "\"extended_x\" or \"extended_v\" for the colvar \""+
                name+"\", but you requested \"extendedLagrangian\".\n");
     }
-    x_reported = xr;
+    x_reported = x_ext;
   } else {
     x_reported = x;
   }
@@ -1875,7 +1999,7 @@ std::istream & colvar::read_restart(std::istream &is)
     }
 
     if (is_enabled(f_cv_extended_Lagrangian)) {
-      v_reported = vr;
+      v_reported = v_ext;
     } else {
       v_reported = v_fdiff;
     }
@@ -1901,8 +2025,8 @@ std::istream & colvar::read_traj(std::istream &is)
     }
 
     if (is_enabled(f_cv_extended_Lagrangian)) {
-      is >> xr;
+      is >> x_ext;
-      x_reported = xr;
+      x_reported = x_ext;
     } else {
       x_reported = x;
     }
@@ -1913,8 +2037,8 @@ std::istream & colvar::read_traj(std::istream &is)
     is >> v_fdiff;
 
     if (is_enabled(f_cv_extended_Lagrangian)) {
-      is >> vr;
+      is >> v_ext;
-      v_reported = vr;
+      v_reported = v_ext;
     } else {
       v_reported = v_fdiff;
     }
@@ -1955,11 +2079,11 @@ std::ostream & colvar::write_restart(std::ostream &os) {
     os << "  extended_x "
        << std::setprecision(cvm::cv_prec)
        << std::setw(cvm::cv_width)
-       << xr << "\n"
+       << x_ext << "\n"
        << "  extended_v "
        << std::setprecision(cvm::cv_prec)
        << std::setw(cvm::cv_width)
-       << vr << "\n";
+       << v_ext << "\n";
   }
 
   os << "}\n\n";
@@ -2190,6 +2314,7 @@ int colvar::calc_acf()
       acf_x_history_p = acf_x_history.begin();
       break;
 
+    case acf_notset:
     default:
       break;
     }
@@ -2222,6 +2347,7 @@ int colvar::calc_acf()
       history_incr(acf_x_history, acf_x_history_p);
       break;
 
+    case acf_notset:
     default:
       break;
     }
@@ -2257,7 +2383,7 @@ void colvar::calc_vel_acf(std::list<colvarvalue> &v_list,
 
 
 void colvar::calc_coor_acf(std::list<colvarvalue> &x_list,
-                            colvarvalue const      &x)
+                           colvarvalue const &x_now)
 {
   // same as above but for coordinates
   if (x_list.size() >= acf_length+acf_offset) {
@@ -2269,7 +2395,7 @@ void colvar::calc_coor_acf(std::list<colvarvalue> &x_list,
 
     *(acf_i++) += x.norm2();
 
-    colvarvalue::inner_opt(x, xs_i, x_list.end(), acf_i);
+    colvarvalue::inner_opt(x_now, xs_i, x_list.end(), acf_i);
 
     acf_nframes++;
   }
@@ -2277,7 +2403,7 @@ void colvar::calc_coor_acf(std::list<colvarvalue> &x_list,
 
 
 void colvar::calc_p2coor_acf(std::list<colvarvalue> &x_list,
-                             colvarvalue const      &x)
+                             colvarvalue const &x_now)
 {
   // same as above but with second order Legendre polynomial instead
   // of just the scalar product
@@ -2291,7 +2417,7 @@ void colvar::calc_p2coor_acf(std::list<colvarvalue> &x_list,
     // value of P2(0) = 1
     *(acf_i++) += 1.0;
 
-    colvarvalue::p2leg_opt(x, xs_i, x_list.end(), acf_i);
+    colvarvalue::p2leg_opt(x_now, xs_i, x_list.end(), acf_i);
 
     acf_nframes++;
   }
@@ -2316,6 +2442,9 @@ int colvar::write_acf(std::ostream &os)
   case acf_p2coor:
     os << "Coordinate (2nd Legendre poly)";
     break;
+  case acf_notset:
+  default:
+    break;
   }
 
   if (acf_colvar_name == name) {
@@ -2420,7 +2549,7 @@ int colvar::calc_runave()
                    << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
                    << runave << " "
                    << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
-                   << std::sqrt(runave_variance) << "\n";
+                   << cvm::sqrt(runave_variance) << "\n";
       }
 
       history_add_value(runave_length, *x_history_p, x);

lib/colvars/colvar.h

@@ -92,7 +92,7 @@ public:
   static std::vector<feature *> cv_features;
 
   /// \brief Implementation of the feature list accessor for colvar
-  virtual const std::vector<feature *> &features()
+  virtual const std::vector<feature *> &features() const
   {
     return cv_features;
   }
@@ -133,7 +133,7 @@ protected:
 
     Here:
     S(x(t)) = x
-    s(t)    = xr
+    s(t)    = x_ext
     DS = Ds = delta
   */
 
@@ -170,13 +170,13 @@ protected:
 
   // Options for extended_lagrangian
   /// Restraint center
-  colvarvalue xr;
+  colvarvalue x_ext;
   /// Previous value of the restraint center;
-  colvarvalue prev_xr;
+  colvarvalue prev_x_ext;
   /// Velocity of the restraint center
-  colvarvalue vr;
+  colvarvalue v_ext;
   /// Previous velocity of the restraint center
-  colvarvalue prev_vr;
+  colvarvalue prev_v_ext;
   /// Mass of the restraint center
   cvm::real ext_mass;
   /// Restraint force constant
@@ -273,6 +273,9 @@ public:
   /// Init output flags
   int init_output_flags(std::string const &conf);
 
+  /// \brief Initialize dependency tree
+  virtual int init_dependencies();
+
 private:
   /// Parse the CVC configuration for all components of a certain type
   template<typename def_class_name> int init_components_type(std::string const &conf,
@@ -373,7 +376,7 @@ protected:
   void update_active_cvc_square_norm();
 
   /// \brief Absolute timestep number when this colvar was last updated
-  int prev_timestep;
+  cvm::step_number prev_timestep;
 
 public:
 
@@ -383,32 +386,32 @@ public:
   /// \brief Return the number of CVC objects with an active flag (as set by update_cvc_flags)
   inline size_t num_active_cvcs() const { return n_active_cvcs; }
 
-  /// \brief Use the internal metrics (as from \link cvc
+  /// \brief Use the internal metrics (as from \link colvar::cvc
   /// \endlink objects) to calculate square distances and gradients
   ///
   /// Handles correctly symmetries and periodic boundary conditions
   cvm::real dist2(colvarvalue const &x1,
                   colvarvalue const &x2) const;
 
-  /// \brief Use the internal metrics (as from \link cvc
+  /// \brief Use the internal metrics (as from \link colvar::cvc
   /// \endlink objects) to calculate square distances and gradients
   ///
   /// Handles correctly symmetries and periodic boundary conditions
   colvarvalue dist2_lgrad(colvarvalue const &x1,
                           colvarvalue const &x2) const;
 
-  /// \brief Use the internal metrics (as from \link cvc
+  /// \brief Use the internal metrics (as from \link colvar::cvc
   /// \endlink objects) to calculate square distances and gradients
   ///
   /// Handles correctly symmetries and periodic boundary conditions
   colvarvalue dist2_rgrad(colvarvalue const &x1,
                           colvarvalue const &x2) const;
 
-  /// \brief Use the internal metrics (as from \link cvc
+  /// \brief Use the internal metrics (as from \link colvar::cvc
   /// \endlink objects) to wrap a value into a standard interval
   ///
   /// Handles correctly symmetries and periodic boundary conditions
-  void wrap(colvarvalue &x) const;
+  void wrap(colvarvalue &x_unwrapped) const;
 
 
   /// Read the analysis tasks
@@ -546,6 +549,7 @@ public:
   class polar_phi;
   class distance_inv;
   class distance_pairs;
+  class dipole_magnitude;
   class angle;
   class dipole_angle;
   class dihedral;
@@ -574,7 +578,7 @@ public:
 
 protected:
 
-  /// \brief Array of \link cvc \endlink objects
+  /// \brief Array of \link colvar::cvc \endlink objects
   std::vector<cvc *> cvcs;
 
   /// \brief Flags to enable or disable cvcs at next colvar evaluation
@@ -619,6 +623,9 @@ public:
   inline size_t n_components() const {
     return cvcs.size();
   }
+
+  /// \brief Get vector of vectors of atom IDs for all atom groups
+  virtual std::vector<std::vector<int> > get_atom_lists();
 };
 
 inline cvm::real const & colvar::force_constant() const
@@ -655,6 +662,8 @@ inline colvarvalue const & colvar::total_force() const
 
 inline void colvar::add_bias_force(colvarvalue const &force)
 {
+  check_enabled(f_cv_gradient,
+                std::string("applying a force to the variable \""+name+"\""));
   if (cvm::debug()) {
     cvm::log("Adding biasing force "+cvm::to_str(force)+" to colvar \""+name+"\".\n");
   }

lib/colvars/colvar_UIestimator.h

@@ -33,24 +33,24 @@ namespace UIestimator {
 
     public:
         n_matrix() {}
-        n_matrix(const std::vector<double> & lowerboundary,   // lowerboundary of x
+        n_matrix(const std::vector<double> & lowerboundary_input,   // lowerboundary of x
-            const std::vector<double> & upperboundary,   // upperboundary of
+            const std::vector<double> & upperboundary_input,   // upperboundary of
-            const std::vector<double> & width,           // width of x
+            const std::vector<double> & width_input,           // width of x
-            const int y_size) {          // size of y, for example, ysize=7, then when x=1, the distribution of y in [-2,4] is considered
+            const int y_size_input) {          // size of y, for example, ysize=7, then when x=1, the distribution of y in [-2,4] is considered
 
             int i;
 
-            this->lowerboundary = lowerboundary;
+            this->lowerboundary = lowerboundary_input;
-            this->upperboundary = upperboundary;
+            this->upperboundary = upperboundary_input;
-            this->width = width;
+            this->width = width_input;
-            this->dimension = lowerboundary.size();
+            this->dimension = lowerboundary_input.size();
-            this->y_size = y_size;     // keep in mind the internal (spare) matrix is stored in diagonal form
+            this->y_size = y_size_input;     // keep in mind the internal (spare) matrix is stored in diagonal form
-            this->y_total_size = int(std::pow(double(y_size), double(dimension)) + EPSILON);
+            this->y_total_size = int(cvm::pow(double(y_size_input), double(dimension)) + EPSILON);
 
             // the range of the matrix is [lowerboundary, upperboundary]
             x_total_size = 1;
             for (i = 0; i < dimension; i++) {
-                x_size.push_back(int((upperboundary[i] - lowerboundary[i]) / width[i] + EPSILON));
+                x_size.push_back(int((upperboundary_input[i] - lowerboundary_input[i]) / width_input[i] + EPSILON));
                 x_total_size *= x_size[i];
             }
 
@@ -89,9 +89,10 @@ namespace UIestimator {
 
         std::vector<int> temp;         // this vector is used in convert_x and convert_y to save computational resource
 
+        int convert_x(const std::vector<double> & x) {       // convert real x value to its interal index
+
             int i, j;
 
-        int convert_x(const std::vector<double> & x) {       // convert real x value to its interal index
             for (i = 0; i < dimension; i++) {
                 temp[i] = int((x[i] - lowerboundary[i]) / width[i] + EPSILON);
             }
@@ -121,7 +122,7 @@ namespace UIestimator {
             int index = 0;
             for (i = 0; i < dimension; i++) {
                 if (i + 1 < dimension)
-                    index += temp[i] * int(std::pow(double(y_size), double(dimension - i - 1)) + EPSILON);
+                    index += temp[i] * int(cvm::pow(double(y_size), double(dimension - i - 1)) + EPSILON);
                 else
                     index += temp[i];
             }
@@ -139,19 +140,19 @@ namespace UIestimator {
 
     public:
         n_vector() {}
-        n_vector(const std::vector<double> & lowerboundary,   // lowerboundary of x
+        n_vector(const std::vector<double> & lowerboundary_input,   // lowerboundary of x
-            const std::vector<double> & upperboundary,   // upperboundary of
+            const std::vector<double> & upperboundary_input,   // upperboundary of
-            const std::vector<double> & width,                // width of x
+            const std::vector<double> & width_input,                // width of x
-            const int y_size,           // size of y, for example, ysize=7, then when x=1, the distribution of y in [-2,4] is considered
+            const int y_size_input,           // size of y, for example, ysize=7, then when x=1, the distribution of y in [-2,4] is considered
             const T & default_value) {         //   the default value of T
 
-            this->width = width;
+            this->width = width_input;
-            this->dimension = lowerboundary.size();
+            this->dimension = lowerboundary_input.size();
 
             x_total_size = 1;
             for (int i = 0; i < dimension; i++) {
-                this->lowerboundary.push_back(lowerboundary[i] - (y_size - 1) / 2 * width[i] - EPSILON);
+                this->lowerboundary.push_back(lowerboundary_input[i] - (y_size_input - 1) / 2 * width_input[i] - EPSILON);
-                this->upperboundary.push_back(upperboundary[i] + (y_size - 1) / 2 * width[i] + EPSILON);
+                this->upperboundary.push_back(upperboundary_input[i] + (y_size_input - 1) / 2 * width_input[i] + EPSILON);
 
                 x_size.push_back(int((this->upperboundary[i] - this->lowerboundary[i]) / this->width[i] + EPSILON));
                 x_total_size *= x_size[i];
@@ -215,26 +216,26 @@ namespace UIestimator {
         UIestimator() {}
 
         //called when (re)start an eabf simulation
-        UIestimator(const std::vector<double> & lowerboundary,
+        UIestimator(const std::vector<double> & lowerboundary_input,
-            const std::vector<double> & upperboundary,
+            const std::vector<double> & upperboundary_input,
-            const std::vector<double> & width,
+            const std::vector<double> & width_input,
-            const std::vector<double> & krestr,                // force constant in eABF
+            const std::vector<double> & krestr_input,                // force constant in eABF
-            const std::string & output_filename,              // the prefix of output files
+            const std::string & output_filename_input,              // the prefix of output files
-            const int output_freq,
+            const int output_freq_input,
-            const bool restart,                              // whether restart from a .count and a .grad file
+            const bool restart_input,                              // whether restart from a .count and a .grad file
-            const std::vector<std::string> & input_filename,   // the prefixes of input files
+            const std::vector<std::string> & input_filename_input,   // the prefixes of input files
-            const double temperature) {
+            const double temperature_input) {
 
             // initialize variables
-            this->lowerboundary = lowerboundary;
+            this->lowerboundary = lowerboundary_input;
-            this->upperboundary = upperboundary;
+            this->upperboundary = upperboundary_input;
-            this->width = width;
+            this->width = width_input;
-            this->krestr = krestr;
+            this->krestr = krestr_input;
-            this->output_filename = output_filename;
+            this->output_filename = output_filename_input;
-            this->output_freq = output_freq;
+            this->output_freq = output_freq_input;
-            this->restart = restart;
+            this->restart = restart_input;
-            this->input_filename = input_filename;
+            this->input_filename = input_filename_input;
-            this->temperature = temperature;
+            this->temperature = temperature_input;
 
             int i, j;
 
@@ -300,7 +301,7 @@ namespace UIestimator {
         ~UIestimator() {}
 
         // called from MD engine every step
-        bool update(const int step, std::vector<double> x, std::vector<double> y) {
+        bool update(cvm::step_number step, std::vector<double> x, std::vector<double> y) {
 
             int i;
 
@@ -431,7 +432,7 @@ namespace UIestimator {
                     loop_flag_y[k] = loop_flag_x[k] - HALF_Y_SIZE * width[k];
                 }
 
-                int j = 0;
+                j = 0;
                 while (j >= 0) {
                     norm += distribution_x_y.get_value(loop_flag_x, loop_flag_y);
                     for (k = 0; k < dimension; k++) {
@@ -672,7 +673,7 @@ namespace UIestimator {
         }
 
         // read input files
-        void read_inputfiles(const std::vector<std::string> input_filename)
+        void read_inputfiles(const std::vector<std::string> filename)
         {
             char sharp;
             double nothing;
@@ -683,11 +684,11 @@ namespace UIestimator {
             std::vector<double> position_temp(dimension, 0);
             std::vector<double> grad_temp(dimension, 0);
             int count_temp = 0;
-            for (i = 0; i < int(input_filename.size()); i++) {
+            for (i = 0; i < int(filename.size()); i++) {
                 int size = 1 , size_temp = 0;
 
-                std::string count_filename = input_filename[i] + ".UI.count";
+                std::string count_filename = filename[i] + ".UI.count";
-                std::string grad_filename = input_filename[i] + ".UI.grad";
+                std::string grad_filename = filename[i] + ".UI.grad";
 
                 std::ifstream count_file(count_filename.c_str(), std::ios::in);
                 std::ifstream grad_file(grad_filename.c_str(), std::ios::in);

lib/colvars/colvaratoms.cpp

@@ -22,7 +22,7 @@ cvm::atom::atom()
   index = -1;
   id = -1;
   mass = 1.0;
-  charge = 1.0;
+  charge = 0.0;
   reset_data();
 }
 
@@ -107,6 +107,8 @@ cvm::atom_group::~atom_group()
     delete fitting_group;
     fitting_group = NULL;
   }
+
+  cvm::main()->unregister_named_atom_group(this);
 }
 
 
@@ -183,10 +185,7 @@ int cvm::atom_group::init()
   // These may be overwritten by parse(), if a name is provided
 
   atoms.clear();
-
+  init_dependencies();
-  // TODO: check with proxy whether atom forces etc are available
-  init_ag_requires();
-
   index = -1;
 
   b_dummy = false;
@@ -207,8 +206,67 @@ int cvm::atom_group::init()
 }
 
 
+int cvm::atom_group::init_dependencies() {
+  size_t i;
+  // Initialize static array once and for all
+  if (features().size() == 0) {
+    for (i = 0; i < f_ag_ntot; i++) {
+      modify_features().push_back(new feature);
+    }
+
+    init_feature(f_ag_active, "active", f_type_dynamic);
+    init_feature(f_ag_center, "translational fit", f_type_static);
+    init_feature(f_ag_rotate, "rotational fit", f_type_static);
+    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);
+    require_feature_self(f_ag_fit_gradients, f_ag_explicit_gradient);
+
+    init_feature(f_ag_atom_forces, "atomic forces", f_type_dynamic);
+
+    // parallel calculation implies that we have at least a scalable center of mass,
+    // but f_ag_scalable is kept as a separate feature to deal with future dependencies
+    init_feature(f_ag_scalable, "scalable group calculation", f_type_static);
+    init_feature(f_ag_scalable_com, "scalable group center of mass calculation", f_type_static);
+    require_feature_self(f_ag_scalable, f_ag_scalable_com);
+
+    // check that everything is initialized
+    for (i = 0; i < colvardeps::f_ag_ntot; i++) {
+      if (is_not_set(i)) {
+        cvm::error("Uninitialized feature " + cvm::to_str(i) + " in " + description);
+      }
+    }
+  }
+
+  // Initialize feature_states for each instance
+  // default as unavailable, not enabled
+  feature_states.reserve(f_ag_ntot);
+  for (i = 0; i < colvardeps::f_ag_ntot; i++) {
+    feature_states.push_back(feature_state(false, false));
+  }
+
+  // Features that are implemented (or not) by all atom groups
+  feature_states[f_ag_active].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)
+  feature_states[f_ag_scalable].available = true;
+  feature_states[f_ag_fit_gradients].available = true;
+  feature_states[f_ag_fitting_group].available = true;
+  feature_states[f_ag_explicit_gradient].available = true;
+
+  return COLVARS_OK;
+}
+
+
 int cvm::atom_group::setup()
 {
+  if (atoms_ids.size() == 0) {
+    atoms_ids.reserve(atoms.size());
+    for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
+      atoms_ids.push_back(ai->id);
+    }
+  }
   for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
     ai->update_mass();
     ai->update_charge();
@@ -237,15 +295,6 @@ void cvm::atom_group::update_total_mass()
 }
 
 
-void cvm::atom_group::reset_mass(std::string &name, int i, int j)
-{
-  update_total_mass();
-  cvm::log("Re-initialized atom group "+name+":"+cvm::to_str(i)+"/"+
-           cvm::to_str(j)+". "+ cvm::to_str(atoms_ids.size())+
-           " atoms: total mass = "+cvm::to_str(total_mass)+".\n");
-}
-
-
 void cvm::atom_group::update_total_charge()
 {
   if (b_dummy) {
@@ -264,6 +313,19 @@ void cvm::atom_group::update_total_charge()
 }
 
 
+void cvm::atom_group::print_properties(std::string const &colvar_name,
+                                       int i, int j)
+{
+  if (cvm::proxy->updated_masses() && cvm::proxy->updated_charges()) {
+    cvm::log("Re-initialized atom group for variable \""+colvar_name+"\":"+
+             cvm::to_str(i)+"/"+
+             cvm::to_str(j)+". "+ cvm::to_str(atoms_ids.size())+
+             " atoms: total mass = "+cvm::to_str(total_mass)+
+             ", total charge = "+cvm::to_str(total_charge)+".\n");
+  }
+}
+
+
 int cvm::atom_group::parse(std::string const &group_conf)
 {
   cvm::log("Initializing atom group \""+key+"\".\n");
@@ -450,10 +512,21 @@ int cvm::atom_group::parse(std::string const &group_conf)
   if (cvm::debug())
     cvm::log("Done initializing atom group \""+key+"\".\n");
 
-  cvm::log("Atom group \""+key+"\" defined, "+
+  {
-            cvm::to_str(atoms_ids.size())+" atoms initialized: total mass = "+
+    std::string init_msg;
-            cvm::to_str(total_mass)+", total charge = "+
+    init_msg.append("Atom group \""+key+"\" defined with "+
-            cvm::to_str(total_charge)+".\n");
+                    cvm::to_str(atoms_ids.size())+" atoms requested");
+    if ((cvm::proxy)->updated_masses()) {
+      init_msg.append(": total mass = "+
+                      cvm::to_str(total_mass));
+      if ((cvm::proxy)->updated_charges()) {
+        init_msg.append(", total charge = "+
+                        cvm::to_str(total_charge));
+      }
+    }
+    init_msg.append(".\n");
+    cvm::log(init_msg);
+  }
 
   if (b_print_atom_ids) {
     cvm::log("Internal definition of the atom group:\n");
@@ -464,7 +537,7 @@ int cvm::atom_group::parse(std::string const &group_conf)
 }
 
 
-int cvm::atom_group::add_atoms_of_group(atom_group const * ag)
+int cvm::atom_group::add_atoms_of_group(atom_group const *ag)
 {
   std::vector<int> const &source_ids = ag->atoms_ids;
 
@@ -696,6 +769,7 @@ int cvm::atom_group::parse_fitting_options(std::string const &group_conf)
           return INPUT_ERROR;
         }
       }
+      enable(f_ag_fitting_group);
     }
 
     atom_group *group_for_fit = fitting_group ? fitting_group : this;
@@ -800,24 +874,24 @@ int cvm::atom_group::create_sorted_ids()
 
   // Sort the internal IDs
   std::list<int> sorted_atoms_ids_list;
-  for (size_t i = 0; i < this->size(); i++) {
+  for (size_t i = 0; i < atoms_ids.size(); i++) {
     sorted_atoms_ids_list.push_back(atoms_ids[i]);
   }
   sorted_atoms_ids_list.sort();
   sorted_atoms_ids_list.unique();
-  if (sorted_atoms_ids_list.size() != this->size()) {
+  if (sorted_atoms_ids_list.size() != atoms_ids.size()) {
     return cvm::error("Error: duplicate atom IDs in atom group? (found " +
                       cvm::to_str(sorted_atoms_ids_list.size()) +
                       " unique atom IDs instead of " +
-                      cvm::to_str(this->size()) + ").\n", BUG_ERROR);
+                      cvm::to_str(atoms_ids.size()) + ").\n", BUG_ERROR);
   }
 
   // Compute map between sorted and unsorted elements
-  sorted_atoms_ids.resize(this->size());
+  sorted_atoms_ids.resize(atoms_ids.size());
-  sorted_atoms_ids_map.resize(this->size());
+  sorted_atoms_ids_map.resize(atoms_ids.size());
   std::list<int>::iterator lsii = sorted_atoms_ids_list.begin();
   size_t ii = 0;
-  for ( ; ii < this->size(); lsii++, ii++) {
+  for ( ; ii < atoms_ids.size(); lsii++, ii++) {
     sorted_atoms_ids[ii] = *lsii;
     size_t const pos = std::find(atoms_ids.begin(), atoms_ids.end(), *lsii) -
       atoms_ids.begin();
@@ -1038,15 +1112,15 @@ int cvm::atom_group::calc_center_of_mass()
 }
 
 
-int cvm::atom_group::calc_dipole(cvm::atom_pos const &com)
+int cvm::atom_group::calc_dipole(cvm::atom_pos const &dipole_center)
 {
   if (b_dummy) {
-    cvm::error("Error: trying to compute the dipole of an empty group.\n", INPUT_ERROR);
+    return cvm::error("Error: trying to compute the dipole "
-    return COLVARS_ERROR;
+                      "of a dummy group.\n", INPUT_ERROR);
   }
   dip.reset();
   for (cvm::atom_const_iter ai = this->begin(); ai != this->end(); ai++) {
-    dip += ai->charge * (ai->pos - com);
+    dip += ai->charge * (ai->pos - dipole_center);
   }
   return COLVARS_OK;
 }
@@ -1056,14 +1130,13 @@ void cvm::atom_group::set_weighted_gradient(cvm::rvector const &grad)
 {
   if (b_dummy) return;
 
-  if (is_enabled(f_ag_scalable)) {
   scalar_com_gradient = grad;
-    return;
-  }
 
+  if (!is_enabled(f_ag_scalable)) {
     for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
       ai->grad = (ai->mass/total_mass) * grad;
     }
+  }
 }
 
 

lib/colvars/colvaratoms.h

@@ -17,7 +17,7 @@
 
 
 /// \brief Stores numeric id, mass and all mutable data for an atom,
-/// mostly used by a \link cvc \endlink
+/// mostly used by a \link colvar::cvc \endlink
 ///
 /// This class may be used to keep atomic data such as id, mass,
 /// position and collective variable derivatives) altogether.
@@ -63,7 +63,7 @@ public:
   /// from the \link colvarvalue \endlink class), which is also the
   /// most frequent case. For more complex types of \link
   /// colvarvalue \endlink objects, atomic gradients should be
-  /// defined within the specific \link cvc \endlink
+  /// defined within the specific \link colvar::cvc \endlink
   /// implementation
   cvm::rvector   grad;
 
@@ -100,13 +100,19 @@ public:
   /// Get the latest value of the mass
   inline void update_mass()
   {
-    mass = (cvm::proxy)->get_atom_mass(index);
+    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()
   {
-    charge = (cvm::proxy)->get_atom_charge(index);
+    colvarproxy *p = cvm::proxy;
+    if (p->updated_charges()) {
+      charge = p->get_atom_charge(index);
+    }
   }
 
   /// Get the current position
@@ -145,7 +151,7 @@ public:
 
 
 /// \brief Group of \link atom \endlink objects, mostly used by a
-/// \link cvc \endlink object to gather all atomic data
+/// \link colvar::cvc \endlink object to gather all atomic data
 class colvarmodule::atom_group
   : public colvarparse, public colvardeps
 {
@@ -174,6 +180,9 @@ public:
   /// \brief Set default values for common flags
   int init();
 
+  /// \brief Initialize dependency tree
+  virtual int init_dependencies();
+
   /// \brief Update data required to calculate cvc's
   int setup();
 
@@ -198,16 +207,16 @@ public:
   /// \brief Remove an atom object from this group
   int remove_atom(cvm::atom_iter ai);
 
-  /// \brief Re-initialize the total mass of a group.
+  /// \brief Print the updated the total mass and charge of a group.
   /// This is needed in case the hosting MD code has an option to
   /// change atom masses after their initialization.
-  void reset_mass(std::string &name, int i, int j);
+  void print_properties(std::string const &colvar_name, int i, int j);
 
   /// \brief Implementation of the feature list for atom group
   static std::vector<feature *> ag_features;
 
   /// \brief Implementation of the feature list accessor for atom group
-  virtual const std::vector<feature *> &features()
+  virtual const std::vector<feature *> &features() const
   {
     return ag_features;
   }
@@ -347,10 +356,14 @@ public:
 
   /// Total mass of the atom group
   cvm::real total_mass;
+
+  /// Update the total mass of the atom group
   void update_total_mass();
 
   /// Total charge of the atom group
   cvm::real total_charge;
+
+  /// Update the total mass of the group
   void update_total_charge();
 
   /// \brief Don't apply any force on this group (use its coordinates
@@ -423,20 +436,32 @@ public:
   /// \brief Calculate the center of mass of the atomic positions, assuming that
   /// they are already pbc-wrapped
   int calc_center_of_mass();
+
 private:
+
   /// \brief Center of mass
   cvm::atom_pos com;
+
   /// \brief The derivative of a scalar variable with respect to the COM
   // TODO for scalable calculations of more complex variables (e.g. rotation),
   // use a colvarvalue of vectors to hold the entire derivative
   cvm::rvector scalar_com_gradient;
+
 public:
-  /// \brief Return the center of mass of the atomic positions
+
+  /// \brief Return the center of mass (COM) of the atomic positions
   inline cvm::atom_pos center_of_mass() const
   {
     return com;
   }
 
+  /// \brief Return previously gradient of scalar variable with respect to the
+  /// COM
+  inline cvm::rvector center_of_mass_scalar_gradient() const
+  {
+    return scalar_com_gradient;
+  }
+
   /// \brief Return a copy of the current atom positions, shifted by a constant vector
   std::vector<cvm::atom_pos> positions_shifted(cvm::rvector const &shift) const;
 
@@ -444,10 +469,15 @@ public:
   std::vector<cvm::rvector> velocities() const;
 
   ///\brief Calculate the dipole of the atom group around the specified center
-  int calc_dipole(cvm::atom_pos const &com);
+  int calc_dipole(cvm::atom_pos const &dipole_center);
+
 private:
+
+  /// Dipole moment of the atom group
   cvm::rvector dip;
+
 public:
+
   ///\brief Return the (previously calculated) dipole of the atom group
   inline cvm::rvector dipole() const
   {

lib/colvars/colvarbias.cpp

@@ -17,15 +17,14 @@
 colvarbias::colvarbias(char const *key)
   : bias_type(to_lower_cppstr(key))
 {
-  init_cvb_requires();
+  description = "uninitialized " + cvm::to_str(key) + " bias";
-
+  init_dependencies();
   rank = 1;
 
   has_data = false;
   b_output_energy = false;
   reset();
-  state_file_step = 0;
+  state_file_step = 0L;
-  description = "uninitialized " + cvm::to_str(key) + " bias";
 }
 
 
@@ -76,6 +75,7 @@ int colvarbias::init(std::string const &conf)
       cvm::error("Error: no collective variables specified.\n", INPUT_ERROR);
       return INPUT_ERROR;
     }
+
   } else {
     cvm::log("Reinitializing bias \""+name+"\".\n");
   }
@@ -98,6 +98,70 @@ int colvarbias::init(std::string const &conf)
 }
 
 
+int colvarbias::init_dependencies() {
+  int i;
+  if (features().size() == 0) {
+    for (i = 0; i < f_cvb_ntot; i++) {
+      modify_features().push_back(new feature);
+    }
+
+    init_feature(f_cvb_active, "active", f_type_dynamic);
+    require_feature_children(f_cvb_active, f_cv_active);
+
+    init_feature(f_cvb_awake, "awake", f_type_static);
+    require_feature_self(f_cvb_awake, f_cvb_active);
+
+    init_feature(f_cvb_apply_force, "apply force", f_type_user);
+    require_feature_children(f_cvb_apply_force, f_cv_gradient);
+
+    init_feature(f_cvb_get_total_force, "obtain total force", f_type_dynamic);
+    require_feature_children(f_cvb_get_total_force, f_cv_total_force);
+
+    init_feature(f_cvb_output_acc_work, "output accumulated work", f_type_user);
+    require_feature_self(f_cvb_output_acc_work, f_cvb_apply_force);
+
+    init_feature(f_cvb_history_dependent, "history-dependent", f_type_static);
+
+    init_feature(f_cvb_time_dependent, "time-dependent", f_type_static);
+
+    init_feature(f_cvb_scalar_variables, "require scalar variables", f_type_static);
+    require_feature_children(f_cvb_scalar_variables, f_cv_scalar);
+
+    init_feature(f_cvb_calc_pmf, "calculate a PMF", f_type_static);
+
+    init_feature(f_cvb_calc_ti_samples, "calculate TI samples", f_type_dynamic);
+    require_feature_self(f_cvb_calc_ti_samples, f_cvb_get_total_force);
+    require_feature_children(f_cvb_calc_ti_samples, f_cv_grid);
+
+    init_feature(f_cvb_write_ti_samples, "write TI samples ", f_type_user);
+    require_feature_self(f_cvb_write_ti_samples, f_cvb_calc_ti_samples);
+
+    init_feature(f_cvb_write_ti_pmf, "write TI PMF", f_type_user);
+    require_feature_self(f_cvb_write_ti_pmf, f_cvb_calc_ti_samples);
+
+    // check that everything is initialized
+    for (i = 0; i < colvardeps::f_cvb_ntot; i++) {
+      if (is_not_set(i)) {
+        cvm::error("Uninitialized feature " + cvm::to_str(i) + " in " + description);
+      }
+    }
+  }
+
+  // Initialize feature_states for each instance
+  feature_states.reserve(f_cvb_ntot);
+  for (i = 0; i < f_cvb_ntot; i++) {
+    feature_states.push_back(feature_state(true, false));
+    // Most features are available, so we set them so
+    // and list exceptions below
+  }
+
+  // only compute TI samples when deriving from colvarbias_ti
+  feature_states[f_cvb_calc_ti_samples].available = false;
+
+  return COLVARS_OK;
+}
+
+
 int colvarbias::reset()
 {
   bias_energy = 0.0;
@@ -217,6 +281,9 @@ int colvarbias::update()
 
 void colvarbias::communicate_forces()
 {
+  if (! is_enabled(f_cvb_apply_force)) {
+    return;
+  }
   size_t i = 0;
   for (i = 0; i < num_variables(); i++) {
     if (cvm::debug()) {
@@ -345,7 +412,8 @@ std::istream & colvarbias::read_state(std::istream &is)
        (set_state_params(conf) != COLVARS_OK) ) {
     cvm::error("Error: in reading state configuration for \""+bias_type+"\" bias \""+
                this->name+"\" at position "+
-               cvm::to_str(is.tellg())+" in stream.\n", INPUT_ERROR);
+               cvm::to_str(static_cast<size_t>(is.tellg()))+
+               " in stream.\n", INPUT_ERROR);
     is.clear();
     is.seekg(start_pos, std::ios::beg);
     is.setstate(std::ios::failbit);
@@ -355,7 +423,8 @@ std::istream & colvarbias::read_state(std::istream &is)
   if (!read_state_data(is)) {
     cvm::error("Error: in reading state data for \""+bias_type+"\" bias \""+
                this->name+"\" at position "+
-               cvm::to_str(is.tellg())+" in stream.\n", INPUT_ERROR);
+               cvm::to_str(static_cast<size_t>(is.tellg()))+
+               " in stream.\n", INPUT_ERROR);
     is.clear();
     is.seekg(start_pos, std::ios::beg);
     is.setstate(std::ios::failbit);
@@ -365,7 +434,8 @@ std::istream & colvarbias::read_state(std::istream &is)
   if (brace != "}") {
     cvm::error("Error: corrupt restart information for \""+bias_type+"\" bias \""+
                this->name+"\": no matching brace at position "+
-               cvm::to_str(is.tellg())+" in stream.\n");
+               cvm::to_str(static_cast<size_t>(is.tellg()))+
+               " in stream.\n");
     is.setstate(std::ios::failbit);
   }
 
@@ -381,7 +451,8 @@ std::istream & colvarbias::read_state_data_key(std::istream &is, char const *key
        !(key_in == to_lower_cppstr(std::string(key))) ) {
     cvm::error("Error: in reading restart configuration for "+
                bias_type+" bias \""+this->name+"\" at position "+
-               cvm::to_str(is.tellg())+" in stream.\n", INPUT_ERROR);
+               cvm::to_str(static_cast<size_t>(is.tellg()))+
+               " in stream.\n", INPUT_ERROR);
     is.clear();
     is.seekg(start_pos, std::ios::beg);
     is.setstate(std::ios::failbit);
@@ -640,7 +711,7 @@ int colvarbias_ti::write_output_files()
       cvm::proxy->close_output_stream(ti_count_file_name);
     }
 
-    std::string const ti_grad_file_name(ti_output_prefix+".ti.grad");
+    std::string const ti_grad_file_name(ti_output_prefix+".ti.force");
     os = cvm::proxy->output_stream(ti_grad_file_name);
     if (os) {
       ti_avg_forces->write_multicol(*os);

lib/colvars/colvarbias.h

@@ -96,6 +96,9 @@ public:
   /// \brief Parse config string and (re)initialize
   virtual int init(std::string const &conf);
 
+  /// \brief Initialize dependency tree
+  virtual int init_dependencies();
+
   /// \brief Set to zero all mutable data
   virtual int reset();
 
@@ -181,7 +184,7 @@ public:
   static std::vector<feature *> cvb_features;
 
   /// \brief Implementation of the feature list accessor for colvarbias
-  virtual const std::vector<feature *> &features()
+  virtual const std::vector<feature *> &features() const
   {
     return cvb_features;
   }
@@ -220,7 +223,7 @@ protected:
   bool                     has_data;
 
   /// \brief Step number read from the last state file
-  size_t                   state_file_step;
+  cvm::step_number         state_file_step;
 
 };
 

lib/colvars/colvarbias_abf.cpp

@@ -191,10 +191,10 @@ int colvarbias_abf::init(std::string const &conf)
       // Projected ABF
       get_keyval(conf, "pABFintegrateFreq", pabf_freq, 0);
       // Parameters for integrating initial (and final) gradient data
-      get_keyval(conf, "integrateInitSteps", integrate_initial_steps, 1e4);
+      get_keyval(conf, "integrateInitMaxIterations", integrate_initial_iterations, 1e4);
       get_keyval(conf, "integrateInitTol", integrate_initial_tol, 1e-6);
       // for updating the integrated PMF on the fly
-      get_keyval(conf, "integrateSteps", integrate_steps, 100);
+      get_keyval(conf, "integrateMaxIterations", integrate_iterations, 100);
       get_keyval(conf, "integrateTol", integrate_tol, 1e-4);
     }
   } else {
@@ -366,10 +366,10 @@ int colvarbias_abf::update()
     if ( b_integrate ) {
       if ( pabf_freq && cvm::step_relative() % pabf_freq == 0 ) {
         cvm::real err;
-        int iter = pmf->integrate(integrate_steps, integrate_tol, err);
+        int iter = pmf->integrate(integrate_iterations, integrate_tol, err);
-        if ( iter == integrate_steps ) {
+        if ( iter == integrate_iterations ) {
           cvm::log("Warning: PMF integration did not converge to " + cvm::to_str(integrate_tol)
-            + " in " + cvm::to_str(integrate_steps)
+            + " in " + cvm::to_str(integrate_iterations)
             + " steps. Residual error: " +  cvm::to_str(err));
         }
         pmf->set_zero_minimum(); // TODO: do this only when necessary
@@ -597,7 +597,7 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
   if (b_integrate) {
     // Do numerical integration (to high precision) and output a PMF
     cvm::real err;
-    pmf->integrate(integrate_initial_steps, integrate_initial_tol, err);
+    pmf->integrate(integrate_initial_iterations, integrate_initial_tol, err);
     pmf->set_zero_minimum();
 
     std::string  pmf_out_name = prefix + ".pmf";
@@ -661,7 +661,7 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
       // Do numerical integration (to high precision) and output a PMF
       cvm::real err;
       czar_pmf->set_div();
-      czar_pmf->integrate(integrate_initial_steps, integrate_initial_tol, err);
+      czar_pmf->integrate(integrate_initial_iterations, integrate_initial_tol, err);
       czar_pmf->set_zero_minimum();
 
       std::string  czar_pmf_out_name = prefix + ".czar.pmf";

lib/colvars/colvarbias_abf.h

@@ -27,9 +27,13 @@ class colvarbias_abf : public colvarbias {
 
 public:
 
+  /// Constructor for ABF bias
   colvarbias_abf(char const *key);
+  /// Initializer for ABF bias
   virtual int init(std::string const &conf);
+  /// Default destructor for ABF bias
   virtual ~colvarbias_abf();
+  /// Per-timestep update of ABF bias
   virtual int update();
 
 private:
@@ -40,11 +44,17 @@ private:
   /// Base filename(s) for reading previous gradient data (replaces data from restart file)
   std::vector<std::string> input_prefix;
 
+  /// Adapt the bias at each time step (as opposed to keeping it constant)?
   bool    update_bias;
+  /// Use normalized definition of PMF for distance functions? (flat at long distances)
+  /// by including the Jacobian term separately of the recorded PMF
   bool    hide_Jacobian;
+  /// Integrate gradients into a PMF on output
   bool    b_integrate;
 
+  /// Number of samples per bin before applying the full biasing force
   size_t  full_samples;
+  /// Number of samples per bin before applying a scaled-down biasing force
   size_t  min_samples;
   /// frequency for updating output files
   int     output_freq;
@@ -52,6 +62,7 @@ private:
   bool    b_history_files;
   /// Write CZAR output file for stratified eABF (.zgrad)
   bool    b_czar_window_file;
+  /// Number of timesteps between recording data in history files (if non-zero)
   size_t  history_freq;
   /// Umbrella Integration estimator of free energy from eABF
   UIestimator::UIestimator eabf_UI;
@@ -63,25 +74,30 @@ private:
   /// Frequency for updating pABF PMF (if zero, pABF is not used)
   int   pabf_freq;
   /// Max number of CG iterations for integrating PMF at startup and for file output
-  int       integrate_initial_steps;
+  int       integrate_initial_iterations;
   /// Tolerance for integrating PMF at startup and for file output
   cvm::real integrate_initial_tol;
   /// Max number of CG iterations for integrating PMF at on-the-fly pABF updates
-  int       integrate_steps;
+  int       integrate_iterations;
   /// Tolerance for integrating PMF at on-the-fly pABF updates
   cvm::real integrate_tol;
 
-  /// Cap the biasing force to be applied?
+  /// Cap the biasing force to be applied? (option maxForce)
   bool                    cap_force;
+  /// Maximum force to be applied
   std::vector<cvm::real>  max_force;
 
-  // Frequency for updating 2D gradients
-  int integrate_freq;
-
   // Internal data and methods
 
-  std::vector<int>  bin, force_bin, z_bin;
+  /// Current bin in sample grid
-  gradient_t    system_force, applied_force;
+  std::vector<int>  bin;
+  /// Current bin in force grid
+  std::vector<int> force_bin;
+  /// Cuurent bin in "actual" coordinate, when running extended Lagrangian dynamics
+  std::vector<int> z_bin;
+
+  /// Measured instantaneous system force
+  gradient_t system_force;
 
   /// n-dim grid of free energy gradients
   colvar_grid_gradient  *gradients;
@@ -118,7 +134,7 @@ private:
   // shared ABF
   bool    shared_on;
   size_t  shared_freq;
-  int     shared_last_step;
+  cvm::step_number shared_last_step;
   // Share between replicas -- may be called independently of update
   virtual int replica_share();
 

lib/colvars/colvarbias_histogram.h

@@ -37,7 +37,7 @@ protected:
   std::string out_name, out_name_dx;
   size_t output_freq;
 
-  /// If one or more of the variables are \link type_vector \endlink, treat them as arrays of this length
+  /// 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;
   /// If colvar_array_size is larger than 1, weigh each one by this number before accumulating the histogram
   std::vector<cvm::real> weights;

lib/colvars/colvarbias_meta.cpp

@@ -11,7 +11,6 @@
 #include <sstream>
 #include <fstream>
 #include <iomanip>
-#include <cmath>
 #include <algorithm>
 
 // used to set the absolute path of a replica file
@@ -39,6 +38,8 @@ colvarbias_meta::colvarbias_meta(char const *key)
   new_hills_begin = hills.end();
   hills_traj_os = NULL;
   replica_hills_os = NULL;
+
+  ebmeta_equil_steps = 0L;
 }
 
 
@@ -61,7 +62,7 @@ int colvarbias_meta::init(std::string const &conf)
     enable(f_cvb_history_dependent);
   }
 
-  get_keyval(conf, "hillWidth", hill_width, std::sqrt(2.0 * PI) / 2.0);
+  get_keyval(conf, "hillWidth", hill_width, cvm::sqrt(2.0 * PI) / 2.0);
   cvm::log("Half-widths of the Gaussian hills (sigma's):\n");
   for (size_t i = 0; i < num_variables(); i++) {
     cvm::log(variables(i)->name+std::string(": ")+
@@ -201,6 +202,7 @@ int colvarbias_meta::init_ebmeta_params(std::string const &conf)
     }
     target_dist = new colvar_grid_scalar();
     target_dist->init_from_colvars(colvars);
+    std::string target_dist_file;
     get_keyval(conf, "targetdistfile", target_dist_file);
     std::ifstream targetdiststream(target_dist_file.c_str());
     target_dist->read_multicol(targetdiststream);
@@ -221,9 +223,9 @@ int colvarbias_meta::init_ebmeta_params(std::string const &conf)
     }
     // normalize target distribution and multiply by effective volume = exp(differential entropy)
     target_dist->multiply_constant(1.0/target_dist->integral());
-    cvm::real volume = std::exp(target_dist->entropy());
+    cvm::real volume = cvm::exp(target_dist->entropy());
     target_dist->multiply_constant(volume);
-    get_keyval(conf, "ebMetaEquilSteps", ebmeta_equil_steps, 0);
+    get_keyval(conf, "ebMetaEquilSteps", ebmeta_equil_steps, ebmeta_equil_steps);
   }
 
   return COLVARS_OK;
@@ -291,7 +293,7 @@ colvarbias_meta::create_hill(colvarbias_meta::hill const &h)
     // need to be computed analytically when the colvar returns
     // off-grid
     cvm::real const min_dist = hills_energy->bin_distance_from_boundaries(h.centers, true);
-    if (min_dist < (3.0 * std::floor(hill_width)) + 1.0) {
+    if (min_dist < (3.0 * cvm::floor(hill_width)) + 1.0) {
       hills_off_grid.push_back(h);
     }
   }
@@ -387,7 +389,7 @@ int colvarbias_meta::update_grid_params()
       // first of all, expand the grids, if specified
       bool changed_grids = false;
       int const min_buffer =
-        (3 * (size_t) std::floor(hill_width)) + 1;
+        (3 * (size_t) cvm::floor(hill_width)) + 1;
 
       std::vector<int>         new_sizes(hills_energy->sizes());
       std::vector<colvarvalue> new_lower_boundaries(hills_energy->lower_boundaries);
@@ -492,9 +494,9 @@ int colvarbias_meta::update_bias()
 
     if (ebmeta) {
       hills_scale *= 1.0/target_dist->value(target_dist->get_colvars_index());
-      if(cvm::step_absolute() <= long(ebmeta_equil_steps)) {
+      if(cvm::step_absolute() <= ebmeta_equil_steps) {
         cvm::real const hills_lambda =
-          (cvm::real(long(ebmeta_equil_steps) - cvm::step_absolute())) /
+          (cvm::real(ebmeta_equil_steps - cvm::step_absolute())) /
           (cvm::real(ebmeta_equil_steps));
         hills_scale = hills_lambda + (1-hills_lambda)*hills_scale;
       }
@@ -508,7 +510,7 @@ int colvarbias_meta::update_bias()
       } else {
         calc_hills(new_hills_begin, hills.end(), hills_energy_sum_here);
       }
-      hills_scale *= std::exp(-1.0*hills_energy_sum_here/(bias_temperature*cvm::boltzmann()));
+      hills_scale *= cvm::exp(-1.0*hills_energy_sum_here/(bias_temperature*cvm::boltzmann()));
     }
 
     switch (comm) {
@@ -710,7 +712,7 @@ void colvarbias_meta::calc_hills(colvarbias_meta::hill_iter      h_first,
       // set it to zero if the exponent is more negative than log(1.0E-05)
       h->value(0.0);
     } else {
-      h->value(std::exp(-0.5*cv_sqdev));
+      h->value(cvm::exp(-0.5*cv_sqdev));
     }
     energy += h->energy();
   }
@@ -904,7 +906,7 @@ void colvarbias_meta::recount_hills_off_grid(colvarbias_meta::hill_iter  h_first
 
   for (hill_iter h = h_first; h != h_last; h++) {
     cvm::real const min_dist = hills_energy->bin_distance_from_boundaries(h->centers, true);
-    if (min_dist < (3.0 * std::floor(hill_width)) + 1.0) {
+    if (min_dist < (3.0 * cvm::floor(hill_width)) + 1.0) {
       hills_off_grid.push_back(*h);
     }
   }
@@ -1427,8 +1429,8 @@ std::istream & colvarbias_meta::read_hill(std::istream &is)
     return is;
   }
 
-  size_t h_it;
+  cvm::step_number h_it;
-  get_keyval(data, "step", h_it, 0, parse_silent);
+  get_keyval(data, "step", h_it, 0L, parse_silent);
   if (h_it <= state_file_step) {
     if (cvm::debug())
       cvm::log("Skipping a hill older than the state file for metadynamics bias \""+
@@ -1457,7 +1459,7 @@ std::istream & colvarbias_meta::read_hill(std::istream &is)
 
   std::vector<cvm::real> h_widths(num_variables());
   get_keyval(data, "widths", h_widths,
-             std::vector<cvm::real>(num_variables(), (std::sqrt(2.0 * PI) / 2.0)),
+             std::vector<cvm::real>(num_variables(), (cvm::sqrt(2.0 * PI) / 2.0)),
              parse_silent);
 
   std::string h_replica = "";
@@ -1482,7 +1484,7 @@ std::istream & colvarbias_meta::read_hill(std::istream &is)
     // be computed analytically
     cvm::real const min_dist =
       hills_energy->bin_distance_from_boundaries((hills.back()).centers, true);
-    if (min_dist < (3.0 * std::floor(hill_width)) + 1.0) {
+    if (min_dist < (3.0 * cvm::floor(hill_width)) + 1.0) {
       hills_off_grid.push_back(hills.back());
     }
   }

lib/colvars/colvarbias_meta.h

@@ -174,12 +174,14 @@ protected:
   /// \brief Biasing temperature in well-tempered metadynamics
   cvm::real  bias_temperature;
 
-  // EBmeta parameters
+  /// Ensemble-biased metadynamics (EBmeta) flag
   bool       ebmeta;
+
+  /// Target distribution for EBmeta
   colvar_grid_scalar* target_dist;
-  std::string target_dist_file;
+
-  cvm::real target_dist_volume;
+  /// Number of equilibration steps for EBmeta
-  size_t ebmeta_equil_steps;
+  cvm::step_number ebmeta_equil_steps;
 
 
   /// \brief Try to read the restart information by allocating new
@@ -285,7 +287,7 @@ public:
   friend class colvarbias_meta;
 
   /// Time step at which this hill was added
-  size_t      it;
+  cvm::step_number it;
 
   /// Identity of the replica who added this hill (only in multiple replica simulations)
   std::string replica;
@@ -325,7 +327,7 @@ public:
   /// weight Weight of the hill \param centers Center of the hill
   /// \param widths Width of the hill around centers \param replica
   /// (optional) Identity of the replica which creates the hill
-  inline hill(size_t                    const &it_in,
+  inline hill(cvm::step_number          const &it_in,
               cvm::real                 const &W_in,
               std::vector<colvarvalue>  const &centers_in,
               std::vector<cvm::real>    const &widths_in,

lib/colvars/colvarbias_restraint.cpp

@@ -7,8 +7,6 @@
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 
-#include <cmath>
-
 #include "colvarmodule.h"
 #include "colvarproxy.h"
 #include "colvarvalue.h"
@@ -179,7 +177,7 @@ int colvarbias_restraint_k::change_configuration(std::string const &conf)
 colvarbias_restraint_moving::colvarbias_restraint_moving(char const *key)
 {
   target_nstages = 0;
-  target_nsteps = 0;
+  target_nsteps = 0L;
   stage = 0;
   acc_work = 0.0;
   b_chg_centers = false;
@@ -241,8 +239,8 @@ int colvarbias_restraint_moving::set_state_params(std::string const &conf)
 {
   if (b_chg_centers || b_chg_force_k) {
     if (target_nstages) {
-      if (!get_keyval(conf, "stage", stage))
+      get_keyval(conf, "stage", stage, stage,
-        cvm::error("Error: current stage is missing from the restart.\n");
+                 colvarparse::parse_restart | colvarparse::parse_required);
     }
   }
   return COLVARS_OK;
@@ -436,13 +434,13 @@ int colvarbias_restraint_centers_moving::set_state_params(std::string const &con
   colvarbias_restraint::set_state_params(conf);
 
   if (b_chg_centers) {
-    //    cvm::log ("Reading the updated restraint centers from the restart.\n");
+    get_keyval(conf, "centers", colvar_centers, colvar_centers,
-    if (!get_keyval(conf, "centers", colvar_centers))
+               colvarparse::parse_restart | colvarparse::parse_required);
-      cvm::error("Error: restraint centers are missing from the restart.\n");
-    if (is_enabled(f_cvb_output_acc_work)) {
-      if (!get_keyval(conf, "accumulatedWork", acc_work))
-        cvm::error("Error: accumulatedWork is missing from the restart.\n");
   }
+
+  if (is_enabled(f_cvb_output_acc_work)) {
+    get_keyval(conf, "accumulatedWork", acc_work, acc_work,
+               colvarparse::parse_restart | colvarparse::parse_required);
   }
 
   return COLVARS_OK;
@@ -563,7 +561,7 @@ int colvarbias_restraint_k_moving::update()
           lambda = 0.0;
         }
         force_k = starting_force_k + (target_force_k - starting_force_k)
-          * std::pow(lambda, force_k_exp);
+          * cvm::pow(lambda, force_k_exp);
           cvm::log("Restraint " + this->name + ", stage " + cvm::to_str(stage)
                   + " : lambda = " + cvm::to_str(lambda)
                   + ", k = " + cvm::to_str(force_k));
@@ -585,7 +583,7 @@ int colvarbias_restraint_k_moving::update()
         for (size_t i = 0; i < num_variables(); i++) {
           dU_dk += d_restraint_potential_dk(i);
         }
-        restraint_FE += force_k_exp * std::pow(lambda, force_k_exp - 1.0)
+        restraint_FE += force_k_exp * cvm::pow(lambda, force_k_exp - 1.0)
           * (target_force_k - starting_force_k) * dU_dk;
       }
 
@@ -608,7 +606,7 @@ int colvarbias_restraint_k_moving::update()
             lambda = cvm::real(stage) / cvm::real(target_nstages);
           }
           force_k = starting_force_k + (target_force_k - starting_force_k)
-            * std::pow(lambda, force_k_exp);
+            * cvm::pow(lambda, force_k_exp);
           cvm::log("Restraint " + this->name + ", stage " + cvm::to_str(stage)
                   + " : lambda = " + cvm::to_str(lambda)
                   + ", k = " + cvm::to_str(force_k));
@@ -622,7 +620,7 @@ int colvarbias_restraint_k_moving::update()
       lambda = cvm::real(cvm::step_absolute()) / cvm::real(target_nsteps);
       cvm::real const force_k_old = force_k;
       force_k = starting_force_k + (target_force_k - starting_force_k)
-        * std::pow(lambda, force_k_exp);
+        * cvm::pow(lambda, force_k_exp);
       force_k_incr = force_k - force_k_old;
     }
   }
@@ -672,13 +670,13 @@ int colvarbias_restraint_k_moving::set_state_params(std::string const &conf)
   colvarbias_restraint::set_state_params(conf);
 
   if (b_chg_force_k) {
-    //    cvm::log ("Reading the updated force constant from the restart.\n");
+    get_keyval(conf, "forceConstant", force_k, force_k,
-    if (!get_keyval(conf, "forceConstant", force_k, force_k))
+               colvarparse::parse_restart | colvarparse::parse_required);
-      cvm::error("Error: force constant is missing from the restart.\n");
-    if (is_enabled(f_cvb_output_acc_work)) {
-      if (!get_keyval(conf, "accumulatedWork", acc_work))
-        cvm::error("Error: accumulatedWork is missing from the restart.\n");
   }
+
+  if (is_enabled(f_cvb_output_acc_work)) {
+    get_keyval(conf, "accumulatedWork", acc_work, acc_work,
+               colvarparse::parse_restart | colvarparse::parse_required);
   }
 
   return COLVARS_OK;
@@ -719,7 +717,8 @@ std::istream & colvarbias_restraint::read_state(std::istream &is)
        (set_state_params(conf) != COLVARS_OK) ) {
     cvm::error("Error: in reading state configuration for \""+bias_type+"\" bias \""+
                this->name+"\" at position "+
-               cvm::to_str(is.tellg())+" in stream.\n", INPUT_ERROR);
+               cvm::to_str(static_cast<size_t>(is.tellg()))+
+               " in stream.\n", INPUT_ERROR);
     is.clear();
     is.seekg(start_pos, std::ios::beg);
     is.setstate(std::ios::failbit);
@@ -729,7 +728,8 @@ std::istream & colvarbias_restraint::read_state(std::istream &is)
   if (!read_state_data(is)) {
     cvm::error("Error: in reading state data for \""+bias_type+"\" bias \""+
                this->name+"\" at position "+
-               cvm::to_str(is.tellg())+" in stream.\n", INPUT_ERROR);
+               cvm::to_str(static_cast<size_t>(is.tellg()))+
+               " in stream.\n", INPUT_ERROR);
     is.clear();
     is.seekg(start_pos, std::ios::beg);
     is.setstate(std::ios::failbit);
@@ -740,7 +740,7 @@ std::istream & colvarbias_restraint::read_state(std::istream &is)
     cvm::log("brace = "+brace+"\n");
     cvm::error("Error: corrupt restart information for \""+bias_type+"\" bias \""+
                this->name+"\": no matching brace at position "+
-               cvm::to_str(is.tellg())+" in stream.\n");
+               cvm::to_str(static_cast<size_t>(is.tellg()))+" in stream.\n");
     is.setstate(std::ios::failbit);
   }
 
@@ -787,11 +787,11 @@ int colvarbias_restraint_harmonic::init(std::string const &conf)
   colvarbias_restraint_k_moving::init(conf);
 
   for (size_t i = 0; i < num_variables(); i++) {
-    if (variables(i)->width != 1.0)
+    cvm::real const w = variables(i)->width;
     cvm::log("The force constant for colvar \""+variables(i)->name+
              "\" will be rescaled to "+
-               cvm::to_str(force_k / (variables(i)->width * variables(i)->width))+
+             cvm::to_str(force_k/(w*w))+
-               " according to the specified width.\n");
+             " according to the specified width ("+cvm::to_str(w)+").\n");
   }
 
   return COLVARS_OK;
@@ -1014,7 +1014,7 @@ int colvarbias_restraint_harmonic_walls::init(std::string const &conf)
                  INPUT_ERROR);
       return INPUT_ERROR;
     }
-    force_k = std::sqrt(lower_wall_k * upper_wall_k);
+    force_k = cvm::sqrt(lower_wall_k * upper_wall_k);
     // transform the two constants to relative values using gemetric mean as ref
     // to preserve force_k if provided as single parameter
     // (allow changing both via force_k)
@@ -1037,25 +1037,21 @@ int colvarbias_restraint_harmonic_walls::init(std::string const &conf)
 
   if (lower_walls.size() > 0) {
     for (i = 0; i < num_variables(); i++) {
-      if (variables(i)->width != 1.0)
+      cvm::real const w = variables(i)->width;
       cvm::log("The lower wall force constant for colvar \""+
-                 variables(i)->name+
+               variables(i)->name+"\" will be rescaled to "+
-                 "\" will be rescaled to "+
+               cvm::to_str(lower_wall_k * force_k / (w*w))+
-                 cvm::to_str(lower_wall_k * force_k /
+               " according to the specified width ("+cvm::to_str(w)+").\n");
-                             (variables(i)->width * variables(i)->width))+
-                 " according to the specified width.\n");
     }
   }
 
   if (upper_walls.size() > 0) {
     for (i = 0; i < num_variables(); i++) {
-      if (variables(i)->width != 1.0)
+      cvm::real const w = variables(i)->width;
       cvm::log("The upper wall force constant for colvar \""+
-                 variables(i)->name+
+               variables(i)->name+"\" will be rescaled to "+
-                 "\" will be rescaled to "+
+               cvm::to_str(upper_wall_k * force_k / (w*w))+
-                 cvm::to_str(upper_wall_k * force_k /
+               " according to the specified width ("+cvm::to_str(w)+").\n");
-                             (variables(i)->width * variables(i)->width))+
-                 " according to the specified width.\n");
     }
   }
 
@@ -1225,11 +1221,11 @@ int colvarbias_restraint_linear::init(std::string const &conf)
                  INPUT_ERROR);
       return INPUT_ERROR;
     }
-    if (variables(i)->width != 1.0)
+    cvm::real const w = variables(i)->width;
     cvm::log("The force constant for colvar \""+variables(i)->name+
              "\" will be rescaled to "+
-               cvm::to_str(force_k / variables(i)->width)+
+             cvm::to_str(force_k / w)+
-               " according to the specified width.\n");
+             " according to the specified width ("+cvm::to_str(w)+").\n");
   }
 
   return COLVARS_OK;
@@ -1367,6 +1363,7 @@ colvarbias_restraint_histogram::colvarbias_restraint_histogram(char const *key)
 int colvarbias_restraint_histogram::init(std::string const &conf)
 {
   colvarbias::init(conf);
+  enable(f_cvb_apply_force);
 
   get_keyval(conf, "lowerBoundary", lower_boundary, lower_boundary);
   get_keyval(conf, "upperBoundary", upper_boundary, upper_boundary);
@@ -1390,7 +1387,7 @@ int colvarbias_restraint_histogram::init(std::string const &conf)
   cvm::real const nbins = (upper_boundary - lower_boundary) / width;
   int const nbins_round = (int)(nbins);
 
-  if (std::fabs(nbins - cvm::real(nbins_round)) > 1.0E-10) {
+  if (cvm::fabs(nbins - cvm::real(nbins_round)) > 1.0E-10) {
     cvm::log("Warning: grid interval ("+
              cvm::to_str(lower_boundary, cvm::cv_width, cvm::cv_prec)+" - "+
              cvm::to_str(upper_boundary, cvm::cv_width, cvm::cv_prec)+
@@ -1440,7 +1437,7 @@ int colvarbias_restraint_histogram::init(std::string const &conf)
     }
   }
   cvm::real const ref_integral = ref_p.sum() * width;
-  if (std::fabs(ref_integral - 1.0) > 1.0e-03) {
+  if (cvm::fabs(ref_integral - 1.0) > 1.0e-03) {
     cvm::log("Reference distribution not normalized, normalizing to unity.\n");
     ref_p /= ref_integral;
   }
@@ -1471,7 +1468,7 @@ int colvarbias_restraint_histogram::update()
     vector_size += variables(icv)->value().size();
   }
 
-  cvm::real const norm = 1.0/(std::sqrt(2.0*PI)*gaussian_width*vector_size);
+  cvm::real const norm = 1.0/(cvm::sqrt(2.0*PI)*gaussian_width*vector_size);
 
   // calculate the histogram
   p.reset();
@@ -1482,7 +1479,7 @@ int colvarbias_restraint_histogram::update()
       size_t igrid;
       for (igrid = 0; igrid < p.size(); igrid++) {
         cvm::real const x_grid = (lower_boundary + (igrid+0.5)*width);
-        p[igrid] += norm * std::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) /
+        p[igrid] += norm * cvm::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) /
                                     (2.0 * gaussian_width * gaussian_width));
       }
     } else if (cv.type() == colvarvalue::type_vector) {
@@ -1492,7 +1489,7 @@ int colvarbias_restraint_histogram::update()
         size_t igrid;
         for (igrid = 0; igrid < p.size(); igrid++) {
           cvm::real const x_grid = (lower_boundary + (igrid+0.5)*width);
-          p[igrid] += norm * std::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) /
+          p[igrid] += norm * cvm::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) /
                                       (2.0 * gaussian_width * gaussian_width));
         }
       }
@@ -1523,7 +1520,7 @@ int colvarbias_restraint_histogram::update()
       for (igrid = 0; igrid < p.size(); igrid++) {
         cvm::real const x_grid = (lower_boundary + (igrid+0.5)*width);
         force += force_k_cv * p_diff[igrid] *
-          norm * std::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) /
+          norm * cvm::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) /
                           (2.0 * gaussian_width * gaussian_width)) *
           (-1.0 * (x_grid - cv_value) / (gaussian_width * gaussian_width));
       }
@@ -1536,7 +1533,7 @@ int colvarbias_restraint_histogram::update()
         for (igrid = 0; igrid < p.size(); igrid++) {
           cvm::real const x_grid = (lower_boundary + (igrid+0.5)*width);
           force += force_k_cv * p_diff[igrid] *
-            norm * std::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) /
+            norm * cvm::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) /
                             (2.0 * gaussian_width * gaussian_width)) *
             (-1.0 * (x_grid - cv_value) / (gaussian_width * gaussian_width));
         }
@@ -1550,7 +1547,7 @@ int colvarbias_restraint_histogram::update()
 }
 
 
-std::ostream & colvarbias_restraint_histogram::write_restart(std::ostream &os)
+int colvarbias_restraint_histogram::write_output_files()
 {
   if (b_write_histogram) {
     std::string file_name(cvm::output_prefix()+"."+this->name+".hist.dat");
@@ -1558,6 +1555,9 @@ std::ostream & colvarbias_restraint_histogram::write_restart(std::ostream &os)
     *os << "# " << cvm::wrap_string(variables(0)->name, cvm::cv_width)
         << "  " << "p(" << cvm::wrap_string(variables(0)->name, cvm::cv_width-3)
         << ")\n";
+
+    os->setf(std::ios::fixed, std::ios::floatfield);
+
     size_t igrid;
     for (igrid = 0; igrid < p.size(); igrid++) {
       cvm::real const x_grid = (lower_boundary + (igrid+1)*width);
@@ -1572,13 +1572,7 @@ std::ostream & colvarbias_restraint_histogram::write_restart(std::ostream &os)
     }
     cvm::proxy->close_output_stream(file_name);
   }
-  return os;
+  return COLVARS_OK;
-}
-
-
-std::istream & colvarbias_restraint_histogram::read_restart(std::istream &is)
-{
-  return is;
 }
 
 

lib/colvars/colvarbias_restraint.h

@@ -132,7 +132,7 @@ protected:
 
   /// \brief Number of steps required to reach the target force constant
   /// or restraint centers
-  long target_nsteps;
+  cvm::step_number target_nsteps;
 
   /// \brief Accumulated work (computed when outputAccumulatedWork == true)
   cvm::real acc_work;
@@ -328,8 +328,7 @@ public:
 
   virtual int update();
 
-  virtual std::istream & read_restart(std::istream &is);
+  virtual int write_output_files();
-  virtual std::ostream & write_restart(std::ostream &os);
   virtual std::ostream & write_traj_label(std::ostream &os);
   virtual std::ostream & write_traj(std::ostream &os);
 

lib/colvars/colvarcomp.cpp

@@ -20,7 +20,8 @@ colvar::cvc::cvc()
     b_periodic(false),
     b_try_scalable(true)
 {
-  init_cvc_requires();
+  description = "uninitialized colvar component";
+  init_dependencies();
   sup_coeff = 1.0;
   period = 0.0;
   wrap_center = 0.0;
@@ -33,7 +34,8 @@ colvar::cvc::cvc(std::string const &conf)
     b_periodic(false),
     b_try_scalable(true)
 {
-  init_cvc_requires();
+  description = "uninitialized colvar component";
+  init_dependencies();
   sup_coeff = 1.0;
   period = 0.0;
   wrap_center = 0.0;
@@ -176,6 +178,100 @@ cvm::atom_group *colvar::cvc::parse_group(std::string const &conf,
 }
 
 
+int colvar::cvc::init_dependencies() {
+  size_t i;
+  // Initialize static array once and for all
+  if (features().size() == 0) {
+    for (i = 0; i < colvardeps::f_cvc_ntot; i++) {
+      modify_features().push_back(new feature);
+    }
+
+    init_feature(f_cvc_active, "active", f_type_dynamic);
+//     The dependency below may become useful if we use dynamic atom groups
+//     require_feature_children(f_cvc_active, f_ag_active);
+
+    init_feature(f_cvc_scalar, "scalar", f_type_static);
+
+    init_feature(f_cvc_gradient, "gradient", f_type_dynamic);
+
+    init_feature(f_cvc_explicit_gradient, "explicit gradient", f_type_static);
+    require_feature_children(f_cvc_explicit_gradient, f_ag_explicit_gradient);
+
+    init_feature(f_cvc_inv_gradient, "inverse gradient", f_type_dynamic);
+    require_feature_self(f_cvc_inv_gradient, f_cvc_gradient);
+
+    init_feature(f_cvc_debug_gradient, "debug gradient", f_type_user);
+    require_feature_self(f_cvc_debug_gradient, f_cvc_gradient);
+    require_feature_self(f_cvc_debug_gradient, f_cvc_explicit_gradient);
+
+    init_feature(f_cvc_Jacobian, "Jacobian derivative", f_type_dynamic);
+    require_feature_self(f_cvc_Jacobian, f_cvc_inv_gradient);
+
+    init_feature(f_cvc_com_based, "depends on group centers of mass", f_type_static);
+
+    init_feature(f_cvc_pbc_minimum_image, "use minimum-image distances with PBCs", f_type_user);
+
+    // Compute total force on first site only to avoid unwanted
+    // coupling to other colvars (see e.g. Ciccotti et al., 2005)
+    init_feature(f_cvc_one_site_total_force, "compute total force from one group", f_type_user);
+    require_feature_self(f_cvc_one_site_total_force, f_cvc_com_based);
+
+    init_feature(f_cvc_scalable, "scalable calculation", f_type_static);
+    require_feature_self(f_cvc_scalable, f_cvc_scalable_com);
+
+    init_feature(f_cvc_scalable_com, "scalable calculation of centers of mass", f_type_static);
+    require_feature_self(f_cvc_scalable_com, f_cvc_com_based);
+
+
+    // TODO only enable this when f_ag_scalable can be turned on for a pre-initialized group
+    // require_feature_children(f_cvc_scalable, f_ag_scalable);
+    // require_feature_children(f_cvc_scalable_com, f_ag_scalable_com);
+
+    // check that everything is initialized
+    for (i = 0; i < colvardeps::f_cvc_ntot; i++) {
+      if (is_not_set(i)) {
+        cvm::error("Uninitialized feature " + cvm::to_str(i) + " in " + description);
+      }
+    }
+  }
+
+  // Initialize feature_states for each instance
+  // default as available, not enabled
+  // except dynamic features which default as unavailable
+  feature_states.reserve(f_cvc_ntot);
+  for (i = 0; i < colvardeps::f_cvc_ntot; i++) {
+    bool avail = is_dynamic(i) ? false : true;
+    feature_states.push_back(feature_state(avail, false));
+  }
+
+  // Features that are implemented by all cvcs by default
+  // Each cvc specifies what other features are available
+  feature_states[f_cvc_active].available = true;
+  feature_states[f_cvc_gradient].available = true;
+
+  // CVCs are enabled from the start - get disabled based on flags
+  enable(f_cvc_active);
+  // feature_states[f_cvc_active].enabled = true;
+
+  // Explicit gradients are implemented in mosts CVCs. Exceptions must be specified explicitly.
+  // feature_states[f_cvc_explicit_gradient].enabled = true;
+  enable(f_cvc_explicit_gradient);
+
+  // Use minimum-image distances by default
+  // feature_states[f_cvc_pbc_minimum_image].enabled = true;
+  enable(f_cvc_pbc_minimum_image);
+
+  // Features that are implemented by default if their requirements are
+  feature_states[f_cvc_one_site_total_force].available = true;
+
+  // Features That are implemented only for certain simulation engine configurations
+  feature_states[f_cvc_scalable_com].available = (cvm::proxy->scalable_group_coms() == COLVARS_OK);
+  feature_states[f_cvc_scalable].available = feature_states[f_cvc_scalable_com].available;
+
+  return COLVARS_OK;
+}
+
+
 int colvar::cvc::setup()
 {
   description = "cvc " + name;
@@ -192,6 +288,7 @@ colvar::cvc::~cvc()
   }
 }
 
+
 void colvar::cvc::read_data()
 {
   size_t ig;
@@ -214,6 +311,66 @@ void colvar::cvc::read_data()
 }
 
 
+std::vector<std::vector<int> > colvar::cvc::get_atom_lists()
+{
+  std::vector<std::vector<int> > lists;
+
+  std::vector<cvm::atom_group *>::iterator agi = atom_groups.begin();
+  for ( ; agi != atom_groups.end(); ++agi) {
+    (*agi)->create_sorted_ids();
+    lists.push_back((*agi)->sorted_ids());
+    if ((*agi)->is_enabled(f_ag_fitting_group) && (*agi)->is_enabled(f_ag_fit_gradients)) {
+      cvm::atom_group &fg = *((*agi)->fitting_group);
+      fg.create_sorted_ids();
+      lists.push_back(fg.sorted_ids());
+    }
+  }
+  return lists;
+}
+
+
+void colvar::cvc::collect_gradients(std::vector<int> const &atom_ids, std::vector<cvm::rvector> &atomic_gradients)
+{
+  // Coefficient: d(a * x^n) = a * n * x^(n-1) * dx
+  cvm::real coeff = sup_coeff * cvm::real(sup_np) *
+    cvm::integer_power(value().real_value, sup_np-1);
+
+  for (size_t j = 0; j < atom_groups.size(); j++) {
+
+    cvm::atom_group &ag = *(atom_groups[j]);
+
+    // If necessary, apply inverse rotation to get atomic
+    // gradient in the laboratory frame
+    if (ag.b_rotate) {
+      cvm::rotation const rot_inv = ag.rot.inverse();
+
+      for (size_t k = 0; k < ag.size(); k++) {
+        size_t a = std::lower_bound(atom_ids.begin(), atom_ids.end(),
+                                    ag[k].id) - atom_ids.begin();
+        atomic_gradients[a] += coeff * rot_inv.rotate(ag[k].grad);
+      }
+
+    } else {
+
+      for (size_t k = 0; k < ag.size(); k++) {
+        size_t a = std::lower_bound(atom_ids.begin(), atom_ids.end(),
+                                    ag[k].id) - atom_ids.begin();
+        atomic_gradients[a] += coeff * ag[k].grad;
+      }
+    }
+    if (ag.is_enabled(f_ag_fitting_group) && ag.is_enabled(f_ag_fit_gradients)) {
+      cvm::atom_group const &fg = *(ag.fitting_group);
+      for (size_t k = 0; k < fg.size(); k++) {
+        size_t a = std::lower_bound(atom_ids.begin(), atom_ids.end(),
+                                    fg[k].id) - atom_ids.begin();
+        // fit gradients are in the unrotated (simulation) frame
+        atomic_gradients[a] += coeff * fg.fit_gradients[k];
+      }
+    }
+  }
+}
+
+
 void colvar::cvc::calc_force_invgrads()
 {
   cvm::error("Error: calculation of inverse gradients is not implemented "
@@ -306,8 +463,8 @@ void colvar::cvc::debug_gradients()
         cvm::log("dx(interp) = "+cvm::to_str(dx_pred,
                               21, 14)+"\n");
         cvm::log("|dx(actual) - dx(interp)|/|dx(actual)| = "+
-                  cvm::to_str(std::fabs(x_1 - x_0 - dx_pred) /
+                  cvm::to_str(cvm::fabs(x_1 - x_0 - dx_pred) /
-                              std::fabs(x_1 - x_0), 12, 5)+"\n");
+                              cvm::fabs(x_1 - x_0), 12, 5)+"\n");
       }
     }
 
@@ -341,8 +498,8 @@ void colvar::cvc::debug_gradients()
           cvm::log("dx(interp) = "+cvm::to_str (dx_pred,
                                 21, 14)+"\n");
           cvm::log ("|dx(actual) - dx(interp)|/|dx(actual)| = "+
-                    cvm::to_str(std::fabs (x_1 - x_0 - dx_pred) /
+                    cvm::to_str(cvm::fabs (x_1 - x_0 - dx_pred) /
-                                std::fabs (x_1 - x_0),
+                                cvm::fabs (x_1 - x_0),
                                 12, 5)+
                     ".\n");
         }
@@ -378,7 +535,7 @@ colvarvalue colvar::cvc::dist2_rgrad(colvarvalue const &x1,
 }
 
 
-void colvar::cvc::wrap(colvarvalue &x) const
+void colvar::cvc::wrap(colvarvalue &x_unwrapped) const
 {
   return;
 }

lib/colvars/colvarcomp.h

@@ -27,12 +27,12 @@
 
 /// \brief Colvar component (base class for collective variables)
 ///
-/// A \link cvc \endlink object (or an object of a
+/// A \link colvar::cvc \endlink object (or an object of a
 /// cvc-derived class) implements the calculation of a collective
 /// variable, its gradients and any other related physical quantities
 /// that depend on microscopic degrees of freedom.
 ///
-/// No restriction is set to what kind of calculation a \link cvc \endlink
+/// No restriction is set to what kind of calculation a \link colvar::cvc \endlink
 /// object performs (usually an analytical function of atomic coordinates).
 /// The only constraints are that: \par
 ///
@@ -42,9 +42,9 @@
 ///   alike, and allows an automatic selection of the applicable algorithms.
 ///
 /// - The object provides an implementation \link apply_force() \endlink to
-///   apply forces to atoms.  Typically, one or more \link cvm::atom_group
+///   apply forces to atoms.  Typically, one or more \link colvarmodule::atom_group
 ///   \endlink objects are used, but this is not a requirement for as long as
-///   the \link cvc \endlink object communicates with the simulation program.
+///   the \link colvar::cvc \endlink object communicates with the simulation program.
 ///
 /// <b> If you wish to implement a new collective variable component, you
 /// should write your own class by inheriting directly from \link
@@ -75,9 +75,9 @@ public:
   /// \brief Description of the type of collective variable
   ///
   /// Normally this string is set by the parent \link colvar \endlink
-  /// object within its constructor, when all \link cvc \endlink
+  /// object within its constructor, when all \link colvar::cvc \endlink
   /// objects are initialized; therefore the main "config string"
-  /// constructor does not need to define it.  If a \link cvc
+  /// constructor does not need to define it.  If a \link colvar::cvc
   /// \endlink is initialized and/or a different constructor is used,
   /// this variable definition should be set within the constructor.
   std::string function_type;
@@ -109,6 +109,9 @@ public:
   /// cvc \endlink
   virtual int init(std::string const &conf);
 
+  /// \brief Initialize dependency tree
+  virtual int init_dependencies();
+
   /// \brief Within the constructor, make a group parse its own
   /// options from the provided configuration string
   /// Returns reference to new group
@@ -122,7 +125,7 @@ public:
   /// \brief After construction, set data related to dependency handling
   int setup();
 
-  /// \brief Default constructor (used when \link cvc \endlink
+  /// \brief Default constructor (used when \link colvar::cvc \endlink
   /// objects are declared within other ones)
   cvc();
 
@@ -133,7 +136,7 @@ public:
   static std::vector<feature *> cvc_features;
 
   /// \brief Implementation of the feature list accessor for colvar
-  virtual const std::vector<feature *> &features()
+  virtual const std::vector<feature *> &features() const
   {
     return cvc_features;
   }
@@ -148,6 +151,9 @@ public:
     cvc_features.clear();
   }
 
+  /// \brief Get vector of vectors of atom IDs for all atom groups
+  virtual std::vector<std::vector<int> > get_atom_lists();
+
   /// \brief Obtain data needed for the calculation for the backend
   virtual void read_data();
 
@@ -164,6 +170,10 @@ public:
   /// \brief Calculate finite-difference gradients alongside the analytical ones, for each Cartesian component
   virtual void debug_gradients();
 
+  /// \brief Calculate atomic gradients and add them to the corresponding item in gradient vector
+  /// May be overridden by CVCs that do not store their gradients in the classic way, see dihedPC
+  virtual void collect_gradients(std::vector<int> const &atom_ids, std::vector<cvm::rvector> &atomic_gradients);
+
   /// \brief Calculate the total force from the system using the
   /// inverse atomic gradients
   virtual void calc_force_invgrads();
@@ -237,7 +247,7 @@ public:
                                   colvarvalue const &x2) const;
 
   /// \brief Wrap value (for periodic/symmetric cvcs)
-  virtual void wrap(colvarvalue &x) const;
+  virtual void wrap(colvarvalue &x_unwrapped) const;
 
   /// \brief Pointers to all atom groups, to let colvars collect info
   /// e.g. atomic gradients
@@ -246,7 +256,7 @@ public:
   /// \brief Store a pointer to new atom group, and list as child for dependencies
   inline void register_atom_group(cvm::atom_group *ag) {
     atom_groups.push_back(ag);
-    add_child((colvardeps *) ag);
+    add_child(ag);
   }
 
   /// \brief Whether or not this CVC will be computed in parallel whenever possible
@@ -415,7 +425,7 @@ public:
   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
                                   colvarvalue const &x2) const;
   /// \brief Redefined to make use of the user-provided period
-  virtual void wrap(colvarvalue &x) const;
+  virtual void wrap(colvarvalue &x_unwrapped) const;
 };
 
 
@@ -474,7 +484,7 @@ public:
   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
                                   colvarvalue const &x2) const;
   /// Redefined to handle the 2*PI periodicity
-  virtual void wrap(colvarvalue &x) const;
+  virtual void wrap(colvarvalue &x_unwrapped) const;
 };
 
 
@@ -559,6 +569,35 @@ public:
 
 
 
+/// \brief Colvar component:  dipole magnitude of a molecule
+class colvar::dipole_magnitude
+  : public colvar::cvc
+{
+protected:
+  /// Dipole atom group
+  cvm::atom_group  *atoms;
+  cvm::atom_pos dipoleV;
+public:
+  /// Initialize by parsing the configuration
+  dipole_magnitude (std::string const &conf);
+  dipole_magnitude (cvm::atom const &a1);
+  dipole_magnitude();
+  virtual inline ~dipole_magnitude() {}
+  virtual void calc_value();
+  virtual void calc_gradients();
+  //virtual void calc_force_invgrads();
+  //virtual void calc_Jacobian_derivative();
+  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;
+};
+
+
+
 /// \brief Colvar component: Radius of gyration of an atom group
 /// (colvarvalue::type_scalar type, range [0:*))
 class colvar::gyration
@@ -818,7 +857,7 @@ public:
   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
                                   colvarvalue const &x2) const;
   /// Redefined to handle the 2*PI periodicity
-  virtual void wrap(colvarvalue &x) const;
+  virtual void wrap(colvarvalue &x_unwrapped) const;
 };
 
 
@@ -1002,7 +1041,7 @@ public:
          cvm::atom const &donor,
          cvm::real r0, int en, int ed);
   h_bond();
-  virtual ~h_bond();
+  virtual ~h_bond() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force(colvarvalue const &force);
@@ -1090,6 +1129,8 @@ public:
   virtual ~alpha_angles();
   void calc_value();
   void calc_gradients();
+  /// Re-implementation of cvc::collect_gradients() to carry over atomic gradients of sub-cvcs
+  void collect_gradients(std::vector<int> const &atom_ids, std::vector<cvm::rvector> &atomic_gradients);
   void apply_force(colvarvalue const &force);
   virtual cvm::real dist2(colvarvalue const &x1,
                           colvarvalue const &x2) const;
@@ -1120,6 +1161,8 @@ public:
   virtual  ~dihedPC();
   void calc_value();
   void calc_gradients();
+  /// Re-implementation of cvc::collect_gradients() to carry over atomic gradients of sub-cvcs
+  void collect_gradients(std::vector<int> const &atom_ids, std::vector<cvm::rvector> &atomic_gradients);
   void apply_force(colvarvalue const &force);
   virtual cvm::real dist2(colvarvalue const &x1,
                           colvarvalue const &x2) const;
@@ -1159,6 +1202,7 @@ public:
 
   orientation(std::string const &conf);
   orientation();
+  virtual int init(std::string const &conf);
   virtual ~orientation() {}
   virtual void calc_value();
   virtual void calc_gradients();
@@ -1183,6 +1227,7 @@ public:
 
   orientation_angle(std::string const &conf);
   orientation_angle();
+  virtual int init(std::string const &conf);
   virtual ~orientation_angle() {}
   virtual void calc_value();
   virtual void calc_gradients();
@@ -1207,6 +1252,7 @@ public:
 
   orientation_proj(std::string const &conf);
   orientation_proj();
+  virtual int init(std::string const &conf);
   virtual ~orientation_proj() {}
   virtual void calc_value();
   virtual void calc_gradients();
@@ -1234,6 +1280,7 @@ public:
 
   tilt(std::string const &conf);
   tilt();
+  virtual int init(std::string const &conf);
   virtual ~tilt() {}
   virtual void calc_value();
   virtual void calc_gradients();
@@ -1261,6 +1308,7 @@ public:
 
   spin_angle(std::string const &conf);
   spin_angle();
+  virtual int init(std::string const &conf);
   virtual ~spin_angle() {}
   virtual void calc_value();
   virtual void calc_gradients();
@@ -1275,7 +1323,7 @@ public:
   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
                                   colvarvalue const &x2) const;
   /// Redefined to handle the 2*PI periodicity
-  virtual void wrap(colvarvalue &x) const;
+  virtual void wrap(colvarvalue &x_unwrapped) const;
 };
 
 

lib/colvars/colvarcomp_angles.cpp

@@ -11,9 +11,6 @@
 #include "colvar.h"
 #include "colvarcomp.h"
 
-#include <cmath>
-
-
 
 colvar::angle::angle(std::string const &conf)
   : cvc(conf)
@@ -77,14 +74,14 @@ void colvar::angle::calc_value()
 
   cvm::real const cos_theta = (r21*r23)/(r21l*r23l);
 
-  x.real_value = (180.0/PI) * std::acos(cos_theta);
+  x.real_value = (180.0/PI) * cvm::acos(cos_theta);
 }
 
 
 void colvar::angle::calc_gradients()
 {
   cvm::real const cos_theta = (r21*r23)/(r21l*r23l);
-  cvm::real const dxdcos = -1.0 / std::sqrt(1.0 - cos_theta*cos_theta);
+  cvm::real const dxdcos = -1.0 / cvm::sqrt(1.0 - cos_theta*cos_theta);
 
   dxdr1 = (180.0/PI) * dxdcos *
     (1.0/r21l) * ( r23/r23l + (-1.0) * cos_theta * r21/r21l );
@@ -126,7 +123,7 @@ void colvar::angle::calc_Jacobian_derivative()
   // det(J) = (2 pi) r^2 * sin(theta)
   // hence Jd = cot(theta)
   const cvm::real theta = x.real_value * PI / 180.0;
-  jd = PI / 180.0 * (theta != 0.0 ? std::cos(theta) / std::sin(theta) : 0.0);
+  jd = PI / 180.0 * (theta != 0.0 ? cvm::cos(theta) / cvm::sin(theta) : 0.0);
 }
 
 
@@ -202,7 +199,7 @@ void colvar::dipole_angle::calc_value()
 
   cvm::real const cos_theta = (r21*r23)/(r21l*r23l);
 
-  x.real_value = (180.0/PI) * std::acos(cos_theta);
+  x.real_value = (180.0/PI) * cvm::acos(cos_theta);
 }
 
 //to be implemented
@@ -212,7 +209,7 @@ void colvar::dipole_angle::calc_value()
 void colvar::dipole_angle::calc_gradients()
 {
   cvm::real const cos_theta = (r21*r23)/(r21l*r23l);
-  cvm::real const dxdcos = -1.0 / std::sqrt(1.0 - cos_theta*cos_theta);
+  cvm::real const dxdcos = -1.0 / cvm::sqrt(1.0 - cos_theta*cos_theta);
 
   dxdr1 = (180.0/PI) * dxdcos *
   (1.0/r21l)* (r23/r23l + (-1.0) * cos_theta * r21/r21l );
@@ -346,7 +343,7 @@ void colvar::dihedral::calc_value()
   cvm::real const cos_phi = n1 * n2;
   cvm::real const sin_phi = n1 * r34 * r23.norm();
 
-  x.real_value = (180.0/PI) * std::atan2(sin_phi, cos_phi);
+  x.real_value = (180.0/PI) * cvm::atan2(sin_phi, cos_phi);
   this->wrap(x);
 }
 
@@ -368,7 +365,7 @@ void colvar::dihedral::calc_gradients()
   rB = 1.0/rB;
   B *= rB;
 
-  if (std::fabs(sin_phi) > 0.1) {
+  if (cvm::fabs(sin_phi) > 0.1) {
     rA = 1.0/rA;
     A *= rA;
     cvm::rvector const dcosdA = rA*(cos_phi*A-B);
@@ -440,8 +437,8 @@ void colvar::dihedral::calc_force_invgrads()
   cvm::real const dot1 = u23 * u12;
   cvm::real const dot4 = u23 * u34;
 
-  cvm::real const fact1 = d12 * std::sqrt(1.0 - dot1 * dot1);
+  cvm::real const fact1 = d12 * cvm::sqrt(1.0 - dot1 * dot1);
-  cvm::real const fact4 = d34 * std::sqrt(1.0 - dot4 * dot4);
+  cvm::real const fact4 = d34 * cvm::sqrt(1.0 - dot4 * dot4);
 
   group1->read_total_forces();
   if (is_enabled(f_cvc_one_site_total_force)) {
@@ -508,19 +505,17 @@ colvarvalue colvar::dihedral::dist2_rgrad(colvarvalue const &x1,
 }
 
 
-void colvar::dihedral::wrap(colvarvalue &x) const
+void colvar::dihedral::wrap(colvarvalue &x_unwrapped) const
 {
-  if ((x.real_value - wrap_center) >= 180.0) {
+  if ((x_unwrapped.real_value - wrap_center) >= 180.0) {
-    x.real_value -= 360.0;
+    x_unwrapped.real_value -= 360.0;
     return;
   }
 
-  if ((x.real_value - wrap_center) < -180.0) {
+  if ((x_unwrapped.real_value - wrap_center) < -180.0) {
-    x.real_value += 360.0;
+    x_unwrapped.real_value += 360.0;
     return;
   }
-
-  return;
 }
 
 
@@ -548,8 +543,8 @@ void colvar::polar_theta::calc_value()
   cvm::rvector pos = atoms->center_of_mass();
   r = atoms->center_of_mass().norm();
   // Internal values of theta and phi are radians
-  theta = (r > 0.) ? std::acos(pos.z / r) : 0.;
+  theta = (r > 0.) ? cvm::acos(pos.z / r) : 0.;
-  phi = std::atan2(pos.y, pos.x);
+  phi = cvm::atan2(pos.y, pos.x);
   x.real_value = (180.0/PI) * theta;
 }
 
@@ -560,9 +555,9 @@ void colvar::polar_theta::calc_gradients()
     atoms->set_weighted_gradient(cvm::rvector(0., 0., 0.));
   else
     atoms->set_weighted_gradient(cvm::rvector(
-      (180.0/PI) *  std::cos(theta) * std::cos(phi) / r,
+      (180.0/PI) *  cvm::cos(theta) * cvm::cos(phi) / r,
-      (180.0/PI) *  std::cos(theta) * std::sin(phi) / r,
+      (180.0/PI) *  cvm::cos(theta) * cvm::sin(phi) / r,
-      (180.0/PI) * -std::sin(theta) / r));
+      (180.0/PI) * -cvm::sin(theta) / r));
 }
 
 
@@ -602,8 +597,8 @@ void colvar::polar_phi::calc_value()
   cvm::rvector pos = atoms->center_of_mass();
   r = atoms->center_of_mass().norm();
   // Internal values of theta and phi are radians
-  theta = (r > 0.) ? std::acos(pos.z / r) : 0.;
+  theta = (r > 0.) ? cvm::acos(pos.z / r) : 0.;
-  phi = std::atan2(pos.y, pos.x);
+  phi = cvm::atan2(pos.y, pos.x);
   x.real_value = (180.0/PI) * phi;
 }
 
@@ -611,8 +606,8 @@ void colvar::polar_phi::calc_value()
 void colvar::polar_phi::calc_gradients()
 {
   atoms->set_weighted_gradient(cvm::rvector(
-    (180.0/PI) * -std::sin(phi) / (r*std::sin(theta)),
+    (180.0/PI) * -cvm::sin(phi) / (r*cvm::sin(theta)),
-    (180.0/PI) *  std::cos(phi) / (r*std::sin(theta)),
+    (180.0/PI) *  cvm::cos(phi) / (r*cvm::sin(theta)),
     0.));
 }
 
@@ -653,15 +648,15 @@ colvarvalue colvar::polar_phi::dist2_rgrad(colvarvalue const &x1,
 }
 
 
-void colvar::polar_phi::wrap(colvarvalue &x) const
+void colvar::polar_phi::wrap(colvarvalue &x_unwrapped) const
 {
-  if ((x.real_value - wrap_center) >= 180.0) {
+  if ((x_unwrapped.real_value - wrap_center) >= 180.0) {
-    x.real_value -= 360.0;
+    x_unwrapped.real_value -= 360.0;
     return;
   }
 
-  if ((x.real_value - wrap_center) < -180.0) {
+  if ((x_unwrapped.real_value - wrap_center) < -180.0) {
-    x.real_value += 360.0;
+    x_unwrapped.real_value += 360.0;
     return;
   }
 

lib/colvars/colvarcomp_coordnums.cpp

@@ -7,8 +7,6 @@
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 
-#include <cmath>
-
 #include "colvarmodule.h"
 #include "colvarparse.h"
 #include "colvaratoms.h"
@@ -102,6 +100,12 @@ colvar::coordnum::coordnum(std::string const &conf)
   group1 = parse_group(conf, "group1");
   group2 = parse_group(conf, "group2");
 
+  if (group1 == NULL || group2 == NULL) {
+    cvm::error("Error: failed to initialize atom groups.\n",
+                INPUT_ERROR);
+    return;
+  }
+
   if (int atom_number = cvm::atom_group::overlap(*group1, *group2)) {
     cvm::error("Error: group1 and group2 share a common atom (number: " +
                cvm::to_str(atom_number) + ")\n", INPUT_ERROR);
@@ -408,12 +412,6 @@ colvar::h_bond::h_bond()
 }
 
 
-colvar::h_bond::~h_bond()
-{
-  delete atom_groups[0];
-}
-
-
 void colvar::h_bond::calc_value()
 {
   int const flags = coordnum::ef_null;
@@ -655,8 +653,6 @@ colvar::groupcoordnum::groupcoordnum()
 
 void colvar::groupcoordnum::calc_value()
 {
-  cvm::rvector const r0_vec(0.0); // TODO enable the flag?
-
   // create fake atoms to hold the com coordinates
   cvm::atom group1_com_atom;
   cvm::atom group2_com_atom;
@@ -680,8 +676,6 @@ void colvar::groupcoordnum::calc_value()
 
 void colvar::groupcoordnum::calc_gradients()
 {
-  cvm::rvector const r0_vec(0.0); // TODO enable the flag?
-
   cvm::atom group1_com_atom;
   cvm::atom group2_com_atom;
   group1_com_atom.pos = group1->center_of_mass();

lib/colvars/colvarcomp_distances.cpp

@@ -7,8 +7,6 @@
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 
-#include <cmath>
-
 #include "colvarmodule.h"
 #include "colvarvalue.h"
 #include "colvarparse.h"
@@ -102,7 +100,7 @@ colvar::distance_vec::distance_vec(std::string const &conf)
 {
   function_type = "distance_vec";
   enable(f_cvc_com_based);
-  enable(f_cvc_implicit_gradient);
+  disable(f_cvc_explicit_gradient);
   x.type(colvarvalue::type_3vector);
 }
 
@@ -112,7 +110,7 @@ colvar::distance_vec::distance_vec()
 {
   function_type = "distance_vec";
   enable(f_cvc_com_based);
-  enable(f_cvc_implicit_gradient);
+  disable(f_cvc_explicit_gradient);
   x.type(colvarvalue::type_3vector);
 }
 
@@ -320,7 +318,7 @@ cvm::real colvar::distance_z::dist2(colvarvalue const &x1,
 {
   cvm::real diff = x1.real_value - x2.real_value;
   if (b_periodic) {
-    cvm::real shift = std::floor(diff/period + 0.5);
+    cvm::real shift = cvm::floor(diff/period + 0.5);
     diff -= shift * period;
   }
   return diff * diff;
@@ -332,7 +330,7 @@ colvarvalue colvar::distance_z::dist2_lgrad(colvarvalue const &x1,
 {
   cvm::real diff = x1.real_value - x2.real_value;
   if (b_periodic) {
-    cvm::real shift = std::floor(diff/period + 0.5);
+    cvm::real shift = cvm::floor(diff/period + 0.5);
     diff -= shift * period;
   }
   return 2.0 * diff;
@@ -344,22 +342,23 @@ colvarvalue colvar::distance_z::dist2_rgrad(colvarvalue const &x1,
 {
   cvm::real diff = x1.real_value - x2.real_value;
   if (b_periodic) {
-    cvm::real shift = std::floor(diff/period + 0.5);
+    cvm::real shift = cvm::floor(diff/period + 0.5);
     diff -= shift * period;
   }
   return (-2.0) * diff;
 }
 
 
-void colvar::distance_z::wrap(colvarvalue &x) const
+void colvar::distance_z::wrap(colvarvalue &x_unwrapped) const
 {
   if (!b_periodic) {
     // don't wrap if the period has not been set
     return;
   }
 
-  cvm::real shift = std::floor((x.real_value - wrap_center) / period + 0.5);
+  cvm::real shift =
-  x.real_value -= shift * period;
+    cvm::floor((x_unwrapped.real_value - wrap_center) / period + 0.5);
+  x_unwrapped.real_value -= shift * period;
   return;
 }
 
@@ -481,7 +480,7 @@ colvar::distance_dir::distance_dir(std::string const &conf)
 {
   function_type = "distance_dir";
   enable(f_cvc_com_based);
-  enable(f_cvc_implicit_gradient);
+  disable(f_cvc_explicit_gradient);
   x.type(colvarvalue::type_unit3vector);
 }
 
@@ -491,7 +490,7 @@ colvar::distance_dir::distance_dir()
 {
   function_type = "distance_dir";
   enable(f_cvc_com_based);
-  enable(f_cvc_implicit_gradient);
+  disable(f_cvc_explicit_gradient);
   x.type(colvarvalue::type_unit3vector);
 }
 
@@ -629,7 +628,7 @@ void colvar::distance_inv::calc_value()
   }
 
   x.real_value *= 1.0 / cvm::real(group1->size() * group2->size());
-  x.real_value = std::pow(x.real_value, -1.0/cvm::real(exponent));
+  x.real_value = cvm::pow(x.real_value, -1.0/cvm::real(exponent));
 
   cvm::real const dxdsum = (-1.0/(cvm::real(exponent))) *
     cvm::integer_power(x.real_value, exponent+1) /
@@ -671,7 +670,7 @@ colvar::distance_pairs::distance_pairs(std::string const &conf)
   group2 = parse_group(conf, "group2");
 
   x.type(colvarvalue::type_vector);
-  enable(f_cvc_implicit_gradient);
+  disable(f_cvc_explicit_gradient);
   x.vector1d_value.resize(group1->size() * group2->size());
 }
 
@@ -679,7 +678,7 @@ colvar::distance_pairs::distance_pairs(std::string const &conf)
 colvar::distance_pairs::distance_pairs()
 {
   function_type = "distance_pairs";
-  enable(f_cvc_implicit_gradient);
+  disable(f_cvc_explicit_gradient);
   x.type(colvarvalue::type_vector);
 }
 
@@ -747,6 +746,63 @@ void colvar::distance_pairs::apply_force(colvarvalue const &force)
 
 
 
+colvar::dipole_magnitude::dipole_magnitude(std::string const &conf)
+  : cvc(conf)
+{
+  function_type = "dipole_magnitude";
+  atoms = parse_group(conf, "atoms");
+  init_total_force_params(conf);
+  x.type(colvarvalue::type_scalar);
+}
+
+
+colvar::dipole_magnitude::dipole_magnitude(cvm::atom const &a1)
+{
+  atoms = new cvm::atom_group(std::vector<cvm::atom>(1, a1));
+  register_atom_group(atoms);
+  x.type(colvarvalue::type_scalar);
+}
+
+
+colvar::dipole_magnitude::dipole_magnitude()
+{
+  function_type = "dipole_magnitude";
+  x.type(colvarvalue::type_scalar);
+}
+
+
+void colvar::dipole_magnitude::calc_value()
+{
+  cvm::atom_pos const atomsCom = atoms->center_of_mass();
+  atoms->calc_dipole(atomsCom);
+  dipoleV = atoms->dipole();
+  x.real_value = dipoleV.norm();
+}
+
+
+void colvar::dipole_magnitude::calc_gradients()
+{
+  cvm::real const aux1 = atoms->total_charge/atoms->total_mass;
+  cvm::atom_pos const dipVunit = dipoleV.unit();
+
+  for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) {
+    ai->grad = (ai->charge - aux1*ai->mass) * dipVunit;
+  }
+}
+
+
+void colvar::dipole_magnitude::apply_force(colvarvalue const &force)
+{
+  if (!atoms->noforce) {
+    atoms->apply_colvar_force(force.real_value);
+  }
+}
+
+
+simple_scalar_dist_functions(dipole_magnitude)
+
+
+
 colvar::gyration::gyration(std::string const &conf)
   : cvc(conf)
 {
@@ -782,7 +838,7 @@ void colvar::gyration::calc_value()
   for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) {
     x.real_value += (ai->pos).norm2();
   }
-  x.real_value = std::sqrt(x.real_value / cvm::real(atoms->size()));
+  x.real_value = cvm::sqrt(x.real_value / cvm::real(atoms->size()));
 }
 
 
@@ -1029,7 +1085,7 @@ void colvar::rmsd::calc_value()
     x.real_value += ((*atoms)[ia].pos - ref_pos[ia]).norm2();
   }
   x.real_value /= cvm::real(atoms->size()); // MSD
-  x.real_value = std::sqrt(x.real_value);
+  x.real_value = cvm::sqrt(x.real_value);
 }
 
 
@@ -1405,7 +1461,7 @@ void colvar::eigenvector::calc_Jacobian_derivative()
     }
   }
 
-  jd.real_value = sum * std::sqrt(eigenvec_invnorm2);
+  jd.real_value = sum * cvm::sqrt(eigenvec_invnorm2);
 }
 
 
@@ -1436,7 +1492,7 @@ colvar::cartesian::cartesian(std::string const &conf)
   }
 
   x.type(colvarvalue::type_vector);
-  enable(f_cvc_implicit_gradient);
+  disable(f_cvc_explicit_gradient);
   x.vector1d_value.resize(atoms->size() * axes.size());
 }
 

lib/colvars/colvarcomp_protein.cpp

@@ -7,8 +7,6 @@
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 
-#include <cmath>
-
 #include "colvarmodule.h"
 #include "colvarvalue.h"
 #include "colvarparse.h"
@@ -27,7 +25,7 @@ colvar::alpha_angles::alpha_angles(std::string const &conf)
     cvm::log("Initializing alpha_angles object.\n");
 
   function_type = "alpha_angles";
-  enable(f_cvc_implicit_gradient);
+  enable(f_cvc_explicit_gradient);
   x.type(colvarvalue::type_scalar);
 
   std::string segment_id;
@@ -118,7 +116,7 @@ colvar::alpha_angles::alpha_angles()
   : cvc()
 {
   function_type = "alpha_angles";
-  enable(f_cvc_implicit_gradient);
+  enable(f_cvc_explicit_gradient);
   x.type(colvarvalue::type_scalar);
 }
 
@@ -133,6 +131,8 @@ colvar::alpha_angles::~alpha_angles()
     delete hb.back();
     hb.pop_back();
   }
+  // Our references to atom groups have become invalid now that children cvcs are deleted
+  atom_groups.clear();
 }
 
 
@@ -191,6 +191,58 @@ void colvar::alpha_angles::calc_gradients()
 }
 
 
+void colvar::alpha_angles::collect_gradients(std::vector<int> const &atom_ids, std::vector<cvm::rvector> &atomic_gradients)
+{
+  cvm::real cvc_coeff = sup_coeff * cvm::real(sup_np) * cvm::integer_power(value().real_value, sup_np-1);
+
+  if (theta.size()) {
+    cvm::real const theta_norm = (1.0-hb_coeff) / cvm::real(theta.size());
+
+    for (size_t i = 0; i < theta.size(); i++) {
+      cvm::real const t = ((theta[i])->value().real_value-theta_ref)/theta_tol;
+      cvm::real const f = ( (1.0 - (t*t)) /
+                            (1.0 - (t*t*t*t)) );
+      cvm::real const dfdt =
+        1.0/(1.0 - (t*t*t*t)) *
+        ( (-2.0 * t) + (-1.0*f)*(-4.0 * (t*t*t)) );
+
+      // Coeficient of this CVC's gradient in the colvar gradient, times coefficient of this
+      // angle's gradient in the CVC's gradient
+      cvm::real const coeff = cvc_coeff * theta_norm * dfdt * (1.0/theta_tol);
+
+      for (size_t j = 0; j < theta[i]->atom_groups.size(); j++) {
+        cvm::atom_group &ag = *(theta[i]->atom_groups[j]);
+        for (size_t k = 0; k < ag.size(); k++) {
+          size_t a = std::lower_bound(atom_ids.begin(), atom_ids.end(),
+                                      ag[k].id) - atom_ids.begin();
+          atomic_gradients[a] += coeff * ag[k].grad;
+        }
+      }
+    }
+  }
+
+  if (hb.size()) {
+
+    cvm::real const hb_norm = hb_coeff / cvm::real(hb.size());
+
+    for (size_t i = 0; i < hb.size(); i++) {
+      // Coeficient of this CVC's gradient in the colvar gradient, times coefficient of this
+      // hbond's gradient in the CVC's gradient
+      cvm::real const coeff = cvc_coeff * 0.5 * hb_norm;
+
+      for (size_t j = 0; j < hb[i]->atom_groups.size(); j++) {
+        cvm::atom_group &ag = *(hb[i]->atom_groups[j]);
+        for (size_t k = 0; k < ag.size(); k++) {
+          size_t a = std::lower_bound(atom_ids.begin(), atom_ids.end(),
+                                      ag[k].id) - atom_ids.begin();
+          atomic_gradients[a] += coeff * ag[k].grad;
+        }
+      }
+    }
+  }
+}
+
+
 void colvar::alpha_angles::apply_force(colvarvalue const &force)
 {
 
@@ -242,7 +294,8 @@ colvar::dihedPC::dihedPC(std::string const &conf)
     cvm::log("Initializing dihedral PC object.\n");
 
   function_type = "dihedPC";
-  enable(f_cvc_implicit_gradient);
+  // Supported through references to atom groups of children cvcs
+  enable(f_cvc_explicit_gradient);
   x.type(colvarvalue::type_scalar);
 
   std::string segment_id;
@@ -372,7 +425,8 @@ colvar::dihedPC::dihedPC()
   : cvc()
 {
   function_type = "dihedPC";
-  enable(f_cvc_implicit_gradient);
+  // Supported through references to atom groups of children cvcs
+  enable(f_cvc_explicit_gradient);
   x.type(colvarvalue::type_scalar);
 }
 
@@ -383,6 +437,8 @@ colvar::dihedPC::~dihedPC()
     delete theta.back();
     theta.pop_back();
   }
+  // Our references to atom groups have become invalid now that children cvcs are deleted
+  atom_groups.clear();
 }
 
 
@@ -392,8 +448,8 @@ void colvar::dihedPC::calc_value()
   for (size_t i = 0; i < theta.size(); i++) {
     theta[i]->calc_value();
     cvm::real const t = (PI / 180.) * theta[i]->value().real_value;
-    x.real_value += coeffs[2*i  ] * std::cos(t)
+    x.real_value += coeffs[2*i  ] * cvm::cos(t)
-                  + coeffs[2*i+1] * std::sin(t);
+                  + coeffs[2*i+1] * cvm::sin(t);
   }
 }
 
@@ -406,12 +462,35 @@ void colvar::dihedPC::calc_gradients()
 }
 
 
+void colvar::dihedPC::collect_gradients(std::vector<int> const &atom_ids, std::vector<cvm::rvector> &atomic_gradients)
+{
+  cvm::real cvc_coeff = sup_coeff * cvm::real(sup_np) * cvm::integer_power(value().real_value, sup_np-1);
+  for (size_t i = 0; i < theta.size(); i++) {
+    cvm::real const t = (PI / 180.) * theta[i]->value().real_value;
+    cvm::real const dcosdt = - (PI / 180.) * cvm::sin(t);
+    cvm::real const dsindt =   (PI / 180.) * cvm::cos(t);
+    // Coeficient of this dihedPC's gradient in the colvar gradient, times coefficient of this
+    // dihedral's gradient in the dihedPC's gradient
+    cvm::real const coeff = cvc_coeff * (coeffs[2*i] * dcosdt + coeffs[2*i+1] * dsindt);
+
+    for (size_t j = 0; j < theta[i]->atom_groups.size(); j++) {
+      cvm::atom_group &ag = *(theta[i]->atom_groups[j]);
+      for (size_t k = 0; k < ag.size(); k++) {
+        size_t a = std::lower_bound(atom_ids.begin(), atom_ids.end(),
+                                    ag[k].id) - atom_ids.begin();
+        atomic_gradients[a] += coeff * ag[k].grad;
+      }
+    }
+  }
+}
+
+
 void colvar::dihedPC::apply_force(colvarvalue const &force)
 {
   for (size_t i = 0; i < theta.size(); i++) {
     cvm::real const t = (PI / 180.) * theta[i]->value().real_value;
-    cvm::real const dcosdt = - (PI / 180.) * std::sin(t);
+    cvm::real const dcosdt = - (PI / 180.) * cvm::sin(t);
-    cvm::real const dsindt =   (PI / 180.) * std::cos(t);
+    cvm::real const dsindt =   (PI / 180.) * cvm::cos(t);
 
     theta[i]->apply_force((coeffs[2*i  ] * dcosdt +
                            coeffs[2*i+1] * dsindt) * force);

lib/colvars/colvarcomp_rotations.cpp

@@ -7,8 +7,6 @@
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 
-#include <cmath>
-
 #include "colvarmodule.h"
 #include "colvarvalue.h"
 #include "colvarparse.h"
@@ -18,21 +16,27 @@
 
 
 colvar::orientation::orientation(std::string const &conf)
-  : cvc(conf)
+  : cvc()
 {
   function_type = "orientation";
-  atoms = parse_group(conf, "atoms");
+  disable(f_cvc_explicit_gradient);
-  enable(f_cvc_implicit_gradient);
   x.type(colvarvalue::type_quaternion);
+  init(conf);
+}
 
+
+int colvar::orientation::init(std::string const &conf)
+{
+  int error_code = cvc::init(conf);
+
+  atoms = parse_group(conf, "atoms");
   ref_pos.reserve(atoms->size());
 
   if (get_keyval(conf, "refPositions", ref_pos, ref_pos)) {
     cvm::log("Using reference positions from input file.\n");
     if (ref_pos.size() != atoms->size()) {
-      cvm::error("Error: reference positions do not "
+      return cvm::error("Error: reference positions do not "
-                        "match the number of requested atoms.\n");
+                        "match the number of requested atoms.\n", INPUT_ERROR);
-      return;
     }
   }
 
@@ -46,9 +50,8 @@ colvar::orientation::orientation(std::string const &conf)
         // use PDB flags if column is provided
         bool found = get_keyval(conf, "refPositionsColValue", file_col_value, 0.0);
         if (found && file_col_value==0.0) {
-          cvm::error("Error: refPositionsColValue, "
+          return cvm::error("Error: refPositionsColValue, "
-                            "if provided, must be non-zero.\n");
+                            "if provided, must be non-zero.\n", INPUT_ERROR);
-          return;
         }
       }
 
@@ -59,9 +62,8 @@ colvar::orientation::orientation(std::string const &conf)
   }
 
   if (!ref_pos.size()) {
-    cvm::error("Error: must define a set of "
+    return cvm::error("Error: must define a set of "
-                      "reference coordinates.\n");
+                      "reference coordinates.\n", INPUT_ERROR);
-    return;
   }
 
 
@@ -85,6 +87,7 @@ colvar::orientation::orientation(std::string const &conf)
     rot.request_group2_gradients(atoms->size());
   }
 
+  return error_code;
 }
 
 
@@ -92,7 +95,7 @@ colvar::orientation::orientation()
   : cvc()
 {
   function_type = "orientation";
-  enable(f_cvc_implicit_gradient);
+  disable(f_cvc_explicit_gradient);
   x.type(colvarvalue::type_quaternion);
 }
 
@@ -158,10 +161,18 @@ colvarvalue colvar::orientation::dist2_rgrad(colvarvalue const &x1,
 
 
 colvar::orientation_angle::orientation_angle(std::string const &conf)
-  : orientation(conf)
+  : orientation()
 {
   function_type = "orientation_angle";
+  enable(f_cvc_explicit_gradient);
   x.type(colvarvalue::type_scalar);
+  init(conf);
+}
+
+
+int colvar::orientation_angle::init(std::string const &conf)
+{
+  return orientation::init(conf);
 }
 
 
@@ -169,6 +180,7 @@ colvar::orientation_angle::orientation_angle()
   : orientation()
 {
   function_type = "orientation_angle";
+  enable(f_cvc_explicit_gradient);
   x.type(colvarvalue::type_scalar);
 }
 
@@ -180,9 +192,9 @@ void colvar::orientation_angle::calc_value()
   rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog));
 
   if ((rot.q).q0 >= 0.0) {
-    x.real_value = (180.0/PI) * 2.0 * std::acos((rot.q).q0);
+    x.real_value = (180.0/PI) * 2.0 * cvm::acos((rot.q).q0);
   } else {
-    x.real_value = (180.0/PI) * 2.0 * std::acos(-1.0 * (rot.q).q0);
+    x.real_value = (180.0/PI) * 2.0 * cvm::acos(-1.0 * (rot.q).q0);
   }
 }
 
@@ -191,7 +203,7 @@ void colvar::orientation_angle::calc_gradients()
 {
   cvm::real const dxdq0 =
     ( ((rot.q).q0 * (rot.q).q0 < 1.0) ?
-      ((180.0 / PI) * (-2.0) / std::sqrt(1.0 - ((rot.q).q0 * (rot.q).q0))) :
+      ((180.0 / PI) * (-2.0) / cvm::sqrt(1.0 - ((rot.q).q0 * (rot.q).q0))) :
       0.0 );
 
   for (size_t ia = 0; ia < atoms->size(); ia++) {
@@ -214,10 +226,18 @@ simple_scalar_dist_functions(orientation_angle)
 
 
 colvar::orientation_proj::orientation_proj(std::string const &conf)
-  : orientation(conf)
+  : orientation()
 {
   function_type = "orientation_proj";
+  enable(f_cvc_explicit_gradient);
   x.type(colvarvalue::type_scalar);
+  init(conf);
+}
+
+
+int colvar::orientation_proj::init(std::string const &conf)
+{
+  return orientation::init(conf);
 }
 
 
@@ -225,6 +245,7 @@ colvar::orientation_proj::orientation_proj()
   : orientation()
 {
   function_type = "orientation_proj";
+  enable(f_cvc_explicit_gradient);
   x.type(colvarvalue::type_scalar);
 }
 
@@ -261,18 +282,28 @@ simple_scalar_dist_functions(orientation_proj)
 
 
 colvar::tilt::tilt(std::string const &conf)
-  : orientation(conf)
+  : orientation()
 {
   function_type = "tilt";
+  enable(f_cvc_explicit_gradient);
+  x.type(colvarvalue::type_scalar);
+  init(conf);
+}
+
+
+int colvar::tilt::init(std::string const &conf)
+{
+  int error_code = COLVARS_OK;
+
+  error_code |= orientation::init(conf);
 
   get_keyval(conf, "axis", axis, cvm::rvector(0.0, 0.0, 1.0));
-
   if (axis.norm2() != 1.0) {
     axis /= axis.norm();
     cvm::log("Normalizing rotation axis to "+cvm::to_str(axis)+".\n");
   }
 
-  x.type(colvarvalue::type_scalar);
+  return error_code;
 }
 
 
@@ -280,6 +311,7 @@ colvar::tilt::tilt()
   : orientation()
 {
   function_type = "tilt";
+  enable(f_cvc_explicit_gradient);
   x.type(colvarvalue::type_scalar);
 }
 
@@ -322,20 +354,30 @@ simple_scalar_dist_functions(tilt)
 
 
 colvar::spin_angle::spin_angle(std::string const &conf)
-  : orientation(conf)
+  : orientation()
 {
   function_type = "spin_angle";
+  period = 360.0;
+  b_periodic = true;
+  enable(f_cvc_explicit_gradient);
+  x.type(colvarvalue::type_scalar);
+  init(conf);
+}
+
+
+int colvar::spin_angle::init(std::string const &conf)
+{
+  int error_code = COLVARS_OK;
+
+  error_code |= orientation::init(conf);
 
   get_keyval(conf, "axis", axis, cvm::rvector(0.0, 0.0, 1.0));
-
   if (axis.norm2() != 1.0) {
     axis /= axis.norm();
     cvm::log("Normalizing rotation axis to "+cvm::to_str(axis)+".\n");
   }
 
-  period = 360.0;
+  return error_code;
-  b_periodic = true;
-  x.type(colvarvalue::type_scalar);
 }
 
 
@@ -345,6 +387,7 @@ colvar::spin_angle::spin_angle()
   function_type = "spin_angle";
   period = 360.0;
   b_periodic = true;
+  enable(f_cvc_explicit_gradient);
   x.type(colvarvalue::type_scalar);
 }
 
@@ -410,15 +453,15 @@ colvarvalue colvar::spin_angle::dist2_rgrad(colvarvalue const &x1,
 }
 
 
-void colvar::spin_angle::wrap(colvarvalue &x) const
+void colvar::spin_angle::wrap(colvarvalue &x_unwrapped) const
 {
-  if ((x.real_value - wrap_center) >= 180.0) {
+  if ((x_unwrapped.real_value - wrap_center) >= 180.0) {
-    x.real_value -= 360.0;
+    x_unwrapped.real_value -= 360.0;
     return;
   }
 
-  if ((x.real_value - wrap_center) < -180.0) {
+  if ((x_unwrapped.real_value - wrap_center) < -180.0) {
-    x.real_value += 360.0;
+    x_unwrapped.real_value += 360.0;
     return;
   }
 

lib/colvars/colvardeps.cpp

@@ -122,12 +122,7 @@ bool colvardeps::get_keyval_feature(colvarparse *cvp,
 
 int colvardeps::enable(int feature_id,
                        bool dry_run /* default: false */,
-                       // dry_run: fail silently, do not enable if available
-                       // flag is passed recursively to deps of this feature
                        bool toplevel /* default: true */)
-// toplevel: false if this is called as part of a chain of dependency resolution
-// this is used to diagnose failed dependencies by displaying the full stack
-// only the toplevel dependency will throw a fatal error
 {
   int res;
   size_t i, j;
@@ -143,8 +138,7 @@ int colvardeps::enable(int feature_id,
 
   if (fs->enabled) {
     if (!(dry_run || toplevel)) {
-      // This is a dependency
+      // This is a dependency: prevent disabling this feature as long
-      // Prevent disabling this feature as long
       // as requirement is enabled
       fs->ref_count++;
       if (cvm::debug())
@@ -173,7 +167,10 @@ int colvardeps::enable(int feature_id,
   if (!toplevel && !is_dynamic(feature_id)) {
     if (!dry_run) {
       cvm::log(feature_type_descr + " feature \"" + f->description
-        + "\" may not be enabled as a dependency in " + description + ".\n");
+        + "\" cannot be enabled automatically in " + description + ".");
+      if (is_user(feature_id)) {
+        cvm::log("Try setting it manually.\n");
+      }
     }
     return COLVARS_ERROR;
   }
@@ -354,6 +351,7 @@ int colvardeps::disable(int feature_id) {
   return COLVARS_OK;
 }
 
+
 int colvardeps::decr_ref_count(int feature_id) {
   int &rc = feature_states[feature_id].ref_count;
   feature *f = features()[feature_id];
@@ -379,324 +377,52 @@ int colvardeps::decr_ref_count(int feature_id) {
   return COLVARS_OK;
 }
 
-void colvardeps::init_feature(int feature_id, const char *description, feature_type type) {
+
-  modify_features()[feature_id]->description = description;
+void colvardeps::init_feature(int feature_id, const char *description_in, feature_type type) {
+  modify_features()[feature_id]->description = description_in;
   modify_features()[feature_id]->type = type;
 }
 
-// Shorthand macros for describing dependencies
-#define f_req_self(f, g) features()[f]->requires_self.push_back(g)
-// This macro ensures that exclusions are symmetric
-#define f_req_exclude(f, g) features()[f]->requires_exclude.push_back(g); \
-                            features()[g]->requires_exclude.push_back(f)
-#define f_req_children(f, g) features()[f]->requires_children.push_back(g)
-#define f_req_alt2(f, g, h) features()[f]->requires_alt.push_back(std::vector<int>(2));\
-  features()[f]->requires_alt.back()[0] = g;                                           \
-  features()[f]->requires_alt.back()[1] = h
-#define f_req_alt3(f, g, h, i) features()[f]->requires_alt.push_back(std::vector<int>(3));\
-  features()[f]->requires_alt.back()[0] = g;                                           \
-  features()[f]->requires_alt.back()[1] = h;                                           \
-  features()[f]->requires_alt.back()[2] = i
-#define f_req_alt4(f, g, h, i, j) features()[f]->requires_alt.push_back(std::vector<int>(4));\
-  features()[f]->requires_alt.back()[0] = g;                                           \
-  features()[f]->requires_alt.back()[1] = h;                                           \
-  features()[f]->requires_alt.back()[2] = i;                                           \
-  features()[f]->requires_alt.back()[3] = j
 
-void colvardeps::init_cvb_requires() {
+// Shorthand functions for describing dependencies
-  int i;
+void colvardeps::require_feature_self(int f, int g) {
-  if (features().size() == 0) {
+  features()[f]->requires_self.push_back(g);
-    for (i = 0; i < f_cvb_ntot; i++) {
-      modify_features().push_back(new feature);
-    }
-
-    init_feature(f_cvb_active, "active", f_type_dynamic);
-    f_req_children(f_cvb_active, f_cv_active);
-
-    init_feature(f_cvb_awake, "awake", f_type_static);
-    f_req_self(f_cvb_awake, f_cvb_active);
-
-    init_feature(f_cvb_apply_force, "apply force", f_type_user);
-    f_req_children(f_cvb_apply_force, f_cv_gradient);
-
-    init_feature(f_cvb_get_total_force, "obtain total force", f_type_dynamic);
-    f_req_children(f_cvb_get_total_force, f_cv_total_force);
-
-    init_feature(f_cvb_output_acc_work, "output accumulated work", f_type_user);
-    f_req_self(f_cvb_output_acc_work, f_cvb_apply_force);
-
-    init_feature(f_cvb_history_dependent, "history-dependent", f_type_static);
-
-    init_feature(f_cvb_time_dependent, "time-dependent", f_type_static);
-
-    init_feature(f_cvb_scalar_variables, "require scalar variables", f_type_static);
-    f_req_children(f_cvb_scalar_variables, f_cv_scalar);
-
-    init_feature(f_cvb_calc_pmf, "calculate a PMF", f_type_static);
-
-    init_feature(f_cvb_calc_ti_samples, "calculate TI samples", f_type_dynamic);
-    f_req_self(f_cvb_calc_ti_samples, f_cvb_get_total_force);
-    f_req_children(f_cvb_calc_ti_samples, f_cv_grid);
-
-    init_feature(f_cvb_write_ti_samples, "write TI samples ", f_type_user);
-    f_req_self(f_cvb_write_ti_samples, f_cvb_calc_ti_samples);
-
-    init_feature(f_cvb_write_ti_pmf, "write TI PMF", f_type_user);
-    f_req_self(f_cvb_write_ti_pmf, f_cvb_calc_ti_samples);
-  }
-
-  // Initialize feature_states for each instance
-  feature_states.reserve(f_cvb_ntot);
-  for (i = 0; i < f_cvb_ntot; i++) {
-    feature_states.push_back(feature_state(true, false));
-    // Most features are available, so we set them so
-    // and list exceptions below
-  }
-
-  // only compute TI samples when deriving from colvarbias_ti
-  feature_states[f_cvb_calc_ti_samples].available = false;
 }
 
 
-void colvardeps::init_cv_requires() {
+// Ensure that exclusions are symmetric
-  size_t i;
+void colvardeps::exclude_feature_self(int f, int g) {
-  if (features().size() == 0) {
+  features()[f]->requires_exclude.push_back(g);
-    for (i = 0; i < f_cv_ntot; i++) {
+  features()[g]->requires_exclude.push_back(f);
-      modify_features().push_back(new feature);
-    }
-
-    init_feature(f_cv_active, "active", f_type_dynamic);
-    // Do not require f_cvc_active in children, as some components may be disabled
-    // Colvars must be either a linear combination, or scalar (and polynomial) or scripted/custom
-    f_req_alt4(f_cv_active, f_cv_scalar, f_cv_linear, f_cv_scripted, f_cv_custom_function);
-
-    init_feature(f_cv_awake, "awake", f_type_static);
-    f_req_self(f_cv_awake, f_cv_active);
-
-    init_feature(f_cv_gradient, "gradient", f_type_dynamic);
-    f_req_children(f_cv_gradient, f_cvc_gradient);
-
-    init_feature(f_cv_collect_gradient, "collect gradient", f_type_dynamic);
-    f_req_self(f_cv_collect_gradient, f_cv_gradient);
-    f_req_self(f_cv_collect_gradient, f_cv_scalar);
-    // The following exlusion could be lifted by implementing the feature
-    f_req_exclude(f_cv_collect_gradient, f_cv_scripted);
-
-    init_feature(f_cv_fdiff_velocity, "velocity from finite differences", f_type_dynamic);
-
-    // System force: either trivial (spring force); through extended Lagrangian, or calculated explicitly
-    init_feature(f_cv_total_force, "total force", f_type_dynamic);
-    f_req_alt2(f_cv_total_force, f_cv_extended_Lagrangian, f_cv_total_force_calc);
-
-    // Deps for explicit total force calculation
-    init_feature(f_cv_total_force_calc, "total force calculation", f_type_dynamic);
-    f_req_self(f_cv_total_force_calc, f_cv_scalar);
-    f_req_self(f_cv_total_force_calc, f_cv_linear);
-    f_req_children(f_cv_total_force_calc, f_cvc_inv_gradient);
-    f_req_self(f_cv_total_force_calc, f_cv_Jacobian);
-
-    init_feature(f_cv_Jacobian, "Jacobian derivative", f_type_dynamic);
-    f_req_self(f_cv_Jacobian, f_cv_scalar);
-    f_req_self(f_cv_Jacobian, f_cv_linear);
-    f_req_children(f_cv_Jacobian, f_cvc_Jacobian);
-
-    init_feature(f_cv_hide_Jacobian, "hide Jacobian force", f_type_user);
-    f_req_self(f_cv_hide_Jacobian, f_cv_Jacobian); // can only hide if calculated
-
-    init_feature(f_cv_extended_Lagrangian, "extended Lagrangian", f_type_user);
-    f_req_self(f_cv_extended_Lagrangian, f_cv_scalar);
-    f_req_self(f_cv_extended_Lagrangian, f_cv_gradient);
-
-    init_feature(f_cv_Langevin, "Langevin dynamics", f_type_user);
-    f_req_self(f_cv_Langevin, f_cv_extended_Lagrangian);
-
-    init_feature(f_cv_linear, "linear", f_type_static);
-
-    init_feature(f_cv_scalar, "scalar", f_type_static);
-
-    init_feature(f_cv_output_energy, "output energy", f_type_user);
-
-    init_feature(f_cv_output_value, "output value", f_type_user);
-
-    init_feature(f_cv_output_velocity, "output velocity", f_type_user);
-    f_req_self(f_cv_output_velocity, f_cv_fdiff_velocity);
-
-    init_feature(f_cv_output_applied_force, "output applied force", f_type_user);
-
-    init_feature(f_cv_output_total_force, "output total force", f_type_user);
-    f_req_self(f_cv_output_total_force, f_cv_total_force);
-
-    init_feature(f_cv_subtract_applied_force, "subtract applied force from total force", f_type_user);
-    f_req_self(f_cv_subtract_applied_force, f_cv_total_force);
-
-    init_feature(f_cv_lower_boundary, "lower boundary", f_type_user);
-    f_req_self(f_cv_lower_boundary, f_cv_scalar);
-
-    init_feature(f_cv_upper_boundary, "upper boundary", f_type_user);
-    f_req_self(f_cv_upper_boundary, f_cv_scalar);
-
-    init_feature(f_cv_grid, "grid", f_type_dynamic);
-    f_req_self(f_cv_grid, f_cv_lower_boundary);
-    f_req_self(f_cv_grid, f_cv_upper_boundary);
-
-    init_feature(f_cv_runave, "running average", f_type_user);
-
-    init_feature(f_cv_corrfunc, "correlation function", f_type_user);
-
-    init_feature(f_cv_scripted, "scripted", f_type_user);
-
-    init_feature(f_cv_custom_function, "custom function", f_type_user);
-    f_req_exclude(f_cv_custom_function, f_cv_scripted);
-
-    init_feature(f_cv_periodic, "periodic", f_type_static);
-    f_req_self(f_cv_periodic, f_cv_scalar);
-    init_feature(f_cv_scalar, "scalar", f_type_static);
-    init_feature(f_cv_linear, "linear", f_type_static);
-    init_feature(f_cv_homogeneous, "homogeneous", f_type_static);
-
-    // because total forces are obtained from the previous time step,
-    // we cannot (currently) have colvar values and total forces for the same timestep
-    init_feature(f_cv_multiple_ts, "multiple timestep colvar");
-    f_req_exclude(f_cv_multiple_ts, f_cv_total_force_calc);
-  }
-
-  // Initialize feature_states for each instance
-  feature_states.reserve(f_cv_ntot);
-  for (i = 0; i < f_cv_ntot; i++) {
-    feature_states.push_back(feature_state(true, false));
-    // Most features are available, so we set them so
-    // and list exceptions below
-   }
-
-  feature_states[f_cv_fdiff_velocity].available =
-    cvm::main()->proxy->simulation_running();
 }
 
 
-void colvardeps::init_cvc_requires() {
+void colvardeps::require_feature_children(int f, int g) {
-  size_t i;
+  features()[f]->requires_children.push_back(g);
-  // Initialize static array once and for all
-  if (features().size() == 0) {
-    for (i = 0; i < colvardeps::f_cvc_ntot; i++) {
-      modify_features().push_back(new feature);
-    }
-
-    init_feature(f_cvc_active, "active", f_type_dynamic);
-//     The dependency below may become useful if we use dynamic atom groups
-//     f_req_children(f_cvc_active, f_ag_active);
-
-    init_feature(f_cvc_scalar, "scalar", f_type_static);
-
-    init_feature(f_cvc_gradient, "gradient", f_type_dynamic);
-
-    init_feature(f_cvc_implicit_gradient, "implicit gradient", f_type_static);
-    f_req_children(f_cvc_implicit_gradient, f_ag_implicit_gradient);
-
-    init_feature(f_cvc_inv_gradient, "inverse gradient", f_type_dynamic);
-    f_req_self(f_cvc_inv_gradient, f_cvc_gradient);
-
-    init_feature(f_cvc_debug_gradient, "debug gradient", f_type_user);
-    f_req_self(f_cvc_debug_gradient, f_cvc_gradient);
-    f_req_exclude(f_cvc_debug_gradient, f_cvc_implicit_gradient);
-
-    init_feature(f_cvc_Jacobian, "Jacobian derivative", f_type_dynamic);
-    f_req_self(f_cvc_Jacobian, f_cvc_inv_gradient);
-
-    init_feature(f_cvc_com_based, "depends on group centers of mass", f_type_static);
-
-    // init_feature(f_cvc_pbc_minimum_image, "use minimum-image distances with PBCs", f_type_user);
-
-    // Compute total force on first site only to avoid unwanted
-    // coupling to other colvars (see e.g. Ciccotti et al., 2005)
-    init_feature(f_cvc_one_site_total_force, "compute total force from one group", f_type_user);
-    f_req_self(f_cvc_one_site_total_force, f_cvc_com_based);
-
-    init_feature(f_cvc_scalable, "scalable calculation", f_type_static);
-    f_req_self(f_cvc_scalable, f_cvc_scalable_com);
-
-    init_feature(f_cvc_scalable_com, "scalable calculation of centers of mass", f_type_static);
-    f_req_self(f_cvc_scalable_com, f_cvc_com_based);
-
-
-    // TODO only enable this when f_ag_scalable can be turned on for a pre-initialized group
-    // f_req_children(f_cvc_scalable, f_ag_scalable);
-    // f_req_children(f_cvc_scalable_com, f_ag_scalable_com);
-  }
-
-  // Initialize feature_states for each instance
-  // default as available, not enabled
-  // except dynamic features which default as unavailable
-  feature_states.reserve(f_cvc_ntot);
-  for (i = 0; i < colvardeps::f_cvc_ntot; i++) {
-    bool avail = is_dynamic(i) ? false : true;
-    feature_states.push_back(feature_state(avail, false));
-  }
-
-  // CVCs are enabled from the start - get disabled based on flags
-  feature_states[f_cvc_active].enabled = true;
-
-  // Features that are implemented by all cvcs by default
-  // Each cvc specifies what other features are available
-  feature_states[f_cvc_active].available = true;
-  feature_states[f_cvc_gradient].available = true;
-
-  // Use minimum-image distances by default
-  feature_states[f_cvc_pbc_minimum_image].enabled = true;
-
-  // Features that are implemented by default if their requirements are
-  feature_states[f_cvc_one_site_total_force].available = true;
-
-  // Features That are implemented only for certain simulation engine configurations
-  feature_states[f_cvc_scalable_com].available = (cvm::proxy->scalable_group_coms() == COLVARS_OK);
-  feature_states[f_cvc_scalable].available = feature_states[f_cvc_scalable_com].available;
 }
 
 
-void colvardeps::init_ag_requires() {
+void colvardeps::require_feature_alt(int f, int g, int h) {
-  size_t i;
+  features()[f]->requires_alt.push_back(std::vector<int>(2));
-  // Initialize static array once and for all
+  features()[f]->requires_alt.back()[0] = g;
-  if (features().size() == 0) {
+  features()[f]->requires_alt.back()[1] = h;
-    for (i = 0; i < f_ag_ntot; i++) {
+}
-      modify_features().push_back(new feature);
-    }
 
-    init_feature(f_ag_active, "active", f_type_dynamic);
-    init_feature(f_ag_center, "translational fit", f_type_static);
-    init_feature(f_ag_rotate, "rotational fit", f_type_static);
-    init_feature(f_ag_fitting_group, "reference positions group", f_type_static);
-    init_feature(f_ag_implicit_gradient, "implicit atom gradient", f_type_dynamic);
-    init_feature(f_ag_fit_gradients, "fit gradients", f_type_user);
-    f_req_exclude(f_ag_fit_gradients, f_ag_implicit_gradient);
 
-    init_feature(f_ag_atom_forces, "atomic forces", f_type_dynamic);
+void colvardeps::require_feature_alt(int f, int g, int h, int i) {
+  features()[f]->requires_alt.push_back(std::vector<int>(3));
+  features()[f]->requires_alt.back()[0] = g;
+  features()[f]->requires_alt.back()[1] = h;
+  features()[f]->requires_alt.back()[2] = i;
+}
 
-    // parallel calculation implies that we have at least a scalable center of mass,
-    // but f_ag_scalable is kept as a separate feature to deal with future dependencies
-    init_feature(f_ag_scalable, "scalable group calculation", f_type_static);
-    init_feature(f_ag_scalable_com, "scalable group center of mass calculation", f_type_static);
-    f_req_self(f_ag_scalable, f_ag_scalable_com);
 
-//     init_feature(f_ag_min_msd_fit, "minimum MSD fit")
+void colvardeps::require_feature_alt(int f, int g, int h, int i, int j) {
-//     f_req_self(f_ag_min_msd_fit, f_ag_center)
+  features()[f]->requires_alt.push_back(std::vector<int>(4));
-//     f_req_self(f_ag_min_msd_fit, f_ag_rotate)
+  features()[f]->requires_alt.back()[0] = g;
-//     f_req_exclude(f_ag_min_msd_fit, f_ag_fitting_group)
+  features()[f]->requires_alt.back()[1] = h;
-  }
+  features()[f]->requires_alt.back()[2] = i;
-
+  features()[f]->requires_alt.back()[3] = j;
-  // Initialize feature_states for each instance
-  // default as unavailable, not enabled
-  feature_states.reserve(f_ag_ntot);
-  for (i = 0; i < colvardeps::f_ag_ntot; i++) {
-    feature_states.push_back(feature_state(false, false));
-  }
-
-  // Features that are implemented (or not) by all atom groups
-  feature_states[f_ag_active].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)
-  feature_states[f_ag_scalable].available = true;
-  feature_states[f_ag_fit_gradients].available = true;
-  feature_states[f_ag_implicit_gradient].available = true;
 }
 
 
@@ -720,7 +446,7 @@ void colvardeps::print_state() {
 void colvardeps::add_child(colvardeps *child) {
 
   children.push_back(child);
-  child->parents.push_back((colvardeps *)this);
+  child->parents.push_back(this);
 
   // Solve dependencies of already enabled parent features
   // in the new child

lib/colvars/colvardeps.h

@@ -23,7 +23,11 @@
 /// 3. Static features are static properties of the object, determined
 ///   programatically at initialization time.
 ///
-/// In all classes, feature 0 is active. When an object is inactivated
+/// The following diagram summarizes the dependency tree at the bias, colvar, and colvarcomp levels.
+/// Isolated and atom group features are not shown to save space.
+/// @image html deps_2019.svg
+///
+/// In all classes, feature 0 is `active`. When an object is inactivated
 /// all its children dependencies are dereferenced (free_children_deps)
 /// While the object is inactive, no dependency solving is done on children
 /// it is done when the object is activated back (restore_children_deps)
@@ -72,7 +76,6 @@ protected:
   /// Unused by lower-level objects (cvcs and atom groups)
   int   time_step_factor;
 
-private:
   /// List of the states of all features
   std::vector<feature_state> feature_states;
 
@@ -89,14 +92,14 @@ public:
   inline int get_time_step_factor() const {return time_step_factor;}
 
   /// Pair a numerical feature ID with a description and type
-  void init_feature(int feature_id, const char *description, feature_type type = f_type_not_set);
+  void init_feature(int feature_id, const char *description, feature_type type);
 
   /// Describes a feature and its dependencies
   /// used in a static array within each subclass
   class feature {
 
   public:
-    feature() {}
+    feature() : type(f_type_not_set) {}
     ~feature() {}
 
     std::string description; // Set by derived object initializer
@@ -126,6 +129,7 @@ public:
     feature_type type;
   };
 
+  inline bool is_not_set(int id) { return features()[id]->type == f_type_not_set; }
   inline bool is_dynamic(int id) { return features()[id]->type == f_type_dynamic; }
   inline bool is_static(int id) { return features()[id]->type == f_type_static; }
   inline bool is_user(int id) { return features()[id]->type == f_type_user; }
@@ -135,7 +139,7 @@ public:
   // with a non-static array
   // Intermediate classes (colvarbias and colvarcomp, which are also base classes)
   // implement this as virtual to allow overriding
-  virtual const std::vector<feature *>&features() = 0;
+  virtual const std::vector<feature *> &features() const = 0;
   virtual std::vector<feature *>&modify_features() = 0;
 
   void add_child(colvardeps *child);
@@ -188,10 +192,12 @@ protected:
 
 public:
 
-  /// enable a feature and recursively solve its dependencies
+  /// Enable a feature and recursively solve its dependencies.
-  /// for proper reference counting, one should not add
+  /// For accurate reference counting, do not add spurious calls to enable()
-  /// spurious calls to enable()
+  /// \param dry_run Recursively test if a feature is available, without enabling it
-  /// dry_run is set to true to recursively test if a feature is available, without enabling it
+  /// \param toplevel False if this is called as part of a chain of dependency resolution.
+  /// This is used to diagnose failed dependencies by displaying the full stack:
+  /// only the toplevel dependency will throw a fatal error.
   int enable(int f, bool dry_run = false, bool toplevel = true);
 
   /// Disable a feature, decrease the reference count of its dependencies
@@ -318,8 +324,8 @@ public:
     f_cvc_active,
     f_cvc_scalar,
     f_cvc_gradient,
-    /// \brief CVC doesn't calculate and store explicit atom gradients
+    /// \brief CVC calculates and stores explicit atom gradients
-    f_cvc_implicit_gradient,
+    f_cvc_explicit_gradient,
     f_cvc_inv_gradient,
     /// \brief If enabled, calc_gradients() will call debug_gradients() for every group needed
     f_cvc_debug_gradient,
@@ -341,7 +347,7 @@ public:
     /// ie. not using refpositionsgroup
 //     f_ag_min_msd_fit,
     /// \brief Does not have explicit atom gradients from parent CVC
-    f_ag_implicit_gradient,
+    f_ag_explicit_gradient,
     f_ag_fit_gradients,
     f_ag_atom_forces,
     f_ag_scalable,
@@ -349,13 +355,39 @@ public:
     f_ag_ntot
   };
 
-  void init_cvb_requires();
+  /// Initialize dependency tree for object of a derived class
-  void init_cv_requires();
+  virtual int init_dependencies() = 0;
-  void init_cvc_requires();
+
-  void init_ag_requires();
+  /// Make feature f require feature g within the same object
+  void require_feature_self(int f, int g);
+
+  /// Make features f and g mutually exclusive within the same object
+  void exclude_feature_self(int f, int g);
+
+  /// Make feature f require feature g within children
+  void require_feature_children(int f, int g);
+
+  /// Make feature f require either g or h within the same object
+  void require_feature_alt(int f, int g, int h);
+
+  /// Make feature f require any of g, h, or i within the same object
+  void require_feature_alt(int f, int g, int h, int i);
+
+  /// Make feature f require any of g, h, i, or j within the same object
+  void require_feature_alt(int f, int g, int h, int i, int j);
 
   /// \brief print all enabled features and those of children, for debugging
   void print_state();
+
+  /// \brief Check that a feature is enabled, raising BUG_ERROR if not
+  inline void check_enabled(int f, std::string const &reason) const
+  {
+    if (! is_enabled(f)) {
+      cvm::error("Error: "+reason+" requires that the feature \""+
+                 features()[f]->description+"\" is active.\n", BUG_ERROR);
+    }
+  }
+
 };
 
 #endif

lib/colvars/colvargrid.cpp

@@ -109,7 +109,7 @@ cvm::real colvar_grid_scalar::entropy() const
 {
   cvm::real sum = 0.0;
   for (size_t i = 0; i < nt; i++) {
-    sum += -1.0 * data[i] * std::log(data[i]);
+    sum += -1.0 * data[i] * cvm::logn(data[i]);
   }
   cvm::real bin_volume = 1.0;
   for (size_t id = 0; id < widths.size(); id++) {

lib/colvars/colvargrid.h

@@ -12,7 +12,6 @@
 
 #include <iostream>
 #include <iomanip>
-#include <cmath>
 
 #include "colvar.h"
 #include "colvarmodule.h"
@@ -53,7 +52,7 @@ protected:
   std::vector<colvar *> cv;
 
   /// Do we request actual value (for extended-system colvars)?
-  std::vector<bool> actual_value;
+  std::vector<bool> use_actual_value;
 
   /// Get the low-level index corresponding to an index
   inline size_t address(std::vector<int> const &ix) const
@@ -136,8 +135,8 @@ public:
   inline void request_actual_value(bool b = true)
   {
     size_t i;
-    for (i = 0; i < actual_value.size(); i++) {
+    for (i = 0; i < use_actual_value.size(); i++) {
-      actual_value[i] = b;
+      use_actual_value[i] = b;
     }
   }
 
@@ -215,7 +214,7 @@ public:
                                          mult(g.mult),
                                          data(),
                                          cv(g.cv),
-                                         actual_value(g.actual_value),
+                                         use_actual_value(g.use_actual_value),
                                          lower_boundaries(g.lower_boundaries),
                                          upper_boundaries(g.upper_boundaries),
                                          periodic(g.periodic),
@@ -290,13 +289,13 @@ public:
       periodic.push_back(cv[i]->periodic_boundaries());
 
       // By default, get reported colvar value (for extended Lagrangian colvars)
-      actual_value.push_back(false);
+      use_actual_value.push_back(false);
 
       // except if a colvar is specified twice in a row
       // then the first instance is the actual value
       // For histograms of extended-system coordinates
       if (i > 0 && cv[i-1] == cv[i]) {
-        actual_value[i-1] = true;
+        use_actual_value[i-1] = true;
       }
 
       if (margin) {
@@ -319,8 +318,7 @@ public:
     return this->setup();
   }
 
-  int init_from_boundaries(T const &t = T(),
+  int init_from_boundaries()
-                           size_t const &mult_i = 1)
   {
     if (cvm::debug()) {
       cvm::log("Configuring grid dimensions from colvars boundaries.\n");
@@ -337,7 +335,7 @@ public:
                           lower_boundaries[i].real_value ) / widths[i];
       int nbins_round = (int)(nbins+0.5);
 
-      if (std::fabs(nbins - cvm::real(nbins_round)) > 1.0E-10) {
+      if (cvm::fabs(nbins - cvm::real(nbins_round)) > 1.0E-10) {
         cvm::log("Warning: grid interval("+
                  cvm::to_str(lower_boundaries[i], cvm::cv_width, cvm::cv_prec)+" - "+
                  cvm::to_str(upper_boundaries[i], cvm::cv_width, cvm::cv_prec)+
@@ -392,20 +390,20 @@ public:
   /// \brief Report the bin corresponding to the current value of variable i
   inline int current_bin_scalar(int const i) const
   {
-    return value_to_bin_scalar(actual_value[i] ? cv[i]->actual_value() : cv[i]->value(), i);
+    return value_to_bin_scalar(use_actual_value[i] ? cv[i]->actual_value() : cv[i]->value(), i);
   }
 
   /// \brief Report the bin corresponding to the current value of variable i
   /// and assign first or last bin if out of boundaries
   inline int current_bin_scalar_bound(int const i) const
   {
-    return value_to_bin_scalar_bound(actual_value[i] ? cv[i]->actual_value() : cv[i]->value(), i);
+    return value_to_bin_scalar_bound(use_actual_value[i] ? cv[i]->actual_value() : cv[i]->value(), i);
   }
 
   /// \brief Report the bin corresponding to the current value of item iv in variable i
   inline int current_bin_scalar(int const i, int const iv) const
   {
-    return value_to_bin_scalar(actual_value[i] ?
+    return value_to_bin_scalar(use_actual_value[i] ?
                                cv[i]->actual_value().vector1d_value[iv] :
                                cv[i]->value().vector1d_value[iv], i);
   }
@@ -414,14 +412,14 @@ public:
   /// the provided value is in
   inline int value_to_bin_scalar(colvarvalue const &value, const int i) const
   {
-    return (int) std::floor( (value.real_value - lower_boundaries[i].real_value) / widths[i] );
+    return (int) cvm::floor( (value.real_value - lower_boundaries[i].real_value) / widths[i] );
   }
 
   /// \brief Use the lower boundary and the width to report which bin
   /// the provided value is in and assign first or last bin if out of boundaries
   inline int value_to_bin_scalar_bound(colvarvalue const &value, const int i) const
   {
-    int bin_index = std::floor( (value.real_value - lower_boundaries[i].real_value) / widths[i] );
+    int bin_index = cvm::floor( (value.real_value - lower_boundaries[i].real_value) / widths[i] );
     if (bin_index < 0) bin_index=0;
     if (bin_index >=int(nx[i])) bin_index=int(nx[i])-1;
     return (int) bin_index;
@@ -432,7 +430,7 @@ public:
                                  colvarvalue const &new_offset,
                                  cvm::real   const &new_width) const
   {
-    return (int) std::floor( (value.real_value - new_offset.real_value) / new_width );
+    return (int) cvm::floor( (value.real_value - new_offset.real_value) / new_width );
   }
 
   /// \brief Use the two boundaries and the width to report the
@@ -611,8 +609,8 @@ public:
 
       if (periodic[i]) continue;
 
-      cvm::real dl = std::sqrt(cv[i]->dist2(values[i], lower_boundaries[i])) / widths[i];
+      cvm::real dl = cvm::sqrt(cv[i]->dist2(values[i], lower_boundaries[i])) / widths[i];
-      cvm::real du = std::sqrt(cv[i]->dist2(values[i], upper_boundaries[i])) / widths[i];
+      cvm::real du = cvm::sqrt(cv[i]->dist2(values[i], upper_boundaries[i])) / widths[i];
 
       if (values[i].real_value < lower_boundaries[i])
         dl *= -1.0;
@@ -841,7 +839,7 @@ public:
 
     if (nd < lower_boundaries.size()) nd = lower_boundaries.size();
 
-    if (! actual_value.size()) actual_value.assign(nd, false);
+    if (! use_actual_value.size()) use_actual_value.assign(nd, false);
     if (! periodic.size()) periodic.assign(nd, false);
     if (! widths.size()) widths.assign(nd, 1.0);
 
@@ -849,7 +847,7 @@ public:
     if (old_nx.size()) {
       for (size_t i = 0; i < nd; i++) {
         if ( (old_nx[i] != nx[i]) ||
-             (std::sqrt(cv[i]->dist2(old_lb[i],
+             (cvm::sqrt(cv[i]->dist2(old_lb[i],
                                      lower_boundaries[i])) > 1.0E-10) ) {
           new_params = true;
         }
@@ -874,11 +872,11 @@ public:
   void check_consistency()
   {
     for (size_t i = 0; i < nd; i++) {
-      if ( (std::sqrt(cv[i]->dist2(cv[i]->lower_boundary,
+      if ( (cvm::sqrt(cv[i]->dist2(cv[i]->lower_boundary,
                                    lower_boundaries[i])) > 1.0E-10) ||
-           (std::sqrt(cv[i]->dist2(cv[i]->upper_boundary,
+           (cvm::sqrt(cv[i]->dist2(cv[i]->upper_boundary,
                                    upper_boundaries[i])) > 1.0E-10) ||
-           (std::sqrt(cv[i]->dist2(cv[i]->width,
+           (cvm::sqrt(cv[i]->dist2(cv[i]->width,
                                    widths[i])) > 1.0E-10) ) {
         cvm::error("Error: restart information for a grid is "
                    "inconsistent with that of its colvars.\n");
@@ -896,11 +894,11 @@ public:
       // we skip dist2(), because periodicities and the like should
       // matter: boundaries should be EXACTLY the same (otherwise,
       // map_grid() should be used)
-      if ( (std::fabs(other_grid.lower_boundaries[i] -
+      if ( (cvm::fabs(other_grid.lower_boundaries[i] -
                       lower_boundaries[i]) > 1.0E-10) ||
-           (std::fabs(other_grid.upper_boundaries[i] -
+           (cvm::fabs(other_grid.upper_boundaries[i] -
                       upper_boundaries[i]) > 1.0E-10) ||
-           (std::fabs(other_grid.widths[i] -
+           (cvm::fabs(other_grid.widths[i] -
                       widths[i]) > 1.0E-10) ||
            (data.size() != other_grid.data.size()) ) {
         cvm::error("Error: inconsistency between "
@@ -1036,11 +1034,11 @@ public:
   std::istream & read_multicol(std::istream &is, bool add = false)
   {
     // Data in the header: nColvars, then for each
-    // xiMin, dXi, nPoints, periodic
+    // xiMin, dXi, nPoints, periodic flag
 
     std::string   hash;
     cvm::real     lower, width, x;
-    size_t        n, periodic;
+    size_t        n, periodic_flag;
     bool          remap;
     std::vector<T>        new_value;
     std::vector<int>      nx_read;
@@ -1053,7 +1051,8 @@ public:
 
     if ( !(is >> hash) || (hash != "#") ) {
       cvm::error("Error reading grid at position "+
-                 cvm::to_str(is.tellg())+" in stream(read \"" + hash + "\")\n");
+                 cvm::to_str(static_cast<size_t>(is.tellg()))+
+                 " in stream(read \"" + hash + "\")\n");
       return is;
     }
 
@@ -1075,15 +1074,16 @@ public:
     for (size_t i = 0; i < nd; i++ ) {
       if ( !(is >> hash) || (hash != "#") ) {
         cvm::error("Error reading grid at position "+
-                   cvm::to_str(is.tellg())+" in stream(read \"" + hash + "\")\n");
+                   cvm::to_str(static_cast<size_t>(is.tellg()))+
+                   " in stream(read \"" + hash + "\")\n");
         return is;
       }
 
-      is >> lower >> width >> nx_read[i] >> periodic;
+      is >> lower >> width >> nx_read[i] >> periodic_flag;
 
 
-      if ( (std::fabs(lower - lower_boundaries[i].real_value) > 1.0e-10) ||
+      if ( (cvm::fabs(lower - lower_boundaries[i].real_value) > 1.0e-10) ||
-           (std::fabs(width - widths[i] ) > 1.0e-10) ||
+           (cvm::fabs(width - widths[i] ) > 1.0e-10) ||
            (nx_read[i] != nx[i]) ) {
         cvm::log("Warning: reading from different grid definition (colvar "
                  + cvm::to_str(i+1) + "); remapping data on new grid.\n");
@@ -1246,7 +1246,7 @@ public:
       if (A0 * A1 == 0) {
         return 0.; // can't handle empty bins
       } else {
-        return (std::log((cvm::real)A1) - std::log((cvm::real)A0))
+        return (cvm::logn((cvm::real)A1) - cvm::logn((cvm::real)A0))
           / (widths[n] * 2.);
       }
     } else if (ix[n] > 0 && ix[n] < nx[n]-1) { // not an edge
@@ -1258,7 +1258,7 @@ public:
       if (A0 * A1 == 0) {
         return 0.; // can't handle empty bins
       } else {
-        return (std::log((cvm::real)A1) - std::log((cvm::real)A0))
+        return (cvm::logn((cvm::real)A1) - cvm::logn((cvm::real)A0))
           / (widths[n] * 2.);
       }
     } else {
@@ -1271,8 +1271,8 @@ public:
       if (A0 * A1 * A2 == 0) {
         return 0.; // can't handle empty bins
       } else {
-        return (-1.5 * std::log((cvm::real)A0) + 2. * std::log((cvm::real)A1)
+        return (-1.5 * cvm::logn((cvm::real)A0) + 2. * cvm::logn((cvm::real)A1)
-          - 0.5 * std::log((cvm::real)A2)) * increment / widths[n];
+          - 0.5 * cvm::logn((cvm::real)A2)) * increment / widths[n];
       }
     }
   }

lib/colvars/colvarmodule.cpp

@@ -28,6 +28,8 @@
 colvarmodule::colvarmodule(colvarproxy *proxy_in)
 {
   depth_s = 0;
+  log_level_ = 10;
+
   cv_traj_os = NULL;
 
   if (proxy == NULL) {
@@ -152,12 +154,12 @@ int colvarmodule::read_config_string(std::string const &config_str)
 {
   cvm::log(cvm::line_marker);
   cvm::log("Reading new configuration:\n");
-  std::istringstream config_s(config_str);
+  std::istringstream new_config_s(config_str);
 
   // strip the comments away
   std::string conf = "";
   std::string line;
-  while (parse->read_config_line(config_s, line)) {
+  while (parse->read_config_line(new_config_s, line)) {
     // Delete lines that contain only white space after removing comments
     if (line.find_first_not_of(colvarparse::white_space) != std::string::npos)
       conf.append(line+"\n");
@@ -255,13 +257,16 @@ int colvarmodule::append_new_config(std::string const &new_conf)
 
 int colvarmodule::parse_global_params(std::string const &conf)
 {
-  colvarmodule *cvm = cvm::main();
+  // TODO document and then echo this keyword
+  parse->get_keyval(conf, "logLevel", log_level_, log_level_,
+                    colvarparse::parse_silent);
 
   {
     std::string index_file_name;
     size_t pos = 0;
     while (parse->key_lookup(conf, "indexFile", &index_file_name, &pos)) {
-      cvm->read_index_file(index_file_name.c_str());
+      cvm::log("# indexFile = \""+index_file_name+"\"\n");
+      read_index_file(index_file_name.c_str());
       index_file_name.clear();
     }
   }
@@ -580,6 +585,24 @@ cvm::atom_group *colvarmodule::atom_group_by_name(std::string const &name)
 }
 
 
+void colvarmodule::register_named_atom_group(atom_group *ag) {
+  named_atom_groups.push_back(ag);
+}
+
+
+void colvarmodule::unregister_named_atom_group(cvm::atom_group *ag)
+{
+  for (std::vector<cvm::atom_group *>::iterator agi = named_atom_groups.begin();
+       agi != named_atom_groups.end();
+       agi++) {
+    if (*agi == ag) {
+      named_atom_groups.erase(agi);
+      break;
+    }
+  }
+}
+
+
 int colvarmodule::change_configuration(std::string const &bias_name,
                                        std::string const &conf)
 {
@@ -1034,9 +1057,6 @@ int colvarmodule::analyze()
     cvm::log("colvarmodule::analyze(), step = "+cvm::to_str(it)+".\n");
   }
 
-  if (cvm::step_relative() == 0)
-    cvm::log("Performing analysis.\n");
-
   // perform colvar-specific analysis
   for (std::vector<colvar *>::iterator cvi = variables_active()->begin();
        cvi != variables_active()->end();
@@ -1089,7 +1109,6 @@ int colvarmodule::end_of_step()
 int colvarmodule::setup()
 {
   if (this->size() == 0) return cvm::get_error();
-  // loop over all components of all colvars to reset masses of all groups
   for (std::vector<colvar *>::iterator cvi = variables()->begin();
        cvi != variables()->end();  cvi++) {
     (*cvi)->setup();
@@ -1248,13 +1267,13 @@ std::istream & colvarmodule::read_restart(std::istream &is)
     std::string restart_conf;
     if (is >> colvarparse::read_block("configuration", restart_conf)) {
       parse->get_keyval(restart_conf, "step",
-                        it_restart, (size_t) 0,
+                        it_restart, static_cast<step_number>(0),
-                        colvarparse::parse_silent);
+                        colvarparse::parse_restart);
         it = it_restart;
       std::string restart_version;
       parse->get_keyval(restart_conf, "version",
                         restart_version, std::string(""),
-                        colvarparse::parse_silent);
+                        colvarparse::parse_restart);
       if (restart_version.size() && (restart_version != std::string(COLVARS_VERSION))) {
         cvm::log("This state file was generated with version "+restart_version+"\n");
       }
@@ -1600,14 +1619,16 @@ std::ostream & colvarmodule::write_traj(std::ostream &os)
 }
 
 
-void cvm::log(std::string const &message)
+void cvm::log(std::string const &message, int min_log_level)
 {
+  if (cvm::log_level() < min_log_level) return;
   // allow logging when the module is not fully initialized
   size_t const d = (cvm::main() != NULL) ? depth() : 0;
-  if (d > 0)
+  if (d > 0) {
     proxy->log((std::string(2*d, ' '))+message);
-  else
+  } else {
     proxy->log(message);
+  }
 }
 
 
@@ -1915,14 +1936,190 @@ int cvm::replica_comm_send(char* msg_data, int msg_len, int dest_rep)
 
 
 
+template<typename T> std::string _to_str(T const &x,
+                                         size_t width, size_t prec)
+{
+  std::ostringstream os;
+  if (width) os.width(width);
+  if (prec) {
+    os.setf(std::ios::scientific, std::ios::floatfield);
+    os.precision(prec);
+  }
+  os << x;
+  return os.str();
+}
+
+
+template<typename T> std::string _to_str_vector(std::vector<T> const &x,
+                                                size_t width, size_t prec)
+{
+  if (!x.size()) return std::string("");
+  std::ostringstream os;
+  if (prec) {
+    os.setf(std::ios::scientific, std::ios::floatfield);
+  }
+  os << "{ ";
+  if (width) os.width(width);
+  if (prec) os.precision(prec);
+  os << x[0];
+  for (size_t i = 1; i < x.size(); i++) {
+    os << ", ";
+    if (width) os.width(width);
+    if (prec) os.precision(prec);
+    os << x[i];
+  }
+  os << " }";
+  return os.str();
+}
+
+
+
+std::string colvarmodule::to_str(std::string const &x)
+{
+  return std::string("\"")+x+std::string("\"");
+}
+
+std::string colvarmodule::to_str(char const *x)
+{
+  return std::string("\"")+std::string(x)+std::string("\"");
+}
+
+std::string colvarmodule::to_str(bool x)
+{
+  return (x ? "on" : "off");
+}
+
+std::string colvarmodule::to_str(int const &x,
+                                 size_t width, size_t prec)
+{
+  return _to_str<int>(x, width, prec);
+}
+
+std::string colvarmodule::to_str(size_t const &x,
+                                 size_t width, size_t prec)
+{
+  return _to_str<size_t>(x, width, prec);
+}
+
+std::string colvarmodule::to_str(long int const &x,
+                                 size_t width, size_t prec)
+{
+  return _to_str<long int>(x, width, prec);
+}
+
+std::string colvarmodule::to_str(step_number const &x,
+                                 size_t width, size_t prec)
+{
+  return _to_str<step_number>(x, width, prec);
+}
+
+std::string colvarmodule::to_str(cvm::real const &x,
+                                 size_t width, size_t prec)
+{
+  return _to_str<cvm::real>(x, width, prec);
+}
+
+std::string colvarmodule::to_str(cvm::rvector const &x,
+                                 size_t width, size_t prec)
+{
+  return _to_str<cvm::rvector>(x, width, prec);
+}
+
+std::string colvarmodule::to_str(cvm::quaternion const &x,
+                                 size_t width, size_t prec)
+{
+  return _to_str<cvm::quaternion>(x, width, prec);
+}
+
+std::string colvarmodule::to_str(colvarvalue const &x,
+                                 size_t width, size_t prec)
+{
+  return _to_str<colvarvalue>(x, width, prec);
+}
+
+std::string colvarmodule::to_str(cvm::vector1d<cvm::real> const &x,
+                                 size_t width, size_t prec)
+{
+  return _to_str< cvm::vector1d<cvm::real> >(x, width, prec);
+}
+
+std::string colvarmodule::to_str(cvm::matrix2d<cvm::real> const &x,
+                                 size_t width, size_t prec)
+{
+  return _to_str< cvm::matrix2d<cvm::real> >(x, width, prec);
+}
+
+
+std::string colvarmodule::to_str(std::vector<int> const &x,
+                                 size_t width, size_t prec)
+{
+  return _to_str_vector<int>(x, width, prec);
+}
+
+std::string colvarmodule::to_str(std::vector<size_t> const &x,
+                                 size_t width, size_t prec)
+{
+  return _to_str_vector<size_t>(x, width, prec);
+}
+
+std::string colvarmodule::to_str(std::vector<long int> const &x,
+                                 size_t width, size_t prec)
+{
+  return _to_str_vector<long int>(x, width, prec);
+}
+
+std::string colvarmodule::to_str(std::vector<cvm::real> const &x,
+                                 size_t width, size_t prec)
+{
+  return _to_str_vector<cvm::real>(x, width, prec);
+}
+
+std::string colvarmodule::to_str(std::vector<cvm::rvector> const &x,
+                                 size_t width, size_t prec)
+{
+  return _to_str_vector<cvm::rvector>(x, width, prec);
+}
+
+std::string colvarmodule::to_str(std::vector<cvm::quaternion> const &x,
+                                 size_t width, size_t prec)
+{
+  return _to_str_vector<cvm::quaternion>(x, width, prec);
+}
+
+std::string colvarmodule::to_str(std::vector<colvarvalue> const &x,
+                                 size_t width, size_t prec)
+{
+  return _to_str_vector<colvarvalue>(x, width, prec);
+}
+
+std::string colvarmodule::to_str(std::vector<std::string> const &x,
+                                 size_t width, size_t prec)
+{
+  return _to_str_vector<std::string>(x, width, prec);
+}
+
+
+std::string cvm::wrap_string(std::string const &s, size_t nchars)
+{
+  if (!s.size()) {
+    return std::string(nchars, ' ');
+  } else {
+    return ( (s.size() <= nchars) ?
+             (s+std::string(nchars-s.size(), ' ')) :
+             (std::string(s, 0, nchars)) );
+  }
+}
+
+
 // shared pointer to the proxy object
 colvarproxy *colvarmodule::proxy = NULL;
 
 // static runtime data
 cvm::real colvarmodule::debug_gradients_step_size = 1.0e-07;
 int       colvarmodule::errorCode = 0;
-long      colvarmodule::it = 0;
+int       colvarmodule::log_level_ = 10;
-long      colvarmodule::it_restart = 0;
+cvm::step_number colvarmodule::it = 0;
+cvm::step_number colvarmodule::it_restart = 0;
 size_t    colvarmodule::restart_out_freq = 0;
 size_t    colvarmodule::cv_traj_freq = 0;
 bool      colvarmodule::use_scripted_forces = false;

lib/colvars/colvarmodule.h

@@ -10,6 +10,8 @@
 #ifndef COLVARMODULE_H
 #define COLVARMODULE_H
 
+#include <cmath>
+
 #include "colvars_version.h"
 
 #ifndef COLVARS_DEBUG
@@ -55,6 +57,7 @@ class colvar;
 class colvarbias;
 class colvarproxy;
 class colvarscript;
+class colvarvalue;
 
 
 /// \brief Collective variables module (main class)
@@ -88,10 +91,16 @@ public:
   // TODO colvarscript should be unaware of colvarmodule's internals
   friend class colvarscript;
 
+  /// Use a 64-bit integer to store the step number
+  typedef long long step_number;
+
   /// Defining an abstract real number allows to switch precision
   typedef  double    real;
 
-  /// Override std::pow with a product for n integer
+
+  // Math functions
+
+  /// Override the STL pow() with a product for n integer
   static inline real integer_power(real const &x, int const n)
   {
     // Original code: math_special.h in LAMMPS
@@ -105,8 +114,68 @@ public:
     return (n > 0) ? yy : 1.0/yy;
   }
 
-  /// Residue identifier
+  /// Reimplemented to work around MS compiler issues
-  typedef  int       residue_id;
+  static inline real pow(real const &x, real const &y)
+  {
+    return ::pow(static_cast<double>(x), static_cast<double>(y));
+  }
+
+  /// Reimplemented to work around MS compiler issues
+  static inline real floor(real const &x)
+  {
+    return ::floor(static_cast<double>(x));
+  }
+
+  /// Reimplemented to work around MS compiler issues
+  static inline real fabs(real const &x)
+  {
+    return ::fabs(static_cast<double>(x));
+  }
+
+  /// Reimplemented to work around MS compiler issues
+  static inline real sqrt(real const &x)
+  {
+    return ::sqrt(static_cast<double>(x));
+  }
+
+  /// Reimplemented to work around MS compiler issues
+  static inline real sin(real const &x)
+  {
+    return ::sin(static_cast<double>(x));
+  }
+
+  /// Reimplemented to work around MS compiler issues
+  static inline real cos(real const &x)
+  {
+    return ::cos(static_cast<double>(x));
+  }
+
+  /// Reimplemented to work around MS compiler issues
+  static inline real acos(real const &x)
+  {
+    return ::acos(static_cast<double>(x));
+  }
+
+  /// Reimplemented to work around MS compiler issues
+  static inline real atan2(real const &x, real const &y)
+  {
+    return ::atan2(static_cast<double>(x), static_cast<double>(y));
+  }
+
+  /// Reimplemented to work around MS compiler issues
+  static inline real exp(real const &x)
+  {
+    return ::exp(static_cast<double>(x));
+  }
+
+  /// Reimplemented to work around MS compiler issues.  Note: log() is
+  /// currently defined as the text logging function, but this can be changed
+  /// at a later time
+  static inline real logn(real const &x)
+  {
+    return ::log(static_cast<double>(x));
+  }
+
 
   class rvector;
   template <class T> class vector1d;
@@ -114,6 +183,10 @@ public:
   class quaternion;
   class rotation;
 
+
+  /// Residue identifier
+  typedef int residue_id;
+
   /// \brief Atom position (different type name from rvector, to make
   /// possible future PBC-transparent implementations)
   typedef rvector atom_pos;
@@ -150,19 +223,19 @@ public:
 
 
   /// Current step number
-  static long it;
+  static step_number it;
   /// Starting step number for this run
-  static long it_restart;
+  static step_number it_restart;
 
   /// Return the current step number from the beginning of this run
-  static inline long step_relative()
+  static inline step_number step_relative()
   {
     return it - it_restart;
   }
 
   /// Return the current step number from the beginning of the whole
   /// calculation
-  static inline long step_absolute()
+  static inline step_number step_absolute()
   {
     return it;
   }
@@ -203,9 +276,10 @@ private:
   std::vector<atom_group *> named_atom_groups;
 public:
   /// Register a named atom group into named_atom_groups
-  inline void register_named_atom_group(atom_group * ag) {
+  void register_named_atom_group(atom_group *ag);
-    named_atom_groups.push_back(ag);
+
-  }
+  /// Remove a named atom group from named_atom_groups
+  void unregister_named_atom_group(atom_group *ag);
 
   /// Array of collective variables
   std::vector<colvar *> *variables();
@@ -254,8 +328,7 @@ public:
   /// \brief How many objects are configured yet?
   size_t size() const;
 
-  /// \brief Constructor \param config_name Configuration file name
+  /// \brief Constructor
-  /// \param restart_name (optional) Restart file name
   colvarmodule(colvarproxy *proxy);
 
   /// Destructor
@@ -265,6 +338,7 @@ public:
   int reset();
 
   /// Open a config file, load its contents, and pass it to config_string()
+  /// \param config_file_name Configuration file name
   int read_config_file(char const *config_file_name);
 
   /// \brief Parse a config string assuming it is a complete configuration
@@ -431,26 +505,92 @@ public:
                 long        traj_read_begin,
                 long        traj_read_end);
 
-  /// Quick conversion of an object to a string
+  /// Convert to string for output purposes
-  template<typename T> static std::string to_str(T const &x,
+  static std::string to_str(char const *s);
-                                                  size_t const &width = 0,
+
-                                                  size_t const &prec = 0);
+  /// Convert to string for output purposes
-  /// Quick conversion of a vector of objects to a string
+  static std::string to_str(std::string const &s);
-  template<typename T> static std::string to_str(std::vector<T> const &x,
+
-                                                  size_t const &width = 0,
+  /// Convert to string for output purposes
-                                                  size_t const &prec = 0);
+  static std::string to_str(bool x);
+
+  /// Convert to string for output purposes
+  static std::string to_str(int const &x,
+                            size_t width = 0, size_t prec = 0);
+
+  /// Convert to string for output purposes
+  static std::string to_str(size_t const &x,
+                            size_t width = 0, size_t prec = 0);
+
+  /// Convert to string for output purposes
+  static std::string to_str(long int const &x,
+                            size_t width = 0, size_t prec = 0);
+
+  /// Convert to string for output purposes
+  static std::string to_str(step_number const &x,
+                            size_t width = 0, size_t prec = 0);
+
+  /// Convert to string for output purposes
+  static std::string to_str(real const &x,
+                            size_t width = 0, size_t prec = 0);
+
+  /// Convert to string for output purposes
+  static std::string to_str(rvector const &x,
+                            size_t width = 0, size_t prec = 0);
+
+  /// Convert to string for output purposes
+  static std::string to_str(quaternion const &x,
+                            size_t width = 0, size_t prec = 0);
+
+  /// Convert to string for output purposes
+  static std::string to_str(colvarvalue const &x,
+                            size_t width = 0, size_t prec = 0);
+
+  /// Convert to string for output purposes
+  static std::string to_str(vector1d<real> const &x,
+                            size_t width = 0, size_t prec = 0);
+
+  /// Convert to string for output purposes
+  static std::string to_str(matrix2d<real> const &x,
+                            size_t width = 0, size_t prec = 0);
+
+
+  /// Convert to string for output purposes
+  static std::string to_str(std::vector<int> const &x,
+                            size_t width = 0, size_t prec = 0);
+
+  /// Convert to string for output purposes
+  static std::string to_str(std::vector<size_t> const &x,
+                            size_t width = 0, size_t prec = 0);
+
+  /// Convert to string for output purposes
+  static std::string to_str(std::vector<long int> const &x,
+                            size_t width = 0, size_t prec = 0);
+
+  /// Convert to string for output purposes
+  static std::string to_str(std::vector<real> const &x,
+                            size_t width = 0, size_t prec = 0);
+
+  /// Convert to string for output purposes
+  static std::string to_str(std::vector<rvector> const &x,
+                            size_t width = 0, size_t prec = 0);
+
+  /// Convert to string for output purposes
+  static std::string to_str(std::vector<quaternion> const &x,
+                            size_t width = 0, size_t prec = 0);
+
+  /// Convert to string for output purposes
+  static std::string to_str(std::vector<colvarvalue> const &x,
+                            size_t width = 0, size_t prec = 0);
+
+  /// Convert to string for output purposes
+  static std::string to_str(std::vector<std::string> const &x,
+                            size_t width = 0, size_t prec = 0);
+
 
   /// Reduce the number of characters in a string
-  static inline std::string wrap_string(std::string const &s,
+  static std::string wrap_string(std::string const &s,
-                                         size_t const &nchars)
+                                 size_t nchars);
-  {
-    if (!s.size())
-      return std::string(nchars, ' ');
-    else
-      return ( (s.size() <= size_t(nchars)) ?
-               (s+std::string(nchars-s.size(), ' ')) :
-               (std::string(s, 0, nchars)) );
-  }
 
   /// Number of characters to represent a time step
   static size_t const it_width;
@@ -485,7 +625,9 @@ public:
   static void request_total_force();
 
   /// Print a message to the main log
-  static void log(std::string const &message);
+  /// \param message Message to print
+  /// \param min_log_level Only print if cvm::log_level() >= min_log_level
+  static void log(std::string const &message, int min_log_level = 10);
 
   /// Print a message to the main log and exit with error code
   static int fatal_error(std::string const &message);
@@ -493,6 +635,50 @@ public:
   /// Print a message to the main log and set global error code
   static int error(std::string const &message, int code = COLVARS_ERROR);
 
+private:
+
+  /// Level of logging requested by the user
+  static int log_level_;
+
+public:
+
+  /// Level of logging requested by the user
+  static inline int log_level()
+  {
+    return log_level_;
+  }
+
+  /// Level at which initialization messages are logged
+  static inline int log_init_messages()
+  {
+    return 1;
+  }
+
+  /// Level at which a keyword's user-provided value is logged
+  static inline int log_user_params()
+  {
+    return 2;
+  }
+
+  /// Level at which a keyword's default value is logged
+  static inline int log_default_params()
+  {
+    return 3;
+  }
+
+  /// Level at which output-file operations are logged
+  static inline int log_output_files()
+  {
+    return 4;
+  }
+
+  /// Level at which input-file operations (configuration, state) are logged
+  static inline int log_input_files()
+  {
+    return 5;
+  }
+
+
   // Replica exchange commands.
   static bool replica_enabled();
   static int replica_index();
@@ -630,41 +816,4 @@ std::ostream & operator << (std::ostream &os, cvm::rvector const &v);
 std::istream & operator >> (std::istream &is, cvm::rvector &v);
 
 
-template<typename T> std::string cvm::to_str(T const &x,
-                                             size_t const &width,
-                                             size_t const &prec) {
-  std::ostringstream os;
-  if (width) os.width(width);
-  if (prec) {
-    os.setf(std::ios::scientific, std::ios::floatfield);
-    os.precision(prec);
-  }
-  os << x;
-  return os.str();
-}
-
-
-template<typename T> std::string cvm::to_str(std::vector<T> const &x,
-                                             size_t const &width,
-                                             size_t const &prec) {
-  if (!x.size()) return std::string("");
-  std::ostringstream os;
-  if (prec) {
-    os.setf(std::ios::scientific, std::ios::floatfield);
-  }
-  os << "{ ";
-  if (width) os.width(width);
-  if (prec) os.precision(prec);
-  os << x[0];
-  for (size_t i = 1; i < x.size(); i++) {
-    os << ", ";
-    if (width) os.width(width);
-    if (prec) os.precision(prec);
-    os << x[i];
-  }
-  os << " }";
-  return os.str();
-}
-
-
 #endif

lib/colvars/colvarparse.cpp

@@ -10,6 +10,7 @@
 
 #include <sstream>
 #include <iostream>
+#include <algorithm>
 
 #include "colvarmodule.h"
 #include "colvarvalue.h"
@@ -20,57 +21,197 @@
 char const * const colvarparse::white_space = " \t";
 
 
-// definition of single-value keyword parsers
+namespace {
 
-template<typename TYPE> bool colvarparse::_get_keyval_scalar_(std::string const &conf,
+  // Avoid having to put the bool assignment in the template :-(
+  void set_bool(void *p, bool val)
+  {
+    bool *v = reinterpret_cast<bool *>(p);
+    *v = val;
+  }
+
+}
+
+
+bool colvarparse::get_key_string_value(std::string const &conf,
+                                       char const *key, std::string &data)
+{
+  bool b_found = false, b_found_any = false;
+  size_t save_pos = 0, found_count = 0;
+
+  do {
+    std::string data_this = "";
+    b_found = key_lookup(conf, key, &data_this, &save_pos);
+    if (b_found) {
+      if (!b_found_any)
+        b_found_any = true;
+      found_count++;
+      data = data_this;
+    }
+  } while (b_found);
+
+  if (found_count > 1) {
+    cvm::error("Error: found more than one instance of \""+
+               std::string(key)+"\".\n", INPUT_ERROR);
+  }
+
+  return b_found_any;
+}
+
+
+template<typename TYPE>
+void colvarparse::mark_key_set_user(std::string const &key_str,
+                                    TYPE const &value,
+                                    Parse_Mode const &parse_mode)
+{
+  key_set_modes[to_lower_cppstr(key_str)] = key_set_user;
+  if (parse_mode & parse_echo) {
+    cvm::log("# "+key_str+" = "+cvm::to_str(value)+"\n",
+             cvm::log_user_params());
+  }
+}
+
+
+template<typename TYPE>
+void colvarparse::mark_key_set_default(std::string const &key_str,
+                                       TYPE const &def_value,
+                                       Parse_Mode const &parse_mode)
+{
+  key_set_modes[to_lower_cppstr(key_str)] = key_set_default;
+  if (parse_mode & parse_echo_default) {
+    cvm::log("# "+key_str+" = "+cvm::to_str(def_value)+
+             " [default]\n", cvm::log_default_params());
+  }
+}
+
+
+void colvarparse::error_key_required(std::string const &key_str,
+                                     Parse_Mode const &parse_mode)
+{
+  if (key_already_set(key_str)) {
+    return;
+  }
+  if (parse_mode & parse_restart) {
+    cvm::error("Error: keyword \""+key_str+
+               "\" is missing from the restart.\n", INPUT_ERROR);
+  } else {
+    cvm::error("Error: keyword \""+key_str+
+               "\" is required.\n", INPUT_ERROR);
+  }
+}
+
+
+template<typename TYPE>
+int colvarparse::_get_keyval_scalar_value_(std::string const &key_str,
+                                           std::string const &data,
+                                           TYPE &value,
+                                           TYPE const &def_value)
+{
+  std::istringstream is(data);
+  size_t value_count = 0;
+  TYPE x(def_value);
+
+  while (is >> x) {
+    value = x;
+    value_count++;
+  }
+
+  if (value_count == 0) {
+    return cvm::error("Error: in parsing \""+
+                      key_str+"\".\n", INPUT_ERROR);
+  }
+
+  if (value_count > 1) {
+    return cvm::error("Error: multiple values "
+                      "are not allowed for keyword \""+
+                      key_str+"\".\n", INPUT_ERROR);
+  }
+
+  return COLVARS_OK;
+}
+
+
+template<>
+int colvarparse::_get_keyval_scalar_value_(std::string const &key_str,
+                                           std::string const &data,
+                                           bool &value,
+                                           bool const &def_value)
+{
+  if ( (data == std::string("on")) ||
+       (data == std::string("yes")) ||
+       (data == std::string("true")) ) {
+    set_bool(reinterpret_cast<void *>(&value), true);
+  } else if ( (data == std::string("off")) ||
+              (data == std::string("no")) ||
+              (data == std::string("false")) ) {
+    set_bool(reinterpret_cast<void *>(&value), false);
+  } else {
+    return cvm::error("Error: boolean values only are allowed "
+                      "for \""+key_str+"\".\n", INPUT_ERROR);
+  }
+  return COLVARS_OK;
+}
+
+
+template<typename TYPE>
+int colvarparse::_get_keyval_scalar_novalue_(std::string const &key_str,
+                                             TYPE &value,
+                                             Parse_Mode const &parse_mode)
+{
+  return cvm::error("Error: improper or missing value "
+                    "for \""+key_str+"\".\n", INPUT_ERROR);
+}
+
+template<>
+int colvarparse::_get_keyval_scalar_novalue_(std::string const &key_str,
+                                             bool &value,
+                                             Parse_Mode const &parse_mode)
+{
+  set_bool(reinterpret_cast<void *>(&value), true);
+  mark_key_set_user<bool>(key_str, value, parse_mode);
+  return COLVARS_OK;
+}
+
+
+template<typename TYPE>
+bool colvarparse::_get_keyval_scalar_(std::string const &conf,
                                       char const *key,
                                       TYPE &value,
                                       TYPE const &def_value,
-                                                              Parse_Mode const parse_mode)
+                                      Parse_Mode const &parse_mode)
 {
+  std::string const key_str(key);
+
   std::string data;
-  bool b_found = false, b_found_any = false;
+  bool const b_found_any = get_key_string_value(conf, key, data);
-  size_t save_pos = 0, found_count = 0;
-
-  do {
-    std::string data_this = "";
-    b_found = key_lookup(conf, key, &data_this, &save_pos);
-    if (b_found) {
-      if (!b_found_any)
-        b_found_any = true;
-      found_count++;
-      data = data_this;
-    }
-  } while (b_found);
-
-  if (found_count > 1) {
-    cvm::error("Error: found more than one instance of \""+
-               std::string(key)+"\".\n", INPUT_ERROR);
-  }
 
   if (data.size()) {
-    std::istringstream is(data);
+
-    TYPE x(def_value);
+    _get_keyval_scalar_value_<TYPE>(key_str, data, value, def_value);
-    if (is >> x) {
+
-      value = x;
+    mark_key_set_user<TYPE>(key_str, value, parse_mode);
-    } else {
+
-      cvm::error("Error: in parsing \""+
+  } else { // No string value
-                 std::string(key)+"\".\n", INPUT_ERROR);
-    }
-    if (parse_mode != parse_silent) {
-      cvm::log("# "+std::string(key)+" = "+
-               cvm::to_str(value)+"\n");
-    }
-  } else {
 
     if (b_found_any) {
-      cvm::error("Error: improper or missing value "
+
-                 "for \""+std::string(key)+"\".\n", INPUT_ERROR);
+      _get_keyval_scalar_novalue_<TYPE>(key_str, value, parse_mode);
+
+    } else {
+
+      if (parse_mode & parse_required) {
+        if (cvm::debug()) {
+          cvm::log("get_keyval, parse_required = "+cvm::to_str(parse_mode & parse_required)+
+                   "\n");
         }
+        error_key_required(key_str, parse_mode);
+        return false;
+      }
+
+      if ( (parse_mode & parse_override) || !(key_already_set(key)) ) {
         value = def_value;
-    if (parse_mode != parse_silent) {
+        mark_key_set_default<TYPE>(key_str, value, parse_mode);
-      cvm::log("# "+std::string(key)+" = "+
+      }
-               cvm::to_str(def_value)+" [default]\n");
     }
   }
 
@@ -78,96 +219,17 @@ template<typename TYPE> bool colvarparse::_get_keyval_scalar_(std::string const
 }
 
 
-bool colvarparse::_get_keyval_scalar_string_(std::string const &conf,
+template<typename TYPE>
-                                             char const *key,
+bool colvarparse::_get_keyval_vector_(std::string const &conf,
-                                             std::string &value,
-                                             std::string const &def_value,
-                                             Parse_Mode const parse_mode)
-{
-  std::string data;
-  bool b_found = false, b_found_any = false;
-  size_t save_pos = 0, found_count = 0;
-
-  do {
-    std::string data_this = "";
-    b_found = key_lookup(conf, key, &data_this, &save_pos);
-    if (b_found) {
-      if (!b_found_any)
-        b_found_any = true;
-      found_count++;
-      data = data_this;
-    }
-  } while (b_found);
-
-  if (found_count > 1) {
-    cvm::error("Error: found more than one instance of \""+
-               std::string(key)+"\".\n", INPUT_ERROR);
-  }
-
-  if (data.size()) {
-    std::istringstream is(data);
-    size_t data_count = 0;
-    std::string x(def_value);
-    while (is >> x) {
-      value = x;
-      data_count++;
-    }
-    if (data_count == 0)
-      cvm::error("Error: in parsing \""+
-                 std::string(key)+"\".\n", INPUT_ERROR);
-    if (data_count > 1) {
-      cvm::error("Error: multiple values "
-                 "are not allowed for keyword \""+
-                 std::string(key)+"\".\n", INPUT_ERROR);
-    }
-    if (parse_mode != parse_silent) {
-      cvm::log("# "+std::string(key)+" = \""+
-               cvm::to_str(value)+"\"\n");
-    }
-  } else {
-
-    if (b_found_any) {
-      cvm::error("Error: improper or missing value "
-                 "for \""+std::string(key)+"\".\n", INPUT_ERROR);
-    }
-    value = def_value;
-    if (parse_mode != parse_silent) {
-      cvm::log("# "+std::string(key)+" = \""+
-               cvm::to_str(def_value)+"\" [default]\n");
-    }
-  }
-
-  return b_found_any;
-}
-
-
-// multiple-value keyword parsers
-
-template<typename TYPE> bool colvarparse::_get_keyval_vector_(std::string const &conf,
                                       char const *key,
                                       std::vector<TYPE> &values,
                                       std::vector<TYPE> const &def_values,
-                                                              Parse_Mode const parse_mode)
+                                      Parse_Mode const &parse_mode)
 {
+  std::string const key_str(key);
+
   std::string data;
-  bool b_found = false, b_found_any = false;
+  bool const b_found_any = get_key_string_value(conf, key, data);
-  size_t save_pos = 0, found_count = 0;
-
-  do {
-    std::string data_this = "";
-    b_found = key_lookup(conf, key, &data_this, &save_pos);
-    if (b_found) {
-      if (!b_found_any)
-        b_found_any = true;
-      found_count++;
-      data = data_this;
-    }
-  } while (b_found);
-
-  if (found_count > 1) {
-    cvm::error("Error: found more than one instance of \""+
-               std::string(key)+"\".\n", INPUT_ERROR);
-  }
 
   if (data.size()) {
     std::istringstream is(data);
@@ -175,10 +237,11 @@ template<typename TYPE> bool colvarparse::_get_keyval_vector_(std::string const
     if (values.size() == 0) {
 
       std::vector<TYPE> x;
-      if (def_values.size())
+      if (def_values.size()) {
         x = def_values;
-      else
+      } else {
         x.assign(1, TYPE());
+      }
 
       for (size_t i = 0;
            ( is >> x[ ((i<x.size()) ? i : x.size()-1) ] );
@@ -195,29 +258,39 @@ template<typename TYPE> bool colvarparse::_get_keyval_vector_(std::string const
           values[i] = x;
         } else {
           cvm::error("Error: in parsing \""+
-                     std::string(key)+"\".\n", INPUT_ERROR);
+                     key_str+"\".\n", INPUT_ERROR);
         }
       }
     }
 
-    if (parse_mode != parse_silent) {
+    mark_key_set_user< std::vector<TYPE> >(key_str, values, parse_mode);
-      cvm::log("# "+std::string(key)+" = "+
-               cvm::to_str(values)+"\n");
-    }
 
   } else {
 
     if (b_found_any) {
       cvm::error("Error: improper or missing values for \""+
-                 std::string(key)+"\".\n", INPUT_ERROR);
+                 key_str+"\".\n", INPUT_ERROR);
+    } else {
+
+      if ((values.size() > 0) && (values.size() != def_values.size())) {
+        cvm::error("Error: the number of default values for \""+
+                   key_str+"\" is different from the number of "
+                   "current values.\n", BUG_ERROR);
       }
 
-    for (size_t i = 0; i < values.size(); i++)
+      if (parse_mode & parse_required) {
-      values[i] = def_values[ (i > def_values.size()) ? 0 : i ];
+        error_key_required(key_str, parse_mode);
+        return false;
+      }
+
+      if ( (parse_mode & parse_override) || !(key_already_set(key)) ) {
+        for (size_t i = 0; i < values.size(); i++) {
+          values[i] = def_values[i];
+        }
+        mark_key_set_default< std::vector<TYPE> >(key_str, def_values,
+                                                  parse_mode);
+      }
 
-    if (parse_mode != parse_silent) {
-      cvm::log("# "+std::string(key)+" = "+
-               cvm::to_str(def_values)+" [default]\n");
     }
   }
 
@@ -255,13 +328,22 @@ bool colvarparse::get_keyval(std::string const &conf,
   return _get_keyval_scalar_<long>(conf, key, value, def_value, parse_mode);
 }
 
+bool colvarparse::get_keyval(std::string const &conf,
+                             char const *key,
+                             cvm::step_number &value,
+                             cvm::step_number const &def_value,
+                             Parse_Mode const parse_mode)
+{
+  return _get_keyval_scalar_<cvm::step_number>(conf, key, value, def_value, parse_mode);
+}
+
 bool colvarparse::get_keyval(std::string const &conf,
                              char const *key,
                              std::string &value,
                              std::string const &def_value,
                              Parse_Mode const parse_mode)
 {
-  return _get_keyval_scalar_string_(conf, key, value, def_value, parse_mode);
+  return _get_keyval_scalar_<std::string>(conf, key, value, def_value, parse_mode);
 }
 
 bool colvarparse::get_keyval(std::string const &conf,
@@ -300,66 +382,13 @@ bool colvarparse::get_keyval(std::string const &conf,
   return _get_keyval_scalar_<colvarvalue>(conf, key, value, def_value, parse_mode);
 }
 
-
 bool colvarparse::get_keyval(std::string const &conf,
                              char const *key,
                              bool &value,
                              bool const &def_value,
                              Parse_Mode const parse_mode)
 {
-  std::string data;
+  return _get_keyval_scalar_<bool>(conf, key, value, def_value, parse_mode);
-  bool b_found = false, b_found_any = false;
-  size_t save_pos = 0, found_count = 0;
-
-  do {
-    std::string data_this = "";
-    b_found = key_lookup(conf, key, &data_this, &save_pos);
-    if (b_found) {
-      if (!b_found_any)
-        b_found_any = true;
-      found_count++;
-      data = data_this;
-    }
-  } while (b_found);
-
-  if (found_count > 1) {
-    cvm::error("Error: found more than one instance of \""+
-               std::string(key)+"\".\n", INPUT_ERROR);
-  }
-
-  if (data.size()) {
-    if ( (data == std::string("on")) ||
-         (data == std::string("yes")) ||
-         (data == std::string("true")) ) {
-      value = true;
-    } else if ( (data == std::string("off")) ||
-                (data == std::string("no")) ||
-                (data == std::string("false")) ) {
-      value = false;
-    } else
-      cvm::error("Error: boolean values only are allowed "
-                 "for \""+std::string(key)+"\".\n", INPUT_ERROR);
-    if (parse_mode != parse_silent) {
-      cvm::log("# "+std::string(key)+" = "+
-               (value ? "on" : "off")+"\n");
-    }
-  } else {
-
-    if (b_found_any) {
-      if (parse_mode != parse_silent) {
-        cvm::log("# "+std::string(key)+" = on\n");
-      }
-      value = true;
-    } else {
-      value = def_value;
-      if (parse_mode != parse_silent) {
-        cvm::log("# "+std::string(key)+" = "+
-                 (def_value ? "on" : "off")+" [default]\n");
-      }
-    }
-  }
-
-  return b_found_any;
 }
 
 
@@ -440,15 +469,27 @@ bool colvarparse::get_keyval(std::string const &conf,
 
 void colvarparse::add_keyword(char const *key)
 {
-  for (std::list<std::string>::iterator ki = allowed_keywords.begin();
+  std::string const key_str_lower(to_lower_cppstr(std::string(key)));
-       ki != allowed_keywords.end(); ki++) {
+
-    if (to_lower_cppstr(std::string(key)) == *ki)
+  if (key_set_modes.find(key_str_lower) != key_set_modes.end()) {
     return;
   }
-  // not found in the list
+
-  //   if (cvm::debug())
+  key_set_modes[key_str_lower] = key_not_set;
-  //     cvm::log("Registering a new keyword, \""+std::string (key)+"\".\n");
+
-  allowed_keywords.push_back(to_lower_cppstr(std::string(key)));
+  allowed_keywords.push_back(key_str_lower);
+}
+
+
+bool colvarparse::key_already_set(std::string const &key_str)
+{
+  std::string const key_str_lower(to_lower_cppstr(key_str));
+
+  if (key_set_modes.find(key_str_lower) == key_set_modes.end()) {
+    return false;
+  }
+
+  return (key_set_modes[key_str_lower] > 0);
 }
 
 
@@ -457,6 +498,10 @@ void colvarparse::strip_values(std::string &conf)
   size_t offset = 0;
   data_begin_pos.sort();
   data_end_pos.sort();
+  std::list<size_t>::iterator data_begin_pos_last = std::unique(data_begin_pos.begin(), data_begin_pos.end());
+  data_begin_pos.erase(data_begin_pos_last, data_begin_pos.end());
+  std::list<size_t>::iterator data_end_pos_last = std::unique(data_end_pos.begin(), data_end_pos.end());
+  data_end_pos.erase(data_end_pos_last, data_end_pos.end());
 
   std::list<size_t>::iterator data_begin = data_begin_pos.begin();
   std::list<size_t>::iterator data_end   = data_end_pos.begin();
@@ -464,28 +509,16 @@ void colvarparse::strip_values(std::string &conf)
   for ( ; (data_begin != data_begin_pos.end()) &&
           (data_end   != data_end_pos.end()) ;
         data_begin++, data_end++) {
-
-    //     std::cerr << "data_begin, data_end "
-    //               << *data_begin << ", " << *data_end
-    //               << "\n";
-
     size_t const nchars = *data_end-*data_begin;
-
-    //     std::cerr << "conf[data_begin:data_end] = \""
-    //               << std::string (conf, *data_begin - offset, nchars)
-    //               << "\"\n";
-
     conf.erase(*data_begin - offset, nchars);
     offset += nchars;
-
-    //     std::cerr << ("Stripped config = \"\n"+conf+"\"\n");
-
   }
 }
 
 
 void colvarparse::clear_keyword_registry()
 {
+  key_set_modes.clear();
   allowed_keywords.clear();
   data_begin_pos.clear();
   data_end_pos.clear();

lib/colvars/colvarparse.h

@@ -12,6 +12,7 @@
 
 #include <cstring>
 #include <string>
+#include <map>
 
 #include "colvarmodule.h"
 #include "colvarvalue.h"
@@ -24,35 +25,9 @@
 /// need to parse input inherit from this
 class colvarparse {
 
-protected:
-
-  /// \brief List of legal keywords for this object: this is updated
-  /// by each call to colvarparse::get_keyval() or
-  /// colvarparse::key_lookup()
-  std::list<std::string> allowed_keywords;
-
-  /// \brief List of delimiters for the values of each keyword in the
-  /// configuration string; all keywords will be stripped of their
-  /// values before the keyword check is performed
-  std::list<size_t>      data_begin_pos;
-
-  /// \brief List of delimiters for the values of each keyword in the
-  /// configuration string; all keywords will be stripped of their
-  /// values before the keyword check is performed
-  std::list<size_t>      data_end_pos;
-
-  /// \brief Add a new valid keyword to the list
-  void add_keyword(char const *key);
-
-  /// \brief Remove all the values from the config string
-  void strip_values(std::string &conf);
-
-  /// \brief Configuration string of the object (includes comments)
-  std::string config_string;
-
 public:
 
-
+  /// Default constructor
   inline colvarparse()
   {
     init();
@@ -88,14 +63,23 @@ public:
 
   /// How a keyword is parsed in a string
   enum Parse_Mode {
-    /// \brief(default) Read the first instance of a keyword (if
+    /// Zero for all flags
-    /// any), report its value, and print a warning when there is more
+    parse_null = 0,
-    /// than one
+    /// Print the value of a keyword if it is given
-    parse_normal,
+    parse_echo = (1<<1),
-    /// \brief Like parse_normal, but don't send any message to the log
+    /// Print the default value of a keyword, if it is NOT given
-    /// (useful e.g. in restart files when such messages are very
+    parse_echo_default = (1<<2),
-    /// numerous and redundant)
+    /// Do not print the keyword
-    parse_silent
+    parse_silent = 0,
+    /// Raise error if the keyword is not provided
+    parse_required = (1<<16),
+    /// Successive calls to get_keyval() will override the previous values
+    /// when the keyword is not given any more
+    parse_override = (1<<17),
+    /// 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)
   };
 
   /// \brief Check that all the keywords within "conf" are in the list
@@ -146,6 +130,11 @@ public:
                   long &value,
                   long const &def_value = 0,
                   Parse_Mode const parse_mode = parse_normal);
+  bool get_keyval(std::string const &conf,
+                  char const *key,
+                  cvm::step_number &value,
+                  cvm::step_number const &def_value = 0,
+                  Parse_Mode const parse_mode = parse_normal);
   bool get_keyval(std::string const &conf,
                   char const *key,
                   std::string &value,
@@ -219,23 +208,57 @@ public:
 
 protected:
 
-  // Templates
+  /// Get the string value of a keyword, and save it for later parsing
-  template<typename TYPE> bool _get_keyval_scalar_(std::string const &conf,
+  bool get_key_string_value(std::string const &conf,
+                            char const *key, std::string &data);
+
+  /// Template for single-value keyword parsers
+  template<typename TYPE>
+  bool _get_keyval_scalar_(std::string const &conf,
                            char const *key,
                            TYPE &value,
                            TYPE const &def_value,
-                                                   Parse_Mode const parse_mode);
+                           Parse_Mode const &parse_mode);
-  bool _get_keyval_scalar_string_(std::string const &conf,
-                                  char const *key,
-                                  std::string &value,
-                                  std::string const &def_value,
-                                  Parse_Mode const parse_mode);
 
-  template<typename TYPE> bool _get_keyval_vector_(std::string const &conf,
+  /// Template for multiple-value keyword parsers
+  template<typename TYPE>
+  bool _get_keyval_vector_(std::string const &conf,
                            char const *key,
                            std::vector<TYPE> &values,
                            std::vector<TYPE> const &def_values,
-                                                   Parse_Mode const parse_mode);
+                           Parse_Mode const &parse_mode);
+
+  /// Extract the value of a variable from a string
+  template<typename TYPE>
+  int _get_keyval_scalar_value_(std::string const &key_str,
+                                std::string const &data,
+                                TYPE &value,
+                                TYPE const &def_value);
+
+  /// Handle the case where the user provides a keyword without value
+  template<typename TYPE>
+  int _get_keyval_scalar_novalue_(std::string const &key_str,
+                                  TYPE &value,
+                                  Parse_Mode const &parse_mode);
+
+  /// Record that the keyword has just been user-defined
+  template<typename TYPE>
+  void mark_key_set_user(std::string const &key_str,
+                         TYPE const &value,
+                         Parse_Mode const &parse_mode);
+
+  /// Record that the keyword has just been set to its default value
+  template<typename TYPE>
+  void mark_key_set_default(std::string const &key_str,
+                            TYPE const &def_value,
+                            Parse_Mode const &parse_mode);
+
+  /// Raise error condition due to the keyword being required!
+  void error_key_required(std::string const &key_str,
+                          Parse_Mode const &parse_mode);
+
+  /// True if the keyword has been set already
+  bool key_already_set(std::string const &key_str);
 
 public:
 
@@ -286,7 +309,7 @@ public:
                   size_t *save_pos = NULL);
 
   /// \brief Reads a configuration line, adds it to config_string, and returns
-  /// the stream \param is Input stream \param s String that will hold the
+  /// the stream \param is Input stream \param line String that will hold the
   /// configuration line, with comments stripped
   std::istream & read_config_line(std::istream &is, std::string &line);
 
@@ -299,7 +322,51 @@ public:
   /// from this position
   static int check_braces(std::string const &conf, size_t const start_pos);
 
+protected:
+
+  /// \brief List of legal keywords for this object: this is updated
+  /// by each call to colvarparse::get_keyval() or
+  /// colvarparse::key_lookup()
+  std::list<std::string> allowed_keywords;
+
+  /// How a keyword has been set
+  enum key_set_mode {
+    key_not_set = 0,
+    key_set_user = 1,
+    key_set_default = 2
+  };
+
+  /// Track which keywords have been already set, and how
+  std::map<std::string, key_set_mode> key_set_modes;
+
+  /// \brief List of delimiters for the values of each keyword in the
+  /// configuration string; all keywords will be stripped of their
+  /// values before the keyword check is performed
+  std::list<size_t>      data_begin_pos;
+
+  /// \brief List of delimiters for the values of each keyword in the
+  /// configuration string; all keywords will be stripped of their
+  /// values before the keyword check is performed
+  std::list<size_t>      data_end_pos;
+
+  /// \brief Add a new valid keyword to the list
+  void add_keyword(char const *key);
+
+  /// \brief Remove all the values from the config string
+  void strip_values(std::string &conf);
+
+  /// \brief Configuration string of the object (includes comments)
+  std::string config_string;
+
 };
 
 
+/// Bitwise OR between two Parse_mode flags
+inline colvarparse::Parse_Mode operator | (colvarparse::Parse_Mode const &mode1,
+                                           colvarparse::Parse_Mode const &mode2)
+{
+  return static_cast<colvarparse::Parse_Mode>(static_cast<int>(mode1) |
+                                              static_cast<int>(mode2));
+}
+
 #endif

lib/colvars/colvarproxy.cpp

@@ -60,7 +60,7 @@ bool colvarproxy_system::total_forces_same_step() const
 
 inline int round_to_integer(cvm::real x)
 {
-  return std::floor(x+0.5);
+  return cvm::floor(x+0.5);
 }
 
 
@@ -129,7 +129,10 @@ cvm::rvector colvarproxy_system::position_distance(cvm::atom_pos const &pos1,
 
 
 
-colvarproxy_atoms::colvarproxy_atoms() {}
+colvarproxy_atoms::colvarproxy_atoms()
+{
+  updated_masses_ = updated_charges_ = false;
+}
 
 
 colvarproxy_atoms::~colvarproxy_atoms()
@@ -544,7 +547,7 @@ int colvarproxy_script::run_colvar_gradient_callback(
 
 colvarproxy_tcl::colvarproxy_tcl()
 {
-  _tcl_interp = NULL;
+  tcl_interp_ = NULL;
 }
 
 
@@ -573,7 +576,7 @@ char const *colvarproxy_tcl::tcl_obj_to_str(unsigned char *obj)
 int colvarproxy_tcl::tcl_run_force_callback()
 {
 #if defined(COLVARS_TCL)
-  Tcl_Interp *const tcl_interp = reinterpret_cast<Tcl_Interp *>(_tcl_interp);
+  Tcl_Interp *const tcl_interp = reinterpret_cast<Tcl_Interp *>(tcl_interp_);
   std::string cmd = std::string("calc_colvar_forces ")
     + cvm::to_str(cvm::step_absolute());
   int err = Tcl_Eval(tcl_interp, cmd.c_str());
@@ -596,7 +599,7 @@ int colvarproxy_tcl::tcl_run_colvar_callback(
 {
 #if defined(COLVARS_TCL)
 
-  Tcl_Interp *const tcl_interp = reinterpret_cast<Tcl_Interp *>(_tcl_interp);
+  Tcl_Interp *const tcl_interp = reinterpret_cast<Tcl_Interp *>(tcl_interp_);
   size_t i;
   std::string cmd = std::string("calc_") + name;
   for (i = 0; i < cvc_values.size(); i++) {
@@ -633,7 +636,7 @@ int colvarproxy_tcl::tcl_run_colvar_gradient_callback(
 {
 #if defined(COLVARS_TCL)
 
-  Tcl_Interp *const tcl_interp = reinterpret_cast<Tcl_Interp *>(_tcl_interp);
+  Tcl_Interp *const tcl_interp = reinterpret_cast<Tcl_Interp *>(tcl_interp_);
   size_t i;
   std::string cmd = std::string("calc_") + name + "_gradient";
   for (i = 0; i < cvc_values.size(); i++) {

lib/colvars/colvarproxy.h

@@ -29,7 +29,7 @@
 ///
 /// To interface to a new MD engine, the simplest solution is to derive a new
 /// class from \link colvarproxy \endlink.  Currently implemented are: \link
-/// colvarproxy_lammps, \endlink, \link colvarproxy_namd, \endlink, \link
+/// colvarproxy_lammps \endlink, \link colvarproxy_namd \endlink, \link
 /// colvarproxy_vmd \endlink.
 
 
@@ -227,11 +227,15 @@ public:
 
   inline std::vector<cvm::real> *modify_atom_masses()
   {
+    // assume that we are requesting masses to change them
+    updated_masses_ = true;
     return &atoms_masses;
   }
 
   inline std::vector<cvm::real> *modify_atom_charges()
   {
+    // assume that we are requesting charges to change them
+    updated_charges_ = true;
     return &atoms_charges;
   }
 
@@ -250,6 +254,18 @@ public:
     return &atoms_new_colvar_forces;
   }
 
+  /// Record whether masses have been updated
+  inline bool updated_masses() const
+  {
+    return updated_masses_;
+  }
+
+  /// Record whether masses have been updated
+  inline bool updated_charges() const
+  {
+    return updated_charges_;
+  }
+
 protected:
 
   /// \brief Array of 0-based integers used to uniquely associate atoms
@@ -268,6 +284,9 @@ protected:
   /// \brief Forces applied from colvars, to be communicated to the MD integrator
   std::vector<cvm::rvector> atoms_new_colvar_forces;
 
+  /// Whether the masses and charges have been updated from the host code
+  bool updated_masses_, updated_charges_;
+
   /// Used by all init_atom() functions: create a slot for an atom not
   /// requested yet; returns the index in the arrays
   int add_atom_slot(int atom_id);
@@ -522,7 +541,7 @@ public:
 protected:
 
   /// Pointer to Tcl interpreter object
-  void *_tcl_interp;
+  void *tcl_interp_;
 
   /// Set Tcl pointers
   virtual void init_tcl_pointers();

lib/colvars/colvars_version.h

@@ -1,5 +1,5 @@
 #ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2018-11-16"
+#define COLVARS_VERSION "2019-04-26"
 // 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

lib/colvars/colvarscript.cpp

@@ -137,6 +137,10 @@ int colvarscript::run(int objc, unsigned char *const objv[])
 
   if (cmd == "update") {
     error_code |= proxy->update_input();
+    if (error_code) {
+      result += "Error updating the Colvars module.\n";
+      return error_code;
+    }
     error_code |= colvars->calc();
     error_code |= proxy->update_output();
     if (error_code) {
@@ -273,6 +277,10 @@ int colvarscript::run(int objc, unsigned char *const objv[])
 
 int colvarscript::proc_colvar(colvar *cv, int objc, unsigned char *const objv[]) {
 
+  if (objc < 3) {
+    result = "Missing arguments";
+    return COLVARSCRIPT_ERROR;
+  }
   std::string const subcmd(obj_to_str(objv[2]));
 
   if (subcmd == "value") {
@@ -323,6 +331,47 @@ int colvarscript::proc_colvar(colvar *cv, int objc, unsigned char *const objv[])
     return COLVARS_OK;
   }
 
+  if (subcmd == "getatomgroups") {
+    std::vector<std::vector<int> > lists = cv->get_atom_lists();
+    std::vector<std::vector<int> >::iterator li = lists.begin();
+
+    for ( ; li != lists.end(); ++li) {
+      result += "{";
+      std::vector<int>::iterator lj = (*li).begin();
+      for ( ; lj != (*li).end(); ++lj) {
+        result += cvm::to_str(*lj);
+        result += " ";
+      }
+      result += "} ";
+    }
+    return COLVARS_OK;
+  }
+
+  if (subcmd == "getatomids") {
+    std::vector<int>::iterator li = cv->atom_ids.begin();
+
+    for ( ; li != cv->atom_ids.end(); ++li) {
+      result += cvm::to_str(*li);
+      result += " ";
+    }
+    return COLVARS_OK;
+  }
+
+  if (subcmd == "getgradients") {
+    std::vector<cvm::rvector>::iterator li = cv->atomic_gradients.begin();
+
+    for ( ; li != cv->atomic_gradients.end(); ++li) {
+      result += "{";
+      int j;
+      for (j = 0; j < 3; ++j) {
+        result += cvm::to_str((*li)[j]);
+        result += " ";
+      }
+      result += "} ";
+    }
+    return COLVARS_OK;
+  }
+
   if (subcmd == "getappliedforce") {
     result = (cv->applied_force()).to_simple_string();
     return COLVARS_OK;
@@ -410,6 +459,10 @@ int colvarscript::proc_colvar(colvar *cv, int objc, unsigned char *const objv[])
 
 int colvarscript::proc_bias(colvarbias *b, int objc, unsigned char *const objv[]) {
 
+  if (objc < 3) {
+    result = "Missing arguments";
+    return COLVARSCRIPT_ERROR;
+  }
   std::string const subcmd(obj_to_str(objv[2]));
 
   if (subcmd == "energy") {

lib/colvars/colvarscript.h

@@ -191,8 +191,9 @@ inline static colvarbias *colvarbias_obj(void *pobj)
 
 #ifdef COLVARSCRIPT_CPP
 #define CVSCRIPT_COMM_FN(COMM,N_ARGS_MIN,N_ARGS_MAX,ARGS,FN_BODY)       \
-  int CVSCRIPT_COMM_FNAME(COMM)(void *pobj,                             \
+  extern "C" int CVSCRIPT_COMM_FNAME(COMM)(void *pobj,                  \
-                                int objc, unsigned char *const objv[])  \
+                                           int objc,                    \
+                                           unsigned char *const objv[]) \
   {                                                                     \
     colvarscript *script = colvarscript_obj();                          \
     script->clear_results();                                            \

lib/colvars/colvartypes.cpp

@@ -19,6 +19,8 @@ bool      colvarmodule::rotation::monitor_crossings = false;
 cvm::real colvarmodule::rotation::crossing_threshold = 1.0E-02;
 
 
+namespace  {
+
 /// Numerical recipes diagonalization
 static int jacobi(cvm::real **a, cvm::real *d, cvm::real **v, int *nrot);
 
@@ -28,6 +30,7 @@ static int eigsrt(cvm::real *d, cvm::real **v);
 /// Transpose the matrix
 static int transpose(cvm::real **v);
 
+}
 
 
 std::string cvm::rvector::to_simple_string() const
@@ -245,13 +248,11 @@ cvm::quaternion::position_derivative_inner(cvm::rvector const &pos,
 // Seok C, Dill KA.  Using quaternions to calculate RMSD.  J Comput
 // Chem. 25(15):1849-57 (2004) DOI: 10.1002/jcc.20110 PubMed: 15376254
 
-void colvarmodule::rotation::build_matrix(std::vector<cvm::atom_pos> const &pos1,
+void colvarmodule::rotation::build_correlation_matrix(
-                                          std::vector<cvm::atom_pos> const &pos2,
+                                        std::vector<cvm::atom_pos> const &pos1,
-                                          cvm::matrix2d<cvm::real>         &S)
+                                        std::vector<cvm::atom_pos> const &pos2)
 {
   // build the correlation matrix
-  C.resize(3, 3);
-  C.reset();
   size_t i;
   for (i = 0; i < pos1.size(); i++) {
     C.xx() += pos1[i].x * pos2[i].x;
@@ -264,7 +265,11 @@ void colvarmodule::rotation::build_matrix(std::vector<cvm::atom_pos> const &pos1
     C.zy() += pos1[i].z * pos2[i].y;
     C.zz() += pos1[i].z * pos2[i].z;
   }
+}
 
+
+void colvarmodule::rotation::compute_overlap_matrix()
+{
   // build the "overlap" matrix, whose eigenvectors are stationary
   // points of the RMSD in the space of rotations
   S[0][0] =    C.xx() + C.yy() + C.zz();
@@ -286,37 +291,38 @@ void colvarmodule::rotation::build_matrix(std::vector<cvm::atom_pos> const &pos1
 }
 
 
-void colvarmodule::rotation::diagonalize_matrix(cvm::matrix2d<cvm::real> &S,
+void colvarmodule::rotation::diagonalize_matrix(
-                                                cvm::vector1d<cvm::real> &S_eigval,
+                                            cvm::matrix2d<cvm::real> &m,
-                                                cvm::matrix2d<cvm::real> &S_eigvec)
+                                            cvm::vector1d<cvm::real> &eigval,
+                                            cvm::matrix2d<cvm::real> &eigvec)
 {
-  S_eigval.resize(4);
+  eigval.resize(4);
-  S_eigval.reset();
+  eigval.reset();
-  S_eigvec.resize(4,4);
+  eigvec.resize(4, 4);
-  S_eigvec.reset();
+  eigvec.reset();
 
   // diagonalize
   int jac_nrot = 0;
-  if (jacobi(S.c_array(), S_eigval.c_array(), S_eigvec.c_array(), &jac_nrot) !=
+  if (jacobi(m.c_array(), eigval.c_array(), eigvec.c_array(), &jac_nrot) !=
       COLVARS_OK) {
     cvm::error("Too many iterations in routine jacobi.\n"
                "This is usually the result of an ill-defined set of atoms for "
                "rotational alignment (RMSD, rotateReference, etc).\n");
   }
-  eigsrt(S_eigval.c_array(), S_eigvec.c_array());
+  eigsrt(eigval.c_array(), eigvec.c_array());
   // jacobi saves eigenvectors by columns
-  transpose(S_eigvec.c_array());
+  transpose(eigvec.c_array());
 
   // normalize eigenvectors
   for (size_t ie = 0; ie < 4; ie++) {
     cvm::real norm2 = 0.0;
     size_t i;
     for (i = 0; i < 4; i++) {
-      norm2 += S_eigvec[ie][i] * S_eigvec[ie][i];
+      norm2 += eigvec[ie][i] * eigvec[ie][i];
     }
-    cvm::real const norm = std::sqrt(norm2);
+    cvm::real const norm = cvm::sqrt(norm2);
     for (i = 0; i < 4; i++) {
-      S_eigvec[ie][i] /= norm;
+      eigvec[ie][i] /= norm;
     }
   }
 }
@@ -324,23 +330,26 @@ void colvarmodule::rotation::diagonalize_matrix(cvm::matrix2d<cvm::real> &S,
 
 // Calculate the rotation, plus its derivatives
 
-void colvarmodule::rotation::calc_optimal_rotation(std::vector<cvm::atom_pos> const &pos1,
+void colvarmodule::rotation::calc_optimal_rotation(
+                                        std::vector<cvm::atom_pos> const &pos1,
                                         std::vector<cvm::atom_pos> const &pos2)
 {
-  S.resize(4,4);
+  C.resize(3, 3);
+  C.reset();
+  build_correlation_matrix(pos1, pos2);
+
+  S.resize(4, 4);
   S.reset();
+  compute_overlap_matrix();
 
-  build_matrix(pos1, pos2, S);
+  S_backup.resize(4, 4);
-
-  S_backup.resize(4,4);
   S_backup = S;
 
   if (b_debug_gradients) {
-    cvm::log("S     = "+cvm::to_str(cvm::to_str(S_backup), cvm::cv_width, cvm::cv_prec)+"\n");
+    cvm::log("S     = "+cvm::to_str(S_backup, cvm::cv_width, cvm::cv_prec)+"\n");
   }
 
   diagonalize_matrix(S, S_eigval, S_eigvec);
-
   // eigenvalues and eigenvectors
   cvm::real const L0 = S_eigval[0];
   cvm::real const L1 = S_eigval[1];
@@ -524,7 +533,7 @@ void colvarmodule::rotation::calc_optimal_rotation(std::vector<cvm::atom_pos> co
                                   dq0_2[3][comp] * colvarmodule::debug_gradients_step_size);
 
         cvm::log(  "|(l_0+dl_0) - l_0^new|/l_0 = "+
-                   cvm::to_str(std::fabs(L0+DL0 - L0_new)/L0, cvm::cv_width, cvm::cv_prec)+
+                   cvm::to_str(cvm::fabs(L0+DL0 - L0_new)/L0, cvm::cv_width, cvm::cv_prec)+
                    ", |(q_0+dq_0) - q_0^new| = "+
                    cvm::to_str((Q0+DQ0 - Q0_new).norm(), cvm::cv_width, cvm::cv_prec)+
                    "\n");
@@ -544,6 +553,9 @@ void colvarmodule::rotation::calc_optimal_rotation(std::vector<cvm::atom_pos> co
 
 #define n 4
 
+
+namespace {
+
 int jacobi(cvm::real **a, cvm::real *d, cvm::real **v, int *nrot)
 {
   int j,iq,ip,i;
@@ -567,7 +579,7 @@ int jacobi(cvm::real **a, cvm::real *d, cvm::real **v, int *nrot)
     sm=0.0;
     for (ip=0;ip<n-1;ip++) {
       for (iq=ip+1;iq<n;iq++)
-        sm += std::fabs(a[ip][iq]);
+        sm += cvm::fabs(a[ip][iq]);
     }
     if (sm == 0.0) {
       return COLVARS_OK;
@@ -578,20 +590,20 @@ int jacobi(cvm::real **a, cvm::real *d, cvm::real **v, int *nrot)
       tresh=0.0;
     for (ip=0;ip<n-1;ip++) {
       for (iq=ip+1;iq<n;iq++) {
-        g=100.0*std::fabs(a[ip][iq]);
+        g=100.0*cvm::fabs(a[ip][iq]);
-        if (i > 4 && (cvm::real)(std::fabs(d[ip])+g) == (cvm::real)std::fabs(d[ip])
+        if (i > 4 && (cvm::real)(cvm::fabs(d[ip])+g) == (cvm::real)cvm::fabs(d[ip])
-            && (cvm::real)(std::fabs(d[iq])+g) == (cvm::real)std::fabs(d[iq]))
+            && (cvm::real)(cvm::fabs(d[iq])+g) == (cvm::real)cvm::fabs(d[iq]))
           a[ip][iq]=0.0;
-        else if (std::fabs(a[ip][iq]) > tresh) {
+        else if (cvm::fabs(a[ip][iq]) > tresh) {
           h=d[iq]-d[ip];
-          if ((cvm::real)(std::fabs(h)+g) == (cvm::real)std::fabs(h))
+          if ((cvm::real)(cvm::fabs(h)+g) == (cvm::real)cvm::fabs(h))
             t=(a[ip][iq])/h;
           else {
             theta=0.5*h/(a[ip][iq]);
-            t=1.0/(std::fabs(theta)+std::sqrt(1.0+theta*theta));
+            t=1.0/(cvm::fabs(theta)+cvm::sqrt(1.0+theta*theta));
             if (theta < 0.0) t = -t;
           }
-          c=1.0/std::sqrt(1+t*t);
+          c=1.0/cvm::sqrt(1+t*t);
           s=t*c;
           tau=s/(1.0+c);
           h=t*a[ip][iq];
@@ -663,5 +675,7 @@ int transpose(cvm::real **v)
   return COLVARS_OK;
 }
 
+}
+
 #undef n
 #undef ROTATE

lib/colvars/colvartypes.h

@@ -10,7 +10,6 @@
 #ifndef COLVARTYPES_H
 #define COLVARTYPES_H
 
-#include <cmath>
 #include <vector>
 
 #include "colvarmodule.h"
@@ -220,7 +219,7 @@ public:
 
   inline cvm::real norm() const
   {
-    return std::sqrt(this->norm2());
+    return cvm::sqrt(this->norm2());
   }
 
   inline cvm::real sum() const
@@ -801,7 +800,7 @@ public:
 
   inline cvm::real norm() const
   {
-    return std::sqrt(this->norm2());
+    return cvm::sqrt(this->norm2());
   }
 
   inline cvm::rvector unit() const
@@ -1008,17 +1007,17 @@ public:
                                     cvm::real theta_in,
                                     cvm::real psi_in)
   {
-    q0 = ( (std::cos(phi_in/2.0)) * (std::cos(theta_in/2.0)) * (std::cos(psi_in/2.0)) +
+    q0 = ( (cvm::cos(phi_in/2.0)) * (cvm::cos(theta_in/2.0)) * (cvm::cos(psi_in/2.0)) +
-           (std::sin(phi_in/2.0)) * (std::sin(theta_in/2.0)) * (std::sin(psi_in/2.0)) );
+           (cvm::sin(phi_in/2.0)) * (cvm::sin(theta_in/2.0)) * (cvm::sin(psi_in/2.0)) );
 
-    q1 = ( (std::sin(phi_in/2.0)) * (std::cos(theta_in/2.0)) * (std::cos(psi_in/2.0)) -
+    q1 = ( (cvm::sin(phi_in/2.0)) * (cvm::cos(theta_in/2.0)) * (cvm::cos(psi_in/2.0)) -
-           (std::cos(phi_in/2.0)) * (std::sin(theta_in/2.0)) * (std::sin(psi_in/2.0)) );
+           (cvm::cos(phi_in/2.0)) * (cvm::sin(theta_in/2.0)) * (cvm::sin(psi_in/2.0)) );
 
-    q2 = ( (std::cos(phi_in/2.0)) * (std::sin(theta_in/2.0)) * (std::cos(psi_in/2.0)) +
+    q2 = ( (cvm::cos(phi_in/2.0)) * (cvm::sin(theta_in/2.0)) * (cvm::cos(psi_in/2.0)) +
-           (std::sin(phi_in/2.0)) * (std::cos(theta_in/2.0)) * (std::sin(psi_in/2.0)) );
+           (cvm::sin(phi_in/2.0)) * (cvm::cos(theta_in/2.0)) * (cvm::sin(psi_in/2.0)) );
 
-    q3 = ( (std::cos(phi_in/2.0)) * (std::cos(theta_in/2.0)) * (std::sin(psi_in/2.0)) -
+    q3 = ( (cvm::cos(phi_in/2.0)) * (cvm::cos(theta_in/2.0)) * (cvm::sin(psi_in/2.0)) -
-           (std::sin(phi_in/2.0)) * (std::sin(theta_in/2.0)) * (std::cos(psi_in/2.0)) );
+           (cvm::sin(phi_in/2.0)) * (cvm::sin(theta_in/2.0)) * (cvm::cos(psi_in/2.0)) );
   }
 
   /// \brief Default constructor
@@ -1115,7 +1114,7 @@ public:
   /// Norm of the quaternion
   inline cvm::real norm() const
   {
-    return std::sqrt(this->norm2());
+    return cvm::sqrt(this->norm2());
   }
 
   /// Return the conjugate quaternion
@@ -1177,7 +1176,7 @@ public:
   }
 
   /// \brief Provides the quaternion product.  \b NOTE: for the inner
-  /// product use: \code h.inner (q); \endcode
+  /// product use: `h.inner (q);`
   friend inline cvm::quaternion operator * (cvm::quaternion const &h,
                                             cvm::quaternion const &q)
   {
@@ -1263,7 +1262,7 @@ public:
     cvm::real const cos_omega = this->q0*Q2.q0 + this->q1*Q2.q1 +
       this->q2*Q2.q2 + this->q3*Q2.q3;
 
-    cvm::real const omega = std::acos( (cos_omega > 1.0) ? 1.0 :
+    cvm::real const omega = cvm::acos( (cos_omega > 1.0) ? 1.0 :
                                        ( (cos_omega < -1.0) ? -1.0 : cos_omega) );
 
     // get the minimum distance: x and -x are the same quaternion
@@ -1278,11 +1277,11 @@ public:
   inline cvm::quaternion dist2_grad(cvm::quaternion const &Q2) const
   {
     cvm::real const cos_omega = this->q0*Q2.q0 + this->q1*Q2.q1 + this->q2*Q2.q2 + this->q3*Q2.q3;
-    cvm::real const omega = std::acos( (cos_omega > 1.0) ? 1.0 :
+    cvm::real const omega = cvm::acos( (cos_omega > 1.0) ? 1.0 :
                                        ( (cos_omega < -1.0) ? -1.0 : cos_omega) );
-    cvm::real const sin_omega = std::sin(omega);
+    cvm::real const sin_omega = cvm::sin(omega);
 
-    if (std::fabs(sin_omega) < 1.0E-14) {
+    if (cvm::fabs(sin_omega) < 1.0E-14) {
       // return a null 4d vector
       return cvm::quaternion(0.0, 0.0, 0.0, 0.0);
     }
@@ -1338,14 +1337,16 @@ public:
   /// \brief Perform gradient tests
   bool b_debug_gradients;
 
-  /// \brief Positions to superimpose: the rotation should brings pos1
+  /// Correlation matrix C (3, 3)
-  /// into pos2
-  std::vector<cvm::atom_pos> pos1, pos2;
-
   cvm::rmatrix C;
 
+  /// Overlap matrix S (4, 4)
   cvm::matrix2d<cvm::real> S;
+
+  /// Eigenvalues of S
   cvm::vector1d<cvm::real> S_eigval;
+
+  /// Eigenvectors of S
   cvm::matrix2d<cvm::real> S_eigvec;
 
   /// Used for debugging gradients
@@ -1404,8 +1405,8 @@ public:
     : b_debug_gradients(false)
   {
     cvm::rvector const axis_n = axis.unit();
-    cvm::real const sina = std::sin(angle/2.0);
+    cvm::real const sina = cvm::sin(angle/2.0);
-    q = cvm::quaternion(std::cos(angle/2.0),
+    q = cvm::quaternion(cvm::cos(angle/2.0),
                         sina * axis_n.x, sina * axis_n.y, sina * axis_n.z);
   }
 
@@ -1437,7 +1438,7 @@ public:
   inline cvm::real spin_angle(cvm::rvector const &axis) const
   {
     cvm::rvector const q_vec = q.get_vector();
-    cvm::real alpha = (180.0/PI) * 2.0 * std::atan2(axis * q_vec, q.q0);
+    cvm::real alpha = (180.0/PI) * 2.0 * cvm::atan2(axis * q_vec, q.q0);
     while (alpha >  180.0) alpha -= 360;
     while (alpha < -180.0) alpha += 360;
     return alpha;
@@ -1473,9 +1474,9 @@ public:
   {
     cvm::rvector const q_vec = q.get_vector();
     cvm::real const alpha =
-      (180.0/PI) * 2.0 * std::atan2(axis * q_vec, q.q0);
+      (180.0/PI) * 2.0 * cvm::atan2(axis * q_vec, q.q0);
 
-    cvm::real const cos_spin_2 = std::cos(alpha * (PI/180.0) * 0.5);
+    cvm::real const cos_spin_2 = cvm::cos(alpha * (PI/180.0) * 0.5);
     cvm::real const cos_theta_2 = ( (cos_spin_2 != 0.0) ?
                                     (q.q0 / cos_spin_2) :
                                     (0.0) );
@@ -1489,7 +1490,7 @@ public:
     cvm::rvector const q_vec = q.get_vector();
     cvm::real const iprod = axis * q_vec;
 
-    cvm::real const cos_spin_2 = std::cos(std::atan2(iprod, q.q0));
+    cvm::real const cos_spin_2 = cvm::cos(cvm::atan2(iprod, q.q0));
 
     if (q.q0 != 0.0)  {
 
@@ -1529,15 +1530,17 @@ protected:
   /// eigenvalue crossing)
   cvm::quaternion q_old;
 
-  /// Build the overlap matrix S (used by calc_optimal_rotation())
+  /// Build the correlation matrix C (used by calc_optimal_rotation())
-  void build_matrix(std::vector<cvm::atom_pos> const &pos1,
+  void build_correlation_matrix(std::vector<cvm::atom_pos> const &pos1,
-                    std::vector<cvm::atom_pos> const &pos2,
+                                std::vector<cvm::atom_pos> const &pos2);
-                    cvm::matrix2d<cvm::real>         &S);
 
-  /// Diagonalize the overlap matrix S (used by calc_optimal_rotation())
+  /// Compute the overlap matrix S (used by calc_optimal_rotation())
-  void diagonalize_matrix(cvm::matrix2d<cvm::real> &S,
+  void compute_overlap_matrix();
-                          cvm::vector1d<cvm::real> &S_eigval,
+
-                          cvm::matrix2d<cvm::real> &S_eigvec);
+  /// Diagonalize a given matrix m (used by calc_optimal_rotation())
+  static void diagonalize_matrix(cvm::matrix2d<cvm::real> &m,
+                                 cvm::vector1d<cvm::real> &eigval,
+                                 cvm::matrix2d<cvm::real> &eigvec);
 };
 
 

lib/colvars/colvarvalue.cpp

@@ -144,10 +144,10 @@ void colvarvalue::apply_constraints()
   case colvarvalue::type_quaternionderiv:
     break;
   case colvarvalue::type_unit3vector:
-    rvector_value /= std::sqrt(rvector_value.norm2());
+    rvector_value /= cvm::sqrt(rvector_value.norm2());
     break;
   case colvarvalue::type_quaternion:
-    quaternion_value /= std::sqrt(quaternion_value.norm2());
+    quaternion_value /= cvm::sqrt(quaternion_value.norm2());
     break;
   case colvarvalue::type_vector:
     if (elem_types.size() > 0) {
@@ -579,7 +579,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 = std::sqrt(1.0 - cos_t*cos_t);
+      cvm::real const sin_t = cvm::sqrt(1.0 - cos_t*cos_t);
       return colvarvalue( 2.0 * sin_t *
                           cvm::rvector((-1.0) * sin_t * v2.x +
                                        cos_t/sin_t * (v1.x - cos_t*v2.x),
@@ -630,7 +630,7 @@ colvarvalue const colvarvalue::interpolate(colvarvalue const &x1,
     break;
   case colvarvalue::type_unit3vector:
   case colvarvalue::type_quaternion:
-    if (interp.norm()/std::sqrt(d2) < 1.0e-6) {
+    if (interp.norm()/cvm::sqrt(d2) < 1.0e-6) {
       cvm::error("Error: interpolation between "+cvm::to_str(x1)+" and "+
                  cvm::to_str(x2)+" with lambda = "+cvm::to_str(lambda)+
                  " is undefined: result = "+cvm::to_str(interp)+"\n",

lib/colvars/colvarvalue.h

@@ -17,22 +17,22 @@
 /// \brief Value of a collective variable: this is a metatype which
 /// can be set at runtime.  By default it is set to be a scalar
 /// number, and can be treated as such in all operations (this is
-/// done by most \link cvc \endlink implementations).
+/// done by most \link colvar::cvc \endlink implementations).
 ///
 /// \link colvarvalue \endlink allows \link colvar \endlink to be
 /// treat different data types.  By default, a \link colvarvalue
 /// \endlink variable is a scalar number.  To use it as
 /// another type, declare and initialize it as
-/// \code colvarvalue x(colvarvalue::type_xxx)\endcode, use \link x.type
+/// `colvarvalue x(colvarvalue::type_xxx)`, use `x.type (colvarvalue::type_xxx)`
-/// (colvarvalue::type_xxx) \endlink at a later stage, or if unset,
+///  at a later stage, or if unset,
-///  assign the type with \code x = y; \endcode, provided y is correctly set.
+///  assign the type with `x = y;`, provided y is correctly set.
 ///
 /// All operators (either unary or binary) on a \link
 /// colvarvalue \endlink object performs one or more checks on the
 /// \link Type \endlink, except when reading from a stream, when there is no way to
-/// detect the \link Type \endlink.  To use  \code is >> x; \endcode x \b MUST
+/// detect the \link Type \endlink.  To use  `is >> x;` x \b MUST
 /// already have a type correcly set up for properly parsing the
-/// stream. No problem of course with the output streams: \code os << x; \endcode
+/// stream. No problem of course with the output streams: `os << x;`
 ///
 /// \em Note \em on \em performance: to avoid type checks in a long array of \link
 /// colvarvalue \endlink objects, use one of the existing "_opt" functions or implement a new one
@@ -159,7 +159,7 @@ public:
   /// \brief If the variable has constraints (e.g. unitvector or
   /// quaternion), transform it to satisfy them; this function needs
   /// to be called only when the \link colvarvalue \endlink
-  /// is calculated outside of \link cvc \endlink objects
+  /// is calculated outside of \link colvar::cvc \endlink objects
   void apply_constraints();
 
   /// Get the current type
@@ -184,7 +184,7 @@ public:
   /// Norm of this colvarvalue
   inline cvm::real norm() const
   {
-    return std::sqrt(this->norm2());
+    return cvm::sqrt(this->norm2());
   }
 
   /// Sum of the components of this colvarvalue (if more than one dimension)
@@ -728,7 +728,7 @@ inline cvm::real colvarvalue::dist2(colvarvalue const &x2) const
   case colvarvalue::type_unit3vector:
   case colvarvalue::type_unit3vectorderiv:
     // angle between (*this) and x2 is the distance
-    return std::acos(this->rvector_value * x2.rvector_value) * std::acos(this->rvector_value * x2.rvector_value);
+    return cvm::acos(this->rvector_value * x2.rvector_value) * cvm::acos(this->rvector_value * x2.rvector_value);
   case colvarvalue::type_quaternion:
   case colvarvalue::type_quaternionderiv:
     // angle between (*this) and x2 is the distance, the quaternion

src/USER-COLVARS/colvarproxy_lammps.cpp

@@ -153,18 +153,29 @@ void colvarproxy_lammps::init(const char *conf_file)
   if (_lmp->update->ntimestep != 0) {
     cvm::log("Setting initial step number from LAMMPS: "+
              cvm::to_str(_lmp->update->ntimestep)+"\n");
-    colvars->it = colvars->it_restart = _lmp->update->ntimestep;
+    colvars->it = colvars->it_restart =
+      static_cast<cvm::step_number>(_lmp->update->ntimestep);
   }
 
   if (cvm::debug()) {
-    log("atoms_ids = "+cvm::to_str(atoms_ids)+"\n");
+    cvm::log("atoms_ids = "+cvm::to_str(atoms_ids)+"\n");
-    log("atoms_ncopies = "+cvm::to_str(atoms_ncopies)+"\n");
+    cvm::log("atoms_ncopies = "+cvm::to_str(atoms_ncopies)+"\n");
-    log("atoms_positions = "+cvm::to_str(atoms_positions)+"\n");
+    cvm::log("atoms_positions = "+cvm::to_str(atoms_positions)+"\n");
-    log(cvm::line_marker);
+    cvm::log(cvm::line_marker);
-    log("Info: done initializing the colvars proxy object.\n");
+    cvm::log("Info: done initializing the colvars proxy object.\n");
   }
 }
 
+void colvarproxy_lammps::add_config_file(const char *conf_file)
+{
+  colvars->read_config_file(conf_file);
+}
+
+void colvarproxy_lammps::add_config_string(const std::string &conf)
+{
+  colvars->read_config_string(conf);
+}
+
 colvarproxy_lammps::~colvarproxy_lammps()
 {
   delete _random;
@@ -185,7 +196,7 @@ int colvarproxy_lammps::setup()
 double colvarproxy_lammps::compute()
 {
   if (cvm::debug()) {
-    log(std::string(cvm::line_marker)+
+    cvm::log(std::string(cvm::line_marker)+
         "colvarproxy_lammps step no. "+
         cvm::to_str(_lmp->update->ntimestep)+" [first - last = "+
         cvm::to_str(_lmp->update->beginstep)+" - "+
@@ -238,20 +249,20 @@ double colvarproxy_lammps::compute()
   bias_energy = 0.0;
 
   if (cvm::debug()) {
-    log("atoms_ids = "+cvm::to_str(atoms_ids)+"\n");
+    cvm::log("atoms_ids = "+cvm::to_str(atoms_ids)+"\n");
-    log("atoms_ncopies = "+cvm::to_str(atoms_ncopies)+"\n");
+    cvm::log("atoms_ncopies = "+cvm::to_str(atoms_ncopies)+"\n");
-    log("atoms_positions = "+cvm::to_str(atoms_positions)+"\n");
+    cvm::log("atoms_positions = "+cvm::to_str(atoms_positions)+"\n");
-    log("atoms_new_colvar_forces = "+cvm::to_str(atoms_new_colvar_forces)+"\n");
+    cvm::log("atoms_new_colvar_forces = "+cvm::to_str(atoms_new_colvar_forces)+"\n");
   }
 
   // call the collective variable module
   colvars->calc();
 
   if (cvm::debug()) {
-    log("atoms_ids = "+cvm::to_str(atoms_ids)+"\n");
+    cvm::log("atoms_ids = "+cvm::to_str(atoms_ids)+"\n");
-    log("atoms_ncopies = "+cvm::to_str(atoms_ncopies)+"\n");
+    cvm::log("atoms_ncopies = "+cvm::to_str(atoms_ncopies)+"\n");
-    log("atoms_positions = "+cvm::to_str(atoms_positions)+"\n");
+    cvm::log("atoms_positions = "+cvm::to_str(atoms_positions)+"\n");
-    log("atoms_new_colvar_forces = "+cvm::to_str(atoms_new_colvar_forces)+"\n");
+    cvm::log("atoms_new_colvar_forces = "+cvm::to_str(atoms_new_colvar_forces)+"\n");
   }
 
   return bias_energy;

src/USER-COLVARS/colvarproxy_lammps.h

@@ -103,6 +103,11 @@ class colvarproxy_lammps : public colvarproxy {
   // Write files expected from Colvars (called by post_run())
   void write_output_files();
 
+  // read additional config from file
+  void add_config_file(char const *config_filename);
+
+  // read additional config from string
+  void add_config_string(const std::string &config);
 
   // implementation of pure methods from base class
  public:

src/USER-COLVARS/colvarproxy_lammps_version.h

@@ -1,5 +1,5 @@
 #ifndef COLVARPROXY_VERSION
-#define COLVARPROXY_VERSION "2018-08-29"
+#define COLVARPROXY_VERSION "2019-04-09"
 // 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

src/USER-COLVARS/fix_colvars.cpp

@@ -480,6 +480,31 @@ void FixColvars::one_time_init()
 
 /* ---------------------------------------------------------------------- */
 
+int FixColvars::modify_param(int narg, char **arg)
+{
+  if (strcmp(arg[0],"configfile") == 0) {
+    if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
+    if (me == 0) {
+      if (! proxy)
+        error->one(FLERR,"Cannot use fix_modify before initialization");
+      proxy->add_config_file(arg[1]);
+    }
+    return 2;
+  } else if (strcmp(arg[0],"config") == 0) {
+    if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
+    if (me == 0) {
+      if (! proxy)
+        error->one(FLERR,"Cannot use fix_modify before initialization");
+      std::string conf(arg[1]);
+      proxy->add_config_string(conf);
+    }
+    return 2;
+  }
+  return 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
 void FixColvars::setup(int vflag)
 {
   const tagint * const tag  = atom->tag;

src/USER-COLVARS/fix_colvars.h

@@ -50,6 +50,7 @@ class FixColvars : public Fix {
   virtual int setmask();
   virtual void init();
   virtual void setup(int);
+  virtual int modify_param(int, char **);
   virtual void min_setup(int vflag) {setup(vflag);};
   virtual void min_post_force(int);
   virtual void post_force(int);