chore: improve the docs of the calculation of fit gradients

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:

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 <bool B_ag_center, bool B_ag_rotate,
             typename main_force_accessor_T, typename fitting_force_accessor_T>
@@ -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 <typename main_force_accessor_T, typename fitting_force_accessor_T>
   void calc_fit_forces(

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) +