fix: try to solve Colvars#87 for non-scala components

We have already had the algorithm for computing the gradients of the
fitting group from the gradients of the main group. In cases that the
gradients of the main group are not available, we can "intercept" the
forces applied to the main group, and then use the same algorithm to
convert the forces applied to the main group with the rotation to the
forces applied to the fitting group.

src/colvaratoms.cpp

@@ -1282,6 +1282,79 @@ void cvm::atom_group::calc_fit_gradients_impl() {
 }
 
 
+template <bool B_ag_center, bool B_ag_rotate>
+std::vector<cvm::rvector> cvm::atom_group::calc_fit_forces_impl(const std::vector<cvm::rvector>& forces_on_main_group) const {
+  const cvm::atom_group *group_for_fit = fitting_group ? fitting_group : this;
+  std::vector<cvm::rvector> forces_on_fitting_group(group_for_fit->size());
+  // the center of geometry contribution to the gradients
+  cvm::rvector atom_grad;
+  // the rotation matrix contribution to the gradients
+  const auto rot_inv = rot.inverse().matrix();
+  // temporary variables for computing and summing derivatives
+  cvm::real sum_dxdq[4] = {0, 0, 0, 0};
+  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 += forces_on_main_group[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, forces_on_main_group[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;
+    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++) {
+    if (B_ag_center) {
+      forces_on_fitting_group[j] = 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
+      forces_on_fitting_group[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];
+    }
+  }
+  return forces_on_fitting_group;
+  // TODO
+}
+
+
+std::vector<cvm::rvector> cvm::atom_group::calc_fit_forces(const std::vector<cvm::rvector>& forces_on_main_group) const {
+  if (cvm::debug())
+    cvm::log("Calculating fit forces.\n");
+
+  if (is_enabled(f_ag_center) && is_enabled(f_ag_rotate))
+    return calc_fit_forces_impl<true, true>(forces_on_main_group);
+  if (is_enabled(f_ag_center) && !is_enabled(f_ag_rotate))
+    return calc_fit_forces_impl<true, false>(forces_on_main_group);
+  if (!is_enabled(f_ag_center) && is_enabled(f_ag_rotate))
+    return calc_fit_forces_impl<false, true>(forces_on_main_group);
+  if (!is_enabled(f_ag_center) && !is_enabled(f_ag_rotate))
+    return calc_fit_forces_impl<false, false>(forces_on_main_group);
+
+  if (cvm::debug())
+    cvm::log("Done calculating fit forces.\n");
+  return std::vector<cvm::rvector>();
+}
+
+
 std::vector<cvm::atom_pos> cvm::atom_group::positions() const
 {
   if (b_dummy) {

src/colvaratoms.h

@@ -497,6 +497,8 @@ class colvarmodule::atom_group
   /// \brief Calculate the derivatives of the fitting transformation
   void calc_fit_gradients();
 
+  std::vector<cvm::rvector> calc_fit_forces(const std::vector<cvm::rvector>& forces_on_main_group) const;
+
 /*! @brief  Actual implementation of `calc_fit_gradients`. 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).
@@ -507,6 +509,8 @@ class colvarmodule::atom_group
  */
   template <bool B_ag_center, bool B_ag_rotate> void calc_fit_gradients_impl();
 
+  template <bool B_ag_center, bool B_ag_rotate> std::vector<cvm::rvector> calc_fit_forces_impl(const std::vector<cvm::rvector>& forces_on_main_group) const;
+
   /// \brief Derivatives of the fitting transformation
   std::vector<cvm::atom_pos> fit_gradients;
 

src/colvarcomp_rotations.cpp

@@ -137,10 +137,30 @@ 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;
+    std::vector<cvm::rvector> main_group_forces;
+    const bool force_on_fitting_group = atoms->fitting_group == nullptr ? false : true;
+    cvm::rmatrix ag_rot;
+    if (force_on_fitting_group) {
+      ag_rot = atoms->rot.inverse().matrix();
+    }
     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];
+      if (force_on_fitting_group) {
+        main_group_forces.push_back(f_ia);
+        (*atoms)[ia].apply_force(ag_rot * f_ia);
+      } else {
+        (*atoms)[ia].apply_force(f_ia);
+      }
+    }
+    if (force_on_fitting_group) {
+      const std::vector<cvm::rvector> fitting_group_forces = atoms->calc_fit_forces(main_group_forces);
+      if (fitting_group_forces.empty()) return;
+      for (size_t ia = 0; ia < atoms->fitting_group->size(); ia++) {
+        (*(atoms->fitting_group))[ia].apply_force(fitting_group_forces[ia]);
       }
     }
   }