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]); } } }