Feature: add gga_grad parameter for noncollinear spin GGA functional

source/source_hamilt/module_xc/CMakeLists.txt

@@ -2,6 +2,7 @@ add_library(
     xc_
     OBJECT
     xc_functional.cpp
+    xc_functional_gga_noncol_sf_builtin.cpp
     xc_pot.cpp
     xc_grad.cpp
     xc_lda_wrap.cpp

source/source_hamilt/module_xc/libxc_abacus.h

@@ -24,30 +24,14 @@ namespace XC_Functional_Libxc
 //-------------------
 
     // sets functional type, which allows combination of LIBXC keyword connected by "+"
-    //        for example, "XC_LDA_X+XC_LDA_C_PZ"
+    //        for example: "XC_LDA_X+XC_LDA_C_PZ"
     extern std::pair<int, std::vector<int>> set_xc_type_libxc(const std::string& xc_func_in);
 
-    /**
-     * @brief instantiate the XC functional by its ID, and set the external parameters if provided.
-     * 
-     * @param func_id libxc ID of functional, see https://libxc.gitlab.io/functionals/ for details
-     * @param xc_polarized 0: unpolarized, 1: spin-polarized
-     * @return std::vector<xc_func_type> 
-     * 
-     * @note the functionality of this method is extended by supporting the user-defined
-     *       external parameters of xc. However, there are several functionals' external 
-     *       parameters are pre-defined in the code, which herein we call those are 
-     *       "in-built" parameters. If the same functional ID is found in both in-built 
-     *       and external parameters, the external parameters will overwrite the in-built ones.
-     *       The external parameters can be passed here by keywords xc_exch_ext and
-     *       xc_corr_ext in the input file. The expected format would be an XC ID 
-     *       followed by a list of parameters.
-     */
     extern std::vector<xc_func_type> init_func(
         const std::vector<int> &func_id,
         const int xc_polarized,
-        const double hybrid_alpha,
-        const double hse_omega);
+        const double hybrid_alpha = 0.0,
+        const double hse_omega = 0.0);
 
     extern void finish_func(std::vector<xc_func_type> &funcs);
 
@@ -65,9 +49,9 @@ namespace XC_Functional_Libxc
         const int nspin,
         const bool domag,
         const bool domag_z,
-        const std::map<int, double>* scaling_factor,
-        const double hybrid_alpha,
-        const double hse_omega);
+        const std::map<int, double>* scaling_factor = nullptr,
+        const double hybrid_alpha = 0.0,
+        const double hse_omega = 0.0);
 
     // for mGGA functional
     extern std::tuple<double, double, ModuleBase::matrix, ModuleBase::matrix> v_xc_meta(
@@ -77,8 +61,8 @@ namespace XC_Functional_Libxc
         const double tpiba,
         const Charge* const chr,
         const int nspin,
-        const double hybrid_alpha,
-        const double hse_omega);
+        const double hybrid_alpha = 0.0,
+        const double hse_omega = 0.0);
 
 
 //-------------------
@@ -157,6 +141,34 @@ namespace XC_Functional_Libxc
         const std::vector<double> &amag,
         const ModuleBase::matrix &v);
 
+    extern std::vector<double> compute_mag_part_nspin4(
+        const std::size_t nrxx,
+        const Charge* const chr);
+
+    extern std::vector<std::vector<ModuleBase::Vector3<double>>> cal_gdr_sf(
+        const int nspin,
+        const std::size_t nrxx,
+        const std::vector<double> &rho,
+        const std::vector<double> &mag_part,
+        const double tpiba,
+        const Charge* const chr);
+
+    extern std::vector<std::vector<double>> cal_dh_sf(
+        const int nspin,
+        const std::size_t nrxx,
+        const std::vector<double> &sgn,
+        const std::vector<std::vector<ModuleBase::Vector3<double>>> &gdr,
+        const std::vector<double> &vsigma,
+        const std::vector<double> &mag_part,
+        const double tpiba,
+        const Charge* const chr);
+
+    extern ModuleBase::matrix convert_v_nspin4_sf(
+        const std::size_t nrxx,
+        const Charge* const chr,
+        const std::vector<double> &mag_part,
+        const ModuleBase::matrix &v);
+
 
 //-------------------
 //  libxc_lda_wrap.cpp
@@ -169,8 +181,8 @@ namespace XC_Functional_Libxc
         double &exc,
         double &vxcup,
         double &vxcdw,
-        const double hybrid_alpha,
-        const double hse_omega);
+        const double hybrid_alpha = 0.0,
+        const double hse_omega = 0.0);
 
 
 //-------------------
@@ -185,8 +197,8 @@ namespace XC_Functional_Libxc
         double &sxc,
         double &v1xc,
         double &v2xc,
-        const double hybrid_alpha,
-        const double hse_omega);
+        const double hybrid_alpha = 0.0,
+        const double hse_omega = 0.0);
 
     // the entire GGA functional, for nspin=2 case
     extern void gcxc_spin_libxc(
@@ -201,8 +213,8 @@ namespace XC_Functional_Libxc
         double &v2xcup,
         double &v2xcdw,
         double &v2xcud,
-        const double hybrid_alpha,
-        const double hse_omega);
+        const double hybrid_alpha = 0.0,
+        const double hse_omega = 0.0);
 
 
 //-------------------
@@ -219,8 +231,8 @@ namespace XC_Functional_Libxc
         double &v1xc,
         double &v2xc,
         double &v3xc,
-        const double &hybrid_alpha,
-        const double &hse_omega);
+        const double &hybrid_alpha = 0.0,
+        const double &hse_omega = 0.0);
 
     extern void tau_xc_spin(
         const std::vector<int> &func_id,
@@ -238,8 +250,8 @@ namespace XC_Functional_Libxc
         double &v2xcud,
         double &v3xcup,
         double &v3xcdw,
-        const double &hybrid_alpha,
-        const double &hse_omega);
+        const double &hybrid_alpha = 0.0,
+        const double &hse_omega = 0.0);
 
 } // namespace XC_Functional_Libxc
 

source/source_hamilt/module_xc/libxc_pot.cpp

@@ -64,9 +64,14 @@ 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;
-    if(1==nspin || 2==nspin_in)
+    std::vector<double> mag_part;
+    if(1==nspin || 2==PARAM.inp.nspin)
     {
         rho = XC_Functional_Libxc::convert_rho(nspin, nrxx, chr);
     }
@@ -75,13 +80,32 @@ 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);
+        }
     }
 
     std::vector<std::vector<ModuleBase::Vector3<double>>> gdr;
     std::vector<double> sigma;
     if(is_gga)
     {
-        gdr = XC_Functional_Libxc::cal_gdr(nspin, nrxx, rho, tpiba, chr);
+        // 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);
+        }
+        else
+        {
+            gdr = XC_Functional_Libxc::cal_gdr(nspin, nrxx, rho, tpiba, chr);
+        }
         sigma = XC_Functional_Libxc::convert_sigma(gdr);
     }
 
@@ -179,7 +203,17 @@ std::tuple<double,double,ModuleBase::matrix> XC_Functional_Libxc::v_xc_libxc(		/
 
     if(4==nspin_in)
     {
-        v = XC_Functional_Libxc::convert_v_nspin4(nrxx, chr, amag, v);
+        // 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);
+        }
+        else
+        {
+            v = XC_Functional_Libxc::convert_v_nspin4(nrxx, chr, amag, v);
+        }
     }
 
     //-------------------------------------------------