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.
HanatoK · Oct 11, 2024 · 38fa427 · 38fa427
1 parent 37fb8df
commit 38fa427
Show file tree

Hide file tree

Showing 3 changed files with 100 additions and 14 deletions.
diff --git a/src/colvaratoms.cpp b/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) {

diff --git a/src/colvaratoms.h b/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;
 

diff --git a/src/colvarcomp_rotations.cpp b/src/colvarcomp_rotations.cpp
@@ -137,21 +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;
-
-    if (atoms->is_enabled(f_ag_rotate)) {
-      auto const rot_inv = 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(rot_inv * (FQ[i] * dq0_2[i]));
-        }
+    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);
+      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);
       }
-    } else {
-      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]);
-        }
+    }
+    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]);
       }
     }
   }