Feature: add gga_grad parameter for noncollinear spin GGA functional

docs/advanced/input_files/input-main.md

@@ -1268,6 +1268,18 @@
   - 4: For the case of noncollinear polarized, nspin will be automatically set to 4 without being specified by the user.
 - **Default**: 1
 
+### gga_grad
+
+- **Type**: Integer
+- **Description**: Method used to evaluate the density gradient entering GGA exchange-correlation terms in noncollinear-spin (nspin=4) calculations. With rho_up/dn = (rho +/- |m|)/2 and m_hat = m/|m|:
+  - 0: original algorithm, identical to the behavior before this keyword was introduced; the default.
+  - 1: collinear approximation, grad(rho_up/dn) = (grad(rho) +/- grad(|m|))/2; only the gradient of the magnetization magnitude is used (equivalent to 0).
+  - 2: projected method. The gradients use the full chain rule, grad(rho_up/dn) = (grad(rho) +/- m_hat . grad(m))/2, but the divergence of h = df/d(grad rho) in the potential is projected onto m_hat: v_mu -= m_hat_mu * div((h_up - h_dn)/2), dropping the (h_up - h_dn) . grad(m_hat_mu) cross terms.
+  - 3: Scalmani-Frisch transformation (G. Scalmani and M. J. Frisch, J. Chem. Theory Comput. 8, 2193 (2012)). Same gradients as 2, but the full divergence is kept: v_mu -= div((h_up - h_dn)/2 * m_hat_mu), retaining all cross terms; the most accurate.
+
+  This parameter only takes effect for nspin=4 with GGA functionals (and magnetic calculation).
+- **Default**: 0
+
 ### smearing_method
 
 - **Type**: String

docs/parameters.yaml

@@ -665,6 +665,19 @@ parameters:
     default_value: "1"
     unit: ""
     availability: ""
+  - name: gga_grad
+    category: Electronic structure
+    type: Integer
+    description: |
+      Method used to evaluate the density gradient entering GGA exchange-correlation terms in noncollinear-spin (nspin=4) calculations. With rho_up/dn = (rho +/- |m|)/2 and m_hat = m/|m|:
+      * 0: original algorithm, identical to the behavior before this keyword was introduced; the default.
+      * 1: collinear approximation, grad(rho_up/dn) = (grad(rho) +/- grad(|m|))/2; only the gradient of the magnetization magnitude is used (equivalent to 0).
+      * 2: projected method. The gradients use the full chain rule, grad(rho_up/dn) = (grad(rho) +/- m_hat . grad(m))/2, but the divergence of h = df/d(grad rho) in the potential is projected onto m_hat: v_mu -= m_hat_mu * div((h_up - h_dn)/2), dropping the (h_up - h_dn) . grad(m_hat_mu) cross terms.
+      * 3: Scalmani-Frisch transformation (G. Scalmani and M. J. Frisch, J. Chem. Theory Comput. 8, 2193 (2012)). Same gradients as 2, but the full divergence is kept: v_mu -= div((h_up - h_dn)/2 * m_hat_mu), retaining all cross terms; the most accurate.
+      This parameter only takes effect for nspin=4 with GGA functionals (and magnetic calculation).
+    default_value: "0"
+    unit: ""
+    availability: ""
   - name: smearing_method
     category: Electronic structure
     type: String

source/source_estate/module_pot/pot_xc.cpp

@@ -54,6 +54,7 @@ void PotXC::cal_v_eff(const Charge*const chg, const UnitCell*const ucell, Module
                                   PARAM.inp.nspin,
                                   PARAM.globalv.domag,
                                   PARAM.globalv.domag_z,
+                                  PARAM.inp.gga_grad,
                                   hybrid_alpha,
                                   hse_omega);
         *(this->etxc_) = std::get<0>(etxc_vtxc_v);

source/source_estate/module_pot/pot_xc_fdm.cpp

@@ -31,6 +31,7 @@ PotXC_FDM::PotXC_FDM(
 							  PARAM.inp.nspin,
 							  PARAM.globalv.domag,
 							  PARAM.globalv.domag_z,
+							  PARAM.inp.gga_grad,
 							  hybrid_alpha,
 							  hse_omega);
 	this->v_xc_0 = std::get<2>(etxc_vtxc_v_0);
