Merge pull request #4556 from Colvars/colvars-update

Update Colvars library to version 2025-04-18

cmake/Modules/Packages/COLVARS.cmake

@@ -26,6 +26,11 @@ if(BUILD_OMP)
   target_link_libraries(colvars PRIVATE OpenMP::OpenMP_CXX)
 endif()
 
+if(BUILD_MPI)
+  target_compile_definitions(colvars PUBLIC -DCOLVARS_MPI)
+  target_link_libraries(colvars PUBLIC MPI::MPI_CXX)
+endif()
+
 if(COLVARS_DEBUG)
   # Need to export the define publicly to be valid in interface code
   target_compile_definitions(colvars PUBLIC -DCOLVARS_DEBUG)

lib/colvars/Makefile.common

@@ -32,6 +32,7 @@ COLVARS_SRCS = \
         colvarbias_histogram_reweight_amd.cpp \
         colvarbias_meta.cpp \
         colvarbias_restraint.cpp \
+        colvarbias_opes.cpp \
         colvarcomp_alchlambda.cpp \
         colvarcomp_angles.cpp \
         colvarcomp_apath.cpp \
@@ -40,6 +41,7 @@ COLVARS_SRCS = \
         colvarcomp_distances.cpp \
         colvarcomp_gpath.cpp \
         colvarcomp_neuralnetwork.cpp \
+        colvarcomp_torchann.cpp \
         colvarcomp_combination.cpp \
         colvarcomp_protein.cpp \
         colvarcomp_rotations.cpp \

lib/colvars/Makefile.deps

@@ -52,6 +52,12 @@ $(COLVARS_OBJ_DIR)colvarbias_restraint.o: colvarbias_restraint.cpp \
  colvarproxy_tcl.h colvarproxy_volmaps.h colvarvalue.h \
  colvarbias_restraint.h colvarbias.h colvar.h colvarparse.h \
  colvarparams.h colvardeps.h
+$(COLVARS_OBJ_DIR)colvarbias_opes.o: colvarbias_opes.cpp \
+ colvarmodule.h colvars_version.h colvarproxy.h colvartypes.h \
+ ../../src/math_eigen_impl.h colvarproxy_io.h colvarproxy_system.h \
+ colvarproxy_tcl.h colvarproxy_volmaps.h colvarvalue.h \
+ colvarbias_opes.h colvarbias.h colvar.h colvarparse.h \
+ colvarparams.h colvardeps.h
 $(COLVARS_OBJ_DIR)colvarcomp_alchlambda.o: colvarcomp_alchlambda.cpp \
  colvarmodule.h colvars_version.h colvarvalue.h colvartypes.h \
  ../../src/math_eigen_impl.h colvar.h colvarparse.h colvarparams.h \
@@ -101,6 +107,11 @@ $(COLVARS_OBJ_DIR)colvarcomp_neuralnetwork.o: \
  colvarproxy.h colvarproxy_io.h colvarproxy_system.h colvarproxy_tcl.h \
  colvarproxy_volmaps.h colvar_geometricpath.h \
  colvar_neuralnetworkcompute.h
+$(COLVARS_OBJ_DIR)colvarcomp_torchann.o: \
+ colvarcomp_torchann.cpp colvarmodule.h colvars_version.h \
+ colvarvalue.h colvartypes.h colvarparse.h colvarparams.h colvar.h \
+ colvardeps.h colvarcomp.h colvarcomp_torchann.h colvaratoms.h colvarproxy.h colvarproxy_io.h \
+ colvarproxy_system.h colvarproxy_tcl.h 
 $(COLVARS_OBJ_DIR)colvarcomp_combination.o: colvarcomp_combination.cpp \
  colvarcomp.h colvarmodule.h colvars_version.h colvaratoms.h \
  colvarproxy.h colvartypes.h ../../src/math_eigen_impl.h colvarproxy_io.h \
@@ -127,7 +138,7 @@ $(COLVARS_OBJ_DIR)colvarcomp_volmaps.o: colvarcomp_volmaps.cpp \
  colvar_geometricpath.h
 $(COLVARS_OBJ_DIR)colvar.o: colvar.cpp colvarmodule.h colvars_version.h \
  colvarvalue.h colvartypes.h ../../src/math_eigen_impl.h colvarparse.h \
- colvarparams.h colvarcomp.h colvaratoms.h colvarproxy.h colvarproxy_io.h \
+ colvarparams.h colvarcomp.h colvarcomp_torchann.h colvaratoms.h colvarproxy.h colvarproxy_io.h \
  colvarproxy_system.h colvarproxy_tcl.h colvarproxy_volmaps.h \
  colvardeps.h colvar.h colvar_geometricpath.h colvarbias.h \
  colvars_memstream.h
@@ -152,7 +163,8 @@ $(COLVARS_OBJ_DIR)colvarmodule.o: colvarmodule.cpp colvarmodule.h \
  colvarbias_histogram_reweight_amd.h colvarbias_meta.h colvarscript.h \
  colvarscript_commands.h colvarscript_commands_colvar.h \
  colvarscript_commands_bias.h colvaratoms.h colvarcomp.h \
- colvar_geometricpath.h colvars_memstream.h colvarmodule_refs.h
+ colvar_geometricpath.h colvars_memstream.h colvarmodule_refs.h \
+ colvarbias_opes.h
 $(COLVARS_OBJ_DIR)colvarparams.o: colvarparams.cpp colvarmodule.h \
  colvars_version.h colvarvalue.h colvartypes.h \
  ../../src/math_eigen_impl.h colvarparams.h

lib/colvars/colvar.cpp

@@ -21,6 +21,7 @@
 #include "colvarbias.h"
 #include "colvars_memstream.h"
 
+#include "colvarcomp_torchann.h"
 
 std::map<std::string, std::function<colvar::cvc *()>> colvar::global_cvc_map =
     std::map<std::string, std::function<colvar::cvc *()>>();
@@ -95,6 +96,12 @@ int colvar::init(std::string const &conf)
   if (error_code != COLVARS_OK) {
     return cvm::get_error();
   }
+#else
+  if (key_lookup(conf, "customFunction")) {
+    return cvm::error(
+        "Error: customFunction keyword is used, but the Lepton library is not available.\n",
+        COLVARS_NOT_IMPLEMENTED);
+  }
 #endif
 
   // Setup colvar as scripted function of components
@@ -175,12 +182,6 @@ int colvar::init(std::string const &conf)
 
   set_enabled(f_cv_scalar, (value().type() == colvarvalue::type_scalar));
 
-  // If using scripted biases, any colvar may receive bias forces
-  // and will need its gradient
-  if (cvm::scripted_forces()) {
-    enable(f_cv_gradient);
-  }
-
   // check for linear combinations
   {
     bool lin = !(is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function));
@@ -311,9 +312,27 @@ int colvar::init(std::string const &conf)
   // Detect if we have a single component that is an alchemical lambda
   if (is_enabled(f_cv_single_cvc) && cvcs[0]->function_type() == "alchLambda") {
     enable(f_cv_external);
+
+    static_cast<colvar::alch_lambda *>(cvcs[0].get())->init_alchemy(time_step_factor);
+  }
+
+  // If using scripted biases, any colvar may receive bias forces
+  if (cvm::scripted_forces()) {
+    enable(f_cv_apply_force);
   }
 
   error_code |= init_extended_Lagrangian(conf);
+
+  // when total atomic forces are obtained from the previous time step,
+  // we cannot (currently) have colvar values and projected total forces for the same timestep
+  // (that would require anticipating the total force request by one timestep)
+  // i.e. the combination of f_cv_total_force_calc and f_cv_multiple_ts requires f_cv_total_force_current_step
+  // Because f_cv_total_force_current_step is static, we can hard-code this, once other features are set
+  // that is f_cv_external and f_cv_extended_Lagrangian
+  if (!is_enabled(f_cv_total_force_current_step)) {
+    exclude_feature_self(f_cv_multiple_ts, f_cv_total_force_calc);
+  }
+
   error_code |= init_output_flags(conf);
 
   // Now that the children are defined we can solve dependencies
