refactor(xc): pass spin controls explicitly to force and stress helpers

Address review feedback by removing PARAM reads from cal_force_cc, stress_cc and stress_gga and updating PW, LCAO, SDFT and OFDFT callers. Name the legacy gga_grad test input explicitly. Production parameter values and numerical formulas are unchanged; no INPUT documentation changes are required.

Validation: Sai DSPRHBM job 1196207 rebuilt abacus_basic_para and four XC test executables. Five CTest entries passed, including serial and MPI2 NCGGA finite differences. CUDA validation remains for CI.

Author: aboys-cb

source/source_hamilt/module_xc/test/test_xc3.cpp

@@ -90,15 +90,16 @@ class XCTest_GRADCORR : public XCTest
 
             XC_Functional::set_xc_type("PBE");
 
+            const int gga_grad = 0;
             double hybrid_alpha = 0.0;
             double hse_omega = 0.0;
-            XC_Functional::gradcorr(et1,vt1,v1,&chr,&rhopw,&ucell,stress1,false,nspin1,domag,domag_z,0, hybrid_alpha, hse_omega);
-            XC_Functional::gradcorr(et1,vt1,v1,&chr,&rhopw,&ucell,stress1,true,nspin1,domag,domag_z,0, hybrid_alpha, hse_omega);
+            XC_Functional::gradcorr(et1,vt1,v1,&chr,&rhopw,&ucell,stress1,false,nspin1,domag,domag_z,gga_grad, hybrid_alpha, hse_omega);
+            XC_Functional::gradcorr(et1,vt1,v1,&chr,&rhopw,&ucell,stress1,true,nspin1,domag,domag_z,gga_grad, hybrid_alpha, hse_omega);
 
-            XC_Functional::gradcorr(et2,vt2,v2,&chr,&rhopw,&ucell,stress2,false,nspin2,domag,domag_z,0, hybrid_alpha, hse_omega);
-            XC_Functional::gradcorr(et2,vt2,v2,&chr,&rhopw,&ucell,stress2,true,nspin2,domag,domag_z,0, hybrid_alpha, hse_omega);
+            XC_Functional::gradcorr(et2,vt2,v2,&chr,&rhopw,&ucell,stress2,false,nspin2,domag,domag_z,gga_grad, hybrid_alpha, hse_omega);
+            XC_Functional::gradcorr(et2,vt2,v2,&chr,&rhopw,&ucell,stress2,true,nspin2,domag,domag_z,gga_grad, hybrid_alpha, hse_omega);
 
-            XC_Functional::gradcorr(et4,vt4,v4,&chr,&rhopw,&ucell,stress4,false,nspin4,domag_true,domag_z,0, hybrid_alpha, hse_omega);
+            XC_Functional::gradcorr(et4,vt4,v4,&chr,&rhopw,&ucell,stress4,false,nspin4,domag_true,domag_z,gga_grad, hybrid_alpha, hse_omega);
         }
 };
 

source/source_lcao/force_stress_lcao.cpp