@@ -69,6 +70,7 @@ void PotXC_FDM::cal_v_eff(
 							  PARAM.inp.nspin,
 							  PARAM.globalv.domag,
 							  PARAM.globalv.domag_z,
+							  PARAM.inp.gga_grad,
 							  hybrid_alpha,
 							  hse_omega);
 	const ModuleBase::matrix &v_xc_01 = std::get<2>(etxc_vtxc_v_01);

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,23 +24,23 @@ 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> 
-     * 
+     * @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 
+     *       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 
+     *       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(
@@ -65,6 +65,7 @@ namespace XC_Functional_Libxc
         const int nspin,
         const bool domag,
         const bool domag_z,
+        const int gga_grad,
         const std::map<int, double>* scaling_factor,
         const double hybrid_alpha,
         const double hse_omega);
@@ -128,6 +129,9 @@ namespace XC_Functional_Libxc
         std::vector<double> exc);
 
     // converting vtxc and v from vrho and vsigma (libxc=>abacus)
+    // use_sf: for nspin=4 magnetic GGA, apply the Scalmani-Frisch
+    // gradient correction instead of the collinear one; gga_grad then
+    // selects the projected (2) or full (3) divergence of h
     extern std::pair<double, ModuleBase::matrix> convert_vtxc_v(
         const xc_func_type &func,
         const int nspin,
@@ -138,7 +142,9 @@ namespace XC_Functional_Libxc
         const std::vector<double> &vrho,
         const std::vector<double> &vsigma,
         const double tpiba,
-        const Charge* const chr);
+        const Charge* const chr,
+        const bool use_sf,
+        const int gga_grad);
 
     // dh for gga v
     extern std::vector<std::vector<double>> cal_dh(
@@ -151,10 +157,42 @@ namespace XC_Functional_Libxc
         const Charge* const chr);
 
     // convert v for NSPIN=4
+    // has_mag: whether the calculation has (noncollinear) magnetization,
+    // i.e. domag || domag_z
     extern ModuleBase::matrix convert_v_nspin4(
         const std::size_t nrxx,
         const Charge* const chr,
         const std::vector<double> &amag,
+        const ModuleBase::matrix &v,
+        const bool has_mag);
+
+    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 int gga_grad,
+        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);
 
 

source/source_hamilt/module_xc/libxc_pot.cpp

@@ -4,7 +4,6 @@
 #include "libxc_abacus.h"
 #include "source_estate/module_charge/charge.h"
 #include "source_base/global_variable.h"
-#include "source_io/module_parameter/parameter.h"
 #include "source_base/parallel_reduce.h"
 #include "source_base/timer.h"
 #include "source_base/tool_title.h"
@@ -25,6 +24,7 @@ std::tuple<double,double,ModuleBase::matrix> XC_Functional_Libxc::v_xc_libxc(		/
         const int nspin_in,
         const bool domag,
         const bool domag_z,
+        const int gga_grad,
         const std::map<int, double>* scaling_factor,
         const double hybrid_alpha,
         const double hse_omega)
@@ -36,6 +36,12 @@ std::tuple<double,double,ModuleBase::matrix> XC_Functional_Libxc::v_xc_libxc(		/
         (nspin_in == 1 || ( nspin_in ==4 && !domag && !domag_z))
         ? 1 : 2;
 
+    // For nspin=4 with noncollinear magnetism, gga_grad=2/3 selects the
+    // Scalmani-Frisch (SF) gradient decomposition; gga_grad=0/1 keeps the
+    // original collinear algorithm.
+    const bool has_mag = domag || domag_z;
+    const bool use_sf = (nspin_in == 4) && has_mag && (gga_grad == 2 || gga_grad == 3);
+
     //----------------------------------------------------------
     // xc_func_type is defined in Libxc package
     // to understand the usage of xc_func_type,
@@ -64,8 +70,13 @@ 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;
     if(1==nspin || 2==nspin_in)
     {
         rho = XC_Functional_Libxc::convert_rho(nspin, nrxx, chr);
@@ -75,13 +86,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=0/1 (else): collinear approximation, uses grad(|m|) only
+        if(use_sf)
+        {
+            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=0/1) for noncollinear magnetism.
+        if(use_sf)
+        {
+            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);
     }
 
@@ -172,14 +202,24 @@ std::tuple<double,double,ModuleBase::matrix> XC_Functional_Libxc::v_xc_libxc(		/
             func, nspin, nrxx,
             sgn, rho, gdr,
             vrho, vsigma,
-            tpiba, chr);
+            tpiba, chr, use_sf, gga_grad);
         vtxc += std::get<0>(vtxc_v) * factor;
         v += std::get<1>(vtxc_v) * factor;
     } // end for( xc_func_type &func : funcs )
 
     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=0/1: standard conversion using |m| decomposition only
+        if(use_sf)
+        {
+            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, has_mag);
+        }
     }
 
     //-------------------------------------------------