docs(gga_grad): add comments explaining SF formulas and branching logic

- libxc_pot.cpp: comments for SF gradient decomposition, gga_grad branching
- libxc_tools.cpp: comment for nspin=4 SF dh path with 4-component potential
- xc_pot.cpp: comment for non-libxc SF builtin dispatch
- xc_functional_gga_noncol_sf_builtin.cpp: header comment block with
  Scalmani-Frisch reference, chain-rule formulas, and step-by-step labels

Author: dyzheng

source/source_hamilt/module_xc/libxc_pot.cpp

@@ -64,6 +64,10 @@ std::tuple<double,double,ModuleBase::matrix> XC_Functional_Libxc::v_xc_libxc(		/
     }();
 
     // converting rho
+    // For nspin=4, the charge density has 4 components:
+    //   rho[0] = total charge, rho[1..3] = magnetization (mx, my, mz)
+    // libxc works with spin-up/spin-down densities:
+    //   rho_up = 0.5*(rho[0] + |m|), rho_dn = 0.5*(rho[0] - |m|)
     std::vector<double> rho;
     std::vector<double> amag;
     std::vector<double> mag_part;
@@ -76,6 +80,11 @@ std::tuple<double,double,ModuleBase::matrix> XC_Functional_Libxc::v_xc_libxc(		/
         std::tuple<std::vector<double>,std::vector<double>> rho_amag = XC_Functional_Libxc::convert_rho_amag_nspin4(nspin, nrxx, chr);
         rho = std::get<0>(std::move(rho_amag));
         amag = std::get<1>(std::move(rho_amag));
+        // gga_grad=2/3: compute magnetization unit vector m_hat = m/|m|
+        // needed for the Scalmani-Frisch (SF) gradient decomposition:
+        //   grad(rho_up)   = 0.5 grad(rho[0]) + 0.5 m_hat_mu * grad(m_mu)
+        //   grad(rho_dn)   = 0.5 grad(rho[0]) - 0.5 m_hat_mu * grad(m_mu)
+        // gga_grad=1 (else): collinear approximation, uses grad(|m|) only
         if((PARAM.inp.gga_grad == 2 || PARAM.inp.gga_grad == 3) && (PARAM.globalv.domag || PARAM.globalv.domag_z))
         {
             mag_part = XC_Functional_Libxc::compute_mag_part_nspin4(nrxx, chr);
@@ -86,6 +95,9 @@ std::tuple<double,double,ModuleBase::matrix> XC_Functional_Libxc::v_xc_libxc(		/
     std::vector<double> sigma;
     if(is_gga)
     {
+        // gga_grad=2/3: use SF method to compute spin-up/spin-down gradients
+        // via the chain-rule decomposition. This is more accurate than the
+        // collinear approximation (gga_grad=1) for noncollinear magnetism.
         if(PARAM.inp.nspin==4 && (PARAM.inp.gga_grad == 2 || PARAM.inp.gga_grad == 3) && (PARAM.globalv.domag || PARAM.globalv.domag_z))
         {
             gdr = XC_Functional_Libxc::cal_gdr_sf(nspin, nrxx, rho, mag_part, tpiba, chr);
@@ -191,6 +203,9 @@ std::tuple<double,double,ModuleBase::matrix> XC_Functional_Libxc::v_xc_libxc(		/
 
     if(4==nspin_in)
     {
+        // gga_grad=2/3: convert libxc spin-up/spin-down potential back to
+        // nspin=4 representation via SF method (v_total, v_mag_hat)
+        // gga_grad=1: standard conversion using |m| decomposition only
         if((PARAM.inp.gga_grad == 2 || PARAM.inp.gga_grad == 3) && (PARAM.globalv.domag || PARAM.globalv.domag_z))
         {
             v = XC_Functional_Libxc::convert_v_nspin4_sf(nrxx, chr, mag_part, v);

source/source_hamilt/module_xc/libxc_tools.cpp

@@ -235,6 +235,11 @@ std::pair<double,ModuleBase::matrix> XC_Functional_Libxc::convert_vtxc_v(
 
 	if(func.info->family == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA)
 	{
+		// For nspin=4 with domag and gga_grad=2/3, use the SF (Scalmani-Frisch)
+		// gradient correction. The potential v has 4 components:
+		//   v(0,:) = total density channel,  v(1..3,:) = magnetization channels
+		// cal_dh_sf returns dh[4][nrxx] in the same 4-component representation.
+		// The vtxc contribution sums dh[is] * rho[is] for each channel.
 		if(PARAM.inp.nspin==4 && (PARAM.inp.gga_grad == 2 || PARAM.inp.gga_grad == 3) && (PARAM.globalv.domag || PARAM.globalv.domag_z))
 		{
 			std::vector<double> mag_part_tmp = XC_Functional_Libxc::compute_mag_part_nspin4(nrxx, chr);

source/source_hamilt/module_xc/xc_functional_gga_noncol_sf_builtin.cpp

@@ -25,6 +25,26 @@ std::tuple<double, double, ModuleBase::matrix> v_xc_ncgga_sf_builtin(
     if (PARAM.inp.nspin != 4 || (!PARAM.globalv.domag && !PARAM.globalv.domag_z))
         throw std::domain_error("v_xc_ncgga_sf_builtin requires NSPIN==4.");
 
+    // ======================================================================
+    // Scalmani-Frisch (SF) builtin for noncollinear GGA (gga_grad=3)
+    //
+    // Reference: Scalmani & Frisch, JCTC 8, 1069 (2012)
+    //
+    // Key formulas:
+    //   rho_up   = 0.5*(rho + |m|),    rho_dn = 0.5*(rho - |m|)
+    //   m_hat    = m / |m|             (magnetization unit vector)
+    //
+    // Gradient decomposition (chain rule via m_hat):
+    //   grad(rho_up) = 0.5 grad(rho) + 0.5 m_hat_mu * grad(m_mu)
+    //   grad(rho_dn) = 0.5 grad(rho) - 0.5 m_hat_mu * grad(m_mu)
+    //
+    // For GGA, the sigma = |grad(rho_up)|^2 etc. are computed from gdr1,gdr2,
+    // then passed to spin-polarized XC functionals (xc_spin, gcx_spin, gcc_spin).
+    // The resulting potential is converted back to nspin=4 representation:
+    //   v_tot = 0.5*(v_up + v_dn)
+    //   v_mag = 0.5*(v_up - v_dn) * m_hat
+    // ======================================================================
+
     ModulePW::PW_Basis* rhopw = chr->rhopw;
     const int npw = rhopw->npw;
     const double e2 = ModuleBase::e2;
@@ -33,6 +53,7 @@ std::tuple<double, double, ModuleBase::matrix> v_xc_ncgga_sf_builtin(
     const double fac = 0.5;
     const bool is_gga = (XC_Functional::get_func_type() == 2 || XC_Functional::get_func_type() == 4);
 
+    // Step 1: decompose charge density into spin-up/spin-down and compute m_hat
     std::vector<double> rhotmp1(nrxx), rhotmp2(nrxx), amag(nrxx);
     std::vector<double> mag_part(3 * nrxx, 0.0);
 
@@ -83,6 +104,9 @@ std::tuple<double, double, ModuleBase::matrix> v_xc_ncgga_sf_builtin(
         }
     }
 
+    // Step 3: LDA contribution (xc_spin) and convert to nspin=4 potential
+    //   v(0,:) = 0.5*(vxc[0] + vxc[1])    -- total density channel
+    //   v(1..3,:) = 0.5*(vxc[0] - vxc[1]) * m_mu / |m|    -- magnetization channels
     double etxc = 0, vtxc = 0;
     ModuleBase::matrix v(PARAM.inp.nspin, nrxx);
 
@@ -114,6 +138,9 @@ std::tuple<double, double, ModuleBase::matrix> v_xc_ncgga_sf_builtin(
         etxc += e2 * exc * arho;
     }
 
+    // Step 4: GGA contribution (gcx_spin + gcc_spin) and divergence correction (grad_dot)
+    //   The gradient correction uses the same SF-decomposed gdr1, gdr2 and converts
+    //   back to nspin=4 via m_hat, matching the LDA conversion above.
     if (is_gga)
     {
         double etxcgc = 0, vtxcgc = 0;

source/source_hamilt/module_xc/xc_pot.cpp

@@ -52,6 +52,10 @@ std::tuple<double, double, ModuleBase::matrix> XC_Functional::v_xc(
 #endif
     }
 
+    // For non-libxc builds: gga_grad=3 uses the SF builtin that computes the
+    // full nspin=4 GGA potential via the Scalmani-Frisch transformation.
+    // This path handles spin-up/spin-down decomposition and gradient
+    // corrections internally, returning (etxc, vtxc, v) directly.
     if (PARAM.inp.nspin == 4 && (PARAM.globalv.domag || PARAM.globalv.domag_z) && PARAM.inp.gga_grad == 3)
     {
         return ModuleXC::NCGGA_SF_Builtin::v_xc_ncgga_sf_builtin(nrxx, ucell->omega, ucell->tpiba, chr);