From d206bbd62507a3160bb3e548dfd1d690a7879cb1 Mon Sep 17 00:00:00 2001 From: HanatoK Date: Thu, 17 Oct 2024 15:27:03 -0500 Subject: [PATCH] chore: improve the docs of the calculation of fit gradients --- doc/doxygen/Doxyfile | 2 +- src/colvaratoms.h | 53 ++++++++++++++++++++++++++++++++++++++++++-- src/colvartypes.h | 36 ++++++++++++++++++++++++++++++ 3 files changed, 88 insertions(+), 3 deletions(-) diff --git a/doc/doxygen/Doxyfile b/doc/doxygen/Doxyfile index 18d638d0e..df85e5377 100644 --- a/doc/doxygen/Doxyfile +++ b/doc/doxygen/Doxyfile @@ -1420,7 +1420,7 @@ FORMULA_TRANSPARENT = YES # The default value is: NO. # This tag requires that the tag GENERATE_HTML is set to YES. -USE_MATHJAX = NO +USE_MATHJAX = YES # When MathJax is enabled you can set the default output format to be used for # the MathJax output. See the MathJax site (see: diff --git a/src/colvaratoms.h b/src/colvaratoms.h index 65a6d4a17..528e849df 100644 --- a/src/colvaratoms.h +++ b/src/colvaratoms.h @@ -264,10 +264,46 @@ class colvarmodule::atom_group public: + /*! @class group_force_object + * @brief A helper class for applying forces on an atom group in a way that + * is aware of the fitting group. NOTE: you are encouraged to use + * get_group_force_object() to get an instance of group_force_object + * instead of constructing directly. + */ class group_force_object { public: + /*! @brief Constructor of group_force_object + * @param ag The pointer to the atom group that forces will be applied on. + */ group_force_object(cvm::atom_group* ag); + /*! @brief Destructor of group_force_object + */ ~group_force_object(); + /*! @brief Apply force to atom i + * @param i The i-th of atom in the atom group. + * @param force The force being added to atom i. + * + * The function can be used as follows, + * @code + * // In your colvar::cvc::apply_force() loop of a component: + * auto ag_force = atoms->get_group_force_object(); + * for (ia = 0; ia < atoms->size(); ia++) { + * const cvm::rvector f = compute_force_on_atom_ia(); + * ag_force.add_atom_force(ia, f); + * } + * @endcode + * There are actually two scenarios under the hood: + * (i) If the atom group does not have a fitting group, then the force is + * added to atom i directly; + * (ii) If the atom group has a fitting group, the force on atom i will just + * be temporary stashed into ag->group_forces. At the end of the loop + * of apply_force(), the destructor ~group_force_object() will be called, + * which then call apply_force_with_fitting_group(). The forces on the + * main group will be rotated back by multiplying ag->group_forces with + * the inverse rotation. The forces on the fitting group (if + * enableFitGradients is on) will be calculated by calling + * calc_fit_forces. + */ void add_atom_force(size_t i, const cvm::rvector& force); private: cvm::atom_group* m_ag; @@ -528,15 +564,25 @@ class colvarmodule::atom_group * @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. + * group forces or gradients acting on the rotated frame. * @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. + * rotated 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. + * + * This function is used to (i) project the gradients of CV with respect to + * rotated main group atoms to fitting group atoms, or (ii) project the forces + * on rotated main group atoms to fitting group atoms, by the following two steps + * (using the goal (ii) for example): + * (1) Loop over the positions of main group atoms and call cvm::quaternion::position_derivative_inner + * to project the forces on rotated main group atoms to the forces on quaternion. + * (2) Loop over the positions of fitting group atoms, compute the gradients of + * \f$\mathbf{q}\f$ with respect to the position of each atom, and then multiply + * that with the force on \f$\mathbf{q}\f$ (chain rule). */ template @@ -555,6 +601,9 @@ class colvarmodule::atom_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. + * + * This function just dispatches the parameters to calc_fit_forces_impl that really + * performs the calculations. */ template void calc_fit_forces( diff --git a/src/colvartypes.h b/src/colvartypes.h index f2ec6eaf0..07a7065e5 100644 --- a/src/colvartypes.h +++ b/src/colvartypes.h @@ -1221,6 +1221,42 @@ class colvarmodule::quaternion { /// \brief Multiply the given vector by the derivative of the given /// (rotated) position with respect to the quaternion + /// \param pos The position \f$\mathbf{x}\f$. + /// \param vec The vector \f$\mathbf{v}\f$. + /// \return A quaternion (see the detailed documentation below). + /// + /// This function is mainly used for projecting the gradients or forces on + /// the rotated atoms to the forces on quaternion. Assume this rotation can + /// be represented as \f$R(\mathbf{q})\f$, + /// where \f$\mathbf{q} := (q_0, q_1, q_2, q_3)\f$ + /// is the current quaternion, the function returns the following new + /// quaternion: + /// \f[ + /// \left(\mathbf{v}^\mathrm{T}\frac{\partial R(\mathbf{q})}{\partial q_0}\mathbf{x}, + /// \mathbf{v}^\mathrm{T}\frac{\partial R(\mathbf{q})}{\partial q_1}\mathbf{x}, + /// \mathbf{v}^\mathrm{T}\frac{\partial R(\mathbf{q})}{\partial q_2}\mathbf{x}, + /// \mathbf{v}^\mathrm{T}\frac{\partial R(\mathbf{q})}{\partial q_3}\mathbf{x}\right) + /// \f] + /// where \f$\mathbf{v}\f$ is usually the gradient of \f$\xi\f$ with respect to + /// the rotated frame \f$\tilde{\mathbf{X}}\f$, + /// \f$\partial \xi / \partial \tilde{\mathbf{X}}\f$, or the force acting on it + /// (\f$\mathbf{F}_{\tilde{\mathbf{X}}}\f$). + /// By using the following loop in pseudo C++ code, + /// either \f$\partial \xi / \partial \tilde{\mathbf{X}}\f$ + /// or \f$\mathbf{F}_{\tilde{\mathbf{X}}}\f$, can be projected to + /// \f$\partial \xi / \partial \mathbf{q}\f$ or \f$\mathbf{F}_q\f$ into `sum_dxdq`: + /// @code + /// cvm::real sum_dxdq[4] = {0, 0, 0, 0}; + /// for (size_t i = 0; i < main_group_size(); ++i) { + /// const cvm::rvector v = grad_or_force_on_rotated_main_group(i); + /// const cvm::rvector x = unrotated_main_group_positions(i); + /// cvm::quaternion const dxdq = position_derivative_inner(x, v); + /// sum_dxdq[0] += dxdq[0]; + /// sum_dxdq[1] += dxdq[1]; + /// sum_dxdq[2] += dxdq[2]; + /// sum_dxdq[3] += dxdq[3]; + /// } + /// @endcode inline cvm::quaternion position_derivative_inner(cvm::rvector const &pos, cvm::rvector const &vec) const { return cvm::quaternion(2.0 * (vec.x * ( q0 * pos.x - q3 * pos.y + q2 * pos.z) +