@@ -921,14 +921,17 @@ void Force_Stress_LCAO<T>::calForcePwPart(UnitCell& ucell,
                                           const pseudopot_cell_vl& locpp,
                                           const Structure_Factor& sf)
 {
+    const auto& xc_input = PARAM.inp;
+    const auto& xc_spin = PARAM.globalv;
     ModuleBase::TITLE("Force_Stress_LCAO", "calForcePwPart");
 #ifdef __CUDA
     if(PARAM.inp.device == "gpu")
     {
         Forces<double, base_device::DEVICE_GPU> f_pw(nat);
         f_pw.cal_force_loc(ucell, fvl_dvl, rhopw, locpp.vloc, chr);
         f_pw.cal_force_ew(ucell, fewalds, rhopw, &sf);
-        f_pw.cal_force_cc(fcc, rhopw, chr, locpp.numeric, ucell);
+        f_pw.cal_force_cc(fcc, rhopw, chr, locpp.numeric, ucell,
+            xc_input.nspin, xc_spin.domag, xc_spin.domag_z, xc_input.gga_grad);
         f_pw.cal_force_scc(fscc, rhopw, vnew, vnew_exist, locpp.numeric, ucell);
     }
     else
@@ -937,7 +940,8 @@ void Force_Stress_LCAO<T>::calForcePwPart(UnitCell& ucell,
         Forces<double, base_device::DEVICE_CPU> f_pw(nat);
         f_pw.cal_force_loc(ucell, fvl_dvl, rhopw, locpp.vloc, chr);
         f_pw.cal_force_ew(ucell, fewalds, rhopw, &sf);
-        f_pw.cal_force_cc(fcc, rhopw, chr, locpp.numeric, ucell);
+        f_pw.cal_force_cc(fcc, rhopw, chr, locpp.numeric, ucell,
+            xc_input.nspin, xc_spin.domag, xc_spin.domag_z, xc_input.gga_grad);
         f_pw.cal_force_scc(fscc, rhopw, vnew, vnew_exist, locpp.numeric, ucell);
     }
 
@@ -1027,6 +1031,8 @@ void Force_Stress_LCAO<T>::calStressPwPart(UnitCell& ucell,
                                            const pseudopot_cell_vl& locpp,
                                            const Structure_Factor& sf)
 {
+    const auto& xc_input = PARAM.inp;
+    const auto& xc_spin = PARAM.globalv;
     ModuleBase::TITLE("Force_Stress_LCAO", "calStressPwPart");
 
     // local pseudopotential stress:
@@ -1039,15 +1045,17 @@ void Force_Stress_LCAO<T>::calStressPwPart(UnitCell& ucell,
     sc_pw.stress_ewa(ucell, sigmaewa, rhopw, 0); // remain problem
 
     // stress due to core correlation.
-    sc_pw.stress_cc(sigmacc, rhopw, ucell, &sf, 0, locpp.numeric, chr);
+    sc_pw.stress_cc(sigmacc, rhopw, ucell, &sf, 0, locpp.numeric, chr,
+        xc_input.nspin, xc_spin.domag, xc_spin.domag_z, xc_input.gga_grad, xc_spin.gamma_only_pw);
 
     // stress due to self-consistent charge.
     for (int i = 0; i < 3; i++)
     {
         sigmaxc(i, i) = -etxc / ucell.omega;
     }
     // Exchange-correlation for PBE
-    sc_pw.stress_gga(ucell, sigmaxc, rhopw, chr);
+    sc_pw.stress_gga(ucell, sigmaxc, rhopw, chr,
+        xc_input.nspin, xc_spin.domag, xc_spin.domag_z, xc_input.gga_grad);
 
     return;
 }

source/source_pw/module_ofdft/of_stress_pw.cpp

@@ -16,6 +16,8 @@ void OF_Stress_PW::cal_stress(ModuleBase::matrix& sigmatot,
                               Structure_Factor* p_sf,
                               K_Vectors* p_kv)
 {
+    const auto& xc_input = PARAM.inp;
+    const auto& xc_spin = PARAM.globalv;
     ModuleBase::TITLE("OF_Stress_PW", "cal_stress");
     ModuleBase::timer::start("OF_Stress_PW", "cal_stress");
 
@@ -74,13 +76,15 @@ void OF_Stress_PW::cal_stress(ModuleBase::matrix& sigmatot,
     {
         sigmaxc(i, i) = -(pelec->f_en.etxc - pelec->f_en.vtxc) / ucell.omega;
     }
-    stress_gga(ucell,sigmaxc, this->rhopw, pelec->charge);
+    stress_gga(ucell,sigmaxc, this->rhopw, pelec->charge,
+        xc_input.nspin, xc_spin.domag, xc_spin.domag_z, xc_input.gga_grad);
 
     // local contribution
     stress_loc(ucell,sigmaloc, this->rhopw, locpp.vloc, p_sf, true, pelec->charge);
 
     // nlcc
-    stress_cc(sigmaxcc, this->rhopw, ucell, p_sf, true, locpp.numeric, pelec->charge);
+    stress_cc(sigmaxcc, this->rhopw, ucell, p_sf, true, locpp.numeric, pelec->charge,
+        xc_input.nspin, xc_spin.domag, xc_spin.domag_z, xc_input.gga_grad, xc_spin.gamma_only_pw);
 
     // vdW term prepared before SCF for this ionic configuration.
     if (vdw_result != nullptr)

source/source_pw/module_pwdft/force_pw.cpp

@@ -41,6 +41,8 @@ void Forces<FPTYPE, Device>::cal_force(UnitCell& ucell,
                                        ModulePW::PW_Basis_K* wfc_basis,
                                        const psi::Psi<std::complex<FPTYPE>, Device>* psi_in)
 {
+    const auto& xc_input = PARAM.inp;
+    const auto& xc_spin = PARAM.globalv;
     ModuleBase::timer::start("Forces", "cal_force");
     ModuleBase::TITLE("Forces", "init");
     this->device = base_device::get_device_type(this->ctx);
@@ -83,7 +85,8 @@ void Forces<FPTYPE, Device>::cal_force(UnitCell& ucell,
     }
 
     // non-linear core correction
-    Forces::cal_force_cc(forcecc, rho_basis, chr, locpp->numeric, ucell);
+    Forces::cal_force_cc(forcecc, rho_basis, chr, locpp->numeric, ucell,
+        xc_input.nspin, xc_spin.domag, xc_spin.domag_z, xc_input.gga_grad);
 
     // force due to core charge
     this->cal_force_scc(forcescc, rho_basis, elec.vnew, elec.vnew_exist, locpp->numeric, ucell);

source/source_pw/module_pwdft/force_pw.h

@@ -77,7 +77,11 @@ class Forces
                       const ModulePW::PW_Basis* const rho_basis,
                       const Charge* const chr,
                       const bool* numeric,
-                      UnitCell& ucell_in);
+                      UnitCell& ucell_in,
+                      const int nspin,
+                      const bool domag,
+                      const bool domag_z,
+                      const int gga_grad);
     /**
      * @brief This routine computes the atomic force of non-local pseudopotential
      *    F^{NL}_i = \sum_{n,k}f_{nk}\sum_I \sum_{lm,l'm'}D_{l,l'}^{I} [

source/source_pw/module_pwdft/force_pw_cc.cpp

@@ -1,7 +1,6 @@
 #include "force_pw.h"
 #include "stress_func.h"
 #include "source_base/parallel_reduce.h"
-#include "source_io/module_parameter/parameter.h"
 // new
 #include "source_base/complexmatrix.h"
 #include "source_base/libm/libm.h"
@@ -32,9 +31,12 @@ void Forces<FPTYPE, Device>::cal_force_cc(ModuleBase::matrix& forcecc,
                                           const ModulePW::PW_Basis* const rho_basis,
                                           const Charge* const chr,
                                           const bool* numeric,
-                                           UnitCell& ucell_in)
+                                          UnitCell& ucell_in,
+                                          const int nspin,
+                                          const bool domag,
+                                          const bool domag_z,
+                                          const int gga_grad)
 {
-    const Parameter& parameters = PARAM;
     ModuleBase::TITLE("Forces", "cal_force_cc");
     // recalculate the exchange-correlation potential.
     ModuleBase::timer::start("Forces", "cal_force_cc");
@@ -54,7 +56,7 @@ void Forces<FPTYPE, Device>::cal_force_cc(ModuleBase::matrix& forcecc,
         return;
     }
 
-    ModuleBase::matrix v(parameters.inp.nspin, rho_basis->nrxx);
+    ModuleBase::matrix v(nspin, rho_basis->nrxx);
 
     const double hybrid_alpha = XC_Functional::get_hybrid_alpha();
 #ifdef __EXX
@@ -67,7 +69,7 @@ void Forces<FPTYPE, Device>::cal_force_cc(ModuleBase::matrix& forcecc,
 #ifdef __LIBXC
         const auto etxc_vtxc_v
             = XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), rho_basis->nrxx, ucell_in.omega, ucell_in.tpiba, chr,
-                                             parameters.inp.nspin, hybrid_alpha, hse_omega);
+                                             nspin, hybrid_alpha, hse_omega);
 
         // etxc = std::get<0>(etxc_vtxc_v);
         // vtxc = std::get<1>(etxc_vtxc_v);
@@ -78,12 +80,12 @@ void Forces<FPTYPE, Device>::cal_force_cc(ModuleBase::matrix& forcecc,
     }
     else
     {
-        unitcell::cal_ux(ucell_in, parameters.inp.nspin);
+        unitcell::cal_ux(ucell_in, nspin);
         const auto etxc_vtxc_v = XC_Functional::v_xc(rho_basis->nrxx, chr, &ucell_in,
-                                              parameters.inp.nspin,
-                                              parameters.globalv.domag,
-                                              parameters.globalv.domag_z,
-                                              parameters.inp.gga_grad,
+                                              nspin,
+                                              domag,
+                                              domag_z,
+                                              gga_grad,
                                               hybrid_alpha,
                                               hse_omega);
 
@@ -94,7 +96,7 @@ void Forces<FPTYPE, Device>::cal_force_cc(ModuleBase::matrix& forcecc,
 
     const ModuleBase::matrix vxc = v;
     std::complex<double>* psiv = new std::complex<double>[rho_basis->nmaxgr];
-    if (parameters.inp.nspin == 1 || parameters.inp.nspin == 4)
+    if (nspin == 1 || nspin == 4)
     {
 #ifdef _OPENMP
 #pragma omp parallel for schedule(static, 1024)
@@ -370,4 +372,4 @@ void Forces<FPTYPE, Device>::deriv_drhoc
 template class Forces<double, base_device::DEVICE_CPU>;
 #if ((defined __CUDA) || (defined __ROCM))
 template class Forces<double, base_device::DEVICE_GPU>;
-#endif
+#endif

source/source_pw/module_pwdft/stress_cc.cpp

@@ -1,7 +1,6 @@
 #include "stress_func.h"
 #include "source_base/parallel_reduce.h"
 #include "source_hamilt/module_xc/xc_functional.h"
-#include "source_io/module_parameter/parameter.h"
 #include "source_base/math_integral.h"
 #include "source_base/timer.h"
 #include "source_cell/cal_ux.h"
@@ -15,19 +14,23 @@
 template <typename FPTYPE, typename Device>
 void Stress_Func<FPTYPE, Device>::stress_cc(ModuleBase::matrix& sigma,
                                             ModulePW::PW_Basis* rho_basis,
-											UnitCell& ucell,
+                                            UnitCell& ucell,
                                             const Structure_Factor* p_sf,
                                             const bool is_pw,
-											const bool *numeric,
-                                            const Charge* const chr)
+                                            const bool *numeric,
+                                            const Charge* const chr,
+                                            const int nspin,
+                                            const bool domag,
+                                            const bool domag_z,
+                                            const int gga_grad,
+                                            const bool gamma_only_pw)
 {
-    const Parameter& parameters = PARAM;
     ModuleBase::TITLE("Stress","stress_cc");
 	ModuleBase::timer::start("Stress","stress_cc");
 
 	FPTYPE fact=1.0;
 
-	if(is_pw&&parameters.globalv.gamma_only_pw)
+	if(is_pw&&gamma_only_pw)
 	{
 		fact = 2.0; //is_pw:PW basis, gamma_only need to FPTYPE.
 	}
@@ -63,7 +66,7 @@ void Stress_Func<FPTYPE, Device>::stress_cc(ModuleBase::matrix& sigma,
 #ifdef __LIBXC
         const auto etxc_vtxc_v
             = XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), rho_basis->nrxx, ucell.omega, ucell.tpiba, chr,
-                                             parameters.inp.nspin, hybrid_alpha, hse_omega);
+                                             nspin, hybrid_alpha, hse_omega);
 
         // etxc = std::get<0>(etxc_vtxc_v);
         // vtxc = std::get<1>(etxc_vtxc_v);
@@ -74,12 +77,12 @@ void Stress_Func<FPTYPE, Device>::stress_cc(ModuleBase::matrix& sigma,
 	}
 	else
 	{
-		unitcell::cal_ux(ucell, parameters.inp.nspin);
+		unitcell::cal_ux(ucell, nspin);
         const auto etxc_vtxc_v = XC_Functional::v_xc(rho_basis->nrxx, chr, &ucell,
-                                              parameters.inp.nspin,
-                                              parameters.globalv.domag,
-                                              parameters.globalv.domag_z,
-                                              parameters.inp.gga_grad,
+                                              nspin,
+                                              domag,
+                                              domag_z,
+                                              gga_grad,
                                               hybrid_alpha,
                                               hse_omega);
         // etxc = std::get<0>(etxc_vtxc_v); // may delete?
@@ -89,7 +92,7 @@ void Stress_Func<FPTYPE, Device>::stress_cc(ModuleBase::matrix& sigma,
 
     std::complex<FPTYPE>* psic = new std::complex<FPTYPE>[rho_basis->nmaxgr];
 
-    if(parameters.inp.nspin==1||parameters.inp.nspin==4)
+    if(nspin==1||nspin==4)
 	{
 #ifdef _OPENMP
 #pragma omp parallel for schedule(static, 1024)

source/source_pw/module_pwdft/stress_func.h

@@ -117,7 +117,12 @@ class Stress_Func
                    const Structure_Factor* p_sf,
                    const bool is_pw,
                    const bool *numeric,
-                   const Charge* const chr); // nonlinear core correction stress in PW or LCAO basis
+                   const Charge* const chr,
+                   const int nspin,
+                   const bool domag,
+                   const bool domag_z,
+                   const int gga_grad,
+                   const bool gamma_only_pw); // nonlinear core correction stress in PW or LCAO basis
 
     void deriv_drhoc(const bool& numeric,
                      const double& omega,
@@ -134,7 +139,11 @@ class Stress_Func
     void stress_gga(const UnitCell& ucell,
                     ModuleBase::matrix& sigma,
                     ModulePW::PW_Basis* rho_basis,
-                    const Charge* const chr); // gga part in both PW and LCAO basis
+                    const Charge* const chr,
+                    const int nspin,
+                    const bool domag,
+                    const bool domag_z,
+                    const int gga_grad); // gga part in both PW and LCAO basis
     void stress_mgga(const UnitCell& ucell,
                      ModuleBase::matrix& sigma,
                      const ModuleBase::matrix& wg,

source/source_pw/module_pwdft/stress_gga.cpp

@@ -1,16 +1,18 @@
 #include "stress_func.h"
 #include "source_base/parallel_reduce.h"
 #include "source_hamilt/module_xc/xc_functional.h"
-#include "source_io/module_parameter/parameter.h"
 
 //calculate the GGA stress correction in PW and LCAO
 template <typename FPTYPE, typename Device>
 void Stress_Func<FPTYPE, Device>::stress_gga(const UnitCell& ucell,
                                              ModuleBase::matrix& sigma,
                                              ModulePW::PW_Basis* rho_basis,
-                                             const Charge* const chr)
+                                             const Charge* const chr,
+                                             const int nspin,
+                                             const bool domag,
+                                             const bool domag_z,
+                                             const int gga_grad)
 {
-    const Parameter& parameters = PARAM;
     ModuleBase::TITLE("Stress","stress_gga");
     ModuleBase::timer::start("Stress","stress_gga");
 
@@ -32,7 +34,7 @@ void Stress_Func<FPTYPE, Device>::stress_gga(const UnitCell& ucell,
     XC_Functional::gradcorr(
         dum1, dum2, dum3, chr, rho_basis, &ucell,
         stress_gga, is_stress,
-        parameters.inp.nspin, parameters.globalv.domag, parameters.globalv.domag_z, parameters.inp.gga_grad,
+        nspin, domag, domag_z, gga_grad,
         hybrid_alpha, hse_omega);
 
     for(int l = 0;l< 3;l++)

source/source_pw/module_pwdft/stress_pw.cpp

@@ -23,6 +23,8 @@ void Stress_PW<FPTYPE, Device>::cal_stress(ModuleBase::matrix& sigmatot,
                                            const General_Exx_Info& exx_info,
                                            const psi::Psi <std::complex<FPTYPE>, Device>* d_psi_in)
 {
+    const auto& xc_input = PARAM.inp;
+    const auto& xc_spin = PARAM.globalv;
     ModuleBase::TITLE("Stress_PW", "cal_stress");
     ModuleBase::timer::start("Stress_PW", "cal_stress");
 
@@ -91,7 +93,8 @@ void Stress_PW<FPTYPE, Device>::cal_stress(ModuleBase::matrix& sigmatot,
     {
         sigmaxc(i, i) = -(pelec->f_en.etxc - pelec->f_en.vtxc) / ucell.omega;
     }
-    this->stress_gga(ucell, sigmaxc, rho_basis, pelec->charge);
+    this->stress_gga(ucell, sigmaxc, rho_basis, pelec->charge,
+        xc_input.nspin, xc_spin.domag, xc_spin.domag_z, xc_input.gga_grad);
     if (XC_Functional::get_ked_flag())
     {
         this->stress_mgga(ucell,
@@ -108,7 +111,8 @@ void Stress_PW<FPTYPE, Device>::cal_stress(ModuleBase::matrix& sigmatot,
     this->stress_loc(ucell, sigmaloc, rho_basis, locpp.vloc, p_sf, 1, pelec->charge);
 
     // nlcc
-    this->stress_cc(sigmaxcc, rho_basis, ucell, p_sf, 1, locpp.numeric, pelec->charge);
+    this->stress_cc(sigmaxcc, rho_basis, ucell, p_sf, 1, locpp.numeric, pelec->charge,
+        xc_input.nspin, xc_spin.domag, xc_spin.domag_z, xc_input.gga_grad, xc_spin.gamma_only_pw);
 
     // nonlocal
     this->stress_nl(sigmanl, this->pelec->wg, this->pelec->ekb, p_sf, p_kv, p_symm, wfc_basis, d_psi_in, nlpp, ucell);