From 9cae417168f8104976b528d45d508c589b676f88 Mon Sep 17 00:00:00 2001 From: HanatoK Date: Wed, 11 Sep 2024 21:57:51 -0500 Subject: [PATCH] fix: try to solve #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 | 73 ++++++++++++++++++++++++++++++++++++ src/colvaratoms.h | 4 ++ src/colvarcomp_rotations.cpp | 37 +++++++++++------- 3 files changed, 100 insertions(+), 14 deletions(-) diff --git a/src/colvaratoms.cpp b/src/colvaratoms.cpp index e15b9301a..f8f6c92f9 100644 --- a/src/colvaratoms.cpp +++ b/src/colvaratoms.cpp @@ -1282,6 +1282,79 @@ void cvm::atom_group::calc_fit_gradients_impl() { } +template +std::vector cvm::atom_group::calc_fit_forces_impl(const std::vector& forces_on_main_group) const { + const cvm::atom_group *group_for_fit = fitting_group ? fitting_group : this; + std::vector 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 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::atom_group::calc_fit_forces(const std::vector& 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(forces_on_main_group); + if (is_enabled(f_ag_center) && !is_enabled(f_ag_rotate)) + return calc_fit_forces_impl(forces_on_main_group); + if (!is_enabled(f_ag_center) && is_enabled(f_ag_rotate)) + return calc_fit_forces_impl(forces_on_main_group); + if (!is_enabled(f_ag_center) && !is_enabled(f_ag_rotate)) + return calc_fit_forces_impl(forces_on_main_group); + + if (cvm::debug()) + cvm::log("Done calculating fit forces.\n"); + return std::vector(); +} + + std::vector cvm::atom_group::positions() const { if (b_dummy) { diff --git a/src/colvaratoms.h b/src/colvaratoms.h index d16ca7bd5..613a36164 100644 --- 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 calc_fit_forces(const std::vector& 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 void calc_fit_gradients_impl(); + template std::vector calc_fit_forces_impl(const std::vector& forces_on_main_group) const; + /// \brief Derivatives of the fitting transformation std::vector fit_gradients; diff --git a/src/colvarcomp_rotations.cpp b/src/colvarcomp_rotations.cpp index 9f401a5f2..f9c265808 100644 --- 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 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 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 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]); } } }