@@ -495,8 +514,6 @@ int colvar::init_grid_parameters(std::string const &conf)
 {
   int error_code = COLVARS_OK;
 
-  colvarmodule *cv = cvm::main();
-
   cvm::real default_width = width;
 
   if (!key_already_set("width")) {
@@ -522,34 +539,68 @@ int colvar::init_grid_parameters(std::string const &conf)
 
   if (is_enabled(f_cv_scalar)) {
 
-    if (is_enabled(f_cv_single_cvc)) {
+    // Record the CVC's intrinsic boundaries, and set them as default values for the user's choice
-      // Get the default boundaries from the component
+    colvarvalue cvc_lower_boundary, cvc_upper_boundary;
+
+    if (is_enabled(f_cv_single_cvc)) { // Get the intrinsic boundaries of the CVC
+
       if (cvcs[0]->is_enabled(f_cvc_lower_boundary)) {
         enable(f_cv_lower_boundary);
         enable(f_cv_hard_lower_boundary);
-        lower_boundary =
+        lower_boundary = cvc_lower_boundary =
           *(reinterpret_cast<colvarvalue const *>(cvcs[0]->get_param_ptr("lowerBoundary")));
       }
+
       if (cvcs[0]->is_enabled(f_cvc_upper_boundary)) {
         enable(f_cv_upper_boundary);
         enable(f_cv_hard_upper_boundary);
-        upper_boundary =
+        upper_boundary = cvc_upper_boundary =
-          *(reinterpret_cast<colvarvalue const *>(cvcs[0]->get_param_ptr("upperBoundary")));
+            *(reinterpret_cast<colvarvalue const *>(cvcs[0]->get_param_ptr("upperBoundary")));
       }
     }
 
     if (get_keyval(conf, "lowerBoundary", lower_boundary, lower_boundary)) {
       enable(f_cv_lower_boundary);
-      // Because this is the user's choice, we cannot assume it is a true
+      if (is_enabled(f_cv_single_cvc) && is_enabled(f_cv_hard_lower_boundary)) {
-      // physical boundary
+        if (cvm::sqrt(dist2(lower_boundary, cvc_lower_boundary))/width > colvar_boundaries_tol)  {
-      disable(f_cv_hard_lower_boundary);
+          // The user choice is different from the CVC's default
+          disable(f_cv_hard_lower_boundary);
+        }
+      }
     }
 
     if (get_keyval(conf, "upperBoundary", upper_boundary, upper_boundary)) {
       enable(f_cv_upper_boundary);
-      disable(f_cv_hard_upper_boundary);
+      if (is_enabled(f_cv_single_cvc) && is_enabled(f_cv_hard_upper_boundary)) {
+        if (cvm::sqrt(dist2(upper_boundary, cvc_upper_boundary))/width > colvar_boundaries_tol)  {
+          disable(f_cv_hard_upper_boundary);
+        }
+      }
     }
 
+    get_keyval_feature(this, conf, "hardLowerBoundary", f_cv_hard_lower_boundary,
+                       is_enabled(f_cv_hard_lower_boundary));
+
+    get_keyval_feature(this, conf, "hardUpperBoundary", f_cv_hard_upper_boundary,
+                       is_enabled(f_cv_hard_upper_boundary));
+
+    get_keyval(conf, "expandBoundaries", expand_boundaries, expand_boundaries);
+
+    error_code |= parse_legacy_wall_params(conf);
+    error_code |= check_grid_parameters();
+  }
+
+  return error_code;
+}
+
+
+int colvar::parse_legacy_wall_params(std::string const &conf)
+{
+  int error_code = COLVARS_OK;
+  colvarmodule *cv = cvm::main();
+
+  if (is_enabled(f_cv_scalar)) {
+
     // Parse legacy wall options and set up a harmonicWalls bias if needed
     cvm::real lower_wall_k = 0.0, upper_wall_k = 0.0;
     cvm::real lower_wall = 0.0, upper_wall = 0.0;
@@ -603,13 +654,14 @@ harmonicWalls {\n\
     }
   }
 
-  get_keyval_feature(this, conf, "hardLowerBoundary", f_cv_hard_lower_boundary,
+  return error_code;
-                     is_enabled(f_cv_hard_lower_boundary));
+}
 
-  get_keyval_feature(this, conf, "hardUpperBoundary", f_cv_hard_upper_boundary,
-                     is_enabled(f_cv_hard_upper_boundary));
 
-  // consistency checks for boundaries and walls
+int colvar::check_grid_parameters()
+{
+  int error_code = COLVARS_OK;
+
   if (is_enabled(f_cv_lower_boundary) && is_enabled(f_cv_upper_boundary)) {
     if (lower_boundary >= upper_boundary) {
       error_code |= cvm::error("Error: the upper boundary, "+
@@ -620,7 +672,6 @@ harmonicWalls {\n\
     }
   }
 
-  get_keyval(conf, "expandBoundaries", expand_boundaries, expand_boundaries);
   if (expand_boundaries && periodic_boundaries()) {
     error_code |= cvm::error("Error: trying to expand boundaries that already "
                              "cover a whole period of a periodic colvar.\n",
@@ -654,14 +705,15 @@ int colvar::init_extended_Lagrangian(std::string const &conf)
     x_ext.type(colvarvalue::type_notset);
     v_ext.type(value());
     fr.type(value());
-    const bool temp_provided = get_keyval(conf, "extendedTemp", temp,
+    const bool temp_provided = get_keyval(conf, "extendedTemp", temp, proxy->target_temperature());
-                                          proxy->target_temperature());
     if (is_enabled(f_cv_external)) {
-      // In the case of an "external" coordinate, there is no coupling potential:
+      // In the case of a driven external parameter in the back-end, there is no coupling potential:
       // only the fictitious mass is meaningful
       get_keyval(conf, "extendedMass", ext_mass);
       // Ensure that the computed restraint energy term is zero
       ext_force_k = 0.0;
+      // Then we need forces from the back-end
+      enable(f_cv_total_force_calc);
     } else {
       // Standard case of coupling to a geometric colvar
       if (temp <= 0.0) { // Then a finite temperature is required
@@ -779,6 +831,7 @@ int colvar::init_components_type(const std::string& conf, const char* def_config
                            &def_conf,
                            &pos) ) {
 
+    cvm::increase_depth();
     cvm::log("Initializing "
              "a new \""+std::string(def_config_key)+"\" component"+
              (cvm::debug() ? ", with configuration:\n"+def_conf
@@ -791,7 +844,6 @@ int colvar::init_components_type(const std::string& conf, const char* def_config
     }
     cvcs.push_back(std::shared_ptr<colvar::cvc>(cvcp));
 
-    cvm::increase_depth();
     int error_code_this = cvcp->init(def_conf);
     if (error_code_this == COLVARS_OK) {
       // Checking for invalid keywords only if the parsing was successful, otherwise any
@@ -851,12 +903,8 @@ void colvar::define_component_types()
   add_component_type<dipole_angle>("dipole angle", "dipoleAngle");
   add_component_type<dihedral>("dihedral", "dihedral");
   add_component_type<h_bond>("hydrogen bond", "hBond");
-
+  add_component_type<alpha_angles>("alpha helix", "alpha");
-  if (proxy->check_atom_name_selections_available() == COLVARS_OK) {
+  add_component_type<dihedPC>("dihedral principal component", "dihedralPC");
-    add_component_type<alpha_angles>("alpha helix", "alpha");
-    add_component_type<dihedPC>("dihedral principal component", "dihedralPC");
-  }
-
   add_component_type<orientation>("orientation", "orientation");
   add_component_type<orientation_angle>("orientation angle", "orientationAngle");
   add_component_type<orientation_proj>("orientation projection", "orientationProj");
@@ -888,6 +936,8 @@ void colvar::define_component_types()
 
   add_component_type<neuralNetwork>("neural network CV for other CVs", "neuralNetwork");
 
+  add_component_type<torchANN>("CV defined by PyTorch artifical neural network models", "torchANN");
+
   if (proxy->check_volmaps_available() == COLVARS_OK) {
     add_component_type<map_total>("total value of atomic map", "mapTotal");
   }
@@ -1098,6 +1148,9 @@ int colvar::init_dependencies() {
     init_feature(f_cv_gradient, "gradient", f_type_dynamic);
     require_feature_children(f_cv_gradient, f_cvc_gradient);
 
+    init_feature(f_cv_apply_force, "apply_force", f_type_dynamic);
+    require_feature_alt(f_cv_apply_force, f_cv_gradient, f_cv_external);
+
     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);
@@ -1116,6 +1169,10 @@ int colvar::init_dependencies() {
     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);
 
+    // If this is active, the total force reported to biases (ABF / TI) is from the current step
+    // therefore it does not include Colvars biases -> it is a "system force"
+    init_feature(f_cv_total_force_current_step, "total_force_current_step", f_type_dynamic);
+
     // 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);
@@ -1134,13 +1191,15 @@ int colvar::init_dependencies() {
 
     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);
+    require_feature_self(f_cv_extended_Lagrangian, f_cv_apply_force);
 
     init_feature(f_cv_Langevin, "Langevin_dynamics", f_type_user);
     require_feature_self(f_cv_Langevin, f_cv_extended_Lagrangian);
 
-    init_feature(f_cv_external, "external", f_type_user);
+    init_feature(f_cv_external, "external_parameter", f_type_static);
     require_feature_self(f_cv_external, f_cv_single_cvc);
+    // External parameters always report the total force for current step
+    require_feature_self(f_cv_external, f_cv_total_force_current_step);
 
     init_feature(f_cv_single_cvc, "single_component", f_type_static);
 
@@ -1201,10 +1260,7 @@ int colvar::init_dependencies() {
     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", 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++) {
@@ -1225,6 +1281,10 @@ int colvar::init_dependencies() {
   feature_states[f_cv_fdiff_velocity].available =
     cvm::main()->proxy->simulation_running();
 
+  // Some back-ends report current total forces for all colvars
+  if (cvm::main()->proxy->total_forces_same_step())
+    enable(f_cv_total_force_current_step);
+
   return COLVARS_OK;
 }
 
@@ -1351,7 +1411,6 @@ int colvar::calc_cvcs(int first_cvc, size_t num_cvcs)
     cvm::log("Calculating colvar \""+this->name+"\", components "+
              cvm::to_str(first_cvc)+" through "+cvm::to_str(first_cvc+num_cvcs)+".\n");
 
-  colvarproxy *proxy = cvm::main()->proxy;
   int error_code = COLVARS_OK;
 
   error_code |= check_cvc_range(first_cvc, num_cvcs);
@@ -1359,7 +1418,7 @@ int colvar::calc_cvcs(int first_cvc, size_t num_cvcs)
     return error_code;
   }
 
-  if ((cvm::step_relative() > 0) && (!proxy->total_forces_same_step())){
+  if ((cvm::step_relative() > 0) && (!is_enabled(f_cv_total_force_current_step))){
     // Use Jacobian derivative from previous timestep
     error_code |= calc_cvc_total_force(first_cvc, num_cvcs);
   }
@@ -1367,7 +1426,7 @@ int colvar::calc_cvcs(int first_cvc, size_t num_cvcs)
   error_code |= calc_cvc_values(first_cvc, num_cvcs);
   error_code |= calc_cvc_gradients(first_cvc, num_cvcs);
   error_code |= calc_cvc_Jacobians(first_cvc, num_cvcs);
-  if (proxy->total_forces_same_step()){
+  if (is_enabled(f_cv_total_force_current_step)){
     // Use Jacobian derivative from this timestep
     error_code |= calc_cvc_total_force(first_cvc, num_cvcs);
   }
@@ -1384,10 +1443,9 @@ int colvar::collect_cvc_data()
   if (cvm::debug())
     cvm::log("Calculating colvar \""+this->name+"\"'s properties.\n");
 
-  colvarproxy *proxy = cvm::main()->proxy;
   int error_code = COLVARS_OK;
 
-  if ((cvm::step_relative() > 0) && (!proxy->total_forces_same_step())){
+  if ((cvm::step_relative() > 0) && (!is_enabled(f_cv_total_force_current_step))){
     // Total force depends on Jacobian derivative from previous timestep
     // collect_cvc_total_forces() uses the previous value of jd
     error_code |= collect_cvc_total_forces();
@@ -1395,7 +1453,7 @@ int colvar::collect_cvc_data()
   error_code |= collect_cvc_values();
   error_code |= collect_cvc_gradients();
   error_code |= collect_cvc_Jacobians();
-  if (proxy->total_forces_same_step()){
+  if (is_enabled(f_cv_total_force_current_step)){
     // Use Jacobian derivative from this timestep
     error_code |= collect_cvc_total_forces();
   }
@@ -1609,22 +1667,20 @@ int colvar::collect_cvc_total_forces()
   if (is_enabled(f_cv_total_force_calc)) {
     ft.reset();
 
-    if (cvm::step_relative() > 0) {
+    for (size_t i = 0; i < cvcs.size();  i++) {
-      // get from the cvcs the total forces from the PREVIOUS step
+      if (!cvcs[i]->is_enabled()) continue;
-      for (size_t i = 0; i < cvcs.size();  i++) {
+          if (cvm::debug())
-        if (!cvcs[i]->is_enabled()) continue;
+          cvm::log("Colvar component no. "+cvm::to_str(i+1)+
-            if (cvm::debug())
+              " within colvar \""+this->name+"\" has total force "+
-            cvm::log("Colvar component no. "+cvm::to_str(i+1)+
+              cvm::to_str((cvcs[i])->total_force(),
-                " within colvar \""+this->name+"\" has total force "+
+              cvm::cv_width, cvm::cv_prec)+".\n");
-                cvm::to_str((cvcs[i])->total_force(),
+      // linear combination is assumed
-                cvm::cv_width, cvm::cv_prec)+".\n");
+      ft += (cvcs[i])->total_force() * (cvcs[i])->sup_coeff / active_cvc_square_norm;
-        // linear combination is assumed
-        ft += (cvcs[i])->total_force() * (cvcs[i])->sup_coeff / active_cvc_square_norm;
-      }
     }
 
     if (!(is_enabled(f_cv_hide_Jacobian) && is_enabled(f_cv_subtract_applied_force))) {
-      // add the Jacobian force to the total force, and don't apply any silent
+      // This is by far the most common case
+      // Add the Jacobian force to the total force, and don't apply any silent
       // correction internally: biases such as colvarbias_abf will handle it
       // If f_cv_hide_Jacobian is enabled, a force of -fj is present in ft due to the
       // Jacobian-compensating force
@@ -1632,6 +1688,10 @@ int colvar::collect_cvc_total_forces()
     }
   }
 
+  if (is_enabled(f_cv_total_force_current_step)) {
+   // Report total force value without waiting for calc_colvar_properties()
+    ft_reported = ft;
+  }
   return COLVARS_OK;
 }
 
@@ -1733,12 +1793,15 @@ int colvar::calc_colvar_properties()
     // But we report values at the beginning of the timestep (value at t=0 on the first timestep)
     x_reported = x_ext;
     v_reported = v_ext;
-    // the "total force" with the extended Lagrangian is
-    // calculated in update_forces_energy() below
 
+    // the "total force" for the extended Lagrangian is calculated in update_forces_energy() below
+    // A future improvement could compute a "system force" here, borrowing a part of update_extended_Lagrangian()
+    // this would change the behavior of eABF with respect to other biases
+    // by enabling f_cv_total_force_current_step, and reducing the total force to a system force
+    // giving the behavior of f_cv_subtract_applied_force - this is correct for WTM-eABF etc.
   } else {
 
-    if (is_enabled(f_cv_subtract_applied_force)) {
+    if (is_enabled(f_cv_subtract_applied_force) && !cvm::proxy->total_forces_same_step()) {
       // correct the total force only if it has been measured
       // TODO add a specific test instead of relying on sq norm
       if (ft.norm2() > 0.0) {
@@ -1825,7 +1888,8 @@ void colvar::update_extended_Lagrangian()
   // Integrate with slow timestep (if time_step_factor != 1)
   cvm::real dt = cvm::dt() * cvm::real(time_step_factor);
 
-  colvarvalue f_ext(fr.type()); // force acting on the extended variable
+  // Force acting on the extended variable
+  colvarvalue f_ext(fr.type());
   f_ext.reset();
 
   if (is_enabled(f_cv_external)) {
@@ -1834,13 +1898,13 @@ void colvar::update_extended_Lagrangian()
     f += fb_actual;
   }
 
-  // fr: bias force on extended variable (without harmonic spring), for output in trajectory
-  fr = f;
-
   // External force has been scaled for an inner-timestep impulse (for the back-end integrator)
   // here we scale it back because this integrator uses only the outer (long) timestep
   f_ext = f / cvm::real(time_step_factor);
 
+  // fr: bias force on extended variable (without harmonic spring), for output in trajectory
+  fr = f_ext;
+
   colvarvalue f_system(fr.type()); // force exterted by the system on the extended DOF
 
   if (is_enabled(f_cv_external)) {
@@ -1863,14 +1927,18 @@ void colvar::update_extended_Lagrangian()
   }
   f_ext += f_system;
 
-  if (is_enabled(f_cv_subtract_applied_force)) {
+  if ( ! is_enabled(f_cv_total_force_current_step)) {
-    // Report a "system" force without the biases on this colvar
+    if (is_enabled(f_cv_subtract_applied_force)) {
-    // that is, just the spring force (or alchemical force)
+      // Report a "system" force without the biases on this colvar
-    ft_reported = f_system;
+      // that is, just the spring force (or alchemical force)
-  } else {
+      ft_reported = f_system;
-    // The total force acting on the extended variable is f_ext
+    } else {
-    // This will be used in the next timestep
+      // The total force acting on the extended variable is f_ext
-    ft_reported = f_ext;
+      // This will be used in the next timestep
+      ft_reported = f_ext;
+    }
+    // Since biases have already been updated, this ft_reported will only be
+    // communicated to biases at the next timestep
   }
 
   // backup in case we need to revert this integration timestep
@@ -2184,12 +2252,10 @@ int colvar::set_cvc_param(std::string const &param_name, void const *new_value)
 bool colvar::periodic_boundaries(colvarvalue const &lb, colvarvalue const &ub) const
 {
   if (period > 0.0) {
-    if ( ((cvm::sqrt(this->dist2(lb, ub))) / this->width)
+    if (((cvm::sqrt(this->dist2(lb, ub))) / this->width) < colvar_boundaries_tol) {
-         < 1.0E-10 ) {
       return true;
     }
   }
-
   return false;
 }
 
@@ -2347,6 +2413,11 @@ int colvar::set_state_params(std::string const &conf)
              cvm::to_str(x)+"\n");
     x_restart = x;
     after_restart = true;
+    // Externally driven cv (e.g. alchemical lambda) is imposed by restart value
+    if (is_enabled(f_cv_external) && is_enabled(f_cv_extended_Lagrangian)) {
+      // Request immediate sync of driven parameter to back-end code
+      cvcs[0]->set_value(x, true);
+    }
   }
 
   if (is_enabled(f_cv_extended_Lagrangian)) {
@@ -2489,8 +2560,14 @@ std::string const colvar::get_state_params() const
   os << "  name " << name << "\n"
      << "  x "
      << std::setprecision(cvm::cv_prec)
-     << std::setw(cvm::cv_width)
+     << std::setw(cvm::cv_width);
-     << x << "\n";
+  if (is_enabled(f_cv_external) && is_enabled(f_cv_extended_Lagrangian)) {
+    // For an external colvar, x is one timestep in the future after integration
+    // write x at beginning of timestep
+    os << x_reported << "\n";
+  } else {
+    os << x << "\n";
+  }
 
   if (is_enabled(f_cv_output_velocity)) {
     os << "  v "

lib/colvars/colvar.h

@@ -263,6 +263,12 @@ public:
   /// Init defaults for grid options
   int init_grid_parameters(std::string const &conf);
 
+  /// Consistency check for the grid paramaters
+  int check_grid_parameters();
+
+  /// Read legacy wall keyword (these are biases now)
+  int parse_legacy_wall_params(std::string const &conf);
+
   /// Init extended Lagrangian parameters
   int init_extended_Lagrangian(std::string const &conf);
 
@@ -633,6 +639,7 @@ public:
   class euler_psi;
   class euler_theta;
   class neuralNetwork;
+  class torchANN;
   class customColvar;
 
   // non-scalar components
@@ -753,7 +760,7 @@ inline colvarvalue const & colvar::total_force() const
 
 inline void colvar::add_bias_force(colvarvalue const &force)
 {
-  check_enabled(f_cv_gradient,
+  check_enabled(f_cv_apply_force,
                 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");
@@ -778,4 +785,10 @@ inline void colvar::reset_bias_force() {
   fb_actual.reset();
 }
 
+
+namespace {
+  // Tolerance parameter to decide when two boundaries coincide
+  constexpr cvm::real colvar_boundaries_tol = 1.0e-10;
+}
+
 #endif

lib/colvars/colvar_rotation_derivative.h

@@ -5,11 +5,21 @@
 #include <type_traits>
 #include <cstring>
 
+#ifndef _noalias
+#if defined(__INTEL_COMPILER) || (defined(__PGI) && !defined(__NVCOMPILER))
+#define _noalias restrict
+#elif defined(__GNUC__) || defined(__INTEL_LLVM_COMPILER) || defined(__NVCOMPILER)
+#define _noalias __restrict
+#else
+#define _noalias
+#endif
+#endif
+
 /// \brief Helper function for loading the ia-th atom in the vector pos to x, y and z (C++11 SFINAE is used)
 template <typename T, typename std::enable_if<std::is_same<T, cvm::atom_pos>::value, bool>::type = true>
 inline void read_atom_coord(
   size_t ia, const std::vector<T>& pos,
-  cvm::real* x, cvm::real* y, cvm::real* z) {
+  cvm::real* _noalias x, cvm::real* _noalias y, cvm::real* _noalias z) {
   *x = pos[ia].x;
   *y = pos[ia].y;
   *z = pos[ia].z;
@@ -18,7 +28,7 @@ inline void read_atom_coord(
 template <typename T, typename std::enable_if<std::is_same<T, cvm::atom>::value, bool>::type = true>
 inline void read_atom_coord(
   size_t ia, const std::vector<T>& pos,
-  cvm::real* x, cvm::real* y, cvm::real* z) {
+  cvm::real* _noalias x, cvm::real* _noalias y, cvm::real* _noalias z) {
   *x = pos[ia].pos.x;
   *y = pos[ia].pos.y;
   *z = pos[ia].pos.z;
@@ -26,9 +36,9 @@ inline void read_atom_coord(
 
 /// \brief Helper enum class for specifying options in rotation_derivative::prepare_derivative
 enum class rotation_derivative_dldq {
-  /// Require the derivative of the leading eigenvalue with respect to the atom coordinats
+  /// Require the derivative of the leading eigenvalue with respect to the atom coordinates
   use_dl = 1 << 0,
-  /// Require the derivative of the leading eigenvector with respect to the atom coordinats
+  /// Require the derivative of the leading eigenvector with respect to the atom coordinates
   use_dq = 1 << 1
 };
 
@@ -327,12 +337,13 @@ struct rotation_derivative {
     *  @param[out] dq0_out The output of derivative of Q
     *  @param[out] ds_out  The output of derivative of overlap matrix S
     */
+  template <bool use_dl, bool use_dq, bool use_ds>
   void calc_derivative_impl(
     const cvm::rvector (&ds)[4][4],
-    cvm::rvector* const dl0_out,
+    cvm::rvector* _noalias const dl0_out,
-    cvm::vector1d<cvm::rvector>* const dq0_out,
+    cvm::vector1d<cvm::rvector>* _noalias const dq0_out,
-    cvm::matrix2d<cvm::rvector>* const ds_out) const {
+    cvm::matrix2d<cvm::rvector>* _noalias const ds_out) const {
-    if (ds_out != nullptr) {
+    if (use_ds) {
       // this code path is for debug_gradients, so not necessary to unroll the loop
       *ds_out = cvm::matrix2d<cvm::rvector>(4, 4);
       for (int i = 0; i < 4; ++i) {
@@ -341,7 +352,7 @@ struct rotation_derivative {
         }
       }
     }
-    if (dl0_out != nullptr) {
+    if (use_dl) {
       /* manually loop unrolling of the following loop:
         dl0_1.reset();
         for (size_t i = 0; i < 4; i++) {
@@ -367,7 +378,7 @@ struct rotation_derivative {
                  tmp_Q0Q0[3][2] * ds[3][2] +
                  tmp_Q0Q0[3][3] * ds[3][3];
     }
-    if (dq0_out != nullptr) {
+    if (use_dq) {
       // we can skip this check if a fixed-size array is used
       if (dq0_out->size() != 4) dq0_out->resize(4);
       /* manually loop unrolling of the following loop:
@@ -462,32 +473,21 @@ struct rotation_derivative {
     *  @param[out] ds_1_out  The output of derivative of overlap matrix S with
     *                        respect to ia-th atom of group 1
     */
+  template <bool use_dl, bool use_dq, bool use_ds>
   void calc_derivative_wrt_group1(
-    size_t ia, cvm::rvector* const dl0_1_out = nullptr,
+    size_t ia, cvm::rvector* _noalias const dl0_1_out = nullptr,
-    cvm::vector1d<cvm::rvector>* const dq0_1_out = nullptr,
+    cvm::vector1d<cvm::rvector>* _noalias const dq0_1_out = nullptr,
-    cvm::matrix2d<cvm::rvector>* const ds_1_out = nullptr) const {
+    cvm::matrix2d<cvm::rvector>* _noalias const ds_1_out = nullptr) const {
-      if (dl0_1_out == nullptr && dq0_1_out == nullptr) return;
+      // if (dl0_1_out == nullptr && dq0_1_out == nullptr) return;
       cvm::real a2x, a2y, a2z;
       // we can get rid of the helper function read_atom_coord if C++17 (constexpr) is available
       read_atom_coord(ia, m_pos2, &a2x, &a2y, &a2z);
-      cvm::rvector ds_1[4][4];
+      const cvm::rvector ds_1[4][4] = {
-      ds_1[0][0].set( a2x,  a2y,  a2z);
+        {{ a2x,  a2y,  a2z}, { 0.0, a2z,  -a2y}, {-a2z,  0.0,  a2x}, { a2y, -a2x,  0.0}},
-      ds_1[1][0].set( 0.0,  a2z, -a2y);
+        {{ 0.0,  a2z, -a2y}, { a2x, -a2y, -a2z}, { a2y,  a2x,  0.0}, { a2z,  0.0,  a2x}},
-      ds_1[0][1] = ds_1[1][0];
+        {{-a2z,  0.0,  a2x}, { a2y,  a2x,  0.0}, {-a2x,  a2y, -a2z}, { 0.0,  a2z,  a2y}},
-      ds_1[2][0].set(-a2z,  0.0,  a2x);
+        {{ a2y, -a2x,  0.0}, { a2z,  0.0,  a2x}, { 0.0,  a2z,  a2y}, {-a2x, -a2y,  a2z}}};
-      ds_1[0][2] = ds_1[2][0];
+      calc_derivative_impl<use_dl, use_dq, use_ds>(ds_1, dl0_1_out, dq0_1_out, ds_1_out);
-      ds_1[3][0].set( a2y, -a2x,  0.0);
-      ds_1[0][3] = ds_1[3][0];
-      ds_1[1][1].set( a2x, -a2y, -a2z);
-      ds_1[2][1].set( a2y,  a2x,  0.0);
-      ds_1[1][2] = ds_1[2][1];
-      ds_1[3][1].set( a2z,  0.0,  a2x);
-      ds_1[1][3] = ds_1[3][1];
-      ds_1[2][2].set(-a2x,  a2y, -a2z);
-      ds_1[3][2].set( 0.0,  a2z,  a2y);
-      ds_1[2][3] = ds_1[3][2];
-      ds_1[3][3].set(-a2x, -a2y,  a2z);
-      calc_derivative_impl(ds_1, dl0_1_out, dq0_1_out, ds_1_out);
     }
   /*! @brief Calculate the derivatives of S, the leading eigenvalue L and
     *         the leading eigenvector Q with respect to `m_pos2`
@@ -499,32 +499,21 @@ struct rotation_derivative {
     *  @param[out] ds_2_out  The output of derivative of overlap matrix S with
     *                        respect to ia-th atom of group 2
     */
+  template <bool use_dl, bool use_dq, bool use_ds>
   void calc_derivative_wrt_group2(
-    size_t ia, cvm::rvector* const dl0_2_out = nullptr,
+    size_t ia, cvm::rvector* _noalias const dl0_2_out = nullptr,
-    cvm::vector1d<cvm::rvector>* const dq0_2_out = nullptr,
+    cvm::vector1d<cvm::rvector>* _noalias const dq0_2_out = nullptr,
-    cvm::matrix2d<cvm::rvector>* const ds_2_out = nullptr) const {
+    cvm::matrix2d<cvm::rvector>* _noalias const ds_2_out = nullptr) const {
-    if (dl0_2_out == nullptr && dq0_2_out == nullptr) return;
+    // if (dl0_2_out == nullptr && dq0_2_out == nullptr) return;
     cvm::real a1x, a1y, a1z;
     // we can get rid of the helper function read_atom_coord if C++17 (constexpr) is available
     read_atom_coord(ia, m_pos1, &a1x, &a1y, &a1z);
-    cvm::rvector ds_2[4][4];
+    const cvm::rvector ds_2[4][4] = {
-    ds_2[0][0].set( a1x,  a1y,  a1z);
+      {{ a1x,  a1y,  a1z}, { 0.0, -a1z,  a1y}, { a1z,  0.0, -a1x}, {-a1y,  a1x,  0.0}},
-    ds_2[1][0].set( 0.0, -a1z,  a1y);
+      {{ 0.0, -a1z,  a1y}, { a1x, -a1y, -a1z}, { a1y,  a1x,  0.0}, { a1z,  0.0,  a1x}},
-    ds_2[0][1] = ds_2[1][0];
+      {{ a1z,  0.0, -a1x}, { a1y,  a1x,  0.0}, {-a1x,  a1y, -a1z}, { 0.0,  a1z,  a1y}},
-    ds_2[2][0].set( a1z,  0.0, -a1x);
+      {{-a1y,  a1x,  0.0}, { a1z,  0.0,  a1x}, { 0.0,  a1z,  a1y}, {-a1x, -a1y,  a1z}}};
-    ds_2[0][2] = ds_2[2][0];
+    calc_derivative_impl<use_dl, use_dq, use_ds>(ds_2, dl0_2_out, dq0_2_out, ds_2_out);
-    ds_2[3][0].set(-a1y,  a1x,  0.0);
-    ds_2[0][3] = ds_2[3][0];
-    ds_2[1][1].set( a1x, -a1y, -a1z);
-    ds_2[2][1].set( a1y,  a1x,  0.0);
-    ds_2[1][2] = ds_2[2][1];
-    ds_2[3][1].set( a1z,  0.0,  a1x);
-    ds_2[1][3] = ds_2[3][1];
-    ds_2[2][2].set(-a1x,  a1y, -a1z);
-    ds_2[3][2].set( 0.0,  a1z,  a1y);
-    ds_2[2][3] = ds_2[3][2];
-    ds_2[3][3].set(-a1x, -a1y,  a1z);
-    calc_derivative_impl(ds_2, dl0_2_out, dq0_2_out, ds_2_out);
   }
 };
 
@@ -585,10 +574,7 @@ void debug_gradients(
   cvm::real S_new_eigval[4];
   cvm::real S_new_eigvec[4][4];
   for (size_t ia = 0; ia < pos2.size(); ++ia) {
-    // cvm::real const &a1x = pos1[ia].x;
+    deriv.template calc_derivative_wrt_group2<true, true, true>(ia, &dl0_2, &dq0_2, &ds_2);
-    // cvm::real const &a1y = pos1[ia].y;
-    // cvm::real const &a1z = pos1[ia].z;
-    deriv.calc_derivative_wrt_group2(ia, &dl0_2, &dq0_2, &ds_2);
     // make an infitesimal move along each cartesian coordinate of
     // this atom, and solve again the eigenvector problem
     for (size_t comp = 0; comp < 3; comp++) {

lib/colvars/colvaratoms.cpp

@@ -673,7 +673,7 @@ int cvm::atom_group::add_atom_numbers(std::string const &numbers_conf)
 }
 
 
-int cvm::atom_group::add_index_group(std::string const &index_group_name)
+int cvm::atom_group::add_index_group(std::string const &index_group_name, bool silent)
 {
   std::vector<std::string> const &index_group_names =
     cvm::main()->index_group_names;
@@ -687,7 +687,10 @@ int cvm::atom_group::add_index_group(std::string const &index_group_name)
   }
 
   if (i_group >= index_group_names.size()) {
-    return cvm::error("Error: could not find index group "+
+    if (silent)
+      return COLVARS_INPUT_ERROR;
+    else
+      return cvm::error("Error: could not find index group "+
                       index_group_name+" among those already provided.\n",
                       COLVARS_INPUT_ERROR);
   }
@@ -1055,6 +1058,14 @@ void cvm::atom_group::calc_apply_roto_translation()
     }
   }
 
+  if (is_enabled(f_ag_fit_gradients) && !b_dummy) {
+    // Save the unrotated frame for fit gradients
+    pos_unrotated.resize(size());
+    for (size_t i = 0; i < size(); ++i) {
+      pos_unrotated[i] = atoms[i].pos;
+    }
+  }
+
   if (is_enabled(f_ag_rotate)) {
     // rotate the group (around the center of geometry if f_ag_center is
     // enabled, around the origin otherwise)
@@ -1217,23 +1228,30 @@ void cvm::atom_group::calc_fit_gradients()
   if (cvm::debug())
     cvm::log("Calculating fit gradients.\n");
 
+  cvm::atom_group *group_for_fit = fitting_group ? fitting_group : this;
+
+  auto accessor_main = [this](size_t i){return atoms[i].grad;};
+  auto accessor_fitting = [&group_for_fit](size_t j, const cvm::rvector& grad){group_for_fit->fit_gradients[j] = grad;};
   if (is_enabled(f_ag_center) && is_enabled(f_ag_rotate))
-    calc_fit_gradients_impl<true, true>();
+    calc_fit_forces_impl<true, true>(accessor_main, accessor_fitting);
   if (is_enabled(f_ag_center) && !is_enabled(f_ag_rotate))
-    calc_fit_gradients_impl<true, false>();
+    calc_fit_forces_impl<true, false>(accessor_main, accessor_fitting);
   if (!is_enabled(f_ag_center) && is_enabled(f_ag_rotate))
-    calc_fit_gradients_impl<false, true>();
+    calc_fit_forces_impl<false, true>(accessor_main, accessor_fitting);
   if (!is_enabled(f_ag_center) && !is_enabled(f_ag_rotate))
-    calc_fit_gradients_impl<false, false>();
+    calc_fit_forces_impl<false, false>(accessor_main, accessor_fitting);
 
   if (cvm::debug())
     cvm::log("Done calculating fit gradients.\n");
 }
 
 
-template <bool B_ag_center, bool B_ag_rotate>
+template <bool B_ag_center, bool B_ag_rotate,
-void cvm::atom_group::calc_fit_gradients_impl() {
+          typename main_force_accessor_T, typename fitting_force_accessor_T>
-  cvm::atom_group *group_for_fit = fitting_group ? fitting_group : this;
+void cvm::atom_group::calc_fit_forces_impl(
+  main_force_accessor_T accessor_main,
+  fitting_force_accessor_T accessor_fitting) const {
+  const cvm::atom_group *group_for_fit = fitting_group ? fitting_group : this;
   // the center of geometry contribution to the gradients
   cvm::rvector atom_grad;
   // the rotation matrix contribution to the gradients
@@ -1243,17 +1261,13 @@ void cvm::atom_group::calc_fit_gradients_impl() {
   cvm::vector1d<cvm::rvector> dq0_1(4);
   // loop 1: iterate over the current atom group
   for (size_t i = 0; i < size(); i++) {
-    cvm::atom_pos pos_orig;
     if (B_ag_center) {
-      atom_grad += atoms[i].grad;
+      atom_grad += accessor_main(i);
-      if (B_ag_rotate) pos_orig = rot_inv * (atoms[i].pos - ref_pos_cog);
-    } else {
-      if (B_ag_rotate) pos_orig = atoms[i].pos;
     }
     if (B_ag_rotate) {
       // calculate \partial(R(q) \vec{x}_i)/\partial q) \cdot \partial\xi/\partial\vec{x}_i
       cvm::quaternion const dxdq =
-        rot.q.position_derivative_inner(pos_orig, atoms[i].grad);
+        rot.q.position_derivative_inner(pos_unrotated[i], accessor_main(i));
       sum_dxdq[0] += dxdq[0];
       sum_dxdq[1] += dxdq[1];
       sum_dxdq[2] += dxdq[2];
@@ -1261,26 +1275,45 @@ void cvm::atom_group::calc_fit_gradients_impl() {
     }
   }
   if (B_ag_center) {
-    if (B_ag_rotate) atom_grad = rot.inverse().matrix() * atom_grad;
+    if (B_ag_rotate) atom_grad = rot_inv * atom_grad;
     atom_grad *= (-1.0)/(cvm::real(group_for_fit->size()));
   }
   // loop 2: iterate over the fitting group
   if (B_ag_rotate) rot_deriv->prepare_derivative(rotation_derivative_dldq::use_dq);
   for (size_t j = 0; j < group_for_fit->size(); j++) {
+    cvm::rvector fitting_force_grad{0, 0, 0};
     if (B_ag_center) {
-      group_for_fit->fit_gradients[j] = atom_grad;
+      fitting_force_grad += atom_grad;
     }
     if (B_ag_rotate) {
-      rot_deriv->calc_derivative_wrt_group1(j, nullptr, &dq0_1);
+      rot_deriv->calc_derivative_wrt_group1<false, true, false>(j, nullptr, &dq0_1);
       // multiply by {\partial q}/\partial\vec{x}_j and add it to the fit gradients
-      group_for_fit->fit_gradients[j] += sum_dxdq[0] * dq0_1[0] +
+      fitting_force_grad += sum_dxdq[0] * dq0_1[0] +
-                                          sum_dxdq[1] * dq0_1[1] +
+                            sum_dxdq[1] * dq0_1[1] +
-                                          sum_dxdq[2] * dq0_1[2] +
+                            sum_dxdq[2] * dq0_1[2] +
-                                          sum_dxdq[3] * dq0_1[3];
+                            sum_dxdq[3] * dq0_1[3];
     }
+    if (cvm::debug()) {
+      cvm::log(cvm::to_str(fitting_force_grad));
+    }
+    accessor_fitting(j, fitting_force_grad);
   }
 }
 
+template <typename main_force_accessor_T, typename fitting_force_accessor_T>
+void cvm::atom_group::calc_fit_forces(
+  main_force_accessor_T accessor_main,
+  fitting_force_accessor_T accessor_fitting) const {
+  if (is_enabled(f_ag_center) && is_enabled(f_ag_rotate))
+    calc_fit_forces_impl<true, true, main_force_accessor_T, fitting_force_accessor_T>(accessor_main, accessor_fitting);
+  if (is_enabled(f_ag_center) && !is_enabled(f_ag_rotate))
+    calc_fit_forces_impl<true, false, main_force_accessor_T, fitting_force_accessor_T>(accessor_main, accessor_fitting);
+  if (!is_enabled(f_ag_center) && is_enabled(f_ag_rotate))
+    calc_fit_forces_impl<false, true, main_force_accessor_T, fitting_force_accessor_T>(accessor_main, accessor_fitting);
+  if (!is_enabled(f_ag_center) && !is_enabled(f_ag_rotate))
+    calc_fit_forces_impl<false, false, main_force_accessor_T, fitting_force_accessor_T>(accessor_main, accessor_fitting);
+}
+
 
 std::vector<cvm::atom_pos> cvm::atom_group::positions() const
 {
@@ -1452,17 +1485,72 @@ void cvm::atom_group::apply_force(cvm::rvector const &force)
     return;
   }
 
-  if (is_enabled(f_ag_rotate)) {
+  auto ag_force = get_group_force_object();
+  for (size_t i = 0; i < size(); ++i) {
+    ag_force.add_atom_force(i, atoms[i].mass / total_mass * force);
+  }
+}
 
-    const auto rot_inv = rot.inverse().matrix();
+cvm::atom_group::group_force_object cvm::atom_group::get_group_force_object() {
-    for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
+  return cvm::atom_group::group_force_object(this);
-      ai->apply_force(rot_inv * ((ai->mass/total_mass) * force));
+}
+
+cvm::atom_group::group_force_object::group_force_object(cvm::atom_group* ag):
+m_ag(ag), m_group_for_fit(m_ag->fitting_group ? m_ag->fitting_group : m_ag),
+m_has_fitting_force(m_ag->is_enabled(f_ag_center) || m_ag->is_enabled(f_ag_rotate)) {
+  if (m_has_fitting_force) {
+    if (m_ag->group_forces.size() != m_ag->size()) {
+      m_ag->group_forces.assign(m_ag->size(), 0);
+    } else {
+      std::fill(m_ag->group_forces.begin(),
+                m_ag->group_forces.end(), 0);
     }
+  }
+}
 
+cvm::atom_group::group_force_object::~group_force_object() {
+  if (m_has_fitting_force) {
+    apply_force_with_fitting_group();
+  }
+}
+
+void cvm::atom_group::group_force_object::add_atom_force(size_t i, const cvm::rvector& force) {
+  if (m_has_fitting_force) {
+    m_ag->group_forces[i] += force;
   } else {
+    // Apply the force directly if we don't use fitting
+    (*m_ag)[i].apply_force(force);
+  }
+}
 
-    for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
+void cvm::atom_group::group_force_object::apply_force_with_fitting_group() {
-      ai->apply_force((ai->mass/total_mass) * force);
+  const cvm::rmatrix rot_inv = m_ag->rot.inverse().matrix();
+  if (cvm::debug()) {
+    cvm::log("Applying force on main group " + m_ag->name + ":\n");
+  }
+  for (size_t ia = 0; ia < m_ag->size(); ++ia) {
+    const cvm::rvector f_ia = rot_inv * m_ag->group_forces[ia];
+    (*m_ag)[ia].apply_force(f_ia);
+    if (cvm::debug()) {
+      cvm::log(cvm::to_str(f_ia));
+    }
+  }
+  // Gradients are only available with scalar components, so for a scalar component,
+  // if f_ag_fit_gradients is disabled, then the forces on the fitting group is not
+  // computed. For a vector component, we can only know the forces on the fitting
+  // group, but checking this flag can mimic results that the users expect (if
+  // "enableFitGradients no" then there is no force on the fitting group).
+  if (!m_ag->b_dummy && m_ag->is_enabled(f_ag_fit_gradients)) {
+    auto accessor_main = [this](size_t i){return m_ag->group_forces[i];};
+    auto accessor_fitting = [this](size_t j, const cvm::rvector& fitting_force){
+      (*(m_group_for_fit))[j].apply_force(fitting_force);
+    };
+    if (cvm::debug()) {
+      cvm::log("Applying force on the fitting group of main group" + m_ag->name + ":\n");
+    }
+    m_ag->calc_fit_forces(accessor_main, accessor_fitting);
+    if (cvm::debug()) {
+      cvm::log("Done applying force on the fitting group of main group" + m_ag->name + ":\n");
     }
   }
 }

lib/colvars/colvaratoms.h

@@ -194,7 +194,7 @@ public:
 
   int add_atom_numbers(std::string const &numbers_conf);
   int add_atoms_of_group(atom_group const * ag);
-  int add_index_group(std::string const &index_group_name);
+  int add_index_group(std::string const &index_group_name, bool silent = false);
   int add_atom_numbers_range(std::string const &range_conf);
   int add_atom_name_residue_range(std::string const &psf_segid,
                                   std::string const &range_conf);
@@ -257,8 +257,63 @@ protected:
   /// \brief Index in the colvarproxy arrays (if the group is scalable)
   int index;
 
+  /// \brief The temporary forces acting on the main group atoms.
+  ///        Currently this is only used for calculating the fitting group forces for
+  ///        non-scalar components.
+  std::vector<cvm::rvector> group_forces;
+
 public:
 
+  /*! @class group_force_object
+   *  @brief A helper class for applying forces on an atom group in a way that
+   *         is aware of the fitting group. NOTE: you are encouraged to use
+   *         get_group_force_object() to get an instance of group_force_object
+   *         instead of constructing directly.
+   */
+  class group_force_object {
+  public:
+    /*! @brief Constructor of group_force_object
+     *  @param ag The pointer to the atom group that forces will be applied on.
+     */
+    group_force_object(cvm::atom_group* ag);
+    /*! @brief Destructor of group_force_object
+     */
+    ~group_force_object();
+    /*! @brief Apply force to atom i
+     *  @param i The i-th of atom in the atom group.
+     *  @param force The force being added to atom i.
+     *
+     * The function can be used as follows,
+     * @code
+     *       // In your colvar::cvc::apply_force() loop of a component:
+     *       auto ag_force = atoms->get_group_force_object();
+     *       for (ia = 0; ia < atoms->size(); ia++) {
+     *         const cvm::rvector f = compute_force_on_atom_ia();
+     *         ag_force.add_atom_force(ia, f);
+     *       }
+     * @endcode
+     * There are actually two scenarios under the hood:
+     * (i) If the atom group does not have a fitting group, then the force is
+     *     added to atom i directly;
+     * (ii) If the atom group has a fitting group, the force on atom i will just
+     *      be temporary stashed into ag->group_forces. At the end of the loop
+     *      of apply_force(), the destructor ~group_force_object() will be called,
+     *      which then call apply_force_with_fitting_group(). The forces on the
+     *      main group will be rotated back by multiplying ag->group_forces with
+     *      the inverse rotation. The forces on the fitting group (if
+     *      enableFitGradients is on) will be calculated by calling
+     *      calc_fit_forces.
+     */
+    void add_atom_force(size_t i, const cvm::rvector& force);
+  private:
+    cvm::atom_group* m_ag;
+    cvm::atom_group* m_group_for_fit;
+    bool m_has_fitting_force;
+    void apply_force_with_fitting_group();
+  };
+
+  group_force_object get_group_force_object();
+
   inline cvm::atom & operator [] (size_t const i)
   {
     return atoms[i];
@@ -423,6 +478,9 @@ private:
   /// \brief Center of geometry before any fitting
   cvm::atom_pos cog_orig;
 
+  /// \brief Unrotated atom positions for fit gradients
+  std::vector<cvm::atom_pos> pos_unrotated;
+
 public:
 
   /// \brief Return the center of geometry of the atomic positions
@@ -497,15 +555,60 @@ public:
   /// \brief Calculate the derivatives of the fitting transformation
   void calc_fit_gradients();
 
-/*! @brief  Actual implementation of `calc_fit_gradients`. The template is
+/*! @brief  Actual implementation of `calc_fit_gradients` and
+ *          `calc_fit_forces`. The template is
  *          used to avoid branching inside the loops in case that the CPU
  *          branch prediction is broken (or further migration to GPU code).
  *  @tparam B_ag_center Centered the reference to origin? This should follow
  *          the value of `is_enabled(f_ag_center)`.
  *  @tparam B_ag_rotate Calculate the optimal rotation? This should follow
  *          the value of `is_enabled(f_ag_rotate)`.
+ *  @tparam main_force_accessor_T The type of accessor of the main
+ *          group forces or gradients acting on the rotated frame.
+ *  @tparam fitting_force_accessor_T The type of accessor of the fitting group
+ *          forces or gradients.
+ *  @param accessor_main The accessor of the main group forces or gradients.
+ *         accessor_main(i) should return the i-th force or gradient of the
+ *         rotated main group.
+ *  @param accessor_fitting The accessor of the fitting group forces or gradients.
+ *         accessor_fitting(j, v) should store/apply the j-th atom gradient or
+ *         force in the fitting group.
+ *
+ *  This function is used to (i) project the gradients of CV with respect to
+ *  rotated main group atoms to fitting group atoms, or (ii) project the forces
+ *  on rotated main group atoms to fitting group atoms, by the following two steps
+ *  (using the goal (ii) for example):
+ *  (1) Loop over the positions of main group atoms and call cvm::quaternion::position_derivative_inner
+ *      to project the forces on rotated main group atoms to the forces on quaternion.
+ *  (2) Loop over the positions of fitting group atoms, compute the gradients of
+ *      \f$\mathbf{q}\f$ with respect to the position of each atom, and then multiply
+ *      that with the force on \f$\mathbf{q}\f$ (chain rule).
  */
-  template <bool B_ag_center, bool B_ag_rotate> void calc_fit_gradients_impl();
+  template <bool B_ag_center, bool B_ag_rotate,
+            typename main_force_accessor_T, typename fitting_force_accessor_T>
+  void calc_fit_forces_impl(
+    main_force_accessor_T accessor_main,
+    fitting_force_accessor_T accessor_fitting) const;
+
+/*! @brief  Calculate or apply the fitting group forces from the main group forces.
+ *  @tparam main_force_accessor_T The type of accessor of the main
+ *          group forces or gradients.
+ *  @tparam fitting_force_accessor_T The type of accessor of the fitting group
+ *          forces or gradients.
+ *  @param accessor_main The accessor of the main group forces or gradients.
+ *         accessor_main(i) should return the i-th force or gradient of the
+ *         main group.
+ *  @param accessor_fitting The accessor of the fitting group forces or gradients.
+ *         accessor_fitting(j, v) should store/apply the j-th atom gradient or
+ *         force in the fitting group.
+ *
+ *  This function just dispatches the parameters to calc_fit_forces_impl that really
+ *  performs the calculations.
+ */
+  template <typename main_force_accessor_T, typename fitting_force_accessor_T>
+  void calc_fit_forces(
+    main_force_accessor_T accessor_main,
+    fitting_force_accessor_T accessor_fitting) const;
 
   /// \brief Derivatives of the fitting transformation
   std::vector<cvm::atom_pos> fit_gradients;

lib/colvars/colvarbias.cpp

@@ -93,6 +93,8 @@ int colvarbias::init(std::string const &conf)
     cvm::log("Reinitializing bias \""+name+"\".\n");
   }
 
+  feature_states[f_cvb_step_zero_data].available = true;
+
   colvar_values.resize(num_variables());
   for (i = 0; i < num_variables(); i++) {
     colvar_values[i].type(colvars[i]->value().type());
@@ -157,7 +159,7 @@ int colvarbias::init_dependencies() {
     init_feature(f_cvb_step_zero_data, "step_zero_data", f_type_user);
 
     init_feature(f_cvb_apply_force, "apply_force", f_type_user);
-    require_feature_children(f_cvb_apply_force, f_cv_gradient);
+    require_feature_children(f_cvb_apply_force, f_cv_apply_force);
 
     init_feature(f_cvb_bypass_ext_lagrangian, "bypass_extended_Lagrangian_coordinates", f_type_user);
 
@@ -199,6 +201,8 @@ int colvarbias::init_dependencies() {
 
     init_feature(f_cvb_extended, "Bias on extended-Lagrangian variables", f_type_static);
 
+    init_feature(f_cvb_smp, "smp_computation", f_type_user);
+
     // check that everything is initialized
     for (i = 0; i < colvardeps::f_cvb_ntot; i++) {
       if (is_not_set(i)) {
@@ -221,8 +225,9 @@ int colvarbias::init_dependencies() {
   // The feature f_cvb_bypass_ext_lagrangian is only implemented by some derived classes
   // (initially, harmonicWalls)
   feature_states[f_cvb_bypass_ext_lagrangian].available = false;
-  // disabled by default; can be changed by derived classes that implement it
+
-  feature_states[f_cvb_bypass_ext_lagrangian].enabled = false;
+  // Most biases cannot currently be processed in parallel over threads
+  feature_states[f_cvb_smp].available = false;
 
   return COLVARS_OK;
 }
@@ -704,7 +709,7 @@ int colvarbias::read_state_string(char const *buffer)
 
 
 std::ostream &colvarbias::write_state_data_key(std::ostream &os, std::string const &key,
-                                               bool header)
+                                               bool header) const
 {
   os << (header ? "\n" : "") << key << (header ? "\n" : " ");
   return os;
@@ -712,7 +717,7 @@ std::ostream &colvarbias::write_state_data_key(std::ostream &os, std::string con
 
 
 cvm::memory_stream &colvarbias::write_state_data_key(cvm::memory_stream &os, std::string const &key,
-                                                     bool /* header */)
+                                                     bool /* header */) const
 {
   os << std::string(key);
   return os;
@@ -792,6 +797,8 @@ int colvarbias_ti::init(std::string const &conf)
 {
   int error_code = COLVARS_OK;
 
+  key_lookup(conf, "grid", &grid_conf);
+
   get_keyval_feature(this, conf, "writeTISamples",
                      f_cvb_write_ti_samples,
                      is_enabled(f_cvb_write_ti_samples));
@@ -800,18 +807,16 @@ int colvarbias_ti::init(std::string const &conf)
                      f_cvb_write_ti_pmf,
                      is_enabled(f_cvb_write_ti_pmf));
 
+  if (is_enabled(f_cvb_write_ti_pmf)) {
+    enable(f_cvb_write_ti_samples);
+  }
+
   if ((num_variables() > 1) && is_enabled(f_cvb_write_ti_pmf)) {
     return cvm::error("Error: only 1-dimensional PMFs can be written "
                       "on the fly.\n"
                       "Consider using writeTISamples instead and "
                       "post-processing the sampled free-energy gradients.\n",
                       COLVARS_NOT_IMPLEMENTED);
-  } else {
-    error_code |= init_grids();
-  }
-
-  if (is_enabled(f_cvb_write_ti_pmf)) {
-    enable(f_cvb_write_ti_samples);
   }
 
   if (is_enabled(f_cvb_calc_ti_samples)) {
@@ -831,6 +836,8 @@ int colvarbias_ti::init(std::string const &conf)
     }
   }
 
+  error_code |= colvarbias_ti::init_grids();
+
   if (is_enabled(f_cvb_write_ti_pmf) || is_enabled(f_cvb_write_ti_samples)) {
     cvm::main()->cite_feature("Internal-forces free energy estimator");
   }
@@ -844,16 +851,15 @@ int colvarbias_ti::init_grids()
   if (is_enabled(f_cvb_calc_ti_samples)) {
     if (!ti_avg_forces) {
       ti_bin.resize(num_variables());
+      ti_bin.assign(ti_bin.size(), -1);
       ti_system_forces.resize(num_variables());
       for (size_t icv = 0; icv < num_variables(); icv++) {
         ti_system_forces[icv].type(variables(icv)->value());
         ti_system_forces[icv].is_derivative();
         ti_system_forces[icv].reset();
       }
-      ti_avg_forces.reset(new colvar_grid_gradient(colvars));
+      ti_count.reset(new colvar_grid_count(colvars, grid_conf));
-      ti_count.reset(new colvar_grid_count(colvars));
+      ti_avg_forces.reset(new colvar_grid_gradient(colvars, ti_count));
-      ti_avg_forces->samples = ti_count;
-      ti_count->has_parent_data = true;
     }
   }
 
@@ -884,8 +890,12 @@ int colvarbias_ti::update_system_forces(std::vector<colvarvalue> const
 
   size_t i;
 
-  if (proxy->total_forces_same_step()) {
+  if (cvm::debug()) {
-    for (i = 0; i < num_variables(); i++) {
+    cvm::log("TI bin for bias \"" + name + "\" = " + cvm::to_str(ti_bin) + ".\n");
+  }
+
+  for (i = 0; i < num_variables(); i++) {
+    if (variables(i)->is_enabled(f_cv_total_force_current_step)) {
       ti_bin[i] = ti_avg_forces->current_bin_scalar(i);
     }
   }
@@ -894,8 +904,10 @@ int colvarbias_ti::update_system_forces(std::vector<colvarvalue> const
   if ((cvm::step_relative() > 0) || proxy->total_forces_same_step()) {
     if (ti_avg_forces->index_ok(ti_bin)) {
       for (i = 0; i < num_variables(); i++) {
-        if (variables(i)->is_enabled(f_cv_subtract_applied_force)) {
+        if (variables(i)->is_enabled(f_cv_subtract_applied_force) ||
+          (cvm::proxy->total_forces_same_step() && !variables(i)->is_enabled(f_cv_external))) {
           // this colvar is already subtracting all applied forces
+          // or the "total force" is really a system force at current step
           ti_system_forces[i] = variables(i)->total_force();
         } else {
           ti_system_forces[i] = variables(i)->total_force() -
@@ -904,14 +916,17 @@ int colvarbias_ti::update_system_forces(std::vector<colvarvalue> const
         }
       }
       if (cvm::step_relative() > 0 || is_enabled(f_cvb_step_zero_data)) {
+        if (cvm::debug()) {
+          cvm::log("Accumulating TI forces for bias \"" + name + "\".\n");
+        }
         ti_avg_forces->acc_value(ti_bin, ti_system_forces);
       }
     }
   }
 
-  if (!proxy->total_forces_same_step()) {
+  for (i = 0; i < num_variables(); i++) {
-    // Set the index for use in the next iteration, when total forces come in
+    if (!variables(i)->is_enabled(f_cv_total_force_current_step)) {
-    for (i = 0; i < num_variables(); i++) {
+      // Set the index for use in the next iteration, when total forces come in
       ti_bin[i] = ti_avg_forces->current_bin_scalar(i);
     }
   }

lib/colvars/colvarbias.h

@@ -174,14 +174,14 @@ public:
   /// \param[in,out] os Output stream
   /// \param[in] key  Keyword labeling the header block
   /// \param[in] header  Whether this is the header of a multi-line segment vs a single line
-  std::ostream &write_state_data_key(std::ostream &os, std::string const &key, bool header = true);
+  std::ostream &write_state_data_key(std::ostream &os, std::string const &key, bool header = true) const;
 
   /// Write a keyword header for a data sequence to an unformatted stream
   /// \param[in,out] os Output stream
   /// \param[in] key  Keyword labeling the header block
   /// \param[in] header  Ignored
   cvm::memory_stream &write_state_data_key(cvm::memory_stream &os, std::string const &key,
-                                           bool header = true);
+                                           bool header = true) const;
 
 private:
 
@@ -358,6 +358,9 @@ protected:
   /// \brief Forces exerted from the system to the associated variables
   std::vector<colvarvalue> ti_system_forces;
 
+  /// Grid configuration parameters (also used by grids in derived classes)
+  std::string grid_conf;
+
   /// Averaged system forces
   std::shared_ptr<colvar_grid_gradient> ti_avg_forces;
 

lib/colvars/colvarbias_abf.cpp

@@ -87,24 +87,25 @@ int colvarbias_abf::init(std::string const &conf)
   get_keyval(conf, "shared", shared_on, false);
   if (shared_on) {
     cvm::main()->cite_feature("Multiple-walker ABF implementation");
-    if ((proxy->replica_enabled() != COLVARS_OK) ||
+    cvm::main()->cite_feature("Updated multiple-walker ABF implementation");
-        (proxy->num_replicas() <= 1)) {
+
-      return cvm::error("Error: shared ABF requires more than one replica.",
+
-                        COLVARS_INPUT_ERROR);
+    // Cannot check this here because the replica communicator is obtained later
-    }
+    // in Gromacs
-    cvm::log("shared ABF will be applied among "+
+
-             cvm::to_str(proxy->num_replicas()) + " replicas.\n");
+    // if ((proxy->check_replicas_enabled() != COLVARS_OK) ||
+    //     (proxy->num_replicas() <= 1)) {
+    //   return cvm::error("Error: shared ABF requires more than one replica.",
+    //                     COLVARS_INPUT_ERROR);
+    // }
+    // cvm::log("shared ABF will be applied among "+
+    //          cvm::to_str(proxy->num_replicas()) + " replicas.\n");
 
     // If shared_freq is not set, we default to output_freq
     get_keyval(conf, "sharedFreq", shared_freq, output_freq);
     if ( shared_freq && output_freq % shared_freq ) {
       return cvm::error("Error: outputFreq must be a multiple of sharedFreq.\n");
     }
-
-    // Allocate these at init time if possible
-    local_samples.reset(new colvar_grid_count(colvars));
-    local_gradients.reset(new colvar_grid_gradient(colvars, local_samples));
-    local_pmf.reset(new integrate_potential(colvars, local_gradients));
   }
 
   // ************* checking the associated colvars *******************
@@ -124,10 +125,17 @@ int colvarbias_abf::init(std::string const &conf)
       colvars[i]->enable(f_cv_hide_Jacobian);
     }
 
-    // If any colvar is extended-system, we need to collect the extended
+    // If any colvar is extended-system (restrained, not driven external param), we are running eABF
-    // system gradient
+    if (colvars[i]->is_enabled(f_cv_extended_Lagrangian)
-    if (colvars[i]->is_enabled(f_cv_extended_Lagrangian))
+        && !colvars[i]->is_enabled(f_cv_external)) {
       enable(f_cvb_extended);
+    }
+
+    if (!colvars[i]->is_enabled(f_cv_total_force_current_step)) {
+      // If any colvar does not have current-step total force, then
+      // we can't do step 0 data
+      provide(f_cvb_step_zero_data, false);
+    }
 
     // Cannot mix and match coarse time steps with ABF because it gives
     // wrong total force averages - total force needs to be averaged over
@@ -181,12 +189,23 @@ int colvarbias_abf::init(std::string const &conf)
     cvm::log("Allocating count and free energy gradient grids.\n");
   }
 
-  samples.reset(new colvar_grid_count(colvars));
+  {
-  gradients.reset(new colvar_grid_gradient(colvars, samples));
+    /// Optional custom configuration string for grid parameters
+    std::string grid_conf;
+    key_lookup(conf, "grid", &grid_conf);
+
+    samples.reset(new colvar_grid_count(colvars, grid_conf));
+  }
+  gradients.reset(new colvar_grid_gradient(colvars, samples)); // Also use samples as template for sizes
 
   gradients->full_samples = full_samples;
   gradients->min_samples = min_samples;
 
+  if (shared_on) {
+    local_samples.reset(new colvar_grid_count(colvars, samples));
+    local_gradients.reset(new colvar_grid_gradient(colvars, local_samples));
+  }
+
   // Data for eABF z-based estimator
   if (is_enabled(f_cvb_extended)) {
     get_keyval(conf, "CZARestimator", b_CZAR_estimator, true);
@@ -198,11 +217,11 @@ int colvarbias_abf::init(std::string const &conf)
                colvarparse::parse_silent);
 
     z_bin.assign(num_variables(), 0);
-    z_samples.reset(new colvar_grid_count(colvars));
+    z_samples.reset(new colvar_grid_count(colvars, samples));
     z_samples->request_actual_value();
     z_gradients.reset(new colvar_grid_gradient(colvars, z_samples));
     z_gradients->request_actual_value();
-    czar_gradients.reset(new colvar_grid_gradient(colvars));
+    czar_gradients.reset(new colvar_grid_gradient(colvars, nullptr, samples));
   }
 
   get_keyval(conf, "integrate", b_integrate, num_variables() <= 3); // Integrate for output if d<=3
@@ -216,6 +235,9 @@ int colvarbias_abf::init(std::string const &conf)
     if (b_CZAR_estimator) {
       czar_pmf.reset(new integrate_potential(colvars, czar_gradients));
     }
+    if (shared_on) {
+      local_pmf.reset(new integrate_potential(colvars, local_gradients));
+    }
     // Parameters for integrating initial (and final) gradient data
     get_keyval(conf, "integrateMaxIterations", integrate_iterations, 10000, colvarparse::parse_silent);
     get_keyval(conf, "integrateTol", integrate_tol, 1e-6, colvarparse::parse_silent);
@@ -228,9 +250,9 @@ int colvarbias_abf::init(std::string const &conf)
   if (b_CZAR_estimator && shared_on && cvm::main()->proxy->replica_index() == 0) {
     // The pointers below are used for outputting CZAR data
     // Allocate grids for collected global data, on replica 0 only
-    global_z_samples.reset(new colvar_grid_count(colvars));
+    global_z_samples.reset(new colvar_grid_count(colvars, samples));
     global_z_gradients.reset(new colvar_grid_gradient(colvars, global_z_samples));
-    global_czar_gradients.reset(new colvar_grid_gradient(colvars));
+    global_czar_gradients.reset(new colvar_grid_gradient(colvars, nullptr, samples));
     global_czar_pmf.reset(new integrate_potential(colvars, global_czar_gradients));
   } else {
     // otherwise they are just aliases for the local CZAR grids
@@ -244,10 +266,10 @@ int colvarbias_abf::init(std::string const &conf)
   // This used to be only if "shared" was defined,
   // but now we allow calling share externally (e.g. from Tcl).
   if (b_CZAR_estimator) {
-    z_samples_in.reset(new colvar_grid_count(colvars));
+    z_samples_in.reset(new colvar_grid_count(colvars, samples));
     z_gradients_in.reset(new colvar_grid_gradient(colvars, z_samples_in));
   }
-  last_samples.reset(new colvar_grid_count(colvars));
+  last_samples.reset(new colvar_grid_count(colvars, samples));
   last_gradients.reset(new colvar_grid_gradient(colvars, last_samples));
   // Any data collected after now is new for shared ABF purposes
   shared_last_step = cvm::step_absolute();
@@ -315,27 +337,36 @@ int colvarbias_abf::update()
   size_t i;
   for (i = 0; i < num_variables(); i++) {
     bin[i] = samples->current_bin_scalar(i);
+    if (colvars[i]->is_enabled(f_cv_total_force_current_step)) {
+      force_bin[i] = bin[i];
+    }
   }
 
-
   // ***********************************************************
   // ******  ABF Part I: update the FE gradient estimate  ******
   // ***********************************************************
 
 
-  if (cvm::proxy->total_forces_same_step()) {
+  // Share data first, so that 2d/3d PMF is refreshed using new data for mw-pABF.
-    // e.g. in LAMMPS, total forces are current
+  // shared_on can be true with shared_freq 0 if we are sharing via script
-    force_bin = bin;
+  if (shared_on && shared_freq &&
+      shared_last_step >= 0 &&                    // we have already collected some data
+      cvm::step_absolute() > shared_last_step &&  // time has passed since the last sharing timestep
+                                                  // (avoid re-sharing at last and first ts of successive run statements)
+      cvm::step_absolute() % shared_freq == 0) {
+    // Share gradients and samples for shared ABF.
+    replica_share();
   }
 
   if (can_accumulate_data() && is_enabled(f_cvb_history_dependent)) {
 
     if (cvm::step_relative() > 0 || cvm::proxy->total_forces_same_step()) {
+      // Note: this will skip step 0 data when available in some cases (extended system),
+      // but not doing so would make the code more complex
       if (samples->index_ok(force_bin)) {
         // Only if requested and within bounds of the grid...
 
-        // get total forces (lagging by 1 timestep) from colvars
+        // get total force and subtract previous ABF force if necessary
-        // and subtract previous ABF force if necessary
         update_system_force();
 
         gradients->acc_force(force_bin, system_force);
@@ -368,21 +399,11 @@ int colvarbias_abf::update()
     }
   }
 
-  if (!(cvm::proxy->total_forces_same_step())) {
+  // In some cases, total forces are stored for next timestep
-    // e.g. in NAMD, total forces will be available for next timestep
+  // hence we store the current colvar bin - this is overwritten on a per-colvar basis
-    // hence we store the current colvar bin
+  // at the top of update()
-    force_bin = bin;
+  force_bin = bin;
-  }
 
-  // Share data after force sample is collected for this time step
-  // shared_on can be true with shared_freq 0 if we are sharing via script
-  if (shared_on && shared_freq &&
-      cvm::step_absolute() > shared_last_step &&  // time has passed since the last sharing timestep
-                                                  // (avoid re-sharing at last and first ts of successive run statements)
-      cvm::step_absolute() % shared_freq == 0) {
-    // Share gradients and samples for shared ABF.
-    replica_share();
-  }
 
   // ******************************************************************
   // ******  ABF Part II: calculate and apply the biasing force  ******
@@ -452,10 +473,13 @@ int colvarbias_abf::update_system_force()
   // System force from atomic forces (or extended Lagrangian if applicable)
 
   for (i = 0; i < num_variables(); i++) {
-    if (colvars[i]->is_enabled(f_cv_subtract_applied_force)) {
+    if (colvars[i]->is_enabled(f_cv_subtract_applied_force)
+      || colvars[i]->is_enabled(f_cv_total_force_current_step)) {
       // this colvar is already subtracting the ABF force
+      // or the "total force" is from current step and cannot possibly contain Colvars biases
       system_force[i] = colvars[i]->total_force().real_value;
     } else {
+      // Subtract previous step's bias force from previous step's total force
       system_force[i] = colvars[i]->total_force().real_value
         - colvar_forces[i].real_value;
     }
@@ -525,7 +549,7 @@ int colvarbias_abf::replica_share() {
 
   colvarproxy *proxy = cvm::main()->proxy;
 
-  if (proxy->replica_enabled() != COLVARS_OK) {
+  if (proxy->check_replicas_enabled() != COLVARS_OK) {
     cvm::error("Error: shared ABF: No replicas.\n");
     return COLVARS_ERROR;
   }
@@ -542,7 +566,7 @@ int colvarbias_abf::replica_share() {
   if (!local_samples) {
     // We arrive here if sharing has just been enabled by a script
     // in which case local arrays have not been initialized yet
-    local_samples.reset(new colvar_grid_count(colvars));
+    local_samples.reset(new colvar_grid_count(colvars, samples));
     local_gradients.reset(new colvar_grid_gradient(colvars, local_samples));
     local_pmf.reset(new integrate_potential(colvars, local_gradients));
   }
@@ -662,9 +686,9 @@ int colvarbias_abf::replica_share_CZAR() {
       // We arrive here if sharing has just been enabled by a script
       // Allocate grids for collective data, on replica 0 only
       // overriding CZAR grids that are equal to local ones by default
-      global_z_samples.reset(new colvar_grid_count(colvars));
+      global_z_samples.reset(new colvar_grid_count(colvars, samples));
       global_z_gradients.reset(new colvar_grid_gradient(colvars, global_z_samples));
-      global_czar_gradients.reset(new colvar_grid_gradient(colvars));
+      global_czar_gradients.reset(new colvar_grid_gradient(colvars, nullptr, samples));
       global_czar_pmf.reset(new integrate_potential(colvars, global_czar_gradients));
     }
 

lib/colvars/colvarbias_histogram.cpp

@@ -98,10 +98,10 @@ int colvarbias_histogram::init(std::string const &conf)
   }
 
   {
-    std::string grid_conf;
+    if (key_lookup(conf, "histogramGrid", &grid_conf) ||
-    if (key_lookup(conf, "histogramGrid", &grid_conf)) {
+        key_lookup(conf, "grid", &grid_conf)) {
       grid->parse_params(grid_conf);
-      grid->check_keywords(grid_conf, "histogramGrid");
+      grid->check_keywords(grid_conf, "grid");
     }
   }
 

lib/colvars/colvarbias_histogram.h

@@ -38,6 +38,7 @@ protected:
 
   /// n-dim histogram
   colvar_grid_scalar *grid;
+  std::string grid_conf;
   std::vector<int> bin;
   std::string out_name, out_name_dx;
 

lib/colvars/colvarbias_histogram_reweight_amd.cpp

@@ -11,43 +11,9 @@
 #include "colvarproxy.h"
 #include "colvars_memstream.h"
 
-colvarbias_reweightaMD::colvarbias_reweightaMD(char const *key)
+colvarbias_reweightaMD::colvarbias_reweightaMD(char const *key) : colvarbias_histogram(key) {}
-  : colvarbias_histogram(key), grid_count(NULL), grid_dV(NULL),
-    grid_dV_square(NULL), pmf_grid_exp_avg(NULL), pmf_grid_cumulant(NULL),
-    grad_grid_exp_avg(NULL), grad_grid_cumulant(NULL)
-{
-}
 
-colvarbias_reweightaMD::~colvarbias_reweightaMD() {
+colvarbias_reweightaMD::~colvarbias_reweightaMD() {}
-  if (grid_dV) {
-    delete grid_dV;
-    grid_dV = NULL;
-  }
-  if (grid_dV_square) {
-    delete grid_dV_square;
-    grid_dV_square = NULL;
-  }
-  if (grid_count) {
-    delete grid_count;
-    grid_count = NULL;
-  }
-  if (pmf_grid_exp_avg) {
-    delete pmf_grid_exp_avg;
-    pmf_grid_exp_avg = NULL;
-  }
-  if (pmf_grid_cumulant) {
-    delete pmf_grid_cumulant;
-    pmf_grid_cumulant = NULL;
-  }
-  if (grad_grid_exp_avg) {
-    delete grad_grid_exp_avg;
-    grad_grid_exp_avg = NULL;
-  }
-  if (grad_grid_cumulant) {
-    delete grad_grid_cumulant;
-    grad_grid_cumulant = NULL;
-  }
-}
 
 int colvarbias_reweightaMD::init(std::string const &conf) {
   if (cvm::proxy->accelMD_enabled() == false) {
@@ -60,21 +26,21 @@ int colvarbias_reweightaMD::init(std::string const &conf) {
   get_keyval(conf, "WritePMFGradients", b_write_gradients, true);
   get_keyval(conf, "historyFreq", history_freq, 0);
   b_history_files = (history_freq > 0);
-  grid_count = new colvar_grid_scalar(colvars);
+  grid_count.reset(new colvar_grid_scalar(colvars, nullptr, false, grid_conf));
   grid_count->request_actual_value();
   grid->request_actual_value();
-  pmf_grid_exp_avg = new colvar_grid_scalar(colvars);
+  pmf_grid_exp_avg.reset(new colvar_grid_scalar(colvars, grid_count));
   if (b_write_gradients) {
-    grad_grid_exp_avg = new colvar_grid_gradient(colvars);
+    grad_grid_exp_avg.reset(new colvar_grid_gradient(colvars, nullptr, grid_count));
   }
   if (b_use_cumulant_expansion) {
-    grid_dV = new colvar_grid_scalar(colvars);
+    grid_dV.reset(new colvar_grid_scalar(colvars, grid_count));
-    grid_dV_square = new colvar_grid_scalar(colvars);
+    grid_dV_square.reset(new colvar_grid_scalar(colvars, grid_count));
-    pmf_grid_cumulant = new colvar_grid_scalar(colvars);
+    pmf_grid_cumulant.reset(new colvar_grid_scalar(colvars, grid_count));
     grid_dV->request_actual_value();
     grid_dV_square->request_actual_value();
     if (b_write_gradients) {
-      grad_grid_cumulant = new colvar_grid_gradient(colvars);
+      grad_grid_cumulant.reset(new colvar_grid_gradient(colvars, nullptr, grid_count));
     }
   }
   previous_bin.assign(num_variables(), -1);
@@ -193,7 +159,7 @@ int colvarbias_reweightaMD::write_exponential_reweighted_pmf(
       pmf_grid_exp_avg->set_value(i, tmp / count);
     }
   }
-  hist_to_pmf(pmf_grid_exp_avg, grid_count);
+  hist_to_pmf(pmf_grid_exp_avg.get(), grid_count.get());
   pmf_grid_exp_avg->write_multicol(pmf_grid_os);
   if (!keep_open) {
     cvm::proxy->close_output_stream(output_pmf);
@@ -231,9 +197,9 @@ int colvarbias_reweightaMD::write_cumulant_expansion_pmf(
   if (!pmf_grid_cumulant_os) {
     return COLVARS_FILE_ERROR;
   }
-  compute_cumulant_expansion_factor(grid_dV, grid_dV_square,
+  compute_cumulant_expansion_factor(grid_dV.get(), grid_dV_square.get(),
-                                    grid_count, pmf_grid_cumulant);
+                                    grid_count.get(), pmf_grid_cumulant.get());
-  hist_to_pmf(pmf_grid_cumulant, grid_count);
+  hist_to_pmf(pmf_grid_cumulant.get(), grid_count.get());
   pmf_grid_cumulant->write_multicol(pmf_grid_cumulant_os);
   if (!keep_open) {
     cvm::proxy->close_output_stream(output_pmf);

lib/colvars/colvarbias_histogram_reweight_amd.h

@@ -68,9 +68,9 @@ protected:
 
   /// Use cumulant expansion to second order?
   bool b_use_cumulant_expansion;
-  colvar_grid_scalar* grid_count;
+  std::shared_ptr<colvar_grid_scalar> grid_count;
-  colvar_grid_scalar* grid_dV;
+  std::unique_ptr<colvar_grid_scalar> grid_dV;
-  colvar_grid_scalar* grid_dV_square;
+  std::unique_ptr<colvar_grid_scalar> grid_dV_square;
 
   /// Number of timesteps between recording data in history files (if non-zero)
   size_t history_freq;
@@ -90,10 +90,10 @@ protected:
 
 private:
   /// temporary grids for evaluating PMFs
-  colvar_grid_scalar  *pmf_grid_exp_avg;
+  std::unique_ptr<colvar_grid_scalar> pmf_grid_exp_avg;
-  colvar_grid_scalar  *pmf_grid_cumulant;
+  std::unique_ptr<colvar_grid_scalar> pmf_grid_cumulant;
-  colvar_grid_gradient *grad_grid_exp_avg;
+  std::unique_ptr<colvar_grid_gradient> grad_grid_exp_avg;
-  colvar_grid_gradient *grad_grid_cumulant;
+  std::unique_ptr<colvar_grid_gradient> grad_grid_cumulant;
 };
 
 #endif // COLVARBIAS_HISTOGRAM_REWEIGHT_AMD

lib/colvars/colvarbias_meta.cpp

@@ -11,27 +11,10 @@
 #include <iomanip>
 #include <algorithm>
 
-// Define function to get the absolute path of a replica file
-#if defined(_WIN32) && !defined(__CYGWIN__)
-#include <direct.h>
-#define GETCWD(BUF, SIZE) ::_getcwd(BUF, SIZE)
-#define PATHSEP "\\"
-#else
-#include <unistd.h>
-#define GETCWD(BUF, SIZE) ::getcwd(BUF, SIZE)
-#define PATHSEP "/"
-#endif
-
-#ifdef __cpp_lib_filesystem
-// When std::filesystem is available, use it
-#include <filesystem>
-#undef GETCWD
-#define GETCWD(BUF, SIZE) (std::filesystem::current_path().string().c_str())
-#endif
-
 #include "colvarmodule.h"
 #include "colvarproxy.h"
 #include "colvar.h"
+#include "colvargrid.h"
 #include "colvarbias_meta.h"
 #include "colvars_memstream.h"
 
@@ -49,8 +32,6 @@ colvarbias_meta::colvarbias_meta(char const *key)
   use_grids = true;
   grids_freq = 0;
   rebin_grids = false;
-  hills_energy = NULL;
-  hills_energy_gradients = NULL;
 
   dump_fes = true;
   keep_hills = false;
@@ -161,9 +142,9 @@ int colvarbias_meta::init(std::string const &conf)
     get_keyval(conf, "keepHills", keep_hills, keep_hills);
     get_keyval(conf, "keepFreeEnergyFiles", dump_fes_save, dump_fes_save);
 
-    if (hills_energy == NULL) {
+    if (!hills_energy) {
-      hills_energy           = new colvar_grid_scalar(colvars);
+      hills_energy.reset(new colvar_grid_scalar(colvars, nullptr, false, grid_conf));
-      hills_energy_gradients = new colvar_grid_gradient(colvars);
+      hills_energy_gradients.reset(new colvar_grid_gradient(colvars, nullptr, hills_energy));
     }
 
   } else {
@@ -209,7 +190,7 @@ int colvarbias_meta::init_replicas_params(std::string const &conf)
 
     get_keyval(conf, "replicaID", replica_id, replica_id);
     if (!replica_id.size()) {
-      if (proxy->replica_enabled() == COLVARS_OK) {
+      if (proxy->check_replicas_enabled() == COLVARS_OK) {
         // Obtain replicaID from the communicator
         replica_id = cvm::to_str(proxy->replica_index());
         cvm::log("Setting replicaID from communication layer: replicaID = "+
@@ -272,7 +253,6 @@ int colvarbias_meta::init_ebmeta_params(std::string const &conf)
 {
   int error_code = COLVARS_OK;
   // for ebmeta
-  target_dist = NULL;
   get_keyval(conf, "ebMeta", ebmeta, false);
   if(ebmeta){
     cvm::main()->cite_feature("Ensemble-biased metadynamics (ebMetaD)");
@@ -283,7 +263,7 @@ int colvarbias_meta::init_ebmeta_params(std::string const &conf)
                                "targetDistFile accordingly.\n",
                                COLVARS_INPUT_ERROR);
     }
-    target_dist = new colvar_grid_scalar();
+    target_dist.reset(new colvar_grid_scalar());
     error_code |= target_dist->init_from_colvars(colvars);
     std::string target_dist_file;
     get_keyval(conf, "targetDistFile", target_dist_file);
@@ -336,33 +316,15 @@ colvarbias_meta::~colvarbias_meta()
 {
   colvarbias_meta::clear_state_data();
   colvarproxy *proxy = cvm::main()->proxy;
-
   proxy->close_output_stream(replica_hills_file);
-
   proxy->close_output_stream(hills_traj_file_name());
-
-  if (target_dist) {
-    delete target_dist;
-    target_dist = NULL;
-  }
 }
 
 
 int colvarbias_meta::clear_state_data()
 {
-  if (hills_energy) {
-    delete hills_energy;
-    hills_energy = NULL;
-  }
-
-  if (hills_energy_gradients) {
-    delete hills_energy_gradients;
-    hills_energy_gradients = NULL;
-  }
-
   hills.clear();
   hills_off_grid.clear();
-
   return COLVARS_OK;
 }
 
@@ -451,8 +413,11 @@ int colvarbias_meta::update()
   error_code |= update_grid_params();
   // add new biasing energy/forces
   error_code |= update_bias();
-  // update grid content to reflect new bias
+
-  error_code |= update_grid_data();
+  if (use_grids) {
+    // update grid content to reflect new bias
+    error_code |= update_grid_data();
+  }
 
   if (comm != single_replica &&
       (cvm::step_absolute() % replica_update_freq) == 0) {
@@ -539,9 +504,9 @@ int colvarbias_meta::update_grid_params()
         // map everything into new grids
 
         colvar_grid_scalar *new_hills_energy =
-          new colvar_grid_scalar(*hills_energy);
+            new colvar_grid_scalar(*hills_energy);
         colvar_grid_gradient *new_hills_energy_gradients =
-          new colvar_grid_gradient(*hills_energy_gradients);
+            new colvar_grid_gradient(*hills_energy_gradients);
 
         // supply new boundaries to the new grids
 
@@ -556,10 +521,8 @@ int colvarbias_meta::update_grid_params()
         new_hills_energy->map_grid(*hills_energy);
         new_hills_energy_gradients->map_grid(*hills_energy_gradients);
 
-        delete hills_energy;
+        hills_energy.reset(new_hills_energy);
-        delete hills_energy_gradients;
+        hills_energy_gradients.reset(new_hills_energy_gradients);
-        hills_energy = new_hills_energy;
-        hills_energy_gradients = new_hills_energy_gradients;
 
         curr_bin = hills_energy->get_colvars_index();
         if (cvm::debug())
@@ -641,8 +604,7 @@ int colvarbias_meta::update_grid_data()
 {
   if ((cvm::step_absolute() % grids_freq) == 0) {
     // map the most recent gaussians to the grids
-    project_hills(new_hills_begin, hills.end(),
+    project_hills(new_hills_begin, hills.end(), hills_energy.get(), hills_energy_gradients.get());
-                  hills_energy,    hills_energy_gradients);
     new_hills_begin = hills.end();
 
     // TODO: we may want to condense all into one replicas array,
@@ -651,8 +613,8 @@ int colvarbias_meta::update_grid_data()
       for (size_t ir = 0; ir < replicas.size(); ir++) {
         replicas[ir]->project_hills(replicas[ir]->new_hills_begin,
                                     replicas[ir]->hills.end(),
-                                    replicas[ir]->hills_energy,
+                                    replicas[ir]->hills_energy.get(),
-                                    replicas[ir]->hills_energy_gradients);
+                                    replicas[ir]->hills_energy_gradients.get());
         replicas[ir]->new_hills_begin = replicas[ir]->hills.end();
       }
     }
@@ -670,11 +632,20 @@ int colvarbias_meta::calc_energy(std::vector<colvarvalue> const *values)
     replicas[ir]->bias_energy = 0.0;
   }
 
-  std::vector<int> const curr_bin = values ?
+  bool index_ok = false;
-    hills_energy->get_colvars_index(*values) :
+  std::vector<int> curr_bin;
-    hills_energy->get_colvars_index();
 
-  if (hills_energy->index_ok(curr_bin)) {
+  if (use_grids) {
+
+    curr_bin = values ?
+      hills_energy->get_colvars_index(*values) :
+      hills_energy->get_colvars_index();
+
+    index_ok = hills_energy->index_ok(curr_bin);
+
+  }
+
+  if ( index_ok ) {
     // index is within the grid: get the energy from there
     for (ir = 0; ir < replicas.size(); ir++) {
 
@@ -723,11 +694,20 @@ int colvarbias_meta::calc_forces(std::vector<colvarvalue> const *values)
     }
   }
 
-  std::vector<int> const curr_bin = values ?
+  bool index_ok = false;
-    hills_energy->get_colvars_index(*values) :
+  std::vector<int> curr_bin;
-    hills_energy->get_colvars_index();
 
-  if (hills_energy->index_ok(curr_bin)) {
+  if (use_grids) {
+
+    curr_bin = values ?
+      hills_energy->get_colvars_index(*values) :
+      hills_energy->get_colvars_index();
+
+    index_ok = hills_energy->index_ok(curr_bin);
+
+  }
+
+  if ( index_ok ) {
     for (ir = 0; ir < replicas.size(); ir++) {
       cvm::real const *f = &(replicas[ir]->hills_energy_gradients->value(curr_bin));
       for (ic = 0; ic < num_variables(); ic++) {
@@ -959,8 +939,7 @@ void colvarbias_meta::project_hills(colvarbias_meta::hill_iter  h_first,
 
 
 void colvarbias_meta::recount_hills_off_grid(colvarbias_meta::hill_iter  h_first,
-                                             colvarbias_meta::hill_iter  h_last,
+                                             colvarbias_meta::hill_iter  h_last)
-                                             colvar_grid_scalar         * /* he */)
 {
   hills_off_grid.clear();
 
@@ -1078,9 +1057,13 @@ int colvarbias_meta::update_replicas_registry()
         (replicas.back())->comm = multiple_replicas;
 
         if (use_grids) {
-          (replicas.back())->hills_energy           = new colvar_grid_scalar(colvars);
+          (replicas.back())
-          (replicas.back())->hills_energy_gradients = new colvar_grid_gradient(colvars);
+              ->hills_energy.reset(new colvar_grid_scalar(colvars, hills_energy));
+          (replicas.back())
+              ->hills_energy_gradients.reset(
+                  new colvar_grid_gradient(colvars, nullptr, hills_energy));
         }
+
         if (is_enabled(f_cvb_calc_ti_samples)) {
           (replicas.back())->enable(f_cvb_calc_ti_samples);
           (replicas.back())->colvarbias_ti::init_grids();
@@ -1336,34 +1319,40 @@ template <typename IST> IST &colvarbias_meta::read_state_data_template_(IST &is)
 {
   if (use_grids) {
 
-    colvar_grid_scalar   *hills_energy_backup = NULL;
+    std::shared_ptr<colvar_grid_scalar> hills_energy_backup;
-    colvar_grid_gradient *hills_energy_gradients_backup = NULL;
+    std::shared_ptr<colvar_grid_gradient> hills_energy_gradients_backup;
 
-    if (has_data) {
+    bool const need_backup = has_data;
+
+    if (need_backup) {
       if (cvm::debug())
-        cvm::log("Backupping grids for metadynamics bias \""+
+        cvm::log("Backing up grids for metadynamics bias \"" + this->name + "\"" +
-                 this->name+"\""+
+                 ((comm != single_replica) ? ", replica \"" + replica_id + "\"" : "") + ".\n");
-                 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n");
+
-      hills_energy_backup           = hills_energy;
+      hills_energy_backup = std::move(hills_energy);
-      hills_energy_gradients_backup = hills_energy_gradients;
+      hills_energy_gradients_backup = std::move(hills_energy_gradients);
-      hills_energy                  = new colvar_grid_scalar(colvars);
+      hills_energy.reset(new colvar_grid_scalar(colvars, hills_energy));
-      hills_energy_gradients        = new colvar_grid_gradient(colvars);
+      hills_energy_gradients.reset(new colvar_grid_gradient(colvars, nullptr, hills_energy));
     }
 
-    read_grid_data_template_<IST, colvar_grid_scalar>(is, "hills_energy", hills_energy,
+    read_grid_data_template_<IST, colvar_grid_scalar>(is, "hills_energy", hills_energy.get(),
-                                                      hills_energy_backup);
+                                                      hills_energy_backup.get());
 
-    read_grid_data_template_<IST, colvar_grid_gradient>(
+    read_grid_data_template_<IST, colvar_grid_gradient>(is, "hills_energy_gradients",
-        is, "hills_energy_gradients", hills_energy_gradients, hills_energy_gradients_backup);
+                                                        hills_energy_gradients.get(),
+                                                        hills_energy_gradients_backup.get());
 
     if (is) {
       cvm::log("  successfully read the biasing potential and its gradients from grids.\n");
-      if (hills_energy_backup != nullptr) {
-        // Now that we have successfully updated the grids, delete the backup copies
-        delete hills_energy_backup;
-        delete hills_energy_gradients_backup;
-      }
     } else {
+      if (need_backup) {
+        if (cvm::debug())
+          cvm::log("Restoring grids from backup for metadynamics bias \"" + this->name + "\"" +
+                   ((comm != single_replica) ? ", replica \"" + replica_id + "\"" : "") + ".\n");
+        // Restoring content from original grid
+        hills_energy->copy_grid(*hills_energy_backup);
+        hills_energy_gradients->copy_grid(*hills_energy_gradients_backup);
+      }
       return is;
     }
   }
@@ -1451,10 +1440,12 @@ void colvarbias_meta::rebin_grids_after_restart()
     // read from the configuration file), and project onto them the
     // grids just read from the restart file
 
-    colvar_grid_scalar   *new_hills_energy =
+    // Create new grids based on the configuration parameters, because reading from the state
-      new colvar_grid_scalar(colvars);
+    // file automatically sets the old parameters
-    colvar_grid_gradient *new_hills_energy_gradients =
+    std::shared_ptr<colvar_grid_scalar> new_hills_energy(
-      new colvar_grid_gradient(colvars);
+        new colvar_grid_scalar(colvars, nullptr, false, grid_conf));
+    std::shared_ptr<colvar_grid_gradient> new_hills_energy_gradients(
+        new colvar_grid_gradient(colvars, nullptr, new_hills_energy));
 
     if (cvm::debug()) {
       std::ostringstream tmp_os;
@@ -1468,9 +1459,9 @@ void colvarbias_meta::rebin_grids_after_restart()
     if (restart_keep_hills && !hills.empty()) {
       // if there are hills, recompute the new grids from them
       cvm::log("Rebinning the energy and forces grids from "+
-               cvm::to_str(hills.size())+" hills (this may take a while)...\n");
+               cvm::to_str(hills.size())+" hills (this may take a bit)...\n");
-      project_hills(hills.begin(), hills.end(),
+      project_hills(hills.begin(), hills.end(), new_hills_energy.get(),
-                    new_hills_energy, new_hills_energy_gradients, true);
+                    new_hills_energy_gradients.get(), true);
       cvm::log("rebinning done.\n");
 
     } else {
@@ -1481,15 +1472,13 @@ void colvarbias_meta::rebin_grids_after_restart()
       new_hills_energy_gradients->map_grid(*hills_energy_gradients);
     }
 
-    delete hills_energy;
+    hills_energy = std::move(new_hills_energy);
-    delete hills_energy_gradients;
+    hills_energy_gradients = std::move(new_hills_energy_gradients);
-    hills_energy = new_hills_energy;
-    hills_energy_gradients = new_hills_energy_gradients;
 
     // assuming that some boundaries have expanded, eliminate those
     // off-grid hills that aren't necessary any more
     if (!hills.empty())
-      recount_hills_off_grid(hills.begin(), hills.end(), hills_energy);
+      recount_hills_off_grid(hills.begin(), hills.end());
   }
 }
 
@@ -1718,29 +1707,17 @@ int colvarbias_meta::setup_output()
 
   if (comm == multiple_replicas) {
 
-    // TODO: one may want to specify the path manually for intricated filesystems?
+    auto const pwd = cvm::main()->proxy->get_current_work_dir();
-    char *pwd = new char[3001];
-    if (GETCWD(pwd, 3000) == nullptr) {
-      if (pwd != nullptr) { //
-        delete[] pwd;
-      }
-      return cvm::error("Error: cannot get the path of the current working directory.\n",
-                        COLVARS_BUG_ERROR);
-    }
-
     replica_list_file =
-      (std::string(pwd)+std::string(PATHSEP)+
+        cvm::main()->proxy->join_paths(pwd, this->name + "." + replica_id + ".files.txt");
-       this->name+"."+replica_id+".files.txt");
     // replica_hills_file and replica_state_file are those written
     // by the current replica; within the mirror biases, they are
     // those by another replica
-    replica_hills_file =
+    replica_hills_file = cvm::main()->proxy->join_paths(
-      (std::string(pwd)+std::string(PATHSEP)+
+        pwd, cvm::output_prefix() + ".colvars." + this->name + "." + replica_id + ".hills");
-       cvm::output_prefix()+".colvars."+this->name+"."+replica_id+".hills");
+
-    replica_state_file =
+    replica_state_file = cvm::main()->proxy->join_paths(
-      (std::string(pwd)+std::string(PATHSEP)+
+        pwd, cvm::output_prefix() + ".colvars." + this->name + "." + replica_id + ".state");
-       cvm::output_prefix()+".colvars."+this->name+"."+replica_id+".state");
-    delete[] pwd;
 
     // now register this replica
 
@@ -1842,7 +1819,7 @@ template <typename OST> OST &colvarbias_meta::write_state_data_template_(OST &os
 
     // this is a very good time to project hills, if you haven't done
     // it already!
-    project_hills(new_hills_begin, hills.end(), hills_energy, hills_energy_gradients);
+    project_hills(new_hills_begin, hills.end(), hills_energy.get(), hills_energy_gradients.get());
     new_hills_begin = hills.end();
 
     // write down the grids to the restart file

lib/colvars/colvarbias_meta.h

@@ -10,12 +10,16 @@
 #ifndef COLVARBIAS_META_H
 #define COLVARBIAS_META_H
 
-#include <vector>
-#include <list>
 #include <iosfwd>
+#include <list>
+#include <memory>
+#include <vector>
 
 #include "colvarbias.h"
-#include "colvargrid.h"
+
+class colvar_grid_scalar;
+class colvar_grid_gradient;
+
 
 
 /// Metadynamics bias (implementation of \link colvarbias \endlink)
@@ -123,8 +127,7 @@ protected:
   hill_iter new_hills_off_grid_begin;
 
   /// Regenerate the hills_off_grid list
-  void recount_hills_off_grid(hill_iter h_first, hill_iter h_last,
+  void recount_hills_off_grid(hill_iter h_first, hill_iter h_last);
-                              colvar_grid_scalar *ge);
 
   template <typename OST> OST &write_hill_template_(OST &os, colvarbias_meta::hill const &h);
 
@@ -211,7 +214,7 @@ protected:
   bool       ebmeta;
 
   /// Target distribution for EBmeta
-  colvar_grid_scalar* target_dist;
+  std::unique_ptr<colvar_grid_scalar> target_dist;
 
   /// Number of equilibration steps for EBmeta
   cvm::step_number ebmeta_equil_steps;
@@ -223,15 +226,14 @@ protected:
   bool       safely_read_restart;
 
   /// Hill energy, cached on a grid
-  colvar_grid_scalar    *hills_energy;
+  std::shared_ptr<colvar_grid_scalar> hills_energy;
 
   /// Hill forces, cached on a grid
-  colvar_grid_gradient  *hills_energy_gradients;
+  std::shared_ptr<colvar_grid_gradient> hills_energy_gradients;
 
-  /// \brief Project the selected hills onto grids
+  /// Project the selected hills onto grids
-  void project_hills(hill_iter h_first, hill_iter h_last,
+  void project_hills(hill_iter h_first, hill_iter h_last, colvar_grid_scalar *ge,
-                      colvar_grid_scalar *ge, colvar_grid_gradient *gf,
+                     colvar_grid_gradient *gf, bool print_progress = false);
-                      bool print_progress = false);
 
 
   // Multiple Replicas variables and functions

lib/colvars/colvarbias_opes.h

@@ -0,0 +1,176 @@
+#ifndef COLVARBIAS_OPES_H
+#define COLVARBIAS_OPES_H
+
+// This code is mainly adapted from the PLUMED opes module, which uses the
+// LGPLv3 license as shown below:
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+   Copyright (c) 2020-2021 of Michele Invernizzi.
+
+   This file is part of the OPES plumed module.
+
+   The OPES plumed module is free software: you can redistribute it and/or modify
+   it under the terms of the GNU Lesser General Public License as published by
+   the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+
+   The OPES plumed module is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public License
+   along with plumed.  If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+
+#include "colvarbias.h"
+
+#include <vector>
+#include <memory>
+
+// OPES_METAD implementation: swiped from OPESmetad.cpp of PLUMED
+class colvarbias_opes: public colvarbias {
+public:
+  /// The Gaussian kernel data structure
+  struct kernel {
+    cvm::real m_height;
+    std::vector<cvm::real> m_center;
+    std::vector<cvm::real> m_sigma;
+    kernel() {}
+    kernel(cvm::real h, const std::vector<cvm::real>& c,
+           const std::vector<cvm::real>& s):
+      m_height(h), m_center(c), m_sigma(s) {}
+  };
+  /// Communication between different replicas
+  enum Communication {
+    /// One replica (default)
+    single_replica,
+    /// Hills added concurrently by several replicas
+    multiple_replicas
+  };
+  /// Constructor
+  colvarbias_opes(char const *key);
+  /// Initializer
+  int init(std::string const &conf) override;
+  /// Per-timestep update
+  int update() override;
+  /// Save the state to a text file for restarting
+  std::ostream &write_state_data(std::ostream &os) override;
+  /// Read the state from a text file for restarting
+  std::istream &read_state_data(std::istream &is) override;
+  /// Save the state to a binary file for restarting
+  cvm::memory_stream &write_state_data(cvm::memory_stream &os) override;
+  /// Read the state from a binary file for restarting
+  cvm::memory_stream &read_state_data(cvm::memory_stream &is) override;
+  /// Write to files at restart steps
+  int write_output_files() override;
+private:
+  int update_opes();
+  int calculate_opes();
+  void save_state();
+  cvm::real getProbAndDerivatives(const std::vector<cvm::real>& cv, std::vector<cvm::real>& der_prob) const;
+  cvm::real evaluateKernel(const kernel& G, const std::vector<cvm::real>& x) const;
+  cvm::real evaluateKernel(const kernel& G, const std::vector<cvm::real>& x, std::vector<cvm::real>& accumulated_derivative, std::vector<cvm::real>& dist) const;
+  void addKernel(const double height, const std::vector<cvm::real>& center, const std::vector<cvm::real>& sigma, const double logweight);
+  void addKernel(const double height, const std::vector<cvm::real>& center, const std::vector<cvm::real>& sigma);
+  size_t getMergeableKernel(const std::vector<cvm::real>& giver_center, const size_t giver_k) const;
+  void mergeKernels(kernel& k1, const kernel& k2) const;
+  void updateNlist(const std::vector<cvm::real>& center);
+  struct traj_line {
+    double rct;
+    double zed;
+    double neff;
+    double work;
+    size_t nker;
+    size_t nlker;
+    size_t nlsteps;
+  };
+  void writeTrajBuffer();
+  void showInfo() const;
+  template <typename OST> OST &write_state_data_template_(OST &os) const;
+  template <typename IST> IST &read_state_data_template_(IST &os);
+  std::string const traj_file_name(const std::string& suffix) const;
+  int collectSampleToPMFGrid();
+  int computePMF();
+  int writePMF(const std::unique_ptr<colvar_grid_scalar>& pmf_grid, const std::string &filename, bool keep_open);
+private:
+  cvm::real m_kbt;
+  cvm::real m_barrier;
+  cvm::real m_biasfactor;
+  cvm::real m_bias_prefactor;
+  cvm::real m_temperature;
+  cvm::step_number m_pace;
+  cvm::step_number m_adaptive_sigma_stride;
+  cvm::step_number m_adaptive_counter;
+  unsigned long long m_counter;
+  cvm::real m_compression_threshold;
+  cvm::real m_compression_threshold2;
+  bool m_adaptive_sigma;
+  bool m_fixed_sigma;
+  bool m_no_zed;
+  // bool m_restart;
+  bool m_nlist;
+  bool m_recursive_merge;
+  std::vector<cvm::real> m_nlist_param;
+  std::vector<cvm::real> m_sigma0;
+  std::vector<cvm::real> m_sigma_min;
+  cvm::real m_epsilon;
+  cvm::real m_sum_weights;
+  cvm::real m_sum_weights2;
+  cvm::real m_cutoff;
+  cvm::real m_cutoff2;
+  cvm::real m_zed;
+  cvm::real m_old_kdenorm;
+  cvm::real m_kdenorm;
+  cvm::real m_val_at_cutoff;
+  cvm::real m_rct;
+  cvm::real m_neff;
+  std::vector<kernel> m_kernels;
+  std::vector<kernel> m_delta_kernels;
+  std::vector<cvm::real> m_av_cv;
+  std::vector<cvm::real> m_av_M2;
+  std::ostringstream m_kernels_output;
+  std::vector<cvm::real> m_nlist_center;
+  std::vector<size_t> m_nlist_index;
+  std::vector<cvm::real> m_nlist_dev2;
+  size_t m_nlist_steps;
+  bool m_nlist_update;
+  bool m_nlist_pace_reset;
+  size_t m_nker;
+  bool m_calc_work;
+  cvm::real m_work;
+  /// Communication between different replicas
+  Communication comm;
+  /// \brief Identifier for this replica
+  std::string            replica_id;
+  size_t m_num_walkers;
+  size_t shared_freq;
+  size_t m_num_threads;
+  size_t m_nlker;
+  // size_t m_state_stride;
+  // std::unordered_map<std::string, std::string> m_kernel_output_components;
+  std::string m_kernels_output_headers;
+  cvm::step_number m_traj_output_frequency;
+  traj_line m_traj_line;
+  std::ostringstream m_traj_oss;
+  bool m_is_first_step;
+  std::vector<cvm::real> m_cv;
+  // For saving states
+  decltype(m_zed) m_saved_zed;
+  decltype(m_sum_weights) m_saved_sum_weights;
+  decltype(m_sum_weights2) m_saved_sum_weights2;
+  decltype(m_kernels) m_saved_kernels;
+  // PMF grid from reweighting
+  bool m_pmf_grid_on;
+  std::vector<colvar*> m_pmf_cvs;
+  std::string grid_conf;
+  std::shared_ptr<colvar_grid_scalar> m_reweight_grid;
+  std::unique_ptr<colvar_grid_scalar> m_pmf_grid;
+  cvm::step_number m_pmf_hist_freq;
+  bool m_pmf_shared; // shared PMF among replicas
+  std::unique_ptr<colvar_grid_scalar> m_global_reweight_grid;
+  std::unique_ptr<colvar_grid_scalar> m_global_pmf_grid;
+  bool m_explore;
+  bool m_inf_biasfactor;
+};
+
+#endif // COLVARBIAS_OPES_H

lib/colvars/colvarcomp.cpp

@@ -261,7 +261,6 @@ int colvar::cvc::init_dependencies() {
     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);
@@ -525,7 +524,7 @@ void colvar::cvc::calc_force_invgrads()
 
 void colvar::cvc::calc_Jacobian_derivative()
 {
-  cvm::error("Error: calculation of inverse gradients is not implemented "
+  cvm::error("Error: calculation of Jacobian derivatives is not implemented "
              "for colvar components of type \""+function_type()+"\".\n",
              COLVARS_NOT_IMPLEMENTED);
 }
@@ -533,8 +532,10 @@ void colvar::cvc::calc_Jacobian_derivative()
 
 void colvar::cvc::calc_fit_gradients()
 {
-  for (size_t ig = 0; ig < atom_groups.size(); ig++) {
+  if (is_enabled(f_cvc_explicit_gradient)) {
-    atom_groups[ig]->calc_fit_gradients();
+    for (size_t ig = 0; ig < atom_groups.size(); ig++) {
+      atom_groups[ig]->calc_fit_gradients();
+    }
   }
 }
 

lib/colvars/colvarcomp.h

@@ -233,8 +233,14 @@ public:
 
   /// Forcibly set value of CVC - useful for driving an external coordinate,
   /// eg. lambda dynamics
-  inline void set_value(colvarvalue const &new_value) {
+  inline void set_value(colvarvalue const &new_value, bool now=false) {
     x = new_value;
+    // Cache value to be communicated to back-end between time steps
+    cvm::proxy->set_alch_lambda(x.real_value);
+    if (now) {
+      // If requested (e.g. upon restarting), sync to back-end
+      cvm::proxy->send_alch_lambda();
+    }
   }
 
 protected:
@@ -1212,9 +1218,11 @@ protected:
   // No atom groups needed
 public:
   alch_lambda();
+  int init_alchemy(int time_step_factor);
   virtual ~alch_lambda() {}
   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);
 };
 

lib/colvars/colvarcomp_alchlambda.cpp

@@ -20,22 +20,46 @@ colvar::alch_lambda::alch_lambda()
 {
   set_function_type("alchLambda");
 
-  disable(f_cvc_explicit_gradient);
+  provide(f_cvc_explicit_gradient, false);
-  disable(f_cvc_gradient);
+  provide(f_cvc_gradient, false); // Cannot apply forces on this CVC
+  provide(f_cvc_collect_atom_ids, false);
+
+  provide(f_cvc_inv_gradient); // Projected force is TI derivative
+  provide(f_cvc_Jacobian);     // Zero
 
   x.type(colvarvalue::type_scalar);
-  // Query initial value from back-end
+
+  // Query initial value from back-end; will be overwritten if restarting from a state file
   cvm::proxy->get_alch_lambda(&x.real_value);
 }
 
 
+int colvar::alch_lambda::init_alchemy(int factor)
+{
+  // We need calculation every time step
+  // default in Tinker-HP and NAMD2, must be enforced in NAMD3
+  // Also checks back-end settings, ie. that alchemy is enabled
+  // (in NAMD3: alchType TI, computeEnergies at the right frequency)
+
+  // Forbid MTS until fully implemented
+  if (factor != 1) {
+    return cvm::error("Error: timeStepFactor > 1 is not yet supported for alchemical variables.");
+  }
+  cvm::proxy->request_alch_energy_freq(factor);
+
+  return COLVARS_OK;
+}
+
+
 void colvar::alch_lambda::calc_value()
 {
-  // Special workflow:
+  // By default, follow external parameter
-  // at the beginning of the timestep we get a force instead of calculating the value
+  // This might get overwritten by driving extended dynamics
+  // (in apply_force() below)
+  cvm::proxy->get_alch_lambda(&x.real_value);
 
   cvm::proxy->get_dE_dlambda(&ft.real_value);
-  ft.real_value *= -1.0; // Energy derivative to force
+  ft.real_value *= -1.0; // Convert energy derivative to force
 
   // Include any force due to bias on Flambda
   ft.real_value += cvm::proxy->indirect_lambda_biasing_force;
@@ -43,19 +67,24 @@ void colvar::alch_lambda::calc_value()
 }
 
 
-void colvar::alch_lambda::calc_gradients()
+void colvar::alch_lambda::calc_force_invgrads()
 {
+  // All the work is done in calc_value()
+}
+
+
+void colvar::alch_lambda::calc_Jacobian_derivative()
+{
+  jd = 0.0;
 }
 
 
 void colvar::alch_lambda::apply_force(colvarvalue const & /* force */)
 {
-  // new value will be cached and sent at end of timestep
+  // Forces, if any, are applied in colvar::update_extended_Lagrangian()
-  cvm::proxy->set_alch_lambda(x.real_value);
 }
 
 
-
 colvar::alch_Flambda::alch_Flambda()
 {
   set_function_type("alch_Flambda");

lib/colvars/colvarcomp_angles.cpp

@@ -267,74 +267,22 @@ void colvar::dihedral::calc_value()
 
 void colvar::dihedral::calc_gradients()
 {
-  cvm::rvector A = cvm::rvector::outer(r12, r23);
+  // Eqs. (27i) ~ (27l) from https://doi.org/10.1002/(SICI)1096-987X(19960715)17:9<1132::AID-JCC5>3.0.CO;2-T.
-  cvm::real   rA = A.norm();
-  cvm::rvector B = cvm::rvector::outer(r23, r34);
-  cvm::real   rB = B.norm();
-  cvm::rvector C = cvm::rvector::outer(r23, A);
-  cvm::real   rC = C.norm();
 
-  cvm::real const cos_phi = (A*B)/(rA*rB);
+  const cvm::rvector A = cvm::rvector::outer(r12, r23);
-  cvm::real const sin_phi = (C*B)/(rC*rB);
+  const cvm::rvector B = cvm::rvector::outer(r23, r34);
+  const cvm::real   nG = r23.norm();
+  const cvm::real   A2 = A.norm2();
+  const cvm::real   B2 = B.norm2();
 
-  cvm::rvector f1, f2, f3;
+  const cvm::real     K = 180.0/PI;
-
+  const cvm::rvector f1 = K * nG / A2 * A;
-  rB = 1.0/rB;
+  const cvm::rvector f2 = K * ((r12 * r23 / (A2 * nG)) * A + (r34 * r23 / (B2 * nG)) * B);
-  B *= rB;
+  const cvm::rvector f3 = K * nG / B2 * B;
-
-  if (cvm::fabs(sin_phi) > 0.1) {
-    rA = 1.0/rA;
-    A *= rA;
-    cvm::rvector const dcosdA = rA*(cos_phi*A-B);
-    cvm::rvector const dcosdB = rB*(cos_phi*B-A);
-    // rA = 1.0;
-
-    cvm::real const K = (1.0/sin_phi) * (180.0/PI);
-
-        f1 = K * cvm::rvector::outer(r23, dcosdA);
-        f3 = K * cvm::rvector::outer(dcosdB, r23);
-        f2 = K * (cvm::rvector::outer(dcosdA, r12)
-                   +  cvm::rvector::outer(r34, dcosdB));
-  }
-  else {
-    rC = 1.0/rC;
-    C *= rC;
-    cvm::rvector const dsindC = rC*(sin_phi*C-B);
-    cvm::rvector const dsindB = rB*(sin_phi*B-C);
-    // rC = 1.0;
-
-    cvm::real    const K = (-1.0/cos_phi) * (180.0/PI);
-
-    f1.x = K*((r23.y*r23.y + r23.z*r23.z)*dsindC.x
-              - r23.x*r23.y*dsindC.y
-              - r23.x*r23.z*dsindC.z);
-    f1.y = K*((r23.z*r23.z + r23.x*r23.x)*dsindC.y
-              - r23.y*r23.z*dsindC.z
-              - r23.y*r23.x*dsindC.x);
-    f1.z = K*((r23.x*r23.x + r23.y*r23.y)*dsindC.z
-              - r23.z*r23.x*dsindC.x
-              - r23.z*r23.y*dsindC.y);
-
-    f3 = cvm::rvector::outer(dsindB, r23);
-    f3 *= K;
-
-    f2.x = K*(-(r23.y*r12.y + r23.z*r12.z)*dsindC.x
-              +(2.0*r23.x*r12.y - r12.x*r23.y)*dsindC.y
-              +(2.0*r23.x*r12.z - r12.x*r23.z)*dsindC.z
-              +dsindB.z*r34.y - dsindB.y*r34.z);
-    f2.y = K*(-(r23.z*r12.z + r23.x*r12.x)*dsindC.y
-              +(2.0*r23.y*r12.z - r12.y*r23.z)*dsindC.z
-              +(2.0*r23.y*r12.x - r12.y*r23.x)*dsindC.x
-              +dsindB.x*r34.z - dsindB.z*r34.x);
-    f2.z = K*(-(r23.x*r12.x + r23.y*r12.y)*dsindC.z
-              +(2.0*r23.z*r12.x - r12.z*r23.x)*dsindC.x
-              +(2.0*r23.z*r12.y - r12.z*r23.y)*dsindC.y
-              +dsindB.y*r34.x - dsindB.x*r34.y);
-  }
 
   group1->set_weighted_gradient(-f1);
-  group2->set_weighted_gradient(-f2 + f1);
+  group2->set_weighted_gradient( f2 + f1);
-  group3->set_weighted_gradient(-f3 + f2);
+  group3->set_weighted_gradient(-f3 - f2);
   group4->set_weighted_gradient(f3);
 }
 

lib/colvars/colvarcomp_distances.cpp

@@ -384,32 +384,30 @@ void colvar::distance_dir::apply_force(colvarvalue const &force)
   cvm::real const iprod = force.rvector_value * x.rvector_value;
   cvm::rvector const force_tang = force.rvector_value - iprod * x.rvector_value;
 
-  if (!group1->noforce)
+  if (!group1->noforce) {
-    group1->apply_force(-1.0 * force_tang);
+    group1->apply_force(-1.0 / dist_v.norm() * force_tang);
-
+  }
-  if (!group2->noforce)
+  if (!group2->noforce) {
-    group2->apply_force(       force_tang);
+    group2->apply_force( 1.0 / dist_v.norm() * force_tang);
+  }
 }
 
 
-cvm::real colvar::distance_dir::dist2(colvarvalue const &x1,
+cvm::real colvar::distance_dir::dist2(colvarvalue const &x1, colvarvalue const &x2) const
-                                      colvarvalue const &x2) const
 {
-  return (x1.rvector_value - x2.rvector_value).norm2();
+  return x1.dist2(x2);
 }
 
 
-colvarvalue colvar::distance_dir::dist2_lgrad(colvarvalue const &x1,
+colvarvalue colvar::distance_dir::dist2_lgrad(colvarvalue const &x1, colvarvalue const &x2) const
-                                              colvarvalue const &x2) const
 {
-  return colvarvalue((x1.rvector_value - x2.rvector_value), colvarvalue::type_unit3vectorderiv);
+  return x1.dist2_grad(x2);
 }
 
 
-colvarvalue colvar::distance_dir::dist2_rgrad(colvarvalue const &x1,
+colvarvalue colvar::distance_dir::dist2_rgrad(colvarvalue const &x1, colvarvalue const &x2) const
-                                              colvarvalue const &x2) const
 {
-  return colvarvalue((x2.rvector_value - x1.rvector_value), colvarvalue::type_unit3vectorderiv);
+  return x2.dist2_grad(x1);
 }
 
 
@@ -1005,7 +1003,7 @@ void colvar::rmsd::calc_Jacobian_derivative()
     for (size_t ia = 0; ia < atoms->size(); ia++) {
 
       // Gradient of optimal quaternion wrt current Cartesian position
-      atoms->rot_deriv->calc_derivative_wrt_group1(ia, nullptr, &dq);
+      atoms->rot_deriv->calc_derivative_wrt_group1<false, true, false>(ia, nullptr, &dq);
 
       g11 = 2.0 * (atoms->rot.q)[1]*dq[1];
       g22 = 2.0 * (atoms->rot.q)[2]*dq[2];
@@ -1304,7 +1302,7 @@ void colvar::eigenvector::calc_Jacobian_derivative()
     // Gradient of optimal quaternion wrt current Cartesian position
     // trick: d(R^-1)/dx = d(R^t)/dx = (dR/dx)^t
     // we can just transpose the derivatives of the direct matrix
-    atoms->rot_deriv->calc_derivative_wrt_group1(ia, nullptr, &dq_1);
+    atoms->rot_deriv->calc_derivative_wrt_group1<false, true, false>(ia, nullptr, &dq_1);
 
     g11 = 2.0 * quat0[1]*dq_1[1];
     g22 = 2.0 * quat0[2]*dq_1[2];
@@ -1403,11 +1401,12 @@ void colvar::cartesian::apply_force(colvarvalue const &force)
   size_t ia, j;
   if (!atoms->noforce) {
     cvm::rvector f;
+    auto ag_force = atoms->get_group_force_object();
     for (ia = 0; ia < atoms->size(); ia++) {
       for (j = 0; j < dim; j++) {
         f[axes[j]] = force.vector1d_value[dim*ia + j];
       }
-      (*atoms)[ia].apply_force(f);
+      ag_force.add_atom_force(ia, f);
     }
   }
 }

lib/colvars/colvarcomp_protein.cpp

@@ -28,34 +28,58 @@ colvar::alpha_angles::alpha_angles()
 int colvar::alpha_angles::init(std::string const &conf)
 {
   int error_code = cvc::init(conf);
+  if (error_code != COLVARS_OK) return error_code;
 
   std::string segment_id;
-  get_keyval(conf, "psfSegID", segment_id, std::string("MAIN"));
-
   std::vector<int> residues;
-  {
+
-    std::string residues_conf = "";
+  bool b_use_index_groups = false;
-    key_lookup(conf, "residueRange", &residues_conf);
+  cvm::atom_group group_CA, group_N, group_O;
+
+  std::string residues_conf = "";
+  std::string prefix;
+
+  // residueRange is mandatory for the topology-based case
+  if (key_lookup(conf, "residueRange", &residues_conf)) {
     if (residues_conf.size()) {
       std::istringstream is(residues_conf);
       int initial, final;
       char dash;
       if ( (is >> initial) && (initial > 0) &&
-           (is >> dash) && (dash == '-') &&
+          (is >> dash) && (dash == '-') &&
-           (is >> final) && (final > 0) ) {
+          (is >> final) && (final > 0) ) {
         for (int rnum = initial; rnum <= final; rnum++) {
           residues.push_back(rnum);
         }
       }
     } else {
-      error_code |=
+      return cvm::error("Error: no residues defined in \"residueRange\".\n", COLVARS_INPUT_ERROR);
-          cvm::error("Error: no residues defined in \"residueRange\".\n", COLVARS_INPUT_ERROR);
     }
-  }
 
-  if (residues.size() < 5) {
+    if (residues.size() < 5) {
-    error_code |= cvm::error("Error: not enough residues defined in \"residueRange\".\n",
+      return cvm::error("Error: not enough residues defined in \"residueRange\".\n", COLVARS_INPUT_ERROR);
-                             COLVARS_INPUT_ERROR);
+    }
+    get_keyval(conf, "psfSegID", segment_id, std::string("MAIN"));
+
+  } else {
+    b_use_index_groups = true;
+    get_keyval(conf, "prefix", prefix, "alpha_");
+
+    // Not all groups are mandatory, parse silently
+    group_CA.add_index_group(prefix + "CA", true);
+    group_N.add_index_group(prefix + "N", true);
+    group_O.add_index_group(prefix + "O", true);
+    int na = group_CA.size();
+    int nn = group_N.size();
+    int no = group_O.size();
+    if ((nn != 0 || no != 0) && (nn != no)) {
+      return cvm::error("Error: If either is provided, atom groups " + prefix + "N and " + prefix + "O must have the same number of atoms.",
+                        COLVARS_INPUT_ERROR);
+    }
+    if (nn != 0 && na != 0 && nn != na) {
+      return cvm::error("Error: If both are provided, atom groups " + prefix + "N and " + prefix + "CA must have the same number of atoms.",
+                        COLVARS_INPUT_ERROR);
+    }
   }
 
   std::string const &sid    = segment_id;
@@ -64,8 +88,7 @@ int colvar::alpha_angles::init(std::string const &conf)
 
   get_keyval(conf, "hBondCoeff", hb_coeff, hb_coeff);
   if ((hb_coeff < 0.0) || (hb_coeff > 1.0)) {
-    error_code |=
+    return cvm::error("Error: hBondCoeff must be defined between 0 and 1.\n", COLVARS_INPUT_ERROR);
-        cvm::error("Error: hBondCoeff must be defined between 0 and 1.\n", COLVARS_INPUT_ERROR);
   }
 
 
@@ -73,14 +96,29 @@ int colvar::alpha_angles::init(std::string const &conf)
   get_keyval(conf, "angleTol", theta_tol, theta_tol);
 
   if (hb_coeff < 1.0) {
-
+    if (b_use_index_groups) {
-    for (size_t i = 0; i < residues.size()-2; i++) {
+      if (group_CA.size() < 5) {
-      theta.push_back(new colvar::angle(cvm::atom(r[i  ], "CA", sid),
+        return cvm::error("Not enough atoms (" + cvm::to_str(group_CA.size()) + ") in index group \"" + prefix + "CA\"",
-                                        cvm::atom(r[i+1], "CA", sid),
+                          COLVARS_INPUT_ERROR);
-                                        cvm::atom(r[i+2], "CA", sid)));
+      }
-      register_atom_group(theta.back()->atom_groups[0]);
+      for (size_t i = 0; i < group_CA.size()-2; i++) {
-      register_atom_group(theta.back()->atom_groups[1]);
+        // Note: the angle constructor constructs copies of the atom objects
-      register_atom_group(theta.back()->atom_groups[2]);
+        theta.push_back(new colvar::angle(group_CA[i],
+                                          group_CA[i+1],
+                                          group_CA[i+2]));
+        register_atom_group(theta.back()->atom_groups[0]);
+        register_atom_group(theta.back()->atom_groups[1]);
+        register_atom_group(theta.back()->atom_groups[2]);
+      }
+    } else {
+      for (size_t i = 0; i < residues.size()-2; i++) {
+        theta.push_back(new colvar::angle(cvm::atom(r[i  ], "CA", sid),
+                                          cvm::atom(r[i+1], "CA", sid),
+                                          cvm::atom(r[i+2], "CA", sid)));
+        register_atom_group(theta.back()->atom_groups[0]);
+        register_atom_group(theta.back()->atom_groups[1]);
+        register_atom_group(theta.back()->atom_groups[2]);
+      }
     }
 
   } else {
@@ -93,14 +131,27 @@ int colvar::alpha_angles::init(std::string const &conf)
     get_keyval(conf, "hBondExpDenom", ed, ed);
 
     if (hb_coeff > 0.0) {
-
+      if (b_use_index_groups) {
-      for (size_t i = 0; i < residues.size()-4; i++) {
+        if (group_N.size() < 5) {
-        hb.push_back(new colvar::h_bond(cvm::atom(r[i  ], "O",  sid),
+          return cvm::error("Not enough atoms (" + cvm::to_str(group_N.size()) + ") in index group \"" + prefix + "N\"",
-                                        cvm::atom(r[i+4], "N",  sid),
+                            COLVARS_INPUT_ERROR);
-                                        r0, en, ed));
+        }
-        register_atom_group(hb.back()->atom_groups[0]);
+        for (size_t i = 0; i < group_N.size()-4; i++) {
+          // Note: we need to call the atom copy constructor here because
+          // the h_bond constructor does not make copies of the provided atoms
+          hb.push_back(new colvar::h_bond(cvm::atom(group_O[i]),
+                                          cvm::atom(group_N[i+4]),
+                                          r0, en, ed));
+          register_atom_group(hb.back()->atom_groups[0]);
+        }
+      } else {
+        for (size_t i = 0; i < residues.size()-4; i++) {
+          hb.push_back(new colvar::h_bond(cvm::atom(r[i  ], "O",  sid),
+                                          cvm::atom(r[i+4], "N",  sid),
+                                          r0, en, ed));
+          register_atom_group(hb.back()->atom_groups[0]);
+        }
       }
-
     } else {
       cvm::log("The hBondCoeff specified will disable the hydrogen bond terms.\n");
     }
@@ -290,41 +341,62 @@ int colvar::dihedPC::init(std::string const &conf)
   if (cvm::debug())
     cvm::log("Initializing dihedral PC object.\n");
 
+  bool b_use_index_groups = false;
   std::string segment_id;
-  get_keyval(conf, "psfSegID", segment_id, std::string("MAIN"));
-
   std::vector<int> residues;
-  {
+  size_t n_residues;
-    std::string residues_conf = "";
+  std::string residues_conf = "";
-    key_lookup(conf, "residueRange", &residues_conf);
+  std::string prefix;
+  cvm::atom_group group_CA, group_N, group_C;
+
+  // residueRange is mandatory for the topology-based case
+  if (key_lookup(conf, "residueRange", &residues_conf)) {
     if (residues_conf.size()) {
       std::istringstream is(residues_conf);
       int initial, final;
       char dash;
       if ( (is >> initial) && (initial > 0) &&
-           (is >> dash) && (dash == '-') &&
+          (is >> dash) && (dash == '-') &&
-           (is >> final) && (final > 0) ) {
+          (is >> final) && (final > 0) ) {
         for (int rnum = initial; rnum <= final; rnum++) {
           residues.push_back(rnum);
         }
       }
     } else {
-      error_code |=
+      return cvm::error("Error: no residues defined in \"residueRange\".\n", COLVARS_INPUT_ERROR);
-          cvm::error("Error: no residues defined in \"residueRange\".\n", COLVARS_INPUT_ERROR);
     }
-  }
+    n_residues = residues.size();
+    get_keyval(conf, "psfSegID", segment_id, std::string("MAIN"));
 
-  if (residues.size() < 2) {
+  } else {
+
+    b_use_index_groups = true;
+    get_keyval(conf, "prefix", prefix, "dihed_");
+
+    // All three groups are required
+    group_CA.add_index_group(prefix + "CA");
+    group_N.add_index_group(prefix + "N");
+    group_C.add_index_group(prefix + "C");
+    int na = group_CA.size();
+    int nn = group_N.size();
+    int nc = group_C.size();
+    if ((nn != na || na != nc)) {
+      return cvm::error("Error: atom groups " + prefix + "N, " + prefix + "CA, and " + prefix +
+                        "C must have the same number of atoms.", COLVARS_INPUT_ERROR);
+    }
+    n_residues = nn;
+  }
+  if (n_residues < 2) {
     error_code |=
-        cvm::error("Error: dihedralPC requires at least two residues.\n", COLVARS_INPUT_ERROR);
+      cvm::error("Error: dihedralPC requires at least two residues.\n", COLVARS_INPUT_ERROR);
   }
 
   std::string const &sid    = segment_id;
   std::vector<int> const &r = residues;
 
   std::string vecFileName;
-  int         vecNumber;
   if (get_keyval(conf, "vectorFile", vecFileName, vecFileName)) {
+    int vecNumber;
     get_keyval(conf, "vectorNumber", vecNumber, 0);
     if (vecNumber < 1) {
       error_code |=
@@ -339,9 +411,8 @@ int colvar::dihedPC::init(std::string const &conf)
     }
 
     // TODO: adapt to different formats by setting this flag
-    bool eigenvectors_as_columns = true;
+    // bool eigenvectors_as_columns = true;
-
+    // if (eigenvectors_as_columns) {
-    if (eigenvectors_as_columns) {
       // Carma-style dPCA file
       std::string line;
       cvm::real c;
@@ -352,9 +423,7 @@ int colvar::dihedPC::init(std::string const &conf)
         for (int i=0; i<vecNumber; i++) ls >> c;
         coeffs.push_back(c);
       }
-    }
+    /* } else { // Uncomment this when different formats are recognized
-/*  TODO Uncomment this when different formats are recognized
-    else {
       // Eigenvectors as lines
       // Skip to the right line
       for (int i = 1; i<vecNumber; i++)
@@ -380,28 +449,42 @@ int colvar::dihedPC::init(std::string const &conf)
     get_keyval(conf, "vector", coeffs, coeffs);
   }
 
-  if ( coeffs.size() != 4 * (residues.size() - 1)) {
+  if ( coeffs.size() != 4 * (n_residues - 1)) {
     error_code |= cvm::error("Error: wrong number of coefficients: " + cvm::to_str(coeffs.size()) +
-                             ". Expected " + cvm::to_str(4 * (residues.size() - 1)) +
+                             ". Expected " + cvm::to_str(4 * (n_residues - 1)) +
                              " (4 coeffs per residue, minus one residue).\n",
                              COLVARS_INPUT_ERROR);
   }
 
-  for (size_t i = 0; i < residues.size()-1; i++) {
+  for (size_t i = 0; i < n_residues-1; i++) {
     // Psi
-    theta.push_back(new colvar::dihedral(cvm::atom(r[i  ], "N", sid),
+    if (b_use_index_groups) {
-                                         cvm::atom(r[i  ], "CA", sid),
+      theta.push_back(new colvar::dihedral( group_N[i],
-                                         cvm::atom(r[i  ], "C", sid),
+                                            group_CA[i],
-                                         cvm::atom(r[i+1], "N", sid)));
+                                            group_C[i],
+                                            group_N[i+1]));
+    } else {
+      theta.push_back(new colvar::dihedral(cvm::atom(r[i  ], "N", sid),
+                                           cvm::atom(r[i  ], "CA", sid),
+                                           cvm::atom(r[i  ], "C", sid),
+                                           cvm::atom(r[i+1], "N", sid)));
+    }
     register_atom_group(theta.back()->atom_groups[0]);
     register_atom_group(theta.back()->atom_groups[1]);
     register_atom_group(theta.back()->atom_groups[2]);
     register_atom_group(theta.back()->atom_groups[3]);
     // Phi (next res)
-    theta.push_back(new colvar::dihedral(cvm::atom(r[i  ], "C", sid),
+    if (b_use_index_groups) {
-                                         cvm::atom(r[i+1], "N", sid),
+      theta.push_back(new colvar::dihedral(group_C[i],
-                                         cvm::atom(r[i+1], "CA", sid),
+                                           group_N[i+1],
-                                         cvm::atom(r[i+1], "C", sid)));
+                                           group_CA[i+1],
+                                           group_C[i+1]));
+    } else {
+      theta.push_back(new colvar::dihedral(cvm::atom(r[i  ], "C", sid),
+                                           cvm::atom(r[i+1], "N", sid),
+                                           cvm::atom(r[i+1], "CA", sid),
+                                           cvm::atom(r[i+1], "C", sid)));
+    }
     register_atom_group(theta.back()->atom_groups[0]);
     register_atom_group(theta.back()->atom_groups[1]);
     register_atom_group(theta.back()->atom_groups[2]);

lib/colvars/colvarcomp_rotations.cpp

@@ -137,11 +137,14 @@ void colvar::orientation::apply_force(colvarvalue const &force)
   if (!atoms->noforce) {
     rot_deriv_impl->prepare_derivative(rotation_derivative_dldq::use_dq);
     cvm::vector1d<cvm::rvector> dq0_2;
+    auto ag_force = atoms->get_group_force_object();
     for (size_t ia = 0; ia < atoms->size(); ia++) {
-      rot_deriv_impl->calc_derivative_wrt_group2(ia, nullptr, &dq0_2);
+      rot_deriv_impl->calc_derivative_wrt_group2<false, true, false>(ia, nullptr, &dq0_2);
-      for (size_t i = 0; i < 4; i++) {
+      const auto f_ia = FQ[0] * dq0_2[0] +
-        (*atoms)[ia].apply_force(FQ[i] * dq0_2[i]);
+                        FQ[1] * dq0_2[1] +
-      }
+                        FQ[2] * dq0_2[2] +
+                        FQ[3] * dq0_2[3];
+      ag_force.add_atom_force(ia, f_ia);
     }
   }
 }
@@ -205,7 +208,7 @@ void colvar::orientation_angle::calc_gradients()
   rot_deriv_impl->prepare_derivative(rotation_derivative_dldq::use_dq);
   cvm::vector1d<cvm::rvector> dq0_2;
   for (size_t ia = 0; ia < atoms->size(); ia++) {
-    rot_deriv_impl->calc_derivative_wrt_group2(ia, nullptr, &dq0_2);
+    rot_deriv_impl->calc_derivative_wrt_group2<false, true, false>(ia, nullptr, &dq0_2);
     (*atoms)[ia].grad = (dxdq0 * dq0_2[0]);
   }
 }
@@ -265,7 +268,7 @@ void colvar::orientation_proj::calc_gradients()
   rot_deriv_impl->prepare_derivative(rotation_derivative_dldq::use_dq);
   cvm::vector1d<cvm::rvector> dq0_2;
   for (size_t ia = 0; ia < atoms->size(); ia++) {
-    rot_deriv_impl->calc_derivative_wrt_group2(ia, nullptr, &dq0_2);
+    rot_deriv_impl->calc_derivative_wrt_group2<false, true, false>(ia, nullptr, &dq0_2);
     (*atoms)[ia].grad = (dxdq0 * dq0_2[0]);
   }
 }
@@ -314,7 +317,7 @@ void colvar::tilt::calc_gradients()
   cvm::vector1d<cvm::rvector> dq0_2;
   for (size_t ia = 0; ia < atoms->size(); ia++) {
     (*atoms)[ia].grad = cvm::rvector(0.0, 0.0, 0.0);
-    rot_deriv_impl->calc_derivative_wrt_group2(ia, nullptr, &dq0_2);
+    rot_deriv_impl->calc_derivative_wrt_group2<false, true, false>(ia, nullptr, &dq0_2);
     for (size_t iq = 0; iq < 4; iq++) {
       (*atoms)[ia].grad += (dxdq[iq] * dq0_2[iq]);
     }
@@ -351,7 +354,7 @@ void colvar::spin_angle::calc_gradients()
   cvm::vector1d<cvm::rvector> dq0_2;
   for (size_t ia = 0; ia < atoms->size(); ia++) {
     (*atoms)[ia].grad = cvm::rvector(0.0, 0.0, 0.0);
-    rot_deriv_impl->calc_derivative_wrt_group2(ia, nullptr, &dq0_2);
+    rot_deriv_impl->calc_derivative_wrt_group2<false, true, false>(ia, nullptr, &dq0_2);
     for (size_t iq = 0; iq < 4; iq++) {
       (*atoms)[ia].grad += (dxdq[iq] * dq0_2[iq]);
     }
@@ -399,7 +402,7 @@ void colvar::euler_phi::calc_gradients()
   rot_deriv_impl->prepare_derivative(rotation_derivative_dldq::use_dq);
   cvm::vector1d<cvm::rvector> dq0_2;
   for (size_t ia = 0; ia < atoms->size(); ia++) {
-    rot_deriv_impl->calc_derivative_wrt_group2(ia, nullptr, &dq0_2);
+    rot_deriv_impl->calc_derivative_wrt_group2<false, true, false>(ia, nullptr, &dq0_2);
     (*atoms)[ia].grad = (dxdq0 * dq0_2[0]) +
                         (dxdq1 * dq0_2[1]) +
                         (dxdq2 * dq0_2[2]) +
@@ -448,7 +451,7 @@ void colvar::euler_psi::calc_gradients()
   rot_deriv_impl->prepare_derivative(rotation_derivative_dldq::use_dq);
   cvm::vector1d<cvm::rvector> dq0_2;
   for (size_t ia = 0; ia < atoms->size(); ia++) {
-    rot_deriv_impl->calc_derivative_wrt_group2(ia, nullptr, &dq0_2);
+    rot_deriv_impl->calc_derivative_wrt_group2<false, true, false>(ia, nullptr, &dq0_2);
     (*atoms)[ia].grad = (dxdq0 * dq0_2[0]) +
                         (dxdq1 * dq0_2[1]) +
                         (dxdq2 * dq0_2[2]) +
@@ -495,7 +498,7 @@ void colvar::euler_theta::calc_gradients()
   rot_deriv_impl->prepare_derivative(rotation_derivative_dldq::use_dq);
   cvm::vector1d<cvm::rvector> dq0_2;
   for (size_t ia = 0; ia < atoms->size(); ia++) {
-    rot_deriv_impl->calc_derivative_wrt_group2(ia, nullptr, &dq0_2);
+    rot_deriv_impl->calc_derivative_wrt_group2<false, true, false>(ia, nullptr, &dq0_2);
     (*atoms)[ia].grad = (dxdq0 * dq0_2[0]) +
                         (dxdq1 * dq0_2[1]) +
                         (dxdq2 * dq0_2[2]) +

lib/colvars/colvarcomp_torchann.cpp

@@ -0,0 +1,233 @@
+// -*- c++ -*-
+
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/Colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
+#include "colvar.h"
+#include "colvarcomp.h"
+#include "colvarmodule.h"
+#include "colvarparse.h"
+#include "colvarvalue.h"
+
+#include "colvarcomp_torchann.h"
+
+
+#ifdef COLVARS_TORCH
+
+colvar::torchANN::torchANN()
+{
+  set_function_type("torchANN");
+  provide(f_cvc_periodic);
+}
+
+colvar::torchANN::~torchANN() {}
+
+
+int colvar::torchANN::init(std::string const &conf) {
+
+  int error_code = linearCombination::init(conf);
+
+  std::string model_file ;
+  get_keyval(conf, "modelFile", model_file, std::string(""));
+  try {
+    nn = torch::jit::load(model_file);
+    nn.to(torch::kCPU);
+    cvm::log("torch model loaded.") ;
+  } catch (const std::exception & e) {
+    return cvm::error("Error: couldn't load libtorch model (see below).\n" + cvm::to_str(e.what()),
+                      COLVARS_INPUT_ERROR);
+  }
+
+  auto const legacy_keyword = get_keyval(conf, "m_output_index", m_output_index, m_output_index);
+  if (legacy_keyword) {
+    cvm::log("Warning: m_output_index is a deprecated keyword, please use output_component instead.\n");
+  }
+  get_keyval(conf, "output_component", m_output_index, m_output_index);
+
+  get_keyval(conf, "doubleInputTensor", use_double_input, use_double_input);
+  //get_keyval(conf, "useGPU", use_gpu, false);
+
+  cvc_indices.resize(cv.size(),0);
+
+  size_t num_inputs = 0;
+  // compute total number of inputs of neural network
+  for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv)
+  {
+    num_inputs += cv[i_cv]->value().size() ;
+    if (i_cv < cv.size() - 1)
+      cvc_indices[i_cv+1] = num_inputs;
+  }
+  cvm::log("Input dimension of model: " + cvm::to_str(num_inputs));
+
+  // initialize the input tensor
+  auto options = torch::TensorOptions().dtype(torch::kFloat32).requires_grad(true);
+
+  /*
+  if (use_gpu) {
+    if (torch::cuda::is_available()) {
+      try {
+        nn.to(torch::kCUDA);
+      } catch(const std::exception & e) {
+        cvm::error("Failed to move model to GPU.");
+        use_gpu = false;
+      }
+    } else {
+      use_gpu = false;
+      cvm::log("GPU not available.");
+    }
+  }
+
+  if (use_gpu) {
+    options = options.device(torch::kCUDA);
+    if (use_double_input) {
+      cvm::log("Data type reset to Float for GPU computation!");
+      use_double_input = false;
+    }
+  }
+  */
+
+  if (use_double_input) {  // set type to double
+    options = options.dtype(torch::kFloat64);
+    nn.to(torch::kFloat64);
+    cvm::log("Model's dtype: kFloat64.");
+  } else {
+    cvm::log("Model's dtype: kFloat32.");
+  }
+
+  input_tensor = torch::zeros({1,(long int) num_inputs}, options);
+
+  try { // test the model
+    std::vector<torch::jit::IValue> inputs={input_tensor};
+    nn_outputs = nn.forward(inputs).toTensor()[0][m_output_index];
+    cvm::log("Evaluating model with zero tensor succeeded.");
+  } catch (const std::exception & e) {
+    error_code |= cvm::error("Error: evaluating model with zero tensor failed (see below).\n" +
+                                 cvm::to_str(e.what()),
+                             COLVARS_INPUT_ERROR);
+  }
+
+  return error_code;
+}
+
+
+void colvar::torchANN::calc_value() {
+
+  for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv)
+    cv[i_cv]->calc_value();
+
+  /*
+  if (use_gpu)
+    input_tensor = input_tensor.to(torch::kCPU);
+  */
+
+  // set input tensor with no_grad
+  {
+    torch::NoGradGuard no_grad;
+    size_t l = 0;
+    for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
+      const colvarvalue& current_cv_value = cv[i_cv]->value();
+      if (current_cv_value.type() == colvarvalue::type_scalar) {
+        input_tensor[0][l++] = cv[i_cv]->sup_coeff * (cvm::pow(current_cv_value.real_value, cv[i_cv]->sup_np));
+      } else {
+        for (size_t j_elem = 0; j_elem < current_cv_value.size(); ++j_elem)
+          input_tensor[0][l++] = cv[i_cv]->sup_coeff * current_cv_value[j_elem];
+      }
+    }
+  }
+
+  /*
+  if (use_gpu)
+    input_tensor = input_tensor.to(torch::kCUDA);
+  */
+
+  std::vector<torch::jit::IValue> inputs={input_tensor};
+
+  // evaluate the value of function
+  nn_outputs = nn.forward(inputs).toTensor()[0][m_output_index];
+
+  input_grad = torch::autograd::grad({nn_outputs}, {input_tensor})[0][0];
+
+  /*
+  if (use_gpu)
+    input_grad = input_grad.to(torch::kCPU);
+  */
+
+  x = nn_outputs.item<double>() ;
+
+  this->wrap(x);
+
+}
+
+void colvar::torchANN::calc_gradients() {
+  for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
+    cv[i_cv]->calc_gradients();
+    if (cv[i_cv]->is_enabled(f_cvc_explicit_gradient)) {
+      const cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv);
+      // get the initial index of this cvc
+      size_t l = cvc_indices[i_cv];
+      for (size_t j_elem = 0; j_elem < cv[i_cv]->value().size(); ++j_elem) {
+        // get derivative of neural network wrt its input
+        const cvm::real factor = input_grad[l+j_elem].item<double>();
+        for (size_t k_ag = 0 ; k_ag < cv[i_cv]->atom_groups.size(); ++k_ag) {
+          for (size_t l_atom = 0; l_atom < (cv[i_cv]->atom_groups)[k_ag]->size(); ++l_atom) {
+            (*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad = factor_polynomial * factor * (*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad;
+          }
+        }
+      }
+    }
+  }
+}
+
+void colvar::torchANN::apply_force(colvarvalue const &force) {
+
+  for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
+    // If this CV uses explicit gradients, then atomic gradients is already calculated
+    // We can apply the force to atom groups directly
+    if (cv[i_cv]->is_enabled(f_cvc_explicit_gradient)) {
+      for (size_t k_ag = 0 ; k_ag < cv[i_cv]->atom_groups.size(); ++k_ag) {
+        (cv[i_cv]->atom_groups)[k_ag]->apply_colvar_force(force.real_value);
+      }
+    } else {
+      const colvarvalue& current_cv_value = cv[i_cv]->value();
+      colvarvalue cv_force(current_cv_value);
+      cv_force.reset();
+      const cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv);
+      // get the initial index of this cvc
+      size_t l = cvc_indices[i_cv];
+      for (size_t j_elem = 0; j_elem < current_cv_value.size(); ++j_elem) {
+        cv_force[j_elem] = factor_polynomial * input_grad[l+j_elem].item<double>() * force.real_value;
+      }
+      cv[i_cv]->apply_force(cv_force);
+    }
+  }
+}
+
+
+#else
+
+colvar::torchANN::torchANN()
+{
+  set_function_type("torchANN");
+}
+
+colvar::torchANN::~torchANN() {}
+
+int colvar::torchANN::init(std::string const &conf) {
+
+  return cvm::error(
+          "torchANN requires the libtorch library, but it is not enabled during compilation.\n"
+          "Please refer to the Compilation Notes section of the Colvars manual for more "
+          "information.\n",
+          COLVARS_NOT_IMPLEMENTED);
+
+}
+
+void colvar::torchANN::calc_value()
+{
+}
+
+#endif

lib/colvars/colvarcomp_torchann.h

@@ -0,0 +1,63 @@
+// -*- c++ -*-
+
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/Colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+//
+#ifndef COLVARCOMP_TORCH_H
+#define COLVARCOMP_TORCH_H
+
+// Declaration of torchann
+
+#include <memory>
+
+#include "colvar.h"
+#include "colvarcomp.h"
+#include "colvarmodule.h"
+
+#ifdef COLVARS_TORCH
+
+#include <torch/torch.h>
+#include <torch/script.h>
+
+class colvar::torchANN
+  : public colvar::linearCombination
+{
+protected:
+    torch::jit::script::Module nn;
+    /// the index of nn output component
+    size_t m_output_index = 0;
+    bool use_double_input = false;
+    //bool use_gpu;
+    // 1d tensor, concatenation of values of sub-cvcs
+    torch::Tensor input_tensor;
+    torch::Tensor nn_outputs;
+    torch::Tensor input_grad;
+    // record the initial index of of sub-cvcs in input_tensor
+    std::vector<int> cvc_indices;
+public:
+    torchANN();
+    virtual ~torchANN();
+    virtual int init(std::string const &conf);
+    virtual void calc_value();
+    virtual void calc_gradients();
+    virtual void apply_force(colvarvalue const &force);
+};
+
+#else
+
+class colvar::torchANN
+  : public colvar::cvc
+{
+public:
+    torchANN();
+    virtual ~torchANN();
+    virtual int init(std::string const &conf);
+    virtual void calc_value();
+};
+#endif // COLVARS_TORCH checking
+
+#endif

lib/colvars/colvardeps.cpp

@@ -92,6 +92,8 @@ void colvardeps::restore_children_deps() {
 
 void colvardeps::provide(int feature_id, bool truefalse) {
   feature_states[feature_id].available = truefalse;
+  // Make sure that we don't leave this feature enabled
+  if (!truefalse) disable(feature_id);
 }
 
 
@@ -123,8 +125,9 @@ bool colvardeps::get_keyval_feature(colvarparse *cvp,
 
 
 int colvardeps::enable(int feature_id,
-                       bool dry_run /* default: false */,
+                       bool dry_run  /* default: false */,
-                       bool toplevel /* default: true */)
+                       bool toplevel /* default: true */,
+                       bool error    /*default: false */)
 {
   int res;
   size_t i, j;
@@ -137,9 +140,12 @@ int colvardeps::enable(int feature_id,
   feature *f = features()[feature_id];
   feature_state *fs = &feature_states[feature_id];
 
+  // dry_run can be true because parent object is not active, yet we are displaying an error message
+  // then error is set to true
+
   if (cvm::debug()) {
     cvm::log("DEPS: " + description +
-      (dry_run ? " testing " : " enabling ") +
+      (dry_run ? " testing " : " enabling ") +  (error ? " [error] " : "") +
       "\"" + f->description +"\"\n");
   }
 
@@ -159,7 +165,7 @@ int colvardeps::enable(int feature_id,
     (is_dynamic(feature_id) ? "Dynamic" : "User-controlled");
 
   if (!fs->available) {
-    if (!dry_run) {
+    if (!dry_run || error) {
       if (toplevel) {
         cvm::error("Error: " + feature_type_descr + " feature unavailable: \""
           + f->description + "\" in " + description + ".\n");
@@ -172,7 +178,7 @@ int colvardeps::enable(int feature_id,
   }
 
   if (!toplevel && !is_dynamic(feature_id)) {
-    if (!dry_run) {
+    if (!dry_run || error) {
       cvm::log(feature_type_descr + " feature \"" + f->description
         + "\" cannot be enabled automatically in " + description + ".\n");
       if (is_user(feature_id)) {
@@ -189,7 +195,7 @@ int colvardeps::enable(int feature_id,
     if (cvm::debug())
       cvm::log(f->description + " requires exclude " + g->description + "\n");
     if (is_enabled(f->requires_exclude[i])) {
-      if (!dry_run) {
+      if (!dry_run || error) {
         cvm::log("Feature \"" + f->description + "\" is incompatible with \""
         + g->description + "\" in " + description + ".\n");
         if (toplevel) {
@@ -204,10 +210,14 @@ int colvardeps::enable(int feature_id,
   for (i=0; i<f->requires_self.size(); i++) {
     if (cvm::debug())
       cvm::log(f->description + " requires self " + features()[f->requires_self[i]]->description + "\n");
-    res = enable(f->requires_self[i], dry_run, false);
+    res = enable(f->requires_self[i], dry_run, false, error);
     if (res != COLVARS_OK) {
-      if (!dry_run) {
+      if (!dry_run || error) {
-        cvm::log("...required by \"" + f->description + "\" in " + description + "\n");
+        if (toplevel) {
+          cvm::log("Cannot enable \"" + f->description + "\" in " + description + "\n");
+        } else {
+          cvm::log("...required by \"" + f->description + "\" in " + description + "\n");
+        }
         if (toplevel) {
           cvm::error("Error: Failed dependency in " + description + ".\n");
         }
@@ -225,11 +235,11 @@ int colvardeps::enable(int feature_id,
       int g = f->requires_alt[i][j];
       if (cvm::debug())
         cvm::log(f->description + " requires alt " + features()[g]->description + "\n");
-      res = enable(g, true, false);  // see if available
+      res = enable(g, true, false, error);  // see if available
       if (res == COLVARS_OK) {
         ok = true;
-        if (!dry_run) {
+        if (!dry_run || error) {
-          enable(g, false, false); // Require again, for real
+          enable(g, false, false, error); // Require again, for real
           fs->alternate_refs.push_back(g); // We remember we enabled this
           // so we can free it if this feature gets disabled
         }
@@ -245,7 +255,7 @@ int colvardeps::enable(int feature_id,
         for (j=0; j<f->requires_alt[i].size(); j++) {
           int g = f->requires_alt[i][j];
           cvm::log(cvm::to_str(j+1) + ". " + features()[g]->description + "\n");
-          enable(g, false, false); // Just for printing error output
+          enable(g, false, false, true); // Just for printing error output
         }
         cvm::decrease_depth();
         cvm::log("-----------------------------------------\n");
@@ -264,10 +274,14 @@ int colvardeps::enable(int feature_id,
   for (i=0; i<f->requires_children.size(); i++) {
     int g = f->requires_children[i];
     for (j=0; j<children.size(); j++) {
-      res = children[j]->enable(g, dry_run || !is_enabled(), false);
+      res = children[j]->enable(g, dry_run || !is_enabled(), false, error);
       if (res != COLVARS_OK) {
-        if (!dry_run) {
+        if (!dry_run || error) {
-          cvm::log("...required by \"" + f->description + "\" in " + description + "\n");
+          if (toplevel) {
+            cvm::log("Cannot enable \"" + f->description + "\" in " + description + "\n");
+          } else {
+            cvm::log("...required by \"" + f->description + "\" in " + description + "\n");
+          }
           if (toplevel) {
             cvm::error("Error: Failed dependency in " + description + ".\n");
           }

lib/colvars/colvardeps.h

@@ -198,7 +198,9 @@ public:
   /// \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);
+  /// \param error Recursively enable, printing error messages along the way
+  /// Necessary when propagating errors across alternate dependencies
+  int enable(int f, bool dry_run = false, bool toplevel = true, bool error = false);
 
   /// Disable a feature, decrease the reference count of its dependencies
   /// and recursively disable them as applicable
@@ -255,6 +257,8 @@ public:
     f_cvb_scale_biasing_force,
     /// \brief whether this bias is applied to one or more ext-Lagrangian colvars
     f_cvb_extended,
+    /// Process this bias's data in parallel over multiple CPU threads
+    f_cvb_smp,
     f_cvb_ntot
   };
 
@@ -263,8 +267,11 @@ public:
     f_cv_active,
     /// \brief Colvar is awake (active on its own accord) this timestep
     f_cv_awake,
-    /// \brief Gradients are calculated and temporarily stored, so
+    /// \brief External force can be applied, either to atoms or to an
-    /// that external forces can be applied
+    /// extended DOF
+    f_cv_apply_force,
+    /// \brief Gradients are calculated and temporarily stored,
+    /// so that external forces can be propagated to atoms
     f_cv_gradient,
     /// \brief Collect atomic gradient data from all cvcs into vector
     /// atomic_gradient
@@ -277,7 +284,10 @@ public:
     /// forces on the inverse gradient
     f_cv_total_force,
     /// \brief Calculate total force from atomic forces
+    /// or get it from the back-end for an external parameter
     f_cv_total_force_calc,
+    /// \brief Total force is that of current time step
+    f_cv_total_force_current_step,
     /// \brief Subtract the applied force from the total force
     f_cv_subtract_applied_force,
     /// \brief Estimate Jacobian derivative
@@ -289,8 +299,10 @@ public:
     /// center with fictitious mass; bias forces will be applied to
     /// the center
     f_cv_extended_Lagrangian,
-    /// \brief An extended variable that sets an external variable in the
+    /// \brief A variable that constrains or follows an external parameter
-    /// back-end (eg. an alchemical coupling parameter for lambda-dynamics)
+    /// in the back-end (eg. an alchemical coupling parameter for lambda-dynamics)
+    /// If extended Lagrangian, then we drive the external parameter
+    /// Otherwise we follow it
     /// Can have a single component
     f_cv_external,
     /// \brief The extended system coordinate undergoes Langevin dynamics

lib/colvars/colvargrid.cpp

@@ -24,15 +24,14 @@ colvar_grid_count::colvar_grid_count()
   mult = 1;
 }
 
-colvar_grid_count::colvar_grid_count(std::vector<int> const &nx_i,
+colvar_grid_count::colvar_grid_count(std::vector<colvar *>  &colvars,
-                                     size_t const &def_count)
+                                     std::string config)
-  : colvar_grid<size_t>(nx_i, def_count, 1)
+  : colvar_grid<size_t>(colvars, 0, 1, false, nullptr, config)
 {}
 
 colvar_grid_count::colvar_grid_count(std::vector<colvar *>  &colvars,
-                                     size_t const &def_count,
+                                     std::shared_ptr<const colvar_grid_params> params)
-                                     bool margin)
+  : colvar_grid<size_t>(colvars, 0, 1, false, params)
-  : colvar_grid<size_t>(colvars, def_count, 1, margin)
 {}
 
 std::string colvar_grid_count::get_state_params() const
@@ -132,13 +131,17 @@ colvar_grid_scalar::colvar_grid_scalar(colvar_grid_scalar const &g)
 {
 }
 
-colvar_grid_scalar::colvar_grid_scalar(std::vector<int> const &nx_i)
+colvar_grid_scalar::colvar_grid_scalar(std::vector<colvar *> &colvars,
-  : colvar_grid<cvm::real>(nx_i, 0.0, 1), samples(NULL)
+                                       std::shared_ptr<const colvar_grid_params> params,
+                                       bool add_extra_bin,
+                                       std::string config)
+  : colvar_grid<cvm::real>(colvars, 0.0, 1, add_extra_bin, params, config), samples(NULL)
 {
 }
 
-colvar_grid_scalar::colvar_grid_scalar(std::vector<colvar *> &colvars, bool margin)
+colvar_grid_scalar::colvar_grid_scalar(std::string const &filename)
-  : colvar_grid<cvm::real>(colvars, 0.0, 1, margin), samples(NULL)
+  : colvar_grid<cvm::real>(filename, 1),
+    samples(nullptr)
 {
 }
 
@@ -330,89 +333,37 @@ cvm::real colvar_grid_scalar::grid_rmsd(colvar_grid_scalar const &other_grid) co
 
 
 colvar_grid_gradient::colvar_grid_gradient()
-  : colvar_grid<cvm::real>(), samples(NULL), full_samples(0), min_samples(0)
+  : colvar_grid<cvm::real>(), samples(NULL)
 {}
 
 
-colvar_grid_gradient::colvar_grid_gradient(std::vector<int> const &nx_i)
+// colvar_grid_gradient::colvar_grid_gradient(std::vector<colvar *> &colvars, std::string config)
-  : colvar_grid<cvm::real>(nx_i, 0.0, nx_i.size()), samples(NULL), full_samples(0), min_samples(0)
+//   : colvar_grid<cvm::real>(colvars, 0.0, colvars.size(), false, nullptr, config), samples(NULL)
-{}
+// {}
 
+// colvar_grid_gradient::colvar_grid_gradient(std::vector<colvar *> &colvars,
+//                                            std::shared_ptr<colvar_grid_count> samples_in)
+//   : colvar_grid<cvm::real>(colvars, 0.0, colvars.size(), false, samples_in), samples(samples_in)
+// {
+//   if (samples_in)
+//     samples_in->has_parent_data = true;
+// }
 
-colvar_grid_gradient::colvar_grid_gradient(std::vector<colvar *> &colvars)
+colvar_grid_gradient::colvar_grid_gradient(std::vector<colvar *> &colvars,
-  : colvar_grid<cvm::real>(colvars, 0.0, colvars.size()), samples(NULL), full_samples(0), min_samples(0)
+                                           std::shared_ptr<colvar_grid_count> samples_in,
-{}
+                                           std::shared_ptr<const colvar_grid_params> params,
-
+                                           std::string config)
-
+  : colvar_grid<cvm::real>(colvars, 0.0, colvars.size(), false, params, config), samples(samples_in)
-colvar_grid_gradient::colvar_grid_gradient(std::vector<colvar *> &colvars, std::shared_ptr<colvar_grid_count> samples_in)
-  : colvar_grid<cvm::real>(colvars, 0.0, colvars.size()), samples(samples_in), full_samples(0), min_samples(0)
 {
-  samples_in->has_parent_data = true;
+  if (samples_in)
+    samples_in->has_parent_data = true;
 }
 
 
-colvar_grid_gradient::colvar_grid_gradient(std::string &filename)
+colvar_grid_gradient::colvar_grid_gradient(std::string const &filename)
-  : colvar_grid<cvm::real>(),
+  : colvar_grid<cvm::real>(filename, 0),
-    samples(NULL)
+    samples(nullptr)
 {
-  std::istream &is = cvm::main()->proxy->input_stream(filename,
-                                                      "gradient file");
-  if (!is) {
-    return;
-  }
-
-  // Data in the header: nColvars, then for each
-  // xiMin, dXi, nPoints, periodic flag
-
-  std::string  hash;
-  size_t i;
-
-  if ( !(is >> hash) || (hash != "#") ) {
-    cvm::error("Error reading grid at position "+
-                cvm::to_str(static_cast<size_t>(is.tellg()))+
-                " in stream(read \"" + hash + "\")\n");
-    return;
-  }
-
-  is >> nd;
-
-  if (nd > 50) {
-    cvm::error("Error: excessive number of dimensions in file \""+
-               filename+"\".  Please ensure that the file is not corrupt.\n",
-               COLVARS_INPUT_ERROR);
-    return;
-  }
-
-  mult = nd;
-  std::vector<cvm::real> lower_in(nd), widths_in(nd);
-  std::vector<int>       nx_in(nd);
-  std::vector<int>       periodic_in(nd);
-
-  for (i = 0; i < nd; i++ ) {
-    if ( !(is >> hash) || (hash != "#") ) {
-      cvm::error("Error reading grid at position "+
-                  cvm::to_str(static_cast<size_t>(is.tellg()))+
-                  " in stream(read \"" + hash + "\")\n");
-      return;
-    }
-
-    is >> lower_in[i] >> widths_in[i] >> nx_in[i] >> periodic_in[i];
-  }
-
-  this->setup(nx_in, 0., mult);
-
-  widths = widths_in;
-
-  for (i = 0; i < nd; i++ ) {
-    lower_boundaries.push_back(colvarvalue(lower_in[i]));
-    periodic.push_back(static_cast<bool>(periodic_in[i]));
-  }
-
-  // Reset the istream for read_multicol, which expects the whole file
-  is.clear();
-  is.seekg(0);
-  read_multicol(is);
-  cvm::main()->proxy->close_input_stream(filename);
 }
 
 std::string colvar_grid_gradient::get_state_params() const
@@ -586,12 +537,13 @@ cvm::real colvar_grid_gradient::grid_rmsd(colvar_grid_gradient const &other_grid
 }
 
 
-integrate_potential::integrate_potential(std::vector<colvar *> &colvars, std::shared_ptr<colvar_grid_gradient> gradients)
+integrate_potential::integrate_potential(std::vector<colvar *> &colvars,
-  : colvar_grid_scalar(colvars, true),
+                                         std::shared_ptr<colvar_grid_gradient> gradients)
+  : colvar_grid_scalar(colvars, gradients, true),
     b_smoothed(false),
     gradients(gradients)
 {
-  // parent class colvar_grid_scalar is constructed with margin option set to true
+  // parent class colvar_grid_scalar is constructed with add_extra_bin option set to true
   // hence PMF grid is wider than gradient grid if non-PBC
 
   if (nd > 1) {

lib/colvars/colvargrid.h

@@ -19,17 +19,13 @@
 #include "colvarparse.h"
 
 
-/// \brief Grid of values of a function of several collective
+/// \brief Unified base class for grid of values of a function of several collective
-/// variables \param T The data type
+/// variables
-///
+class colvar_grid_params  {
-/// Only scalar colvars supported so far: vector colvars are treated as arrays
-template <class T> class colvar_grid : public colvarparse {
-
-  //protected:
-public: // TODO create accessors for these after all instantiations work
 
+public:
   /// Number of dimensions
-  size_t nd;
+  size_t nd = 0;
 
   /// Number of points along each dimension
   std::vector<int> nx;
@@ -37,6 +33,27 @@ public: // TODO create accessors for these after all instantiations work
   /// Cumulative number of points along each dimension
   std::vector<int> nxc;
 
+  /// Lower boundaries of the colvars in this grid
+  std::vector<colvarvalue>  lower_boundaries;
+
+  /// Upper boundaries of the colvars in this grid
+  std::vector<colvarvalue>  upper_boundaries;
+
+  /// Widths of the colvars in this grid
+  std::vector<cvm::real>    widths;
+};
+
+
+/// \brief Grid of values of a function of several collective
+/// variables \param T The data type
+///
+/// Only scalar colvars supported so far: vector colvars are treated as arrays
+/// All common, type-independent members are collected in the base class colvar_grid_base
+template <class T> class colvar_grid : public colvar_grid_params, public colvarparse {
+
+  //protected:
+public: // TODO create accessors for these after all instantiations work
+
   /// \brief Multiplicity of each datum (allow the binning of
   /// non-scalar types such as atomic gradients)
   size_t mult;
@@ -73,13 +90,6 @@ public: // TODO create accessors for these after all instantiations work
   }
 
 public:
-
-  /// Lower boundaries of the colvars in this grid
-  std::vector<colvarvalue> lower_boundaries;
-
-  /// Upper boundaries of the colvars in this grid
-  std::vector<colvarvalue> upper_boundaries;
-
   /// Whether some colvars are periodic
   std::vector<bool>        periodic;
 
@@ -89,9 +99,6 @@ public:
   /// Whether some colvars have hard upper boundaries
   std::vector<bool>        hard_upper_boundaries;
 
-  /// Widths of the colvars in this grid
-  std::vector<cvm::real>   widths;
-
   /// True if this is a count grid related to another grid of data
   bool has_parent_data;
 
@@ -218,19 +225,15 @@ public:
   /// \brief "Almost copy-constructor": only copies configuration
   /// parameters from another grid, but doesn't reallocate stuff;
   /// setup() must be called after that;
-  colvar_grid(colvar_grid<T> const &g) : colvarparse(),
+  colvar_grid(colvar_grid<T> const &g) : colvar_grid_params(colvar_grid_params(g)),
-                                         nd(g.nd),
+                                         colvarparse(),
-                                         nx(g.nx),
                                          mult(g.mult),
                                          data(),
                                          cv(g.cv),
                                          use_actual_value(g.use_actual_value),
-                                         lower_boundaries(g.lower_boundaries),
-                                         upper_boundaries(g.upper_boundaries),
                                          periodic(g.periodic),
                                          hard_lower_boundaries(g.hard_lower_boundaries),
                                          hard_upper_boundaries(g.hard_upper_boundaries),
-                                         widths(g.widths),
                                          has_parent_data(false),
                                          has_data(false)
   {}
@@ -247,22 +250,31 @@ public:
     this->setup(nx_i, t, mult_i);
   }
 
-  /// \brief Constructor from a vector of colvars
+  /// \brief Constructor from a vector of colvars or an optional grid config string
   /// \param add_extra_bin requests that non-periodic dimensions are extended
   /// by 1 bin to accommodate the integral (PMF) of another gridded quantity (gradient)
   colvar_grid(std::vector<colvar *> const &colvars,
               T const &t = T(),
               size_t mult_i = 1,
-              bool add_extra_bin = false)
+              bool add_extra_bin = false,
+              std::shared_ptr<const colvar_grid_params> params = nullptr,
+              std::string config = std::string())
     : has_parent_data(false), has_data(false)
   {
     (void) t;
-    this->init_from_colvars(colvars, mult_i, add_extra_bin);
+    this->init_from_colvars(colvars, mult_i, add_extra_bin, params, config);
   }
 
+  /// \brief Constructor from a multicol file
+  /// \param filename multicol file containing data to be read
+  /// \param multi_i multiplicity of the data - if 0, assume gradient multiplicity (mult = nd)
+  colvar_grid(std::string const &filename, size_t mult_i = 1);
+
   int init_from_colvars(std::vector<colvar *> const &colvars,
                         size_t mult_i = 1,
-                        bool add_extra_bin = false)
+                        bool add_extra_bin = false,
+                        std::shared_ptr<const colvar_grid_params> params = nullptr,
+                        std::string config = std::string())
   {
     if (cvm::debug()) {
       cvm::log("Reading grid configuration from collective variables.\n");
@@ -279,8 +291,7 @@ public:
                " collective variables, multiplicity = "+cvm::to_str(mult_i)+".\n");
     }
 
-    for (i =  0; i < cv.size(); i++) {
+    for (i =  0; i < nd; i++) {
-
       if (cv[i]->value().type() != colvarvalue::type_scalar) {
         cvm::error("Colvar grids can only be automatically "
                    "constructed for scalar variables.  "
@@ -298,7 +309,6 @@ public:
       widths.push_back(cv[i]->width);
       hard_lower_boundaries.push_back(cv[i]->is_enabled(colvardeps::f_cv_hard_lower_boundary));
       hard_upper_boundaries.push_back(cv[i]->is_enabled(colvardeps::f_cv_hard_upper_boundary));
-      periodic.push_back(cv[i]->periodic_boundaries());
 
       // By default, get reported colvar value (for extended Lagrangian colvars)
       use_actual_value.push_back(false);
@@ -310,22 +320,55 @@ public:
         use_actual_value[i-1] = true;
       }
 
+      // This needs to work if the boundaries are undefined in the colvars
+      lower_boundaries.push_back(cv[i]->lower_boundary);
+      upper_boundaries.push_back(cv[i]->upper_boundary);
+    }
+
+    // Replace widths and boundaries with optional custom configuration
+    if (!config.empty()) {
+      this->parse_params(config);
+      this->check_keywords(config, "grid");
+
+      if (params) {
+        cvm::error("Error: init_from_colvars was passed both a grid config and a template grid.", COLVARS_BUG_ERROR);
+        return COLVARS_BUG_ERROR;
+      }
+    } else if (params) {
+      // Match grid sizes with template
+
+      if (params->nd != nd) {
+        cvm::error("Trying to initialize grid from template with wrong dimension (" +
+                    cvm::to_str(params->nd) + " instead of " +
+                    cvm::to_str(this->nd) + ").");
+        return COLVARS_ERROR;
+      }
+
+      widths =params->widths;
+      lower_boundaries =params->lower_boundaries;
+      upper_boundaries =params->upper_boundaries;
+    }
+
+    // Only now can we determine periodicity
+    for (i =  0; i < nd; i++) {
+      periodic.push_back(cv[i]->periodic_boundaries(lower_boundaries[i].real_value,
+                                                    upper_boundaries[i].real_value));
+
       if (add_extra_bin) {
+        // Shift the grid by half the bin width (values at edges instead of center of bins)
+        lower_boundaries[i] -= 0.5 * widths[i];
+
         if (periodic[i]) {
-          // Shift the grid by half the bin width (values at edges instead of center of bins)
+          // Just shift
-          lower_boundaries.push_back(cv[i]->lower_boundary.real_value - 0.5 * widths[i]);
+          upper_boundaries[i] -= 0.5 * widths[i];
-          upper_boundaries.push_back(cv[i]->upper_boundary.real_value - 0.5 * widths[i]);
         } else {
-          // Make this grid larger by one bin width
+          // Widen grid by one bin width
-          lower_boundaries.push_back(cv[i]->lower_boundary.real_value - 0.5 * widths[i]);
+          upper_boundaries[i] += 0.5 * widths[i];
-          upper_boundaries.push_back(cv[i]->upper_boundary.real_value + 0.5 * widths[i]);
         }
-      } else {
-        lower_boundaries.push_back(cv[i]->lower_boundary);
-        upper_boundaries.push_back(cv[i]->upper_boundary);
       }
     }
 
+    // Reset grid sizes based on widths and boundaries
     this->init_from_boundaries();
     return this->setup();
   }
@@ -966,14 +1009,12 @@ public:
   virtual ~colvar_grid_count()
   {}
 
-  /// Constructor
+  /// Constructor from a vector of colvars or a config string
-  colvar_grid_count(std::vector<int> const &nx_i,
-                    size_t const           &def_count = 0);
-
-  /// Constructor from a vector of colvars
   colvar_grid_count(std::vector<colvar *>  &colvars,
-                    size_t const           &def_count = 0,
+                    std::shared_ptr<const colvar_grid_params> params = nullptr);
-                    bool                   add_extra_bin = false);
+
+  colvar_grid_count(std::vector<colvar *>  &colvars,
+                    std::string            config);
 
   /// Increment the counter at given position
   inline void incr_count(std::vector<int> const &ix)
@@ -1255,12 +1296,14 @@ public:
   /// Destructor
   virtual ~colvar_grid_scalar();
 
-  /// Constructor from specific sizes arrays
-  colvar_grid_scalar(std::vector<int> const &nx_i);
-
   /// Constructor from a vector of colvars
   colvar_grid_scalar(std::vector<colvar *> &colvars,
-                     bool add_extra_bin = false);
+                     std::shared_ptr<const colvar_grid_params> params = nullptr,
+                     bool add_extra_bin = false,
+                     std::string config = std::string());
+
+  /// Constructor from a multicol file
+  colvar_grid_scalar(std::string const &filename);
 
   /// Accumulate the value
   inline void acc_value(std::vector<int> const &ix,
@@ -1334,8 +1377,8 @@ public:
 
   /// \brief Return the gradient of the scalar field from finite differences
   /// Input coordinates are those of gradient grid, shifted wrt scalar grid
-  /// Should not be called on edges of scalar grid, provided the latter has margins
+  /// Should not be called on edges of scalar grid, provided the latter has
-  /// wrt gradient grid
+  /// margins (extra bins) wrt gradient grid
   inline void vector_gradient_finite_diff( const std::vector<int> &ix0, std::vector<cvm::real> &grad)
   {
     cvm::real A0, A1;
@@ -1566,17 +1609,21 @@ public:
   virtual ~colvar_grid_gradient()
   {}
 
-  /// Constructor from specific sizes arrays
+  // /// Constructor from specific sizes arrays
-  colvar_grid_gradient(std::vector<int> const &nx_i);
+  // colvar_grid_gradient(std::vector<int> const &nx_i);
 
-  /// Constructor from a vector of colvars
+  // /// Constructor from a vector of colvars
-  colvar_grid_gradient(std::vector<colvar *>  &colvars);
+  // colvar_grid_gradient(std::vector<colvar *>  &colvars,
+  //                      std::string config = std::string());
 
   /// Constructor from a multicol file
-  colvar_grid_gradient(std::string &filename);
+  colvar_grid_gradient(std::string const &filename);
 
   /// Constructor from a vector of colvars and a pointer to the count grid
-  colvar_grid_gradient(std::vector<colvar *> &colvars, std::shared_ptr<colvar_grid_count> samples_in);
+  colvar_grid_gradient(std::vector<colvar *> &colvars,
+                       std::shared_ptr<colvar_grid_count> samples_in = nullptr,
+                       std::shared_ptr<const colvar_grid_params> params = nullptr,
+                       std::string config = std::string());
 
   /// Parameters for smoothing data with low sampling
   int full_samples;
@@ -1829,7 +1876,8 @@ class integrate_potential : public colvar_grid_scalar
   {}
 
   /// Constructor from a vector of colvars + gradient grid
-  integrate_potential(std::vector<colvar *> &colvars, std::shared_ptr<colvar_grid_gradient> gradients);
+  integrate_potential(std::vector<colvar *> &colvars,
+                      std::shared_ptr<colvar_grid_gradient> gradients);
 
   /// Constructor from a gradient grid (for processing grid files without a Colvars config)
   integrate_potential(std::shared_ptr<colvar_grid_gradient> gradients);

lib/colvars/colvargrid_def.h

@@ -22,6 +22,62 @@
 #include "colvars_memstream.h"
 
 
+template <class T>
+colvar_grid<T>::colvar_grid(std::string const &filename, size_t mult_i)
+{
+std::istream &is = cvm::main()->proxy->input_stream(filename, "multicol grid file");
+if (!is) {
+  return;
+}
+
+// Data in the header: nColvars, then for each
+// xiMin, dXi, nPoints, periodic flag
+
+std::string  hash;
+size_t i;
+
+if ( !(is >> hash) || (hash != "#") ) {
+  cvm::error("Error reading grid at position "+
+              cvm::to_str(static_cast<size_t>(is.tellg()))+
+              " in stream(read \"" + hash + "\")\n");
+  return;
+}
+
+is >> nd;
+mult = (mult_i == 0) ? nd : mult_i;
+
+std::vector<cvm::real> lower_in(nd), widths_in(nd);
+std::vector<int>       nx_in(nd);
+std::vector<int>       periodic_in(nd);
+
+for (i = 0; i < nd; i++ ) {
+  if ( !(is >> hash) || (hash != "#") ) {
+    cvm::error("Error reading grid at position "+
+                cvm::to_str(static_cast<size_t>(is.tellg()))+
+                " in stream(read \"" + hash + "\")\n");
+    return;
+  }
+
+  is >> lower_in[i] >> widths_in[i] >> nx_in[i] >> periodic_in[i];
+}
+
+this->setup(nx_in, 0., mult);
+
+widths = widths_in;
+
+for (i = 0; i < nd; i++ ) {
+  lower_boundaries.push_back(colvarvalue(lower_in[i]));
+  periodic.push_back(static_cast<bool>(periodic_in[i]));
+}
+
+// Reset the istream for read_multicol, which expects the whole file
+is.clear();
+is.seekg(0);
+read_multicol(is);
+cvm::main()->proxy->close_input_stream(filename);
+}
+
+
 template <class T, class IST> IST &read_restart_template_(colvar_grid<T> &g, IST &is)
 {
   auto const start_pos = is.tellg();
@@ -203,14 +259,16 @@ template <class T> int colvar_grid<T>::parse_params(std::string const &conf,
                           lower_boundaries, lower_boundaries, colvarparse::parse_silent);
   colvarparse::get_keyval(conf, "upper_boundaries",
                           upper_boundaries, upper_boundaries, colvarparse::parse_silent);
+  // plural form is used in state file
+  colvarparse::get_keyval(conf, "widths", widths, widths, colvarparse::parse_silent);
 
   // camel case keywords are used in config file
-  colvarparse::get_keyval(conf, "lowerBoundaries",
+  colvarparse::get_keyval(conf, "lowerBoundary",
                           lower_boundaries, lower_boundaries, parse_mode);
-  colvarparse::get_keyval(conf, "upperBoundaries",
+  colvarparse::get_keyval(conf, "upperBoundary",
                           upper_boundaries, upper_boundaries, parse_mode);
 
-  colvarparse::get_keyval(conf, "widths", widths, widths, parse_mode);
+  colvarparse::get_keyval(conf, "width", widths, widths, parse_mode);
 
   // only used in state file
   colvarparse::get_keyval(conf, "sizes", nx, nx, colvarparse::parse_silent);

lib/colvars/colvarmodule.cpp

@@ -24,6 +24,7 @@
 #include "colvarbias_histogram_reweight_amd.h"
 #include "colvarbias_meta.h"
 #include "colvarbias_restraint.h"
+#include "colvarbias_opes.h"
 #include "colvarscript.h"
 #include "colvaratoms.h"
 #include "colvarcomp.h"
@@ -109,23 +110,23 @@ colvarmodule::colvarmodule(colvarproxy *proxy_in)
            "  https://doi.org/10.1080/00268976.2013.813594\n"
            "as well as all other papers listed below for individual features used.\n");
 
-#if (__cplusplus >= 201103L)
-  cvm::log("This version was built with the C++11 standard or higher.\n");
-#else
-  cvm::log("This version was built without the C++11 standard: some features are disabled.\n"
-    "Please see the following link for details:\n"
-    "  https://colvars.github.io/README-c++11.html\n");
-#endif
-
   cvm::log("Summary of compile-time features available in this build:\n");
 
-  if (proxy->check_smp_enabled() == COLVARS_NOT_IMPLEMENTED) {
+  std::string cxx_lang_msg("  - C++ language version: " + cvm::to_str(__cplusplus));
-    cvm::log("  - SMP parallelism: not available\n");
+#if defined(_WIN32) && !defined(__CYGWIN__)
+  cxx_lang_msg += std::string(" (warning: may not be accurate for this build)");
+#endif
+  cxx_lang_msg += std::string("\n");
+  cvm::log(cxx_lang_msg);
+
+  if (proxy->check_replicas_enabled() == COLVARS_NOT_IMPLEMENTED) {
+    cvm::log("  - Multiple replicas: not available\n");
   } else {
-    if (proxy->check_smp_enabled() == COLVARS_OK) {
+    if (proxy->check_replicas_enabled() == COLVARS_OK) {
-      cvm::log("  - SMP parallelism: enabled (num. threads = " + to_str(proxy->smp_num_threads()) + ")\n");
+      cvm::log("  - Multiple replicas: enabled (replica number " +
+               to_str(proxy->replica_index() + 1) + " of " + to_str(proxy->num_replicas()) + ")\n");
     } else {
-      cvm::log("  - SMP parallelism: available, but not enabled\n");
+      cvm::log("  - Multiple replicas: available, but not (yet) enabled\n");
     }
   }
 
@@ -201,6 +202,20 @@ std::vector<int> *colvarmodule::variables_active_smp_items()
 }
 
 
+int colvarmodule::calc_component_smp(int i)
+{
+  colvar *x = (*(variables_active_smp()))[i];
+  int x_item = (*(variables_active_smp_items()))[i];
+  if (cvm::debug()) {
+    cvm::log("Thread "+cvm::to_str(proxy->smp_thread_id())+"/"+
+             cvm::to_str(proxy->smp_num_threads())+
+             ": calc_component_smp(), i = "+cvm::to_str(i)+", cv = "+
+             x->name+", cvc = "+cvm::to_str(x_item)+"\n");
+  }
+  return x->calc_cvcs(x_item, 1);
+}
+
+
 std::vector<colvarbias *> *colvarmodule::biases_active()
 {
   return &(biases_active_);
@@ -387,8 +402,26 @@ int colvarmodule::parse_global_params(std::string const &conf)
     }
   }
 
-  if (parse->get_keyval(conf, "smp", proxy->b_smp_active, proxy->b_smp_active)) {
+  std::string smp;
-    if (proxy->b_smp_active == false) {
+  if (parse->get_keyval(conf, "smp", smp, "cvcs")) {
+    if (smp == "cvcs" || smp == "on" || smp == "yes") {
+      if (proxy->set_smp_mode(colvarproxy_smp::smp_mode_t::cvcs) != COLVARS_OK) {
+        cvm::error("Colvars component-based parallelism is not implemented.\n");
+        return COLVARS_INPUT_ERROR;
+      } else {
+        cvm::log("SMP parallelism will be applied to Colvars components.\n");
+        cvm::log("  - SMP parallelism: enabled (num. threads = " + to_str(proxy->smp_num_threads()) + ")\n");
+      }
+    } else if (smp == "inner_loop") {
+      if (proxy->set_smp_mode(colvarproxy_smp::smp_mode_t::inner_loop) != COLVARS_OK) {
+        cvm::error("SMP parallelism inside the calculation of Colvars components is not implemented.\n");
+        return COLVARS_INPUT_ERROR;
+      } else {
+        cvm::log("SMP parallelism will be applied inside the Colvars components.\n");
+      cvm::log("  - SMP parallelism: enabled (num. threads = " + to_str(proxy->smp_num_threads()) + ")\n");
+      }
+    } else {
+      proxy->set_smp_mode(colvarproxy_smp::smp_mode_t::none);
       cvm::log("SMP parallelism has been disabled.\n");
     }
   }
@@ -589,6 +622,9 @@ int colvarmodule::parse_biases(std::string const &conf)
   /// initialize reweightaMD instances
   parse_biases_type<colvarbias_reweightaMD>(conf, "reweightaMD");
 
+  /// initialize OPES instances
+  parse_biases_type<colvarbias_opes>(conf, "opes_metad");
+
   if (use_scripted_forces) {
     cvm::log(cvm::line_marker);
     cvm::increase_depth();
@@ -922,7 +958,7 @@ int colvarmodule::calc_colvars()
   }
 
   // if SMP support is available, split up the work
-  if (proxy->check_smp_enabled() == COLVARS_OK) {
+  if (proxy->get_smp_mode() == colvarproxy_smp::smp_mode_t::cvcs) {
 
     // first, calculate how much work (currently, how many active CVCs) each colvar has
 
@@ -948,8 +984,10 @@ int colvarmodule::calc_colvars()
     }
     cvm::decrease_depth();
 
-    // calculate colvar components in parallel
+    // calculate active colvar components in parallel
-    error_code |= proxy->smp_colvars_loop();
+    error_code |= proxy->smp_loop(variables_active_smp()->size(), [](int i) {
+        return cvm::main()->calc_component_smp(i);
+      });
 
     cvm::increase_depth();
     for (cvi = variables_active()->begin(); cvi != variables_active()->end(); cvi++) {
@@ -1013,7 +1051,7 @@ int colvarmodule::calc_biases()
   }
 
   // If SMP support is available, split up the work (unless biases need to use main thread's memory)
-  if (proxy->check_smp_enabled() == COLVARS_OK && !biases_need_main_thread) {
+  if (proxy->get_smp_mode() == colvarproxy::smp_mode_t::cvcs && !biases_need_main_thread) {
 
     if (use_scripted_forces && !scripting_after_biases) {
       // calculate biases and scripted forces in parallel
@@ -1097,7 +1135,7 @@ int colvarmodule::update_colvar_forces()
     cvm::log("Communicating forces from the colvars to the atoms.\n");
   cvm::increase_depth();
   for (cvi = variables_active()->begin(); cvi != variables_active()->end(); cvi++) {
-    if ((*cvi)->is_enabled(colvardeps::f_cv_gradient)) {
+    if ((*cvi)->is_enabled(colvardeps::f_cv_apply_force)) {
       (*cvi)->communicate_forces();
       if (cvm::get_error()) {
         return COLVARS_ERROR;
@@ -1986,7 +2024,7 @@ size_t & colvarmodule::depth()
 {
   // NOTE: do not call log() or error() here, to avoid recursion
   colvarmodule *cv = cvm::main();
-  if (proxy->check_smp_enabled() == COLVARS_OK) {
+  if (proxy->get_smp_mode() == colvarproxy::smp_mode_t::cvcs) {
     int const nt = proxy->smp_num_threads();
     if (int(cv->depth_v.size()) != nt) {
       proxy->smp_lock();

lib/colvars/colvarmodule.h

@@ -18,6 +18,11 @@
 #define COLVARS_DEBUG false
 #endif
 
+#if defined(__FAST_MATH__)
+// NOTE: This is used for fixing https://github.com/Colvars/colvars/issues/767
+#define COLVARS_BOUNDED_INV_TRIGONOMETRIC_FUNC
+#endif
+
 /*! \mainpage Main page
 This is the Developer's documentation for the Collective Variables module (Colvars).
 
@@ -147,17 +152,44 @@ public:
     return ::cos(static_cast<double>(x));
   }
 
-  /// Reimplemented to work around MS compiler issues
+#ifndef PI
-  static inline real asin(real const &x)
+#define PI   3.14159265358979323846
-  {
+#endif
-    return ::asin(static_cast<double>(x));
+#ifndef PI_2
-  }
+#define PI_2 1.57079632679489661923
+#endif
 
-  /// Reimplemented to work around MS compiler issues
+/// Reimplemented to work around compiler issues; return hard-coded values for boundary conditions
-  static inline real acos(real const &x)
+static inline real asin(real const &x)
-  {
+{
+#ifdef COLVARS_BOUNDED_INV_TRIGONOMETRIC_FUNC
+    if (x <= -1.0) {
+        return -PI_2;
+    } else if (x >= 1.0) {
+        return PI_2;
+    } else {
+        return ::asin(static_cast<double>(x));
+    }
+#else
+    return ::asin(static_cast<double>(x));
+#endif
+}
+
+/// Reimplemented to work around compiler issues; return hard-coded values for boundary conditions
+static inline real acos(real const &x)
+{
+#ifdef COLVARS_BOUNDED_INV_TRIGONOMETRIC_FUNC
+    if (x <= -1.0) {
+        return PI;
+    } else if (x >= 1.0) {
+        return 0.0;
+    } else {
+        return ::acos(static_cast<double>(x));
+    }
+#else
     return ::acos(static_cast<double>(x));
-  }
+#endif
+}
 
   /// Reimplemented to work around MS compiler issues
   static inline real atan2(real const &x, real const &y)
@@ -307,6 +339,9 @@ public:
   /// Indexes of the items to calculate for each colvar
   std::vector<int> *variables_active_smp_items();
 
+  /// Calculate the value of the specified component (to be called in a SMP loop)
+  int calc_component_smp(int i);
+
   /// Array of collective variable biases
   std::vector<colvarbias *> biases;
 

lib/colvars/colvarmodule_refs.h

@@ -129,6 +129,23 @@
     "  url = {https://doi.org/10.1002/jcc.26075}\n"
     "}\n";
 
+  paper_count_[std::string("Fiorin2024")] = 0;
+  paper_url_[std::string("Fiorin2024")] = "https://doi.org/10.1021/acs.jpcb.4c05604";
+  paper_bibtex_[std::string("Fiorin2024")] =
+    "\n"
+    "@article{Fiorin2024,\n"
+    "  author = {Fiorin, Giacomo and Marinelli, Fabrizio and Forrest, Lucy R. and Chen, Haochuan and Chipot, Christophe and Kohlmeyer, Axel and Santuz, Hubert and H{\\'e}nin, J{\\'e}rôme},\n"
+    "  title = {Expanded Functionality and Portability for the Colvars Library},\n"
+    "  journal = {J. Phys. Chem. {B}},\n"
+    "  volume = {128},\n"
+    "  number = {45},\n"
+    "  pages = {11108--11123},\n"
+    "  year = {2024},\n"
+    "  doi = {10.1021/acs.jpcb.4c05604},\n"
+    "  pmid = 39501453,\n"
+    "  url = { https://doi.org/10.1021/acs.jpcb.4c05604}\n"
+    "}\n";
+
   paper_count_[std::string("Fu2016")] = 0;
   paper_url_[std::string("Fu2016")] = "https://doi.org/10.1021/acs.jctc.6b00447";
   paper_bibtex_[std::string("Fu2016")] =
@@ -227,6 +244,20 @@
     "  url = {https://doi.org/10.1016/0263-7855(96)00018-5}\n"
     "}\n";
 
+  paper_count_[std::string("Lagardere2023")] = 0;
+  paper_url_[std::string("Lagardere2023")] = "https://arxiv.org/abs/2307.08006";
+  paper_bibtex_[std::string("Lagardere2023")] =
+    "\n"
+    "@misc{Lagardere2023,\n"
+    "      title={Lambda-ABF: Simplified, Accurate and Cost-effective Alchemical Free Energy Computations},\n"
+    "      author={Louis Lagard\\`ere and Lise Maurin and Olivier Adjoua and Krystel El Hage and Pierre Monmarch\\'e and Jean-Philip Piquemal and J\\'er\\^ome H\\'enin},\n"
+    "      year={2023},\n"
+    "      eprint={2307.08006},\n"
+    "      archivePrefix={arXiv},\n"
+    "      primaryClass={physics.chem-ph},\n"
+    "      url = {https://arxiv.org/abs/2307.08006}\n"
+    "}\n";
+
   paper_count_[std::string("Lesage2017")] = 0;
   paper_url_[std::string("Lesage2017")] = "https://doi.org/10.1021/acs.jpcb.6b10055";
   paper_bibtex_[std::string("Lesage2017")] =
@@ -344,6 +375,45 @@
     "  url = {https://doi.org/10.1021/ct500320c}\n"
     "}\n";
 
+  paper_count_[std::string("Invernizzi2020")] = 0;
+  paper_url_[std::string("Invernizzi2020")] = "https://pubs.acs.org/doi/10.1021/acs.jpclett.0c00497";
+  paper_bibtex_[std::string("Invernizzi2020")] =
+    "\n"
+    "@article{Invernizzi2020,\n"
+    "	title = {Rethinking {Metadynamics}: {From} {Bias} {Potentials} to {Probability} {Distributions}},\n"
+    "	volume = {11},\n"
+    "	issn = {1948-7185, 1948-7185},\n"
+    "	shorttitle = {Rethinking {Metadynamics}},\n"
+    "	url = {https://pubs.acs.org/doi/10.1021/acs.jpclett.0c00497},\n"
+    "	doi = {10.1021/acs.jpclett.0c00497},\n"
+    "	number = {7},\n"
+    "	urldate = {2020-09-30},\n"
+    "	journal = {J. Phys. Chem. Lett.},\n"
+    "	author = {Invernizzi, Michele and Parrinello, Michele},\n"
+    "	month = apr,\n"
+    "	year = {2020},\n"
+    "	pages = {2731--2736},\n"
+    "}\n";
+
+  paper_count_[std::string("Invernizzi2022")] = 0;
+  paper_url_[std::string("Invernizzi2022")] = "https://doi.org/10.1021/acs.jctc.2c00152";
+  paper_bibtex_[std::string("Invernizzi2022")] =
+    "\n"
+    "@article{Invernizzi2022,\n"
+    "	title = {Exploration vs {Convergence} {Speed} in {Adaptive}-{Bias} {Enhanced} {Sampling}},\n"
+    "	volume = {18},\n"
+    "	issn = {1549-9618},\n"
+    "	url = {https://doi.org/10.1021/acs.jctc.2c00152},\n"
+    "	doi = {10.1021/acs.jctc.2c00152},\n"
+    "	number = {6},\n"
+    "	urldate = {2024-07-02},\n"
+    "	journal = {J. Chem. Theory Comput.},\n"
+    "	author = {Invernizzi, Michele and Parrinello, Michele},\n"
+    "	month = jun,\n"
+    "	year = {2022},\n"
+    "	pages = {3988--3996},\n"
+    "}\n";
+
   paper_count_[std::string("n/a")] = 0;
   paper_url_[std::string("n/a")] = "";
   paper_bibtex_[std::string("n/a")] = "";
@@ -489,6 +559,42 @@
   feature_count_[std::string("Multi-Map collective variables")] = 0;
   feature_paper_map_[std::string("Multi-Map collective variables")] = "Fiorin2020";
 
+  feature_count_[std::string("Colvars-GROMACS interface")] = 0;
+  feature_paper_map_[std::string("Colvars-GROMACS interface")] = "Fiorin2024";
+
+  feature_count_[std::string("gspath colvar component")] = 0;
+  feature_paper_map_[std::string("gspath colvar component")] = "Fiorin2024";
+
+  feature_count_[std::string("gzpath colvar component")] = 0;
+  feature_paper_map_[std::string("gzpath colvar component")] = "Fiorin2024";
+
+  feature_count_[std::string("linearCombination colvar component")] = 0;
+  feature_paper_map_[std::string("linearCombination colvar component")] = "Fiorin2024";
+
+  feature_count_[std::string("gspathCV colvar component")] = 0;
+  feature_paper_map_[std::string("gspathCV colvar component")] = "Fiorin2024";
+
+  feature_count_[std::string("gzpathCV colvar component")] = 0;
+  feature_paper_map_[std::string("gzpathCV colvar component")] = "Fiorin2024";
+
+  feature_count_[std::string("aspathCV colvar component")] = 0;
+  feature_paper_map_[std::string("aspathCV colvar component")] = "Fiorin2024";
+
+  feature_count_[std::string("azpathCV colvar component")] = 0;
+  feature_paper_map_[std::string("azpathCV colvar component")] = "Fiorin2024";
+
+  feature_count_[std::string("Custom functions (Lepton)")] = 0;
+  feature_paper_map_[std::string("Custom functions (Lepton)")] = "Fiorin2024";
+
+  feature_count_[std::string("Scripted functions (Tcl)")] = 0;
+  feature_paper_map_[std::string("Scripted functions (Tcl)")] = "Fiorin2024";
+
+  feature_count_[std::string("ABMD bias")] = 0;
+  feature_paper_map_[std::string("ABMD bias")] = "Fiorin2024";
+
+  feature_count_[std::string("Updated multiple-walker ABF implementation")] = 0;
+  feature_paper_map_[std::string("Updated multiple-walker ABF implementation")] = "Fiorin2024";
+
   feature_count_[std::string("Umbrella-integration eABF estimator")] = 0;
   feature_paper_map_[std::string("Umbrella-integration eABF estimator")] = "Fu2016";
 
@@ -525,6 +631,15 @@
   feature_count_[std::string("VMD engine")] = 0;
   feature_paper_map_[std::string("VMD engine")] = "Humphrey1996";
 
+  feature_count_[std::string("alchLambda colvar component")] = 0;
+  feature_paper_map_[std::string("alchLambda colvar component")] = "Lagardere2023";
+
+  feature_count_[std::string("alchFLambda colvar component")] = 0;
+  feature_paper_map_[std::string("alchFLambda colvar component")] = "Lagardere2023";
+
+  feature_count_[std::string("Tinker-HP interface")] = 0;
+  feature_paper_map_[std::string("Tinker-HP interface")] = "Lagardere2023";
+
   feature_count_[std::string("eABF implementation")] = 0;
   feature_paper_map_[std::string("eABF implementation")] = "Lesage2017";
 
@@ -555,38 +670,14 @@
   feature_count_[std::string("ALB colvar bias implementation")] = 0;
   feature_paper_map_[std::string("ALB colvar bias implementation")] = "White2014";
 
-  feature_count_[std::string("Colvars-GROMACS interface")] = 0;
+  feature_count_[std::string("OPES")] = 0;
-  feature_paper_map_[std::string("Colvars-GROMACS interface")] = "n/a";
+  feature_paper_map_[std::string("OPES")] = "Invernizzi2020";
 
-  feature_count_[std::string("gspath colvar component")] = 0;
+  feature_count_[std::string("OPES explore or adaptive kernels")] = 0;
-  feature_paper_map_[std::string("gspath colvar component")] = "n/a";
+  feature_paper_map_[std::string("OPES explore or adaptive kernels")] = "Invernizzi2022";
-
-  feature_count_[std::string("gzpath colvar component")] = 0;
-  feature_paper_map_[std::string("gzpath colvar component")] = "n/a";
-
-  feature_count_[std::string("linearCombination colvar component")] = 0;
-  feature_paper_map_[std::string("linearCombination colvar component")] = "n/a";
-
-  feature_count_[std::string("gspathCV colvar component")] = 0;
-  feature_paper_map_[std::string("gspathCV colvar component")] = "n/a";
-
-  feature_count_[std::string("gzpathCV colvar component")] = 0;
-  feature_paper_map_[std::string("gzpathCV colvar component")] = "n/a";
-
-  feature_count_[std::string("aspathCV colvar component")] = 0;
-  feature_paper_map_[std::string("aspathCV colvar component")] = "n/a";
-
-  feature_count_[std::string("azpathCV colvar component")] = 0;
-  feature_paper_map_[std::string("azpathCV colvar component")] = "n/a";
 
   feature_count_[std::string("coordNum pairlist")] = 0;
   feature_paper_map_[std::string("coordNum pairlist")] = "n/a";
 
-  feature_count_[std::string("Custom functions (Lepton)")] = 0;
+  feature_count_[std::string("torchANN colvar component")] = 0;
-  feature_paper_map_[std::string("Custom functions (Lepton)")] = "n/a";
+  feature_paper_map_[std::string("torchANN colvar component")] = "n/a";
-
-  feature_count_[std::string("Scripted functions (Tcl)")] = 0;
-  feature_paper_map_[std::string("Scripted functions (Tcl)")] = "n/a";
-
-  feature_count_[std::string("ABMD bias")] = 0;
-  feature_paper_map_[std::string("ABMD bias")] = "n/a";

lib/colvars/colvarparse.cpp

@@ -592,7 +592,7 @@ int colvarparse::check_keywords(std::string &conf, char const *key)
 {
   if (cvm::debug())
     cvm::log("Configuration string for \""+std::string(key)+
-             "\": \"\n"+conf+"\".\n");
+             "\":\n\""+conf+"\".\n");
 
   strip_values(conf);
   // after stripping, the config string has either empty lines, or
@@ -833,7 +833,8 @@ bool colvarparse::key_lookup(std::string const &conf,
                                        data_end) + 1;
     }
 
-    if (data != NULL) {
+    // data_end < data_begin means that the data or block contains only whitespace
+    if (data != NULL && data_end > data_begin) {
       data->append(line, data_begin, (data_end-data_begin));
 
       if (cvm::debug()) {

lib/colvars/colvarproxy.cpp

@@ -243,7 +243,7 @@ void colvarproxy_atom_groups::compute_max_atom_groups_applied_force()
 
 colvarproxy_smp::colvarproxy_smp()
 {
-  b_smp_active = true; // May be disabled by user option
+  smp_mode = smp_mode_t::cvcs; // May be disabled by user option
   omp_lock_state = NULL;
 #if defined(_OPENMP)
   if (omp_get_thread_num() == 0) {
@@ -265,41 +265,45 @@ colvarproxy_smp::~colvarproxy_smp()
 #endif
 }
 
-
+colvarproxy::smp_mode_t colvarproxy_smp::get_smp_mode() const {
-int colvarproxy_smp::check_smp_enabled()
-{
 #if defined(_OPENMP)
-  if (b_smp_active) {
+  return smp_mode;
-    return COLVARS_OK;
-  }
-  return COLVARS_ERROR;
 #else
-  return COLVARS_NOT_IMPLEMENTED;
+  return colvarproxy::smp_mode_t::none;
+#endif
+}
+
+int colvarproxy_smp::set_smp_mode(smp_mode_t mode) {
+#if defined(_OPENMP)
+  smp_mode = mode;
+  return COLVARS_OK;
+#else
+  if (mode != colvarproxy::smp_mode_t::none) {
+    return COLVARS_NOT_IMPLEMENTED;
+  } else {
+    smp_mode = colvarproxy::smp_mode_t::none;
+  }
+  return COLVARS_OK;
 #endif
 }
 
 
-int colvarproxy_smp::smp_colvars_loop()
+int colvarproxy_smp::smp_loop(int n_items, std::function<int (int)> const &worker)
 {
+  int error_code = COLVARS_OK;
 #if defined(_OPENMP)
-  colvarmodule *cv = cvm::main();
+  cvm::increase_depth();
-  colvarproxy *proxy = cv->proxy;
 #pragma omp parallel for
-  for (int i = 0; i < static_cast<int>(cv->variables_active_smp()->size()); i++) {
+  for (int i = 0; i < n_items; i++) {
-    colvar *x = (*(cv->variables_active_smp()))[i];
+    int const retcode = worker(i);
-    int x_item = (*(cv->variables_active_smp_items()))[i];
+#pragma omp atomic
-    if (cvm::debug()) {
+    error_code |= retcode;
-      cvm::log("["+cvm::to_str(proxy->smp_thread_id())+"/"+
-               cvm::to_str(proxy->smp_num_threads())+
-               "]: calc_colvars_items_smp(), i = "+cvm::to_str(i)+", cv = "+
-               x->name+", cvc = "+cvm::to_str(x_item)+"\n");
-    }
-    x->calc_cvcs(x_item, 1);
   }
-  return cvm::get_error();
+  cvm::decrease_depth();
 #else
-  return COLVARS_NOT_IMPLEMENTED;
+  error_code |= COLVARS_NOT_IMPLEMENTED;
 #endif
+  return error_code;
 }
 
 
@@ -470,8 +474,8 @@ colvarproxy::~colvarproxy()
 
 bool colvarproxy::io_available()
 {
-  return (check_smp_enabled() == COLVARS_OK && smp_thread_id() == 0) ||
+  return ((get_smp_mode() != smp_mode_t::none) && smp_thread_id() == 0) ||
-    (check_smp_enabled() != COLVARS_OK);
+    (get_smp_mode() == smp_mode_t::none);
 }
 
 

lib/colvars/colvarproxy.h

@@ -10,9 +10,12 @@
 #ifndef COLVARPROXY_H
 #define COLVARPROXY_H
 
+#include <functional>
+
 #include "colvarmodule.h"
 #include "colvartypes.h"
 #include "colvarproxy_io.h"
+#include "colvarproxy_replicas.h"
 #include "colvarproxy_system.h"
 #include "colvarproxy_tcl.h"
 #include "colvarproxy_volmaps.h"
@@ -447,21 +450,22 @@ class colvarproxy_smp {
 
 public:
 
+  enum class smp_mode_t {cvcs, inner_loop, none};
+
   /// Constructor
   colvarproxy_smp();
 
   /// Destructor
   virtual ~colvarproxy_smp();
 
-  /// Whether threaded parallelization should be used (TODO: make this a
+  /// Get the current SMP mode
-  /// cvm::deps feature)
+  virtual smp_mode_t get_smp_mode() const;
-  bool b_smp_active;
 
-  /// Whether threaded parallelization is available (TODO: make this a cvm::deps feature)
+  /// Set the current SMP mode
-  virtual int check_smp_enabled();
+  virtual int set_smp_mode(smp_mode_t mode);
 
-  /// Distribute calculation of colvars (and their components) across threads
+  /// Distribute computation over threads using OpenMP, unless overridden in the backend (e.g. NAMD)
-  virtual int smp_colvars_loop();
+  virtual int smp_loop(int n_items, std::function<int (int)> const &worker);
 
   /// Distribute calculation of biases across threads
   virtual int smp_biases_loop();
@@ -488,38 +492,10 @@ protected:
 
   /// Lock state for OpenMP
   omp_lock_t *omp_lock_state;
-};
-
-
-/// \brief Methods for multiple-replica communication
-class colvarproxy_replicas {
-
-public:
-
-  /// Constructor
-  colvarproxy_replicas();
-
-  /// Destructor
-  virtual ~colvarproxy_replicas();
-
-  /// \brief Indicate if multi-replica support is available and active
-  virtual int replica_enabled();
-
-  /// \brief Index of this replica
-  virtual int replica_index();
-
-  /// \brief Total number of replicas
-  virtual int num_replicas();
-
-  /// \brief Synchronize replica with others
-  virtual void replica_comm_barrier();
-
-  /// \brief Receive data from other replica
-  virtual int replica_comm_recv(char* msg_data, int buf_len, int src_rep);
-
-  /// \brief Send data to other replica
-  virtual int replica_comm_send(char* msg_data, int msg_len, int dest_rep);
 
+  /// Whether threaded parallelization should be used (TODO: make this a
+  /// cvm::deps feature)
+  smp_mode_t smp_mode;
 };
 
 

lib/colvars/colvarproxy_io.cpp

@@ -7,10 +7,28 @@
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 
+
+#if defined(_WIN32) && !defined(__CYGWIN__)
+
 // Using access() to check if a file exists (until we can assume C++14/17)
-#if !defined(_WIN32) || defined(__CYGWIN__)
+#include <direct.h>
-#include <unistd.h>
+
+#if defined(__has_include)
+# if __has_include(<filesystem>)
+#  include <filesystem> // MSVC only defines __cpp_lib_filesystem after include
+# endif
 #endif
+
+#else
+
+#include <unistd.h>
+
+#ifdef __cpp_lib_filesystem
+#include <filesystem>
+#endif
+
+#endif
+
 #if defined(_WIN32)
 #include <io.h>
 #endif
@@ -64,6 +82,53 @@ int colvarproxy_io::set_frame(long int)
 }
 
 
+std::string colvarproxy_io::get_current_work_dir() const
+{
+#ifdef __cpp_lib_filesystem
+
+  return std::filesystem::current_path().string();
+
+#else
+
+  // Legacy code
+  size_t constexpr buf_size = 3001;
+  char buf[buf_size];
+
+#if defined(_WIN32) && !defined(__CYGWIN__)
+  char *getcwd_result = ::_getcwd(buf, buf_size);
+#else
+  char *getcwd_result = ::getcwd(buf, buf_size);
+#endif
+
+  if (getcwd_result == nullptr) {
+    cvm::error("Error: cannot read the current working directory.\n", COLVARS_INPUT_ERROR);
+    return std::string("");
+  }
+
+  return std::string(getcwd_result);
+#endif
+}
+
+
+std::string colvarproxy_io::join_paths(std::string const &path1, std::string const &path2) const
+{
+#ifdef __cpp_lib_filesystem
+
+  return (std::filesystem::path(path1) / std::filesystem::path(path2)).string();
+
+#else
+
+  // Legacy code
+#if defined(_WIN32) && !defined(__CYGWIN__)
+  return (path1 + "\\" + path2);
+#else
+  return (path1 + "/" + path2);
+#endif
+
+#endif
+}
+
+
 int colvarproxy_io::backup_file(char const *filename)
 {
   // Simplified version of NAMD_file_exists()

lib/colvars/colvarproxy_io.h

@@ -38,6 +38,12 @@ public:
   // Returns error code
   virtual int set_frame(long int);
 
+  /// Get the current working directory of this process
+  std::string get_current_work_dir() const;
+
+  /// Join two paths using the operating system's path separation
+  std::string join_paths(std::string const &path1, std::string const &path2) const;
+
   /// \brief Rename the given file, before overwriting it
   virtual int backup_file(char const *filename);
 

lib/colvars/colvarproxy_replicas.cpp

@@ -7,50 +7,103 @@
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 
+
 #include "colvarmodule.h"
-#include "colvarproxy.h"
+#include "colvarproxy_replicas.h"
 
 
-colvarproxy_replicas::colvarproxy_replicas() {}
+colvarproxy_replicas::colvarproxy_replicas()
+{
+#ifdef COLVARS_MPI
+  replicas_mpi_comm = MPI_COMM_NULL;
+#endif
+}
 
 
 colvarproxy_replicas::~colvarproxy_replicas() {}
 
 
-int colvarproxy_replicas::replica_enabled()
+void colvarproxy_replicas::set_replicas_mpi_communicator(replicas_mpi_comm_t comm)
 {
+  replicas_mpi_comm = comm;
+#ifdef COLVARS_MPI
+  if (comm != MPI_COMM_NULL) {
+    MPI_Comm_rank(comm, &replicas_mpi_rank);
+    MPI_Comm_size(comm, &replicas_mpi_num);
+    cvm::log("Enabling multiple replicas: this is replica number " +
+             cvm::to_str(replica_index() + 1) + " of " + cvm::to_str(num_replicas()) + ".\n");
+  }
+#endif
+}
+
+
+int colvarproxy_replicas::check_replicas_enabled()
+{
+#ifdef COLVARS_MPI
+  if (replicas_mpi_comm != MPI_COMM_NULL) {
+    return num_replicas() > 1 ? COLVARS_OK : COLVARS_ERROR;
+  }
+  return COLVARS_ERROR;
+#else
   return COLVARS_NOT_IMPLEMENTED;
+#endif
 }
 
 
 int colvarproxy_replicas::replica_index()
 {
-  return 0;
+  return replicas_mpi_rank;
 }
 
 
 int colvarproxy_replicas::num_replicas()
 {
-  return 1;
+  return replicas_mpi_num;
 }
 
 
-void colvarproxy_replicas::replica_comm_barrier() {}
+void colvarproxy_replicas::replica_comm_barrier()
-
-
-int colvarproxy_replicas::replica_comm_recv(char* /* msg_data */,
-                                            int /* buf_len */,
-                                            int /* src_rep */)
 {
-  return COLVARS_NOT_IMPLEMENTED;
+#ifdef COLVARS_MPI
+  MPI_Barrier(replicas_mpi_comm);
+#endif
 }
 
 
-int colvarproxy_replicas::replica_comm_send(char* /* msg_data */,
+int colvarproxy_replicas::replica_comm_recv(char *buffer, int buffer_length, int source_rank)
-                                            int /* msg_len */,
-                                            int /* dest_rep */)
 {
+#ifdef COLVARS_MPI
+  MPI_Status status;
+  int retval = MPI_Recv(buffer, buffer_length, MPI_CHAR, source_rank, 0, replicas_mpi_comm, &status);
+  if (retval == MPI_SUCCESS) {
+    MPI_Get_count(&status, MPI_CHAR, &retval);
+  } else {
+    retval = 0;
+  }
+  return retval;
+#else
+  (void)buffer;
+  (void)buffer_length;
+  (void)source_rank;
   return COLVARS_NOT_IMPLEMENTED;
+#endif
 }
 
 
+int colvarproxy_replicas::replica_comm_send(char *buffer, int buffer_length, int destination_rank)
+{
+#ifdef COLVARS_MPI
+  int retval = MPI_Send(buffer, buffer_length, MPI_CHAR, destination_rank, 0, replicas_mpi_comm);
+  if (retval == MPI_SUCCESS) {
+    retval = buffer_length;
+  } else {
+    retval = 0;
+  }
+  return retval;
+#else
+  (void)buffer;
+  (void)buffer_length;
+  (void)destination_rank;
+  return COLVARS_NOT_IMPLEMENTED;
+#endif
+}

lib/colvars/colvarproxy_replicas.h

@@ -0,0 +1,66 @@
+// -*- c++ -*-
+
+// This file is part of the Collective Variables module (Colvars).
+// The original version of Colvars and its updates are located at:
+// https://github.com/Colvars/colvars
+// Please update all Colvars source files before making any changes.
+// If you wish to distribute your changes, please submit them to the
+// Colvars repository at GitHub.
+
+#ifndef COLVARPROXY_REPLICAS_H
+#define COLVARPROXY_REPLICAS_H
+
+
+#ifdef COLVARS_MPI
+#include <mpi.h>
+typedef MPI_Comm replicas_mpi_comm_t;
+#else
+typedef void * replicas_mpi_comm_t;
+#endif
+
+
+/// \brief Methods for multiple-replica communication
+class colvarproxy_replicas {
+
+public:
+
+  /// Constructor
+  colvarproxy_replicas();
+
+  /// Destructor
+  virtual ~colvarproxy_replicas();
+
+  /// Set the multiple replicas communicator
+  virtual void set_replicas_mpi_communicator(replicas_mpi_comm_t comm);
+
+  /// Indicate if multi-replica support is available and active
+  virtual int check_replicas_enabled();
+
+  /// Index of this replica
+  virtual int replica_index();
+
+  /// Total number of replicas
+  virtual int num_replicas();
+
+  /// Synchronize replica with others
+  virtual void replica_comm_barrier();
+
+  /// Receive data from other replica
+  virtual int replica_comm_recv(char* msg_data, int buf_len, int src_rep);
+
+  /// Send data to other replica
+  virtual int replica_comm_send(char* msg_data, int msg_len, int dest_rep);
+
+protected:
+
+  /// MPI communicator containint 1 root proc from each world
+  replicas_mpi_comm_t replicas_mpi_comm;
+
+  /// Index (rank) of this replica in the MPI implementation
+  int replicas_mpi_rank = 0;
+
+  /// Number of replicas in the MPI implementation
+  int replicas_mpi_num = 1;
+};
+
+#endif

lib/colvars/colvarproxy_system.h

@@ -94,6 +94,7 @@ public:
   virtual bool total_forces_enabled() const;
 
   /// Are total forces from the current step available?
+  /// in which case they are really system forces
   virtual bool total_forces_same_step() const;
 
   /// Get the molecule ID when called in VMD; raise error otherwise
@@ -109,6 +110,11 @@ public:
   /// Send cached value of alchemical lambda parameter to back-end (if available)
   virtual int send_alch_lambda();
 
+  /// Request energy computation every freq steps (necessary for NAMD3, not all back-ends)
+  virtual int request_alch_energy_freq(int const freq) {
+    return COLVARS_OK;
+  }
+
   /// Get energy derivative with respect to lambda (if available)
   virtual int get_dE_dlambda(cvm::real* dE_dlambda);
 

lib/colvars/colvars_memstream.h

@@ -108,6 +108,9 @@ public:
   /// Ignore formatting operators
   inline void setf(decltype(std::ios::fmtflags(0)), decltype(std::ios::floatfield)) {}
 
+  /// Ignore formatting operators
+  inline void setf(decltype(std::ios::fmtflags(0))) {}
+
   /// Ignore formatting operators
   inline void flags(decltype(std::ios::fmtflags(0))) {}
 

lib/colvars/colvars_version.h

@@ -1,3 +1,3 @@
 #ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2024-06-04"
+#define COLVARS_VERSION "2025-04-30"
 #endif

lib/colvars/colvarscript_commands.h

@@ -541,6 +541,15 @@ CVSCRIPT(cv_printframe,
          return COLVARS_OK;
          )
 
+CVSCRIPT(cv_patchversion,
+         "Get the Colvars patch version number (used for bugfixes only)\n"
+         "version : string - Colvars version",
+         0, 0,
+         "",
+         script->set_result_int(cvm::main()->patch_version_number());
+         return COLVARS_OK;
+         )
+
 CVSCRIPT(cv_printframelabels,
          "Return the labels that would be written to colvars.traj\n"
          "Labels : string - The labels",
@@ -656,7 +665,7 @@ CVSCRIPT(cv_update,
          )
 
 CVSCRIPT(cv_version,
-         "Get the Colvars Module version string\n"
+         "Get the Colvars version string\n"
          "version : string - Colvars version",
          0, 0,
          "",
@@ -665,7 +674,7 @@ CVSCRIPT(cv_version,
          )
 
 // This guard allows compiling colvar and bias function bodies in their
-// respecitve files instead of colvarscript_commands.o
+// respective files instead of colvarscript_commands.o
 #ifndef COLVARSCRIPT_COMMANDS_GLOBAL
 #include "colvarscript_commands_colvar.h"
 #include "colvarscript_commands_bias.h"

lib/colvars/colvarscript_commands_colvar.h

@@ -23,6 +23,7 @@ CVSCRIPT(colvar_addforce,
            script->add_error_msg("addforce : error parsing force value");
            return COLVARSCRIPT_ERROR;
          }
+         this_colvar->enable(colvardeps::f_cv_apply_force);
          this_colvar->add_bias_force(force);
          script->set_result_colvarvalue(force);
          return COLVARS_OK;

lib/colvars/colvartypes.cpp

@@ -137,71 +137,6 @@ std::istream & operator >> (std::istream &is, colvarmodule::quaternion &q)
 }
 
 
-cvm::quaternion
-cvm::quaternion::position_derivative_inner(cvm::rvector const &pos,
-                                            cvm::rvector const &vec) const
-{
-  cvm::quaternion result(0.0, 0.0, 0.0, 0.0);
-
-
-  result.q0 =   2.0 * pos.x * q0 * vec.x
-               +2.0 * pos.y * q0 * vec.y
-               +2.0 * pos.z * q0 * vec.z
-
-               -2.0 * pos.y * q3 * vec.x
-               +2.0 * pos.z * q2 * vec.x
-
-               +2.0 * pos.x * q3 * vec.y
-               -2.0 * pos.z * q1 * vec.y
-
-               -2.0 * pos.x * q2 * vec.z
-               +2.0 * pos.y * q1 * vec.z;
-
-
-  result.q1 =  +2.0 * pos.x * q1 * vec.x
-               -2.0 * pos.y * q1 * vec.y
-               -2.0 * pos.z * q1 * vec.z
-
-               +2.0 * pos.y * q2 * vec.x
-               +2.0 * pos.z * q3 * vec.x
-
-               +2.0 * pos.x * q2 * vec.y
-               -2.0 * pos.z * q0 * vec.y
-
-               +2.0 * pos.x * q3 * vec.z
-               +2.0 * pos.y * q0 * vec.z;
-
-
-  result.q2 =  -2.0 * pos.x * q2 * vec.x
-               +2.0 * pos.y * q2 * vec.y
-               -2.0 * pos.z * q2 * vec.z
-
-               +2.0 * pos.y * q1 * vec.x
-               +2.0 * pos.z * q0 * vec.x
-
-               +2.0 * pos.x * q1 * vec.y
-               +2.0 * pos.z * q3 * vec.y
-
-               -2.0 * pos.x * q0 * vec.z
-               +2.0 * pos.y * q3 * vec.z;
-
-
-  result.q3 =  -2.0 * pos.x * q3 * vec.x
-               -2.0 * pos.y * q3 * vec.y
-               +2.0 * pos.z * q3 * vec.z
-
-               -2.0 * pos.y * q0 * vec.x
-               +2.0 * pos.z * q1 * vec.x
-
-               +2.0 * pos.x * q0 * vec.y
-               +2.0 * pos.z * q2 * vec.y
-
-               +2.0 * pos.x * q1 * vec.z
-               +2.0 * pos.y * q2 * vec.z;
-
-  return result;
-}
-
 #ifdef COLVARS_LAMMPS
 namespace {
   inline void *new_Jacobi_solver(int size) {
@@ -336,7 +271,7 @@ void colvarmodule::rotation::compute_overlap_matrix()
 #ifndef COLVARS_LAMMPS
 namespace NR {
 
-void diagonalize_matrix(cvm::real m[4][4],
+int diagonalize_matrix(cvm::real m[4][4],
                         cvm::real eigval[4],
                         cvm::real eigvec[4][4])
 {
@@ -347,9 +282,7 @@ void diagonalize_matrix(cvm::real m[4][4],
   int jac_nrot = 0;
   if (NR_Jacobi::jacobi(m, eigval, eigvec, &jac_nrot) !=
       COLVARS_OK) {
-    cvm::error("Too many iterations in jacobi diagonalization.\n"
+    return COLVARS_ERROR;
-               "This is usually the result of an ill-defined set of atoms for "
-               "rotational alignment (RMSD, rotateReference, etc).\n");
   }
   NR_Jacobi::eigsrt(eigval, eigvec);
   // jacobi saves eigenvectors by columns
@@ -367,6 +300,7 @@ void diagonalize_matrix(cvm::real m[4][4],
       eigvec[ie][i] /= norm;
     }
   }
+  return COLVARS_OK;
 }
 
 }
@@ -429,14 +363,25 @@ void colvarmodule::rotation::calc_optimal_rotation_impl() {
                                        cvm::real[4][4]> *>(jacobi);
 
   int ierror = ecalc->Diagonalize(S, S_eigval, S_eigvec);
+#else
+  int ierror = NR::diagonalize_matrix(S, S_eigval, S_eigvec);
+#endif
   if (ierror) {
+    cvm::log("Failed to diagonalize the following overlapping matrix:\n");
+    for (size_t i = 0; i < 4; ++i) {
+      for (size_t j = 0; j < 4; ++j) {
+        cvm::log(cvm::to_str(S[i][j]) + " ");
+      }
+      cvm::log("\n");
+    }
+    cvm::log("The corresponding correlation matrix is:\n");
+    cvm::log(" " + cvm::to_str(C.xx) + " " + cvm::to_str(C.xy) + " " + cvm::to_str(C.xz));
+    cvm::log(" " + cvm::to_str(C.yx) + " " + cvm::to_str(C.yy) + " " + cvm::to_str(C.yz));
+    cvm::log(" " + cvm::to_str(C.zx) + " " + cvm::to_str(C.zy) + " " + cvm::to_str(C.zz) + "\n");
     cvm::error("Too many iterations in jacobi diagonalization.\n"
                "This is usually the result of an ill-defined set of atoms for "
                "rotational alignment (RMSD, rotateReference, etc).\n");
   }
-#else
-  NR::diagonalize_matrix(S, S_eigval, S_eigvec);
-#endif
   q = cvm::quaternion{S_eigvec[0][0], S_eigvec[0][1], S_eigvec[0][2], S_eigvec[0][3]};
 
   if (cvm::rotation::monitor_crossings) {

lib/colvars/colvartypes.h

@@ -20,10 +20,6 @@
 
 #include "colvarmodule.h"
 
-#ifndef PI
-#define PI 3.14159265358979323846
-#endif
-
 // ----------------------------------------------------------------------
 /// Linear algebra functions and data types used in the collective
 /// variables implemented so far
@@ -1221,8 +1217,57 @@ public:
 
   /// \brief Multiply the given vector by the derivative of the given
   /// (rotated) position with respect to the quaternion
-  cvm::quaternion position_derivative_inner(cvm::rvector const &pos,
+  /// \param pos The position \f$\mathbf{x}\f$.
-                                            cvm::rvector const &vec) const;
+  /// \param vec The vector \f$\mathbf{v}\f$.
+  /// \return A quaternion (see the detailed documentation below).
+  ///
+  /// This function is mainly used for projecting the gradients or forces on
+  /// the rotated atoms to the forces on quaternion. Assume this rotation can
+  /// be represented as \f$R(\mathbf{q})\f$,
+  /// where \f$\mathbf{q} := (q_0, q_1, q_2, q_3)\f$
+  /// is the current quaternion, the function returns the following new
+  /// quaternion:
+  /// \f[
+  /// \left(\mathbf{v}^\mathrm{T}\frac{\partial R(\mathbf{q})}{\partial q_0}\mathbf{x},
+  ///       \mathbf{v}^\mathrm{T}\frac{\partial R(\mathbf{q})}{\partial q_1}\mathbf{x},
+  ///       \mathbf{v}^\mathrm{T}\frac{\partial R(\mathbf{q})}{\partial q_2}\mathbf{x},
+  ///       \mathbf{v}^\mathrm{T}\frac{\partial R(\mathbf{q})}{\partial q_3}\mathbf{x}\right)
+  /// \f]
+  /// where \f$\mathbf{v}\f$ is usually the gradient of \f$\xi\f$ with respect to
+  /// the rotated frame \f$\tilde{\mathbf{X}}\f$,
+  /// \f$\partial \xi / \partial \tilde{\mathbf{X}}\f$, or the force acting on it
+  /// (\f$\mathbf{F}_{\tilde{\mathbf{X}}}\f$).
+  /// By using the following loop in pseudo C++ code,
+  /// either \f$\partial \xi / \partial \tilde{\mathbf{X}}\f$
+  /// or \f$\mathbf{F}_{\tilde{\mathbf{X}}}\f$, can be projected to
+  /// \f$\partial \xi / \partial \mathbf{q}\f$ or \f$\mathbf{F}_q\f$ into `sum_dxdq`:
+  /// @code
+  /// cvm::real sum_dxdq[4] = {0, 0, 0, 0};
+  /// for (size_t i = 0; i < main_group_size(); ++i) {
+  ///   const cvm::rvector v = grad_or_force_on_rotated_main_group(i);
+  ///   const cvm::rvector x = unrotated_main_group_positions(i);
+  ///   cvm::quaternion const dxdq = position_derivative_inner(x, v);
+  ///   sum_dxdq[0] += dxdq[0];
+  ///   sum_dxdq[1] += dxdq[1];
+  ///   sum_dxdq[2] += dxdq[2];
+  ///   sum_dxdq[3] += dxdq[3];
+  /// }
+  /// @endcode
+  inline cvm::quaternion position_derivative_inner(cvm::rvector const &pos,
+                                            cvm::rvector const &vec) const {
+    return cvm::quaternion(2.0 * (vec.x * ( q0 * pos.x - q3 * pos.y + q2 * pos.z) +
+                                  vec.y * ( q3 * pos.x + q0 * pos.y - q1 * pos.z) +
+                                  vec.z * (-q2 * pos.x + q1 * pos.y + q0 * pos.z)),
+                           2.0 * (vec.x * ( q1 * pos.x + q2 * pos.y + q3 * pos.z) +
+                                  vec.y * ( q2 * pos.x - q1 * pos.y - q0 * pos.z) +
+                                  vec.z * ( q3 * pos.x + q0 * pos.y - q1 * pos.z)),
+                           2.0 * (vec.x * (-q2 * pos.x + q1 * pos.y + q0 * pos.z) +
+                                  vec.y * ( q1 * pos.x + q2 * pos.y + q3 * pos.z) +
+                                  vec.z * (-q0 * pos.x + q3 * pos.y - q2 * pos.z)),
+                           2.0 * (vec.x * (-q3 * pos.x - q0 * pos.y + q1 * pos.z) +
+                                  vec.y * ( q0 * pos.x - q3 * pos.y + q2 * pos.z) +
+                                  vec.z * ( q1 * pos.x + q2 * pos.y + q3 * pos.z)));
+  }
 
 
   /// \brief Return the cosine between the orientation frame
@@ -1301,7 +1346,7 @@ public:
 
 #ifndef COLVARS_LAMMPS
 namespace NR {
-void diagonalize_matrix(cvm::real m[4][4],
+int diagonalize_matrix(cvm::real m[4][4],
                         cvm::real eigval[4],
                         cvm::real eigvec[4][4]);
 }

lib/colvars/colvarvalue.cpp

@@ -153,29 +153,6 @@ std::string const colvarvalue::type_keyword(Type t)
 }
 
 
-size_t colvarvalue::num_df(Type t)
-{
-  switch (t) {
-  case colvarvalue::type_notset:
-  default:
-    return 0; break;
-  case colvarvalue::type_scalar:
-    return 1; break;
-  case colvarvalue::type_3vector:
-    return 3; break;
-  case colvarvalue::type_unit3vector:
-  case colvarvalue::type_unit3vectorderiv:
-    return 2; break;
-  case colvarvalue::type_quaternion:
-  case colvarvalue::type_quaternionderiv:
-    return 3; break;
-  case colvarvalue::type_vector:
-    // the size of a vector is unknown without its object
-    return 0; break;
-  }
-}
-
-
 size_t colvarvalue::num_dimensions(Type t)
 {
   switch (t) {
@@ -591,34 +568,132 @@ cvm::real operator * (colvarvalue const &x1,
 }
 
 
+cvm::real colvarvalue::norm2() const
+{
+  switch (value_type) {
+  case colvarvalue::type_scalar:
+    return (this->real_value)*(this->real_value);
+  case colvarvalue::type_3vector:
+  case colvarvalue::type_unit3vector:
+  case colvarvalue::type_unit3vectorderiv:
+    return (this->rvector_value).norm2();
+  case colvarvalue::type_quaternion:
+  case colvarvalue::type_quaternionderiv:
+    return (this->quaternion_value).norm2();
+  case colvarvalue::type_vector:
+    if (elem_types.size() > 0) {
+      // if we have information about non-scalar types, use it
+      cvm::real result = 0.0;
+      size_t i;
+      for (i = 0; i < elem_types.size(); i++) {
+        result += (this->get_elem(i)).norm2();
+      }
+      return result;
+    } else {
+      return vector1d_value.norm2();
+    }
+    break;
+  case colvarvalue::type_notset:
+  default:
+    return 0.0;
+  }
+}
+
+
+cvm::real colvarvalue::sum() const
+{
+  switch (value_type) {
+  case colvarvalue::type_scalar:
+    return (this->real_value);
+  case colvarvalue::type_3vector:
+  case colvarvalue::type_unit3vector:
+  case colvarvalue::type_unit3vectorderiv:
+    return (this->rvector_value).x + (this->rvector_value).y +
+      (this->rvector_value).z;
+  case colvarvalue::type_quaternion:
+  case colvarvalue::type_quaternionderiv:
+    return (this->quaternion_value).q0 + (this->quaternion_value).q1 +
+      (this->quaternion_value).q2 + (this->quaternion_value).q3;
+  case colvarvalue::type_vector:
+    return (this->vector1d_value).sum();
+  case colvarvalue::type_notset:
+  default:
+    return 0.0;
+  }
+}
+
+
+cvm::real colvarvalue::dist2(colvarvalue const &x2) const
+{
+  colvarvalue::check_types(*this, x2);
+
+  switch (this->type()) {
+  case colvarvalue::type_scalar:
+    return (this->real_value - x2.real_value) * (this->real_value - x2.real_value);
+  case colvarvalue::type_3vector:
+    return (this->rvector_value - x2.rvector_value).norm2();
+  case colvarvalue::type_unit3vector: {
+    cvm::rvector const &v1 = this->rvector_value;
+    cvm::rvector const &v2 = x2.rvector_value;
+    cvm::real const theta = cvm::acos(v1 * v2);
+    return theta * theta;
+  }
+  case colvarvalue::type_quaternion:
+    // angle between (*this) and x2 is the distance, the quaternion
+    // object has it implemented internally
+    return this->quaternion_value.dist2(x2.quaternion_value);
+  case colvarvalue::type_vector:
+    return (this->vector1d_value - x2.vector1d_value).norm2();
+  case colvarvalue::type_unit3vectorderiv:
+  case colvarvalue::type_quaternionderiv:
+    cvm::error("Error: computing a squared-distance between two variables of type \"" +
+                   type_desc(this->type()) + "\", for which it is not defined.\n",
+               COLVARS_BUG_ERROR);
+  case colvarvalue::type_notset:
+  default:
+    this->undef_op();
+    return 0.0;
+  };
+
+  return 0.0;
+}
+
+
 colvarvalue colvarvalue::dist2_grad(colvarvalue const &x2) const
 {
   colvarvalue::check_types(*this, x2);
 
+  // Compute derivative with respect to (*this)
+
   switch (this->value_type) {
   case colvarvalue::type_scalar:
     return 2.0 * (this->real_value - x2.real_value);
   case colvarvalue::type_3vector:
     return 2.0 * (this->rvector_value - x2.rvector_value);
-  case colvarvalue::type_unit3vector:
+  case colvarvalue::type_unit3vector: {
-  case colvarvalue::type_unit3vectorderiv:
+    cvm::rvector const &v1 = this->rvector_value;
-    {
+    cvm::rvector const &v2 = x2.rvector_value;
-      cvm::rvector const &v1 = this->rvector_value;
+    cvm::real const cos_t = v1 * v2;
-      cvm::rvector const &v2 = x2.rvector_value;
+    return colvarvalue(2.0 * cvm::acos(cos_t) * -1.0 / cvm::sqrt(1.0 - cos_t * cos_t) * v2,
-      cvm::real const cos_t = v1 * v2;
+                       colvarvalue::type_unit3vectorderiv);
-      return colvarvalue(2.0 * (cos_t * v1 - v2), colvarvalue::type_unit3vectorderiv);
+  }
-    }
   case colvarvalue::type_quaternion:
-  case colvarvalue::type_quaternionderiv:
     return this->quaternion_value.dist2_grad(x2.quaternion_value);
   case colvarvalue::type_vector:
     return colvarvalue(2.0 * (this->vector1d_value - x2.vector1d_value), colvarvalue::type_vector);
     break;
+  case colvarvalue::type_unit3vectorderiv:
+  case colvarvalue::type_quaternionderiv:
+    cvm::error("Error: computing a squared-distance gradient between two variables of type \"" +
+                   type_desc(this->type()) + "\", for which it is not defined.\n",
+               COLVARS_BUG_ERROR);
   case colvarvalue::type_notset:
   default:
     this->undef_op();
     return colvarvalue(colvarvalue::type_notset);
   };
+
+  return colvarvalue(colvarvalue::type_notset);
 }
 
 

lib/colvars/colvarvalue.h

@@ -109,9 +109,6 @@ public:
   /// User keywords for specifying value types in the configuration
   static std::string const type_keyword(Type t);
 
-  /// Number of degrees of freedom for each supported type
-  static size_t num_df(Type t);
-
   /// Number of dimensions for each supported type (used to allocate vector1d_value)
   static size_t num_dimensions(Type t);
 
@@ -671,87 +668,4 @@ inline cvm::vector1d<cvm::real> const colvarvalue::as_vector() const
 }
 
 
-inline cvm::real colvarvalue::norm2() const
-{
-  switch (value_type) {
-  case colvarvalue::type_scalar:
-    return (this->real_value)*(this->real_value);
-  case colvarvalue::type_3vector:
-  case colvarvalue::type_unit3vector:
-  case colvarvalue::type_unit3vectorderiv:
-    return (this->rvector_value).norm2();
-  case colvarvalue::type_quaternion:
-  case colvarvalue::type_quaternionderiv:
-    return (this->quaternion_value).norm2();
-  case colvarvalue::type_vector:
-    if (elem_types.size() > 0) {
-      // if we have information about non-scalar types, use it
-      cvm::real result = 0.0;
-      size_t i;
-      for (i = 0; i < elem_types.size(); i++) {
-        result += (this->get_elem(i)).norm2();
-      }
-      return result;
-    } else {
-      return vector1d_value.norm2();
-    }
-    break;
-  case colvarvalue::type_notset:
-  default:
-    return 0.0;
-  }
-}
-
-
-inline cvm::real colvarvalue::sum() const
-{
-  switch (value_type) {
-  case colvarvalue::type_scalar:
-    return (this->real_value);
-  case colvarvalue::type_3vector:
-  case colvarvalue::type_unit3vector:
-  case colvarvalue::type_unit3vectorderiv:
-    return (this->rvector_value).x + (this->rvector_value).y +
-      (this->rvector_value).z;
-  case colvarvalue::type_quaternion:
-  case colvarvalue::type_quaternionderiv:
-    return (this->quaternion_value).q0 + (this->quaternion_value).q1 +
-      (this->quaternion_value).q2 + (this->quaternion_value).q3;
-  case colvarvalue::type_vector:
-    return (this->vector1d_value).sum();
-  case colvarvalue::type_notset:
-  default:
-    return 0.0;
-  }
-}
-
-
-inline cvm::real colvarvalue::dist2(colvarvalue const &x2) const
-{
-  colvarvalue::check_types(*this, x2);
-
-  switch (this->type()) {
-  case colvarvalue::type_scalar:
-    return (this->real_value - x2.real_value)*(this->real_value - x2.real_value);
-  case colvarvalue::type_3vector:
-    return (this->rvector_value - x2.rvector_value).norm2();
-  case colvarvalue::type_unit3vector:
-  case colvarvalue::type_unit3vectorderiv:
-    // angle between (*this) and x2 is the distance
-    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
-    // object has it implemented internally
-    return this->quaternion_value.dist2(x2.quaternion_value);
-  case colvarvalue::type_vector:
-    return (this->vector1d_value - x2.vector1d_value).norm2();
-  case colvarvalue::type_notset:
-  default:
-    this->undef_op();
-    return 0.0;
-  };
-}
-
-
 #endif

src/COLVARS/colvarproxy_lammps.cpp

@@ -33,12 +33,9 @@ colvarproxy_lammps::colvarproxy_lammps(LAMMPS_NS::LAMMPS *lmp)  : _lmp(lmp), _ra
   previous_step = -1;
   do_exit = false;
 
-  inter_me = 0;
-  inter_num = 1;
   bias_energy = 0.0;
 
   engine_ready_ = false;
-  inter_comm = MPI_COMM_NULL;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -83,19 +80,6 @@ void colvarproxy_lammps::set_random_seed(int seed)
   _random = new LAMMPS_NS::RanPark(_lmp, seed);
 }
 
-/* ---------------------------------------------------------------------- */
-
-void colvarproxy_lammps::set_replicas_communicator(MPI_Comm root2root)
-{
-  inter_comm = root2root;
-
-  // initialize multi-replica support, if available
-  if (replica_enabled() == COLVARS_OK) {
-    MPI_Comm_rank(inter_comm, &inter_me);
-    MPI_Comm_size(inter_comm, &inter_num);
-  }
-}
-
 /* ----------------------------------------------------------------------
    re-initialize data where needed
 ------------------------------------------------------------------------- */
@@ -255,63 +239,7 @@ int colvarproxy_lammps::set_unit_system(std::string const &units_in, bool /*chec
   return COLVARS_OK;
 }
 
-/* ----------------------------------------------------------------------
-   multi-replica support
-------------------------------------------------------------------------- */
 
-int colvarproxy_lammps::replica_enabled()
-{
-  return (inter_comm != MPI_COMM_NULL) ? COLVARS_OK : COLVARS_NOT_IMPLEMENTED;
-}
-
-/* ---------------------------------------------------------------------- */
-
-int colvarproxy_lammps::replica_index()
-{
-  return inter_me;
-}
-
-/* ---------------------------------------------------------------------- */
-
-int colvarproxy_lammps::num_replicas()
-{
-  return inter_num;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void colvarproxy_lammps::replica_comm_barrier()
-{
-  MPI_Barrier(inter_comm);
-}
-
-/* ---------------------------------------------------------------------- */
-
-int colvarproxy_lammps::replica_comm_recv(char* msg_data, int buf_len, int src_rep)
-{
-  MPI_Status status;
-  int retval;
-
-  retval = MPI_Recv(msg_data,buf_len,MPI_CHAR,src_rep,0,inter_comm,&status);
-  if (retval == MPI_SUCCESS) {
-    MPI_Get_count(&status, MPI_CHAR, &retval);
-  } else retval = 0;
-  return retval;
-}
-
-/* ---------------------------------------------------------------------- */
-
-int colvarproxy_lammps::replica_comm_send(char* msg_data, int msg_len, int dest_rep)
-{
-  int retval;
-  retval = MPI_Send(msg_data,msg_len,MPI_CHAR,dest_rep,0,inter_comm);
-  if (retval == MPI_SUCCESS) {
-    retval = msg_len;
-  } else retval = 0;
-  return retval;
-}
-
-/* ---------------------------------------------------------------------- */
 
 int colvarproxy_lammps::check_atom_id(int atom_number)
 {

src/COLVARS/colvarproxy_lammps.h

@@ -45,9 +45,6 @@ class colvarproxy_lammps : public colvarproxy {
 
   std::vector<int> atoms_types;
 
-  MPI_Comm inter_comm;        // MPI comm with 1 root proc from each world
-  int inter_me, inter_num;    // rank for the inter replica comm
-
  public:
   friend class cvm::atom;
 
@@ -59,9 +56,6 @@ class colvarproxy_lammps : public colvarproxy {
   /// Set the internal seed used by \link rand_gaussian() \endlink
   void set_random_seed(int seed);
 
-  /// Set the multiple replicas communicator
-  void set_replicas_communicator(MPI_Comm root2root);
-
   int setup() override;
 
   // disable default and copy constructor
@@ -72,7 +66,8 @@ class colvarproxy_lammps : public colvarproxy {
   // methods for lammps to move data or trigger actions in the proxy
  public:
   bool total_forces_enabled() const override { return total_force_requested; };
-  bool total_forces_same_step() const override { return true; };
+  // Total forces are saved at end of step, only processed at the next step
+  bool total_forces_same_step() const override { return false; };
   bool want_exit() const { return do_exit; };
 
   // perform colvars computation. returns biasing energy
@@ -102,14 +97,6 @@ class colvarproxy_lammps : public colvarproxy {
   int check_atom_id(int atom_number) override;
 
   inline std::vector<int> *modify_atom_types() { return &atoms_types; }
-
-  int replica_enabled() override;
-  int replica_index() override;
-  int num_replicas() override;
-
-  void replica_comm_barrier() override;
-  int replica_comm_recv(char *msg_data, int buf_len, int src_rep) override;
-  int replica_comm_send(char *msg_data, int msg_len, int dest_rep) override;
 };
 
 #endif

src/COLVARS/colvarproxy_lammps_version.h

@@ -1,3 +1,3 @@
 #ifndef COLVARPROXY_VERSION
-#define COLVARPROXY_VERSION "2024-07-05"
+#define COLVARPROXY_VERSION "2025-03-31"
 #endif

src/COLVARS/fix_colvars.cpp

@@ -267,6 +267,7 @@ void FixColvars::init()
   if (init_flag) return;
   init_flag = 1;
 
+#if defined(COLVARS_MPI)
   if (universe->nworlds > 1) {
     // create inter root communicator
     int color = 1;
@@ -275,9 +276,10 @@ void FixColvars::init()
     }
     MPI_Comm_split(universe->uworld, color, universe->iworld, &root2root);
     if (me == 0) {
-      proxy->set_replicas_communicator(root2root);
+      proxy->set_replicas_mpi_communicator(root2root);
     }
   }
+#endif
 }