diff --git a/lib/colvars/colvaratoms.cpp b/lib/colvars/colvaratoms.cpp
index e15b9301a12..745cb032816 100644
--- a/lib/colvars/colvaratoms.cpp
+++ b/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,15 +1252,15 @@ 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];
@@ -1261,26 +1268,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);
       // 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");
     }
   }
 }
diff --git a/lib/colvars/colvaratoms.h b/lib/colvars/colvaratoms.h
index d16ca7bd567..1493046e033 100644
--- a/lib/colvars/colvaratoms.h
+++ b/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;
diff --git a/lib/colvars/colvarbias_meta.cpp b/lib/colvars/colvarbias_meta.cpp
index 905cd17883b..004dcbb2e63 100644
--- a/lib/colvars/colvarbias_meta.cpp
+++ b/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
 
diff --git a/lib/colvars/colvarcomp_distances.cpp b/lib/colvars/colvarcomp_distances.cpp
index 319190c385a..1c7267f1e5e 100644
--- a/lib/colvars/colvarcomp_distances.cpp
+++ b/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)
-    group1->apply_force(-1.0 * force_tang);
-
-  if (!group2->noforce)
-    group2->apply_force(       force_tang);
+  if (!group1->noforce) {
+    group1->apply_force(-1.0 / dist_v.norm() * force_tang);
+  }
+  if (!group2->noforce) {
+    group2->apply_force( 1.0 / dist_v.norm() * force_tang);
+  }
 }
 
 
-cvm::real colvar::distance_dir::dist2(colvarvalue const &x1,
-                                      colvarvalue const &x2) const
+cvm::real colvar::distance_dir::dist2(colvarvalue const &x1, 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 const &x2) const
+colvarvalue colvar::distance_dir::dist2_lgrad(colvarvalue const &x1, 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 const &x2) const
+colvarvalue colvar::distance_dir::dist2_rgrad(colvarvalue const &x1, colvarvalue const &x2) const
 {
-  return colvarvalue((x2.rvector_value - x1.rvector_value), colvarvalue::type_unit3vectorderiv);
+  return x2.dist2_grad(x1);
 }
 
 
@@ -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);
     }
   }
 }
diff --git a/lib/colvars/colvarcomp_rotations.cpp b/lib/colvars/colvarcomp_rotations.cpp
index a04ace851a0..12dd90e33be 100644
--- a/lib/colvars/colvarcomp_rotations.cpp
+++ b/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);
-      for (size_t i = 0; i < 4; i++) {
-        (*atoms)[ia].apply_force(FQ[i] * dq0_2[i]);
-      }
+      const auto f_ia = FQ[0] * dq0_2[0] +
+                        FQ[1] * dq0_2[1] +
+                        FQ[2] * dq0_2[2] +
+                        FQ[3] * dq0_2[3];
+      ag_force.add_atom_force(ia, f_ia);
     }
   }
 }
diff --git a/lib/colvars/colvargrid.cpp b/lib/colvars/colvargrid.cpp
index 11693a75874..36f83b31ae3 100644
--- a/lib/colvars/colvargrid.cpp
+++ b/lib/colvars/colvargrid.cpp
@@ -617,7 +617,7 @@ integrate_potential::integrate_potential(std::vector<colvar *> &colvars, std::sh
 }
 
 
-integrate_potential::integrate_potential(std::shared_ptr<colvar_grid_gradient> gradients)
+integrate_potential::integrate_potential(colvar_grid_gradient * gradients)
   : b_smoothed(false),
     gradients(gradients)
 {
diff --git a/lib/colvars/colvargrid.h b/lib/colvars/colvargrid.h
index 4cbbb109614..0f168af07fb 100644
--- a/lib/colvars/colvargrid.h
+++ b/lib/colvars/colvargrid.h
@@ -1832,7 +1832,7 @@ class integrate_potential : public colvar_grid_scalar
   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);
+  integrate_potential(colvar_grid_gradient * gradients);
 
   /// \brief Calculate potential from divergence (in 2D); return number of steps
   int integrate(const int itmax, const cvm::real & tol, cvm::real & err, bool verbose = true);
diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h
index fa84b1ad75d..cf551e3d7c3 100644
--- a/lib/colvars/colvarmodule.h
+++ b/lib/colvars/colvarmodule.h
@@ -84,7 +84,7 @@ class colvarmodule {
   int version_int = 0;
 
   /// Patch version number (non-zero in patch releases of other packages)
-  int patch_version_int = 0;
+  int patch_version_int = 4;
 
 public:
 
diff --git a/lib/colvars/colvarproxy_io.cpp b/lib/colvars/colvarproxy_io.cpp
index 4cfdfeec26f..08666d138e6 100644
--- a/lib/colvars/colvarproxy_io.cpp
+++ b/lib/colvars/colvarproxy_io.cpp
@@ -8,13 +8,20 @@
 // Colvars repository at GitHub.
 
 // Using access() to check if a file exists (until we can assume C++14/17)
-#if !defined(_WIN32) || defined(__CYGWIN__)
+#if defined(_WIN32) && !defined(__CYGWIN__)
+#include <direct.h>
+#else
 #include <unistd.h>
 #endif
+
 #if defined(_WIN32)
 #include <io.h>
 #endif
 
+#ifdef __cpp_lib_filesystem
+#include <filesystem>
+#endif
+
 #include <cerrno>
 #include <cstdio>
 
@@ -64,6 +71,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()
diff --git a/lib/colvars/colvarproxy_io.h b/lib/colvars/colvarproxy_io.h
index 726f915c970..eaf750366d6 100644
--- a/lib/colvars/colvarproxy_io.h
+++ b/lib/colvars/colvarproxy_io.h
@@ -38,6 +38,12 @@ class colvarproxy_io {
   // 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);
 
diff --git a/lib/colvars/colvarvalue.cpp b/lib/colvars/colvarvalue.cpp
index 3b8077d2e77..0a7223d9e77 100644
--- a/lib/colvars/colvarvalue.cpp
+++ b/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,6 +568,97 @@ 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);
@@ -600,25 +668,30 @@ colvarvalue colvarvalue::dist2_grad(colvarvalue const &x2) const
     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_unit3vectorderiv:
-    {
-      cvm::rvector const &v1 = this->rvector_value;
-      cvm::rvector const &v2 = x2.rvector_value;
-      cvm::real const cos_t = v1 * v2;
-      return colvarvalue(2.0 * (cos_t * v1 - v2), colvarvalue::type_unit3vectorderiv);
-    }
+  case colvarvalue::type_unit3vector: {
+    cvm::rvector const &v1 = this->rvector_value;
+    cvm::rvector const &v2 = x2.rvector_value;
+    cvm::real const cos_t = v1 * v2;
+    return colvarvalue(2.0 * std::acos(cos_t) * -1.0 / cvm::sqrt(1.0 - cos_t * cos_t) * 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);
 }
 
 
diff --git a/lib/colvars/colvarvalue.h b/lib/colvars/colvarvalue.h
index e8a6a849d35..61f1bf718b0 100644
--- a/lib/colvars/colvarvalue.h
+++ b/lib/colvars/colvarvalue.h
@@ -109,9 +109,6 @@ class colvarvalue {
   /// 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