feat(xc): add gga_grad 1/2 with variationally consistent mode 2

docs/advanced/input_files/input-main.md

@@ -95,6 +95,7 @@
     - [pseudo\_rcut](#pseudo_rcut)
     - [pseudo\_mesh](#pseudo_mesh)
     - [nspin](#nspin)
+    - [gga\_grad](#gga_grad)
     - [smearing\_method](#smearing_method)
     - [smearing\_sigma](#smearing_sigma)
     - [smearing\_sigma\_temp](#smearing_sigma_temp)
@@ -1353,6 +1354,17 @@
   - 4: Noncollinear or spin-orbit calculations. Set nspin to 4 explicitly when noncolin or lspinorb is enabled.
 - **Default**: 1
 
+### gga_grad
+
+- **Type**: Integer
+- **Description**: Selects the local spin mapping for LDA/GGA functionals in magnetic nspin=4 calculations.
+  - 0: preserves the original algorithm (default).
+  - 1: uses the local magnetization magnitude instead of the global quantization axis in the built-in GGA gradient correction. For LIBXC functionals, 0 and 1 are equivalent.
+  - 2: uses a C2-regularized magnetization magnitude with eta = 1e-3 in atomic density units. The spin densities are (abs(n + rho_core) +/- min(S_eta(m), abs(n + rho_core)))/2. GGA gradients are the local-map Jacobian applied to the FFT gradients of the four density channels. The potential reverses this same discrete energy graph, including the radial Hessian and density/sigma clipping branches; the GGA stress uses the corresponding metric derivative.
+  For r = |m| and x = r/eta, S_eta = eta*x^3*(3*x^2 - 8*x + 6) for r < eta, and S_eta = r otherwise. The regularization is part of the functional definition, including its first and second derivatives.
+  Mode 2 also uses this local map for the LDA contribution. Other spin configurations retain their existing behavior.
+- **Default**: 0
+
 ### smearing_method
 
 - **Type**: String

docs/parameters.yaml

@@ -712,6 +712,19 @@ parameters:
     default_value: "1"
     unit: ""
     availability: ""
+  - name: gga_grad
+    category: Electronic structure
+    type: Integer
+    description: |
+      Selects the local spin mapping for LDA/GGA functionals in magnetic nspin=4 calculations.
+      * 0: preserves the original algorithm (default).
+      * 1: uses the local magnetization magnitude instead of the global quantization axis in the built-in GGA gradient correction. For LIBXC functionals, 0 and 1 are equivalent.
+      * 2: uses a C2-regularized magnetization magnitude with eta = 1e-3 in atomic density units. The spin densities are (abs(n + rho_core) +/- min(S_eta(m), abs(n + rho_core)))/2. GGA gradients are the local-map Jacobian applied to the FFT gradients of the four density channels. The potential reverses this same discrete energy graph, including the radial Hessian and density/sigma clipping branches; the GGA stress uses the corresponding metric derivative.
+      For r = |m| and x = r/eta, S_eta = eta*x^3*(3*x^2 - 8*x + 6) for r < eta, and S_eta = r otherwise. The regularization is part of the functional definition, including its first and second derivatives.
+      Mode 2 also uses this local map for the LDA contribution. Other spin configurations retain their existing behavior.
+    default_value: "0"
+    unit: ""
+    availability: ""
   - name: smearing_method
     category: Electronic structure
     type: String

source/Makefile.Objects

@@ -589,6 +589,8 @@ OBJS_SYMMETRY=symm_other.o\
     symmetry.o\
 
 OBJS_XC=xc_functional.o\
+    xc_functional_ncgga_sf.o\
+    xc_ncgga_radial.o\
     xc_functional_op.o\
     xc_pot.o\
     xc_grad.o\

source/source_estate/module_pot/pot_xc.cpp

@@ -14,9 +14,11 @@ namespace elecstate
 
 void PotXC::cal_v_eff(const Charge*const chg, const UnitCell*const ucell, ModuleBase::matrix& v_eff)
 {
+    const Parameter& parameters = PARAM;
     ModuleBase::TITLE("PotXC", "cal_veff");
     ModuleBase::timer::start("PotXC", "cal_veff");
     const int nrxx_current = chg->nrxx;
+    const int nspin = parameters.inp.nspin;
 
     //----------------------------------------------------------
     //  calculate the exchange-correlation potential
@@ -33,7 +35,7 @@ void PotXC::cal_v_eff(const Charge*const chg, const UnitCell*const ucell, Module
 #endif
         const std::tuple<double, double, ModuleBase::matrix, ModuleBase::matrix> etxc_vtxc_v
             = XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), nrxx_current, ucell->omega, ucell->tpiba, chg,
-                                             PARAM.inp.nspin, hybrid_alpha, hse_omega);
+                                             nspin, hybrid_alpha, hse_omega);
         *(this->etxc_) = std::get<0>(etxc_vtxc_v);
         *(this->vtxc_) = std::get<1>(etxc_vtxc_v);
         v_eff += std::get<2>(etxc_vtxc_v);
@@ -52,9 +54,10 @@ void PotXC::cal_v_eff(const Charge*const chg, const UnitCell*const ucell, Module
 #endif
         const std::tuple<double, double, ModuleBase::matrix> etxc_vtxc_v
             = XC_Functional::v_xc(nrxx_current, chg, ucell,
-                                  PARAM.inp.nspin,
-                                  PARAM.globalv.domag,
-                                  PARAM.globalv.domag_z,
+                                  nspin,
+                                  parameters.globalv.domag,
+                                  parameters.globalv.domag_z,
+                                  parameters.inp.gga_grad,
                                   hybrid_alpha,
                                   hse_omega);
         *(this->etxc_) = std::get<0>(etxc_vtxc_v);

source/source_estate/module_pot/pot_xc_fdm.cpp

@@ -16,6 +16,7 @@ PotXC_FDM::PotXC_FDM(
 	const UnitCell*const ucell)
 	: chg_0(chg_0_in)
 {
+    const Parameter& parameters = PARAM;
 	this->rho_basis_ = rho_basis_in;
 	this->dynamic_mode = true;
 	this->fixed_mode = false;
@@ -28,9 +29,10 @@ PotXC_FDM::PotXC_FDM(
 #endif
 	const std::tuple<double, double, ModuleBase::matrix> etxc_vtxc_v_0
 		= XC_Functional::v_xc(this->chg_0->nrxx, this->chg_0, ucell,
-							  PARAM.inp.nspin,
-							  PARAM.globalv.domag,
-							  PARAM.globalv.domag_z,
+							  parameters.inp.nspin,
+							  parameters.globalv.domag,
+							  parameters.globalv.domag_z,
+							  parameters.inp.gga_grad,
 							  hybrid_alpha,
 							  hse_omega);
 	this->v_xc_0 = std::get<2>(etxc_vtxc_v_0);
@@ -41,6 +43,7 @@ void PotXC_FDM::cal_v_eff(
 	const UnitCell*const ucell,
 	ModuleBase::matrix& v_eff)
 {
+    const Parameter& parameters = PARAM;
 	ModuleBase::TITLE("PotXC_FDM", "cal_veff");
 	ModuleBase::timer::start("PotXC_FDM", "cal_veff");
 
@@ -66,9 +69,10 @@ void PotXC_FDM::cal_v_eff(
 #endif
 	const std::tuple<double, double, ModuleBase::matrix> etxc_vtxc_v_01
 		= XC_Functional::v_xc(chg_01.nrxx, &chg_01, ucell,
-							  PARAM.inp.nspin,
-							  PARAM.globalv.domag,
-							  PARAM.globalv.domag_z,
+							  parameters.inp.nspin,
+							  parameters.globalv.domag,
+							  parameters.globalv.domag_z,
+							  parameters.inp.gga_grad,
 							  hybrid_alpha,
 							  hse_omega);
 	const ModuleBase::matrix &v_xc_01 = std::get<2>(etxc_vtxc_v_01);

source/source_hamilt/module_surchem/test/CMakeLists.txt

@@ -29,6 +29,7 @@ AddTest(
   LIBS parameter planewave device base container
   SOURCES cal_vcav_test.cpp  ../cal_vcav.cpp ../surchem.cpp
   ../../module_xc/xc_grad.cpp ../../module_xc/xc_grad_prepare.cpp ../../module_xc/xc_grad_kernel.cpp ../../module_xc/xc_grad_assemble.cpp ../../module_xc/xc_grad_wfc.cpp ../../module_xc/xc_grad_utils.cpp ../../module_xc/xc_functional.cpp
+  ../../module_xc/xc_functional_ncgga_sf.cpp ../../module_xc/xc_ncgga_radial.cpp
   ../../module_xc/xc_lda_wrap.cpp ../../module_xc/xc_gga_wrap.cpp
   ../../module_xc/libxc_setup.cpp
   ../../module_xc/libxc_pot.cpp
@@ -44,6 +45,7 @@ AddTest(
   LIBS parameter planewave device base container
   SOURCES cal_vel_test.cpp  ../cal_vel.cpp ../surchem.cpp ../cal_epsilon.cpp ../minimize_cg.cpp
   ../../module_xc/xc_grad.cpp ../../module_xc/xc_grad_prepare.cpp ../../module_xc/xc_grad_kernel.cpp ../../module_xc/xc_grad_assemble.cpp ../../module_xc/xc_grad_wfc.cpp ../../module_xc/xc_grad_utils.cpp ../../module_xc/xc_functional.cpp
+  ../../module_xc/xc_functional_ncgga_sf.cpp ../../module_xc/xc_ncgga_radial.cpp
   ../../module_xc/xc_lda_wrap.cpp ../../module_xc/xc_gga_wrap.cpp
   ../../module_xc/libxc_setup.cpp
   ../../module_xc/libxc_pot.cpp

source/source_hamilt/module_xc/CMakeLists.txt

@@ -2,6 +2,8 @@ add_library(
     xc_
     OBJECT
     xc_functional.cpp
+    xc_functional_ncgga_sf.cpp
+    xc_ncgga_radial.cpp
     xc_pot.cpp
     xc_grad.cpp
     xc_grad_prepare.cpp

source/source_hamilt/module_xc/libxc_abacus.h

@@ -5,10 +5,12 @@
 
 #include "source_base/matrix.h"
 #include "source_base/vector3.h"
+#include "xc_ncgga_radial.h"
 
 #include <xc.h>
 #include <xc_funcs.h>
 
+#include <array>
 #include <tuple>
 #include <vector>
 
@@ -19,28 +21,48 @@ class Charge;
 
 namespace XC_Functional_Libxc
 {
+    struct LibxcWeightedDerivatives
+    {
+        double energy_sum;
+        std::vector<double> drho;
+        std::vector<double> dsigma;
+    };
+
+    // Complete forward data for the gga_grad=2 noncollinear Libxc graph:
+    //   rho_s = N_s(x),
+    //   g_s   = sum_A (d N_s / d x_A) G_h x_A.
+    // Keeping the local map and all input gradients together lets the reverse
+    // use the exact same branch choices and radial Hessian as the forward.
+    struct NclSfDiscreteData
+    {
+        std::vector<ModuleXC::NcggaSpinMapPoint> spin_map;
+        std::vector<double> rho;
+        std::vector<std::vector<ModuleBase::Vector3<double>>> spin_gradient;
+        std::array<std::vector<ModuleBase::Vector3<double>>, 3> grad_m;
+    };
+
 //-------------------
 //  libxc_setup.cpp
 //-------------------
 
     // 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,10 +87,24 @@ 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);
 
+    // Reciprocal-metric derivative of the exact gga_grad=2 Libxc energy
+    // graph. The returned lower-triangular tensor is the unnormalized local
+    // grid sum; Stress_Func performs the pool reduction and divides by nxyz.
+    extern void gradcorr_ncgga_sf_libxc(
+        const std::vector<int>& func_id,
+        const std::size_t nrxx,
+        const double tpiba,
+        const Charge* const chr,
+        const std::map<int, double>* scaling_factor,
+        const double hybrid_alpha,
+        const double hse_omega,
+        std::vector<double>& stress_gga);
+
     // for mGGA functional
     extern std::tuple<double, double, ModuleBase::matrix, ModuleBase::matrix> v_xc_meta(
         const std::vector<int> &func_id,
@@ -97,6 +133,25 @@ namespace XC_Functional_Libxc
         const std::size_t nrxx,
         const Charge* const chr);
 
+    // Build the exact gga_grad=2 local spin map and, when requested, its
+    // projected FFT-gradient graph.  LDA-only callers set need_gradient=false.
+    extern NclSfDiscreteData make_ncl_sf_discrete_data(
+        const std::size_t nrxx,
+        const double tpiba,
+        const Charge* const chr,
+        const bool need_gradient);
+
+    // Reverse one aggregate of all scaled Libxc components.  The returned
+    // potential is already in (n,mx,my,mz) representation.  An empty dsigma
+    // selects the LDA-only local reverse and performs no FFT divergence.
+    extern ModuleBase::matrix reverse_ncl_sf_discrete(
+        const std::size_t nrxx,
+        const NclSfDiscreteData& data,
+        const std::vector<double>& drho,
+        const std::vector<double>& dsigma,
+        const double tpiba,
+        const Charge* const chr);
+
     // calculating grho
     extern std::vector<std::vector<ModuleBase::Vector3<double>>> cal_gdr(
         const int nspin,
@@ -137,7 +192,21 @@ namespace XC_Functional_Libxc
         const std::vector<double> &rho,
         std::vector<double> exc);
 
-    // converting vtxc and v from vrho and vsigma (libxc=>abacus)
+    // Reverse the density and sigma sanitizers for the weighted energy
+    // accumulated by ABACUS. The result is in Hartree units and excludes the
+    // real-space grid weight and ModuleBase::e2.
+    extern LibxcWeightedDerivatives make_libxc_weighted_derivatives(
+        const xc_func_type &func,
+        const int nspin,
+        const std::size_t nrxx,
+        const std::vector<double> &sgn,
+        const std::vector<double> &rho,
+        const std::vector<double> &sigma,
+        const std::vector<double> &exc,
+        const std::vector<double> &vrho,
+        const std::vector<double> &vsigma);
+
+    // Convert collinear LibXC derivatives to the potential.
     extern std::pair<double, ModuleBase::matrix> convert_vtxc_v(
         const xc_func_type &func,
         const int nspin,
@@ -161,14 +230,14 @@ 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 domag,
-        const bool domag_z);
-
+        const bool has_mag);
 
 //-------------------
 //  libxc_lda_wrap.cpp