Fixes for the Colvars library

- 759 min_image fix; Saves long runs from crashes;
  Colvars/colvars#759 (@PolyachenkoYA)

- 728 Fix undefined behavior when getting the current working directory from std::filesystem
  Colvars/colvars#728 (@giacomofiorin)

- 724 Fix gradients and metric functions of distanceDir
  Colvars/colvars#724 (@giacomofiorin)

- 715 Add missing rotation in orientation component
  Colvars/colvars#715 (@giacomofiorin)

- 713 fix: try to solve lammps#87 for non-scala components
  Colvars/colvars#713 (@HanatoK)

- 706 BUGFIX for Segmentation fault in colvarbias_meta::calc_energy() with useGrids off
  Colvars/colvars#706 (@alphataubio)

- 701 Do not try accessing LAMMPS proxy object before allocating it
  Colvars/colvars#701 (@giacomofiorin)

Authors: @alphataubio, @giacomofiorin, @HanatoK, @PolyachenkoYA

lib/colvars/colvaratoms.cpp

@@ -1217,23 +1217,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>
-void cvm::atom_group::calc_fit_gradients_impl() {
-  cvm::atom_group *group_for_fit = fitting_group ? fitting_group : this;
+template <bool B_ag_center, bool B_ag_rotate,
+          typename main_force_accessor_T, typename fitting_force_accessor_T>
+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
@@ -1245,42 +1252,61 @@ void cvm::atom_group::calc_fit_gradients_impl() {
   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) pos_orig = rot_inv * 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_orig, accessor_main(i));
       sum_dxdq[0] += dxdq[0];
       sum_dxdq[1] += dxdq[1];
       sum_dxdq[2] += dxdq[2];
       sum_dxdq[3] += dxdq[3];
     }
   }
   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);
       // 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] +
-                                          sum_dxdq[1] * dq0_1[1] +
-                                          sum_dxdq[2] * dq0_1[2] +
-                                          sum_dxdq[3] * dq0_1[3];
+      fitting_force_grad += sum_dxdq[0] * dq0_1[0] +
+                            sum_dxdq[1] * dq0_1[1] +
+                            sum_dxdq[2] * dq0_1[2] +
+                            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 +1478,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();
-    for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
-      ai->apply_force(rot_inv * ((ai->mass/total_mass) * force));
+cvm::atom_group::group_force_object cvm::atom_group::get_group_force_object() {
+  return cvm::atom_group::group_force_object(this);
+}
+
+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++) {
-      ai->apply_force((ai->mass/total_mass) * force);
+void cvm::atom_group::group_force_object::apply_force_with_fitting_group() {
+  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->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

@@ -257,8 +257,27 @@ class colvarmodule::atom_group
   /// \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 {
+  public:
+    group_force_object(cvm::atom_group* ag);
+    ~group_force_object();
+    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];
@@ -497,15 +516,47 @@ class colvarmodule::atom_group
   /// \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.
+ *  @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.
+ */
+  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.
  */
-  template <bool B_ag_center, bool B_ag_rotate> void calc_fit_gradients_impl();
+  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_meta.cpp

@@ -11,24 +11,6 @@
 #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"
@@ -451,8 +433,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) {
@@ -670,11 +655,20 @@ int colvarbias_meta::calc_energy(std::vector<colvarvalue> const *values)
     replicas[ir]->bias_energy = 0.0;
   }
 
-  std::vector<int> const curr_bin = values ?
-    hills_energy->get_colvars_index(*values) :
-    hills_energy->get_colvars_index();
+  bool index_ok = false;
+  std::vector<int> 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 (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 +717,20 @@ int colvarbias_meta::calc_forces(std::vector<colvarvalue> const *values)
     }
   }
 
-  std::vector<int> const curr_bin = values ?
-    hills_energy->get_colvars_index(*values) :
-    hills_energy->get_colvars_index();
+  bool index_ok = false;
+  std::vector<int> curr_bin;
 
-  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++) {
@@ -1718,29 +1721,17 @@ int colvarbias_meta::setup_output()
 
   if (comm == multiple_replicas) {
 
-    // TODO: one may want to specify the path manually for intricated filesystems?
-    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);
-    }
-
+    auto const pwd = cvm::main()->proxy->get_current_work_dir();
     replica_list_file =
-      (std::string(pwd)+std::string(PATHSEP)+
-       this->name+"."+replica_id+".files.txt");
+        cvm::main()->proxy->join_paths(pwd, 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 =
-      (std::string(pwd)+std::string(PATHSEP)+
-       cvm::output_prefix()+".colvars."+this->name+"."+replica_id+".hills");
-    replica_state_file =
-      (std::string(pwd)+std::string(PATHSEP)+
-       cvm::output_prefix()+".colvars."+this->name+"."+replica_id+".state");
-    delete[] pwd;
+    replica_hills_file = cvm::main()->proxy->join_paths(
+        pwd, cvm::output_prefix() + ".colvars." + this->name + "." + replica_id + ".hills");
+
+    replica_state_file = cvm::main()->proxy->join_paths(
+        pwd, cvm::output_prefix() + ".colvars." + this->name + "." + replica_id + ".state");
 
     // now register this replica