diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 428f38fd01d..5fe3c57bdcc 100644 --- a/docs/advanced/input_files/input-main.md +++ b/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 diff --git a/docs/parameters.yaml b/docs/parameters.yaml index a0f76236e86..139e37c192d 100644 --- a/docs/parameters.yaml +++ b/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 diff --git a/source/Makefile.Objects b/source/Makefile.Objects index c47a23f4686..53049cffae0 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -593,6 +593,8 @@ OBJS_SYMMETRY=symm_other.o\ symm_rotation_k.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\ diff --git a/source/source_esolver/esolver_dfpt_pw.cpp b/source/source_esolver/esolver_dfpt_pw.cpp index e125b1d431c..4be0eb2db8b 100644 --- a/source/source_esolver/esolver_dfpt_pw.cpp +++ b/source/source_esolver/esolver_dfpt_pw.cpp @@ -37,12 +37,16 @@ namespace { class XC_First_Order_FDM : public ModuleDFPT::XC_First_Order { public: - XC_First_Order_FDM(ModulePW::PW_Basis* rho_basis, + XC_First_Order_FDM(const int nspin, + const bool domag, + const bool domag_z, + const int gga_grad, + ModulePW::PW_Basis* rho_basis, const Charge* chg0, const UnitCell* ucell) : ucell_(ucell) { - fdm_ = new elecstate::PotXC_FDM(rho_basis, chg0, ucell); + fdm_ = new elecstate::PotXC_FDM(nspin, domag, domag_z, gga_grad, rho_basis, chg0, ucell); chg1_ = new Charge(); chg1_->set_rhopw(rho_basis); chg1_->allocate(chg0->nspin, false); @@ -291,7 +295,8 @@ void ESolver_DFPT_PW::init_dfpt(UnitCell& ucell) } // first-order XC kernel adapter around the converged ground-state density - xc_adapter_ = new XC_First_Order_FDM(this->pw_rho, this->pelec->charge, &ucell); + xc_adapter_ = new XC_First_Order_FDM(this->inp_->nspin, PARAM.globalv.domag, PARAM.globalv.domag_z, this->inp_->gga_grad, + this->pw_rho, this->pelec->charge, &ucell); if (getenv("DFPT_VKB") != nullptr) { diff --git a/source/source_esolver/esolver_double_xc.cpp b/source/source_esolver/esolver_double_xc.cpp index 85fab14ab1f..2b66c8665db 100644 --- a/source/source_esolver/esolver_double_xc.cpp +++ b/source/source_esolver/esolver_double_xc.cpp @@ -398,8 +398,10 @@ void ESolver_DoubleXC::cal_force(BaseCell& basecell, ModuleBase::matrix& this->deepks.dpks_out_type = "base"; // for deepks method - FSCalcConfig fs_cfg{this->inp_->nspin, this->inp_->nbands, this->inp_->t_in_h, - this->inp_->sc_mag_switch, this->inp_->device}; + const FSCalcConfig fs_cfg{this->inp_->nspin, this->inp_->nbands, this->inp_->t_in_h, + this->inp_->sc_mag_switch, this->inp_->device, + PARAM.globalv.domag, PARAM.globalv.domag_z, this->inp_->gga_grad, + PARAM.globalv.gamma_only_pw}; fsl.getForceStress(ucell, this->get_vdw_result(), diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index e6747b42aa1..962ec47c77a 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -22,8 +22,8 @@ #include "source_estate/module_charge/chgmixing.h" // use charge mixing, mohan add 20251006 #include "source_estate/module_dm/init_dm.h" // init dm from electronic wave functions #include "source_io/module_restart/restart.h" // GlobalC::restart for load_exx_flag -#include "source_io/module_ctrl/ctrl_runner_lcao.h" // use ctrl_runner_lcao() -#include "source_io/module_ctrl/ctrl_iter_lcao.h" // use ctrl_iter_lcao() +#include "source_io/module_ctrl/ctrl_runner_lcao.h" // use ctrl_runner_lcao() +#include "source_io/module_ctrl/ctrl_iter_lcao.h" // use ctrl_iter_lcao() #include "source_io/module_ctrl/ctrl_scf_lcao.h" // use ctrl_scf_lcao() #include "source_io/module_output/print_info.h" #include "source_lcao/rho_tau_lcao.h" // mohan add 20251024 @@ -267,8 +267,10 @@ void ESolver_KS_LCAO::cal_force(BaseCell& basecell, ModuleBase::matrix& deepks.dpks_out_type = "tot"; // for deepks method - FSCalcConfig fs_cfg{this->inp_->nspin, this->inp_->nbands, this->inp_->t_in_h, - this->inp_->sc_mag_switch, this->inp_->device}; + const FSCalcConfig fs_cfg{this->inp_->nspin, this->inp_->nbands, this->inp_->t_in_h, + this->inp_->sc_mag_switch, this->inp_->device, + PARAM.globalv.domag, PARAM.globalv.domag_z, this->inp_->gga_grad, + PARAM.globalv.gamma_only_pw}; fsl.getForceStress(ucell, this->get_vdw_result(), this->inp_->cal_force, this->inp_->cal_stress, this->inp_->test_force, this->inp_->test_stress, @@ -329,7 +331,7 @@ void ESolver_KS_LCAO::after_all_runners(BaseCell& basecell) } ModuleIO::ctrl_runner_lcao(ucell, - *this->inp_, this->kv, this->pelec, this->dmat, this->pv, this->Pgrid, + *this->inp_, this->kv, this->pelec, this->dmat, this->pv, this->Pgrid, this->gd, this->psi, this->chr, hamilt_lcao, this->two_center_bundle_, this->orb_, this->pw_rho, this->pw_rhod, @@ -354,8 +356,8 @@ void ESolver_KS_LCAO::iter_init(UnitCell& ucell, const int istep, const // call iter_init() of ESolver_KS ESolver_KS::iter_init(ucell, istep, iter); - module_charge::chgmixing_ks_lcao(iter, this->p_chgmix, *this->dftu_, - this->dmat.dm->get_DMR_pointer(1)->get_nnr(), *this->inp_); + module_charge::chgmixing_ks_lcao(iter, this->p_chgmix, *this->dftu_, + this->dmat.dm->get_DMR_pointer(1)->get_nnr(), *this->inp_); if (iter == 1) { @@ -386,7 +388,7 @@ void ESolver_KS_LCAO::iter_init(UnitCell& ucell, const int istep, const { // the following steps are only needed in the first outer exx loop exx_two_level_step - = exx_info_.info_ri.real_number ? + = exx_info_.info_ri.real_number ? this->exx_nao.exd->two_level_step : this->exx_nao.exc->two_level_step; } #endif @@ -483,7 +485,7 @@ void ESolver_KS_LCAO::hamilt2rho_single(UnitCell& ucell, int istep, int GlobalV::MY_RANK); // the eigensolvers only ever ask the Hamiltonian for H(k) and S(k) hamilt::HamiltHSMatrix hs(static_cast*>(this->p_hamilt)); - hsolver_lcao_obj.solve(hs, this->psi[0], this->pelec, *this->dmat.dm, + hsolver_lcao_obj.solve(hs, this->psi[0], this->pelec, *this->dmat.dm, this->chr, this->inp_->nspin, skip_charge); } else @@ -549,7 +551,7 @@ void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& // call iter_finish() of ESolver_KS, where band gap is printed, // eig and occ are printed, magnetization is calculated, - // charge mixing is performed, potential is updated, + // charge mixing is performed, potential is updated, // HF and kS energies are computed, meta-GGA, Jason and restart ESolver_KS::iter_finish(ucell, istep, iter, conv_esolver); const bool precision_switched = this->gint_precision_controller_.update_after_iteration(this->drho, this->scf_thr); @@ -573,8 +575,8 @@ void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& // control the output related to the finished iteration ModuleIO::ctrl_iter_lcao(ucell, *this->inp_, this->kv, this->pelec, *this->dmat.dm, - this->pv, this->gd, this->psi, this->chr, this->p_chgmix, - hamilt_lcao, this->orb_, this->deepks, + this->pv, this->gd, this->psi, this->chr, this->p_chgmix, + hamilt_lcao, this->orb_, this->deepks, this->exx_nao, this->exx_info_, iter, istep, conv_esolver, this->scf_ene_thr); } diff --git a/source/source_esolver/esolver_ks_pw.cpp b/source/source_esolver/esolver_ks_pw.cpp index cce89b61005..3cf3391bf7d 100644 --- a/source/source_esolver/esolver_ks_pw.cpp +++ b/source/source_esolver/esolver_ks_pw.cpp @@ -379,7 +379,9 @@ void ESolver_KS_PW::cal_force(BaseCell& basecell, ModuleBase::matrix& this->stp.update_psi_d(); // Calculate forces - ff.cal_force(ucell, + ff.cal_force(this->inp_->nspin, PARAM.globalv.domag, PARAM.globalv.domag_z, this->inp_->gga_grad, + this->inp_->dft_plus_u || this->inp_->sc_mag_switch, + ucell, force, this->get_vdw_result(), *this->pelec, @@ -406,7 +408,9 @@ void ESolver_KS_PW::cal_stress(BaseCell& basecell, ModuleBase::matrix // mohan add 2025-10-12 this->stp.update_psi_d(); - ss.cal_stress(stress, + ss.cal_stress(this->inp_->nspin, PARAM.globalv.domag, PARAM.globalv.domag_z, + this->inp_->gga_grad, PARAM.globalv.gamma_only_pw, + stress, ucell, this->get_vdw_result(), *this->dftu_, diff --git a/source/source_esolver/esolver_of.cpp b/source/source_esolver/esolver_of.cpp index cfdf1a56e25..fbd9eeae9f9 100644 --- a/source/source_esolver/esolver_of.cpp +++ b/source/source_esolver/esolver_of.cpp @@ -561,7 +561,9 @@ void ESolver_OF::cal_force(BaseCell& basecell, ModuleBase::matrix& force) // here nullptr is for DFT+U, which may cause bugs, mohan note 2025-11-07 // solvent can be used? mohan ask 2025-11-07 - ff.cal_force(ucell, force, this->get_vdw_result(), *pelec, this->pw_rho, &ucell.symm, &sf, + ff.cal_force(this->inp_->nspin, PARAM.globalv.domag, PARAM.globalv.domag_z, this->inp_->gga_grad, + this->inp_->dft_plus_u || this->inp_->sc_mag_switch, + ucell, force, this->get_vdw_result(), *pelec, this->pw_rho, &ucell.symm, &sf, this->solvent, nullptr, &this->locpp); } @@ -581,6 +583,8 @@ void ESolver_OF::cal_stress(BaseCell& basecell, ModuleBase::matrix& stress) this->pphi_, this->pw_rho, kinetic_stress_); // kinetic stress OF_Stress_PW ss(this->pelec, this->pw_rho); - ss.cal_stress(stress, kinetic_stress_, ucell, this->get_vdw_result(), &ucell.symm, this->locpp, &sf, &kv); + ss.cal_stress(this->inp_->nspin, PARAM.globalv.domag, PARAM.globalv.domag_z, + this->inp_->gga_grad, PARAM.globalv.gamma_only_pw, + stress, kinetic_stress_, ucell, this->get_vdw_result(), &ucell.symm, this->locpp, &sf, &kv); } } // namespace ModuleESolver diff --git a/source/source_esolver/esolver_sdft_pw.cpp b/source/source_esolver/esolver_sdft_pw.cpp index 218bda84637..59a6958cdd1 100644 --- a/source/source_esolver/esolver_sdft_pw.cpp +++ b/source/source_esolver/esolver_sdft_pw.cpp @@ -225,7 +225,8 @@ void ESolver_SDFT_PW::cal_force(BaseCell& basecell, ModuleBase::matri Sto_Forces ff(ucell.nat); - ff.cal_stoforce(force, + ff.cal_stoforce(this->inp_->nspin, PARAM.globalv.domag, PARAM.globalv.domag_z, this->inp_->gga_grad, + force, *this->pelec, this->pw_rho, &ucell.symm, @@ -246,7 +247,9 @@ void ESolver_SDFT_PW::cal_stress(BaseCell& basecell, ModuleBase::matr UnitCell& ucell = static_cast(basecell); Sto_Stress_PW ss; - ss.cal_stress(stress, + ss.cal_stress(this->inp_->nspin, PARAM.globalv.domag, PARAM.globalv.domag_z, + this->inp_->gga_grad, PARAM.globalv.gamma_only_pw, + stress, *this->pelec, this->pw_rho, &ucell.symm, diff --git a/source/source_estate/module_pot/pot_xc.cpp b/source/source_estate/module_pot/pot_xc.cpp index eb311cc9225..a4e2061699c 100644 --- a/source/source_estate/module_pot/pot_xc.cpp +++ b/source/source_estate/module_pot/pot_xc.cpp @@ -3,7 +3,6 @@ #include "source_base/timer.h" #include "source_base/constants.h" #include "source_hamilt/module_xc/xc_functional.h" -#include "source_io/module_parameter/parameter.h" #ifdef __LIBXC #include "source_hamilt/module_xc/libxc_abacus.h" @@ -17,6 +16,7 @@ void PotXC::cal_v_eff(const Charge*const chg, const UnitCell*const ucell, Module ModuleBase::TITLE("PotXC", "cal_veff"); ModuleBase::timer::start("PotXC", "cal_veff"); const int nrxx_current = chg->nrxx; + const int nspin = chg->nspin; //---------------------------------------------------------- // calculate the exchange-correlation potential @@ -33,7 +33,7 @@ void PotXC::cal_v_eff(const Charge*const chg, const UnitCell*const ucell, Module #endif const std::tuple 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 +52,10 @@ void PotXC::cal_v_eff(const Charge*const chg, const UnitCell*const ucell, Module #endif const std::tuple etxc_vtxc_v = XC_Functional::v_xc(nrxx_current, chg, ucell, - PARAM.inp.nspin, - PARAM.globalv.domag, - PARAM.globalv.domag_z, + nspin, + domag_, + domag_z_, + gga_grad_, hybrid_alpha, hse_omega); *(this->etxc_) = std::get<0>(etxc_vtxc_v); diff --git a/source/source_estate/module_pot/pot_xc.h b/source/source_estate/module_pot/pot_xc.h index 47243c633d7..50737e3ec95 100644 --- a/source/source_estate/module_pot/pot_xc.h +++ b/source/source_estate/module_pot/pot_xc.h @@ -11,11 +11,15 @@ class PotXC : public PotBase public: // constructor for exchange-correlation potential // meta-GGA should input matrix of kinetic potential, it is optional - PotXC(const ModulePW::PW_Basis* rho_basis_in, + PotXC(const bool domag, + const bool domag_z, + const int gga_grad, + const ModulePW::PW_Basis* rho_basis_in, double* etxc_in, double* vtxc_in, ModuleBase::matrix* vofk_in = nullptr) - : etxc_(etxc_in), vtxc_(vtxc_in), vofk(vofk_in) + : vofk(vofk_in), etxc_(etxc_in), vtxc_(vtxc_in), + domag_(domag), domag_z_(domag_z), gga_grad_(gga_grad) { this->rho_basis_ = rho_basis_in; this->dynamic_mode = true; @@ -27,8 +31,12 @@ class PotXC : public PotBase ModuleBase::matrix* vofk = nullptr; double* etxc_ = nullptr; double* vtxc_ = nullptr; + private: + const bool domag_; + const bool domag_z_; + const int gga_grad_; }; } // namespace elecstate -#endif \ No newline at end of file +#endif diff --git a/source/source_estate/module_pot/pot_xc_fdm.cpp b/source/source_estate/module_pot/pot_xc_fdm.cpp index 03349fe4e7e..4281553e1f9 100644 --- a/source/source_estate/module_pot/pot_xc_fdm.cpp +++ b/source/source_estate/module_pot/pot_xc_fdm.cpp @@ -5,16 +5,19 @@ #include "pot_xc_fdm.h" #include "source_hamilt/module_xc/xc_functional.h" -#include "source_io/module_parameter/parameter.h" namespace elecstate { PotXC_FDM::PotXC_FDM( + const int nspin, + const bool domag, + const bool domag_z, + const int gga_grad, const ModulePW::PW_Basis* rho_basis_in, const Charge*const chg_0_in, const UnitCell*const ucell) - : chg_0(chg_0_in) + : chg_0(chg_0_in), nspin_(nspin), domag_(domag), domag_z_(domag_z), gga_grad_(gga_grad) { this->rho_basis_ = rho_basis_in; this->dynamic_mode = true; @@ -28,9 +31,10 @@ PotXC_FDM::PotXC_FDM( #endif const std::tuple 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, + nspin_, + domag_, + domag_z_, + gga_grad_, hybrid_alpha, hse_omega); this->v_xc_0 = std::get<2>(etxc_vtxc_v_0); @@ -66,9 +70,10 @@ void PotXC_FDM::cal_v_eff( #endif const std::tuple etxc_vtxc_v_01 = XC_Functional::v_xc(chg_01.nrxx, &chg_01, ucell, - PARAM.inp.nspin, - PARAM.globalv.domag, - PARAM.globalv.domag_z, + nspin_, + domag_, + domag_z_, + gga_grad_, hybrid_alpha, hse_omega); const ModuleBase::matrix &v_xc_01 = std::get<2>(etxc_vtxc_v_01); diff --git a/source/source_estate/module_pot/pot_xc_fdm.h b/source/source_estate/module_pot/pot_xc_fdm.h index eb1cd880576..2ee16b591e9 100644 --- a/source/source_estate/module_pot/pot_xc_fdm.h +++ b/source/source_estate/module_pot/pot_xc_fdm.h @@ -16,6 +16,10 @@ class PotXC_FDM : public PotBase public: PotXC_FDM( + const int nspin, + const bool domag, + const bool domag_z, + const int gga_grad, const ModulePW::PW_Basis* rho_basis_in, const Charge*const chg_0_in, const UnitCell*const ucell); @@ -27,8 +31,13 @@ class PotXC_FDM : public PotBase const Charge*const chg_0 = nullptr; ModuleBase::matrix v_xc_0; +private: + const int nspin_; + const bool domag_; + const bool domag_z_; + const int gga_grad_; }; } // namespace elecstate -#endif \ No newline at end of file +#endif diff --git a/source/source_estate/module_pot/potential_types.cpp b/source/source_estate/module_pot/potential_types.cpp index 2aca1d5075d..4aa14ba6283 100644 --- a/source/source_estate/module_pot/potential_types.cpp +++ b/source/source_estate/module_pot/potential_types.cpp @@ -33,7 +33,7 @@ PotBase* Potential::get_pot_type(const std::string& pot_type) } else if (pot_type == "xc") { - return new PotXC(this->rho_basis_, this->etxc_, this->vtxc_, &(this->vofk_eff)); + return new PotXC(PARAM.globalv.domag, PARAM.globalv.domag_z, PARAM.inp.gga_grad, this->rho_basis_, this->etxc_, this->vtxc_, &(this->vofk_eff)); } else if (pot_type == "surchem") { diff --git a/source/source_hamilt/module_surchem/test/CMakeLists.txt b/source/source_hamilt/module_surchem/test/CMakeLists.txt index 50c14955cee..b843b7cb6e1 100644 --- a/source/source_hamilt/module_surchem/test/CMakeLists.txt +++ b/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 diff --git a/source/source_hamilt/module_xc/CMakeLists.txt b/source/source_hamilt/module_xc/CMakeLists.txt index 0b6309a5a24..495c5ae1997 100644 --- a/source/source_hamilt/module_xc/CMakeLists.txt +++ b/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 diff --git a/source/source_hamilt/module_xc/libxc_abacus.h b/source/source_hamilt/module_xc/libxc_abacus.h index 087a3577fc5..022cfc11122 100644 --- a/source/source_hamilt/module_xc/libxc_abacus.h +++ b/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 #include +#include #include #include @@ -19,28 +21,48 @@ class Charge; namespace XC_Functional_Libxc { + struct LibxcWeightedDerivatives + { + double energy_sum; + std::vector drho; + std::vector 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 spin_map; + std::vector rho; + std::vector>> spin_gradient; + std::array>, 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> 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 - * + * @return std::vector + * * @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 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* 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& func_id, + const std::size_t nrxx, + const double tpiba, + const Charge* const chr, + const std::map* scaling_factor, + const double hybrid_alpha, + const double hse_omega, + std::vector& stress_gga); + // for mGGA functional extern std::tuple v_xc_meta( const std::vector &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& drho, + const std::vector& dsigma, + const double tpiba, + const Charge* const chr); + // calculating grho extern std::vector>> cal_gdr( const int nspin, @@ -137,7 +192,21 @@ namespace XC_Functional_Libxc const std::vector &rho, std::vector 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 &sgn, + const std::vector &rho, + const std::vector &sigma, + const std::vector &exc, + const std::vector &vrho, + const std::vector &vsigma); + + // Convert collinear LibXC derivatives to the potential. extern std::pair 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 &amag, const ModuleBase::matrix &v, - const bool domag, - const bool domag_z); - + const bool has_mag); //------------------- // libxc_lda_wrap.cpp diff --git a/source/source_hamilt/module_xc/libxc_pot.cpp b/source/source_hamilt/module_xc/libxc_pot.cpp index 12351a45868..9eb5e96505c 100644 --- a/source/source_hamilt/module_xc/libxc_pot.cpp +++ b/source/source_hamilt/module_xc/libxc_pot.cpp @@ -15,25 +15,195 @@ #include #include -std::tuple XC_Functional_Libxc::v_xc_libxc( // Peize Lin update for nspin==4 at 2023.01.14 - const std::vector &func_id, - const int &nrxx, // number of real-space grid - const double &omega, // volume of cell - const double tpiba, - const Charge* const chr, - const int nspin_in, - const bool domag, - const bool domag_z, - const std::map* scaling_factor, - const double hybrid_alpha, - const double hse_omega) +void XC_Functional_Libxc::gradcorr_ncgga_sf_libxc(const std::vector& func_id, + const std::size_t nrxx, + const double tpiba, + const Charge* const chr, + const std::map* scaling_factor, + const double hybrid_alpha, + const double hse_omega, + std::vector& stress_gga) { - ModuleBase::TITLE("XC_Functional_Libxc","v_xc_libxc"); - ModuleBase::timer::start("XC_Functional_Libxc","v_xc_libxc"); + constexpr int nspin = 2; + stress_gga.assign(9, 0.0); - const int nspin = - (nspin_in == 1 || ( nspin_in ==4 && !domag && !domag_z)) - ? 1 : 2; + std::vector funcs = XC_Functional_Libxc::init_func(func_id, XC_POLARIZED, hybrid_alpha, hse_omega); + bool has_gga = false; + for (const xc_func_type& func: funcs) + { + has_gga = has_gga || func.info->family == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA; + } + if (!has_gga) + { + XC_Functional_Libxc::finish_func(funcs); + return; + } + + // This is the same forward graph used by v_xc_libxc: the local spin map, + // its projected FFT gradients, and the sigma invariants are constructed + // once and shared by all Libxc components. + const XC_Functional_Libxc::NclSfDiscreteData sf_data + = XC_Functional_Libxc::make_ncl_sf_discrete_data(nrxx, tpiba, chr, true); + const std::vector& rho = sf_data.rho; + const std::vector sigma = XC_Functional_Libxc::convert_sigma(sf_data.spin_gradient); + std::vector aggregate_dsigma(3 * nrxx, 0.0); + + for (xc_func_type& func: funcs) + { + if (func.info->family != XC_FAMILY_GGA && func.info->family != XC_FAMILY_HYB_GGA) + { + continue; + } + + constexpr double rho_threshold = 1.0e-6; + constexpr double grho_threshold = 1.0e-10; + xc_func_set_dens_threshold(&func, rho_threshold); + const std::vector sgn + = XC_Functional_Libxc::cal_sgn(rho_threshold, grho_threshold, func, nspin, nrxx, rho, sigma); + std::vector exc(nrxx); + std::vector vrho(nspin * nrxx); + std::vector vsigma(3 * nrxx); + constexpr int nr_batch_size = 1024; +#ifdef _OPENMP +#pragma omp parallel for schedule(static, nr_batch_size) +#endif + for (int ir_start = 0; ir_start < static_cast(nrxx); ir_start += nr_batch_size) + { + const int ir_end = std::min(ir_start + nr_batch_size, static_cast(nrxx)); + const int nrxx_thread = ir_end - ir_start; + xc_gga_exc_vxc(&func, + nrxx_thread, + rho.data() + ir_start * nspin, + sigma.data() + ir_start * 3, + exc.data() + ir_start, + vrho.data() + ir_start * nspin, + vsigma.data() + ir_start * 3); + } + + double factor = 1.0; + if (scaling_factor != nullptr) + { + const std::map::const_iterator entry = scaling_factor->find(func.info->number); + if (entry != scaling_factor->end()) + { + factor = entry->second; + } + } + const XC_Functional_Libxc::LibxcWeightedDerivatives weighted + = XC_Functional_Libxc::make_libxc_weighted_derivatives(func, + nspin, + nrxx, + sgn, + rho, + sigma, + exc, + vrho, + vsigma); + for (std::size_t index = 0; index < aggregate_dsigma.size(); ++index) + { + aggregate_dsigma[index] += factor * weighted.dsigma[index]; + } + } + +// For g_s=sum_A J_sA G_h x_A, a reciprocal deformation changes G_h but +// not the pointwise map J. Therefore the metric derivative is exactly +// sum_s h_s,l g_s,m, with h_s=dE/dg_s built from the sanitizer-reversed, +// component-scaled aggregate above. +#ifdef _OPENMP +#pragma omp parallel + { + std::vector local_stress(9, 0.0); +#pragma omp for schedule(static, 512) + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + const std::size_t sigma_index = 3 * ir; + const ModuleBase::Vector3& grad_up = sf_data.spin_gradient[0][ir]; + const ModuleBase::Vector3& grad_down = sf_data.spin_gradient[1][ir]; + const ModuleBase::Vector3 h_up + = ModuleBase::e2 + * (2.0 * aggregate_dsigma[sigma_index] * grad_up + aggregate_dsigma[sigma_index + 1] * grad_down); + const ModuleBase::Vector3 h_down + = ModuleBase::e2 + * (2.0 * aggregate_dsigma[sigma_index + 2] * grad_down + aggregate_dsigma[sigma_index + 1] * grad_up); + const double grad_up_component[3] = {grad_up.x, grad_up.y, grad_up.z}; + const double grad_down_component[3] = {grad_down.x, grad_down.y, grad_down.z}; + const double h_up_component[3] = {h_up.x, h_up.y, h_up.z}; + const double h_down_component[3] = {h_down.x, h_down.y, h_down.z}; + for (int l = 0; l < 3; ++l) + { + for (int m = 0; m <= l; ++m) + { + local_stress[l * 3 + m] + += h_up_component[l] * grad_up_component[m] + h_down_component[l] * grad_down_component[m]; + } + } + } +#pragma omp critical(libxc_ncgga_stress_reduce) + { + for (int l = 0; l < 3; ++l) + { + for (int m = 0; m <= l; ++m) + { + stress_gga[l * 3 + m] += local_stress[l * 3 + m]; + } + } + } + } +#else + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + const std::size_t sigma_index = 3 * ir; + const ModuleBase::Vector3& grad_up = sf_data.spin_gradient[0][ir]; + const ModuleBase::Vector3& grad_down = sf_data.spin_gradient[1][ir]; + const ModuleBase::Vector3 h_up + = ModuleBase::e2 + * (2.0 * aggregate_dsigma[sigma_index] * grad_up + aggregate_dsigma[sigma_index + 1] * grad_down); + const ModuleBase::Vector3 h_down + = ModuleBase::e2 + * (2.0 * aggregate_dsigma[sigma_index + 2] * grad_down + aggregate_dsigma[sigma_index + 1] * grad_up); + const double grad_up_component[3] = {grad_up.x, grad_up.y, grad_up.z}; + const double grad_down_component[3] = {grad_down.x, grad_down.y, grad_down.z}; + const double h_up_component[3] = {h_up.x, h_up.y, h_up.z}; + const double h_down_component[3] = {h_down.x, h_down.y, h_down.z}; + for (int l = 0; l < 3; ++l) + { + for (int m = 0; m <= l; ++m) + { + stress_gga[l * 3 + m] + += h_up_component[l] * grad_up_component[m] + h_down_component[l] * grad_down_component[m]; + } + } + } +#endif + + XC_Functional_Libxc::finish_func(funcs); +} + +std::tuple XC_Functional_Libxc::v_xc_libxc( // Peize Lin update for nspin==4 at + // 2023.01.14 + const std::vector& func_id, + const int& nrxx, // number of real-space grid + const double& omega, // volume of cell + const double tpiba, + const Charge* const chr, + const int nspin_in, + const bool domag, + const bool domag_z, + const int gga_grad, + const std::map* scaling_factor, + const double hybrid_alpha, + const double hse_omega) +{ + ModuleBase::TITLE("XC_Functional_Libxc", "v_xc_libxc"); + ModuleBase::timer::start("XC_Functional_Libxc", "v_xc_libxc"); + + const int nspin = (nspin_in == 1 || (nspin_in == 4 && !domag && !domag_z)) ? 1 : 2; + + // For nspin=4 with noncollinear magnetism, gga_grad=2 selects the + // regularized projected local-collinear graph; gga_grad=0/1 keeps the + // original collinear algorithm. + const bool has_mag = domag || domag_z; + const bool use_lca = (nspin_in == 4) && has_mag && gga_grad == 2; //---------------------------------------------------------- // xc_func_type is defined in Libxc package @@ -43,52 +213,81 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( / //---------------------------------------------------------- std::vector funcs = XC_Functional_Libxc::init_func( - /* func_id = */ func_id, - /* xc_polarized = */ (1==nspin) ? XC_UNPOLARIZED : XC_POLARIZED, + /* func_id = */ func_id, + /* xc_polarized = */ (1 == nspin) ? XC_UNPOLARIZED : XC_POLARIZED, /* hybrid_alpha = */ hybrid_alpha, /* hse_omega = */ hse_omega); - const bool is_gga = [&funcs]() - { - for( xc_func_type &func : funcs ) + const bool is_gga = [&funcs]() { + for (xc_func_type& func: funcs) { - switch( func.info->family ) + switch (func.info->family) { - case XC_FAMILY_GGA: - case XC_FAMILY_HYB_GGA: - return true; + case XC_FAMILY_GGA: + case XC_FAMILY_HYB_GGA: + return true; } } return false; }(); // 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 rho; std::vector amag; - if(1==nspin || 2==nspin_in) + XC_Functional_Libxc::NclSfDiscreteData sf_data; + if (1 == nspin || 2 == nspin_in) { rho = XC_Functional_Libxc::convert_rho(nspin, nrxx, chr); } + else if (use_lca) + { + // gga_grad=2 uses one complete local map for both the Libxc density + // input and the projected FFT-gradient graph. LDA-only functionals + // need the same local map but do not pay for gradients. + sf_data = XC_Functional_Libxc::make_ncl_sf_discrete_data(nrxx, tpiba, chr, is_gga); + rho = sf_data.rho; + } else { - std::tuple,std::vector> rho_amag = XC_Functional_Libxc::convert_rho_amag_nspin4(nspin, nrxx, chr); + std::tuple, std::vector> 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)); } std::vector>> gdr; std::vector sigma; - if(is_gga) + if (is_gga) { - gdr = XC_Functional_Libxc::cal_gdr(nspin, nrxx, rho, tpiba, chr); + if (use_lca) + { + gdr = sf_data.spin_gradient; + } + else + gdr = XC_Functional_Libxc::cal_gdr(nspin, nrxx, rho, tpiba, chr); + sigma = XC_Functional_Libxc::convert_sigma(gdr); } double etxc = 0.0; double vtxc = 0.0; - ModuleBase::matrix v(nspin,nrxx); + ModuleBase::matrix v(use_lca ? 4 : nspin, nrxx); + XC_Functional_Libxc::LibxcWeightedDerivatives sf_weighted; + if (use_lca) + { + sf_weighted.energy_sum = 0.0; + sf_weighted.drho.assign(nrxx * nspin, 0.0); + if (is_gga) + { + sf_weighted.dsigma.assign(nrxx * 3, 0.0); + } + } - for( xc_func_type &func : funcs ) + for (xc_func_type& func: funcs) { // jiyy add for threshold constexpr double rho_threshold = 1E-6; @@ -97,97 +296,146 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( / xc_func_set_dens_threshold(&func, rho_threshold); // sgn for threshold mask - const std::vector sgn = XC_Functional_Libxc::cal_sgn(rho_threshold, grho_threshold, func, nspin, nrxx, rho, sigma); + const std::vector sgn + = XC_Functional_Libxc::cal_sgn(rho_threshold, grho_threshold, func, nspin, nrxx, rho, sigma); - std::vector exc ( nrxx ); - std::vector vrho ( nrxx * nspin ); - std::vector vsigma( nrxx * ((1==nspin)?1:3) ); + std::vector exc(nrxx); + std::vector vrho(nrxx * nspin); + std::vector vsigma(nrxx * ((1 == nspin) ? 1 : 3)); - ModuleBase::timer::start("Libxc","xc_lda/gga_exc_vxc"); - switch( func.info->family ) + ModuleBase::timer::start("Libxc", "xc_lda/gga_exc_vxc"); + switch (func.info->family) { - case XC_FAMILY_LDA: + case XC_FAMILY_LDA: { + constexpr int nr_batch_size = 1024; +#ifdef _OPENMP +#pragma omp parallel for schedule(static, nr_batch_size) +#endif + for (int ir_start = 0; ir_start < nrxx; ir_start += nr_batch_size) { - constexpr int nr_batch_size = 1024; - #ifdef _OPENMP - #pragma omp parallel for schedule(static, nr_batch_size) - #endif - for( int ir_start = 0; ir_start < nrxx; ir_start += nr_batch_size ) - { - const int ir_end = std::min(ir_start + nr_batch_size, nrxx); - const int nrxx_thread = ir_end - ir_start; - xc_lda_exc_vxc( - &func, - nrxx_thread, - rho.data() + ir_start * nspin, - exc.data() + ir_start, - vrho.data() + ir_start * nspin ); - } - break; + const int ir_end = std::min(ir_start + nr_batch_size, nrxx); + const int nrxx_thread = ir_end - ir_start; + xc_lda_exc_vxc(&func, + nrxx_thread, + rho.data() + ir_start * nspin, + exc.data() + ir_start, + vrho.data() + ir_start * nspin); } - case XC_FAMILY_GGA: - case XC_FAMILY_HYB_GGA: - { - constexpr int nr_batch_size = 1024; - #ifdef _OPENMP - #pragma omp parallel for schedule(static, nr_batch_size) - #endif - for( int ir_start = 0; ir_start < nrxx; ir_start += nr_batch_size ) - { - const int ir_end = std::min(ir_start + nr_batch_size, nrxx); - const int nrxx_thread = ir_end - ir_start; - xc_gga_exc_vxc( - &func, - nrxx_thread, - rho.data() + ir_start * nspin, - sigma.data() + ir_start * ((1==nspin)?1:3), - exc.data() + ir_start, - vrho.data() + ir_start * nspin, - vsigma.data() + ir_start * ((1==nspin)?1:3) ); - } - break; - } - default: + break; + } + case XC_FAMILY_GGA: + case XC_FAMILY_HYB_GGA: { + constexpr int nr_batch_size = 1024; +#ifdef _OPENMP +#pragma omp parallel for schedule(static, nr_batch_size) +#endif + for (int ir_start = 0; ir_start < nrxx; ir_start += nr_batch_size) { - throw std::domain_error("func.info->family ="+std::to_string(func.info->family) - +" unfinished in "+std::string(__FILE__)+" line "+std::to_string(__LINE__)); - + const int ir_end = std::min(ir_start + nr_batch_size, nrxx); + const int nrxx_thread = ir_end - ir_start; + xc_gga_exc_vxc(&func, + nrxx_thread, + rho.data() + ir_start * nspin, + sigma.data() + ir_start * ((1 == nspin) ? 1 : 3), + exc.data() + ir_start, + vrho.data() + ir_start * nspin, + vsigma.data() + ir_start * ((1 == nspin) ? 1 : 3)); } + break; + } + default: { + throw std::domain_error("func.info->family =" + std::to_string(func.info->family) + " unfinished in " + + std::string(__FILE__) + " line " + std::to_string(__LINE__)); + } } - ModuleBase::timer::end("Libxc","xc_lda/gga_exc_vxc"); + ModuleBase::timer::end("Libxc", "xc_lda/gga_exc_vxc"); // added by jghan, 2024-10-10 double factor = 1.0; - if( scaling_factor ) + if (scaling_factor) { auto pair_factor = scaling_factor->find(func.info->number); - if( pair_factor != scaling_factor->end() ) - { factor = pair_factor->second; } + if (pair_factor != scaling_factor->end()) + { + factor = pair_factor->second; + } } - // time factor is added by jghan, 2024-10-10 + // Keep the established energy accumulation and reduction order. In + // gga_grad=2, reverse every sanitizer now, apply the component scaling, + // and aggregate before traversing the shared projected graph once. etxc += XC_Functional_Libxc::convert_etxc(nspin, nrxx, sgn, rho, exc) * factor; - const std::pair vtxc_v = XC_Functional_Libxc::convert_vtxc_v( - func, nspin, nrxx, - sgn, rho, gdr, - vrho, vsigma, - tpiba, chr); - vtxc += std::get<0>(vtxc_v) * factor; - v += std::get<1>(vtxc_v) * factor; + if (use_lca) + { + const XC_Functional_Libxc::LibxcWeightedDerivatives weighted + = XC_Functional_Libxc::make_libxc_weighted_derivatives(func, + nspin, + nrxx, + sgn, + rho, + sigma, + exc, + vrho, + vsigma); + for (std::size_t index = 0; index < sf_weighted.drho.size(); ++index) + { + sf_weighted.drho[index] += factor * weighted.drho[index]; + } + for (std::size_t index = 0; index < weighted.dsigma.size(); ++index) + { + sf_weighted.dsigma[index] += factor * weighted.dsigma[index]; + } + } + else + { + const std::pair vtxc_v + = XC_Functional_Libxc::convert_vtxc_v(func, nspin, nrxx, sgn, rho, gdr, vrho, vsigma, tpiba, chr); + 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) + if (use_lca) { - v = XC_Functional_Libxc::convert_v_nspin4(nrxx, chr, amag, v, domag, domag_z); + v = XC_Functional_Libxc::reverse_ncl_sf_discrete(nrxx, + sf_data, + sf_weighted.drho, + sf_weighted.dsigma, + tpiba, + chr); } - //------------------------------------------------- - // for MPI, reduce the exchange-correlation energy - //------------------------------------------------- - #ifdef __MPI + if (4 == nspin_in && !use_lca) + { + v = XC_Functional_Libxc::convert_v_nspin4(nrxx, chr, amag, v, has_mag); + } + + if (use_lca) + { + // Define vtxc from the potential that this routine actually returns. + // The nonlinear core density belongs to the XC energy graph, but the + // electronic variational density here is the four-channel valence + // density stored in chr->rho. + vtxc = 0.0; +#ifdef _OPENMP +#pragma omp parallel for collapse(2) reduction(+ : vtxc) schedule(static, 256) +#endif + for (int channel = 0; channel < 4; ++channel) + { + for (int ir = 0; ir < nrxx; ++ir) + { + vtxc += v(channel, ir) * chr->rho[channel][ir]; + } + } + } + +//------------------------------------------------- +// for MPI, reduce the exchange-correlation energy +//------------------------------------------------- +#ifdef __MPI Parallel_Reduce::reduce_pool(etxc); Parallel_Reduce::reduce_pool(vtxc); - #endif +#endif etxc *= omega / chr->rhopw->nxyz; vtxc *= omega / chr->rhopw->nxyz; @@ -198,7 +446,6 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( / return std::make_tuple( etxc, vtxc, std::move(v) ); } - //the interface to libxc xc_mgga_exc_vxc(xc_func,n,rho,grho,laplrho,tau,e,v1,v2,v3,v4) //xc_func : LIBXC data type, contains information on xc functional //n: size of array, nspin*nnr @@ -238,7 +485,7 @@ std::tuple XC_Functional_Li // https://www.tddft.org/programs/libxc/manual/libxc-5.1.x/ //---------------------------------------------------------- std::vector funcs = XC_Functional_Libxc::init_func( - /* func_id = */ func_id, + /* func_id = */ func_id, /* xc_polarized = */ (1==nspin) ? XC_UNPOLARIZED:XC_POLARIZED, /* hybrid_alpha = */ hybrid_alpha, /* hse_omega = */ hse_omega); diff --git a/source/source_hamilt/module_xc/libxc_tools.cpp b/source/source_hamilt/module_xc/libxc_tools.cpp index 5ed04eff8db..1579b4df3b8 100644 --- a/source/source_hamilt/module_xc/libxc_tools.cpp +++ b/source/source_hamilt/module_xc/libxc_tools.cpp @@ -4,6 +4,10 @@ #include "xc_functional.h" #include "source_estate/module_charge/charge.h" +#include +#include +#include + // converting rho (abacus=>libxc) std::vector XC_Functional_Libxc::convert_rho( const int nspin, @@ -31,10 +35,7 @@ XC_Functional_Libxc::convert_rho_amag_nspin4( const std::size_t nrxx, const Charge* const chr) { - // `nspin` is the number of LibXC spin channels, not the physical nspin of the - // system. This nspin4 branch is only entered for a magnetized nspin==4 system, - // which LibXC always evaluates in the polarized (2-channel) representation. - assert(nspin==2); + assert(nspin==2); // nspin here is the collapsed spin dimension for libxc (up/down) std::vector rho(nrxx*nspin); std::vector amag(nrxx); #ifdef _OPENMP @@ -53,6 +54,192 @@ XC_Functional_Libxc::convert_rho_amag_nspin4( return std::make_tuple(std::move(rho), std::move(amag)); } +XC_Functional_Libxc::NclSfDiscreteData +XC_Functional_Libxc::make_ncl_sf_discrete_data( + const std::size_t nrxx, + const double tpiba, + const Charge* const chr, + const bool need_gradient) +{ + constexpr int nspin = 2; + NclSfDiscreteData data; + data.spin_map.resize(nrxx); + data.rho.resize(nrxx * nspin); + + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + const std::array magnetization + = {{chr->rho[1][ir], chr->rho[2][ir], chr->rho[3][ir]}}; + const ModuleXC::NcggaRadialPoint radial + = ModuleXC::make_ncgga_radial_point( + magnetization, ModuleXC::ncgga_lca_radial_eta()); + data.spin_map[ir] = ModuleXC::make_ncgga_spin_map_point( + chr->rho[0][ir] + chr->rho_core[ir], radial); + data.rho[ir * nspin] = data.spin_map[ir].spin_density[0]; + data.rho[ir * nspin + 1] = data.spin_map[ir].spin_density[1]; + } + + if (!need_gradient) + { + return data; + } + + std::vector> grad_total(nrxx); + std::vector real_field(nrxx); + std::vector> reciprocal(chr->rhopw->npw); + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + real_field[ir] = chr->rho[0][ir] + chr->rho_core[ir]; + } + chr->rhopw->real2recip(real_field.data(), reciprocal.data()); + XC_Functional::grad_rho( + reciprocal.data(), grad_total.data(), chr->rhopw, tpiba); + + for (int mu = 0; mu < 3; ++mu) + { + data.grad_m[mu].resize(nrxx); + chr->rhopw->real2recip(chr->rho[mu + 1], reciprocal.data()); + XC_Functional::grad_rho( + reciprocal.data(), data.grad_m[mu].data(), chr->rhopw, tpiba); + } + + data.spin_gradient.resize(nspin); + for (int spin = 0; spin < nspin; ++spin) + { + data.spin_gradient[spin].resize(nrxx); + } + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 512) + #endif + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + for (int spin = 0; spin < nspin; ++spin) + { + ModuleBase::Vector3 gradient + = data.spin_map[ir].jacobian(spin, 0) * grad_total[ir]; + for (int mu = 0; mu < 3; ++mu) + { + gradient += data.spin_map[ir].jacobian(spin, mu + 1) + * data.grad_m[mu][ir]; + } + data.spin_gradient[spin][ir] = gradient; + } + } + return data; +} + +ModuleBase::matrix XC_Functional_Libxc::reverse_ncl_sf_discrete( + const std::size_t nrxx, + const NclSfDiscreteData& data, + const std::vector& drho, + const std::vector& dsigma, + const double tpiba, + const Charge* const chr) +{ + constexpr int nspin = 2; + constexpr int nchannel = 4; + assert(data.spin_map.size() == nrxx); + assert(data.rho.size() == nrxx * nspin); + assert(drho.size() == nrxx * nspin); + + ModuleBase::matrix potential(nchannel, nrxx); + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 512) + #endif + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + for (int channel = 0; channel < nchannel; ++channel) + { + potential(channel, ir) + = ModuleBase::e2 + * (data.spin_map[ir].jacobian(0, channel) + * drho[ir * nspin] + + data.spin_map[ir].jacobian(1, channel) + * drho[ir * nspin + 1]); + } + } + + if (dsigma.empty()) + { + return potential; + } + + assert(dsigma.size() == nrxx * 3); + assert(data.spin_gradient.size() == nspin); + for (int spin = 0; spin < nspin; ++spin) + { + assert(data.spin_gradient[spin].size() == nrxx); + } + for (int mu = 0; mu < 3; ++mu) + { + assert(data.grad_m[mu].size() == nrxx); + } + + std::vector> h_up(nrxx), h_down(nrxx); + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 512) + #endif + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + const std::size_t sigma_index = 3 * ir; + h_up[ir] + = ModuleBase::e2 + * (2.0 * dsigma[sigma_index] * data.spin_gradient[0][ir] + + dsigma[sigma_index + 1] * data.spin_gradient[1][ir]); + h_down[ir] + = ModuleBase::e2 + * (2.0 * dsigma[sigma_index + 2] * data.spin_gradient[1][ir] + + dsigma[sigma_index + 1] * data.spin_gradient[0][ir]); + } + + std::vector> flux(nrxx); + std::vector divergence(nrxx); + for (int channel = 0; channel < nchannel; ++channel) + { + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 512) + #endif + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + flux[ir] + = data.spin_map[ir].jacobian(0, channel) * h_up[ir] + + data.spin_map[ir].jacobian(1, channel) * h_down[ir]; + } + XC_Functional::grad_dot( + flux.data(), divergence.data(), chr->rhopw, tpiba); + + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 512) + #endif + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + potential(channel, ir) -= divergence[ir]; + if (channel == 0 || data.spin_map[ir].saturated) + { + continue; + } + + const ModuleBase::Vector3 spin_flux + = 0.5 * (h_up[ir] - h_down[ir]); + double local_response = 0.0; + for (int nu = 0; nu < 3; ++nu) + { + local_response + += data.spin_map[ir].radial.jacobian(nu, channel - 1) + * (spin_flux * data.grad_m[nu][ir]); + } + potential(channel, ir) += local_response; + } + } + return potential; +} + // calculating grho std::vector>> XC_Functional_Libxc::cal_gdr( @@ -213,6 +400,129 @@ double XC_Functional_Libxc::convert_etxc( return etxc; } +XC_Functional_Libxc::LibxcWeightedDerivatives +XC_Functional_Libxc::make_libxc_weighted_derivatives( + const xc_func_type &func, + const int nspin, + const std::size_t nrxx, + const std::vector &sgn, + const std::vector &rho, + const std::vector &sigma, + const std::vector &exc, + const std::vector &vrho, + const std::vector &vsigma) +{ + assert(nspin == 1 || nspin == 2); + assert(sgn.size() == nrxx * nspin); + assert(rho.size() == nrxx * nspin); + assert(exc.size() == nrxx); + assert(vrho.size() == nrxx * nspin); + assert(func.nspin == nspin); + + const bool is_gga + = func.info->family == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA; + const std::size_t nsigma = nspin == 1 ? 1 : 3; + if (is_gga) + { + assert(sigma.size() == nrxx * nsigma); + assert(vsigma.size() == nrxx * nsigma); + } + + LibxcWeightedDerivatives weighted; + weighted.energy_sum = 0.0; + weighted.drho.assign(nrxx * nspin, 0.0); + if (is_gga) + { + weighted.dsigma.assign(nrxx * nsigma, 0.0); + } + + const double density_floor = func.dens_threshold; + const double sigma_floor = func.sigma_threshold * func.sigma_threshold; + double energy_sum = 0.0; + #ifdef _OPENMP + #pragma omp parallel for reduction(+:energy_sum) schedule(static, 512) + #endif + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + double raw_density_sum = 0.0; + double sanitized_density_sum = 0.0; + double energy_weight = 0.0; + for (int is = 0; is < nspin; ++is) + { + const std::size_t index = ir * nspin + is; + raw_density_sum += rho[index]; + sanitized_density_sum += std::max(density_floor, rho[index]); + energy_weight += sgn[index] * rho[index]; + } + + // Libxc leaves all outputs zero below the total-density threshold. + // Inside that branch, the ABACUS weighted energy is locally constant. + if (raw_density_sum < density_floor) + { + continue; + } + + // ABACUS accumulates M*eps while Libxc differentiates Y*eps after + // y_s=max(T,rho_s). Hence d eps/d y_s=(vrho_s-eps)/Y. + energy_sum += energy_weight * exc[ir]; + const double libxc_weight = energy_weight / sanitized_density_sum; + for (int is = 0; is < nspin; ++is) + { + const std::size_t index = ir * nspin + is; + const double floor_jacobian = rho[index] > density_floor ? 1.0 : 0.0; + weighted.drho[index] + = sgn[index] * exc[ir] + + libxc_weight * floor_jacobian * (vrho[index] - exc[ir]); + } + + if (!is_gga) + { + continue; + } + + if (nspin == 1) + { + const double floor_jacobian = sigma[ir] > sigma_floor ? 1.0 : 0.0; + weighted.dsigma[ir] = libxc_weight * floor_jacobian * vsigma[ir]; + continue; + } + + const std::size_t sigma_index = 3 * ir; + const double sigma_uu = sigma[sigma_index]; + const double sigma_ud = sigma[sigma_index + 1]; + const double sigma_dd = sigma[sigma_index + 2]; + const double jacobian_uu = sigma_uu > sigma_floor ? 1.0 : 0.0; + const double jacobian_dd = sigma_dd > sigma_floor ? 1.0 : 0.0; + const double sanitized_uu = std::max(sigma_floor, sigma_uu); + const double sanitized_dd = std::max(sigma_floor, sigma_dd); + const double cross_limit = 0.5 * (sanitized_uu + sanitized_dd); + + double cross_to_diagonal = 0.0; + double cross_jacobian = 1.0; + if (sigma_ud < -cross_limit) + { + cross_to_diagonal = -0.5; + cross_jacobian = 0.0; + } + else if (sigma_ud > cross_limit) + { + cross_to_diagonal = 0.5; + cross_jacobian = 0.0; + } + + weighted.dsigma[sigma_index] + = libxc_weight * jacobian_uu + * (vsigma[sigma_index] + cross_to_diagonal * vsigma[sigma_index + 1]); + weighted.dsigma[sigma_index + 1] + = libxc_weight * cross_jacobian * vsigma[sigma_index + 1]; + weighted.dsigma[sigma_index + 2] + = libxc_weight * jacobian_dd + * (vsigma[sigma_index + 2] + cross_to_diagonal * vsigma[sigma_index + 1]); + } + weighted.energy_sum = energy_sum; + return weighted; +} + // converting vtxc and v from vrho and vsigma (libxc=>abacus) std::pair XC_Functional_Libxc::convert_vtxc_v( const xc_func_type &func, @@ -267,7 +577,6 @@ std::pair XC_Functional_Libxc::convert_vtxc_v( return std::make_pair(vtxc, std::move(v)); } - // dh for gga v std::vector> XC_Functional_Libxc::cal_dh( const int nspin, @@ -317,32 +626,30 @@ std::vector> XC_Functional_Libxc::cal_dh( return dh; } - // convert v for NSPIN=4 ModuleBase::matrix XC_Functional_Libxc::convert_v_nspin4( const std::size_t nrxx, const Charge* const chr, const std::vector &amag, const ModuleBase::matrix &v, - const bool domag, - const bool domag_z) + const bool has_mag) { //assert(nrxx>0); - constexpr int nspin = 4; + constexpr int nspin4 = 4; constexpr double vanishing_charge = 1.0e-10; - ModuleBase::matrix v_nspin4(nspin, nrxx); + ModuleBase::matrix v_nspin4(nspin4, nrxx); for( int ir=0; ir vanishing_charge ) { const double vs = 0.5 * (v(0,ir)-v(1,ir)); - for(int ipol=1; ipolrho[ipol][ir] / amag[ir]; } diff --git a/source/source_hamilt/module_xc/test/CMakeLists.txt b/source/source_hamilt/module_xc/test/CMakeLists.txt index 51fe559c966..a07e4eb4d7e 100644 --- a/source/source_hamilt/module_xc/test/CMakeLists.txt +++ b/source/source_hamilt/module_xc/test/CMakeLists.txt @@ -56,14 +56,14 @@ AddTest( ../libxc_lda_wrap.cpp ../libxc_gga_wrap.cpp ../libxc_mgga_wrap.cpp - ../xc_gga_corr.cpp ../xc_lda_corr.cpp - ../xc_gga_exch.cpp ../xc_lda_exch.cpp ../xc_hcth.cpp + ../xc_gga_corr.cpp ../xc_lda_corr.cpp + ../xc_gga_exch.cpp ../xc_lda_exch.cpp ../xc_hcth.cpp ) AddTest( TARGET MODULE_HAMILT_XCTest_VXC LIBS parameter MPI::MPI_CXX Libxc::xc psi device container - SOURCES test_xc5.cpp ../xc_grad.cpp ../xc_grad_prepare.cpp ../xc_grad_kernel.cpp ../xc_grad_assemble.cpp ../xc_grad_wfc.cpp ../xc_grad_utils.cpp ../xc_functional.cpp + SOURCES test_xc5.cpp test_libxc_tools.cpp ../xc_grad.cpp ../xc_grad_prepare.cpp ../xc_grad_kernel.cpp ../xc_grad_assemble.cpp ../xc_grad_wfc.cpp ../xc_grad_utils.cpp ../xc_functional.cpp ../xc_lda_wrap.cpp ../xc_gga_wrap.cpp ../libxc_setup.cpp ../libxc_lda_wrap.cpp @@ -72,6 +72,8 @@ AddTest( ../xc_gga_corr.cpp ../xc_lda_corr.cpp ../xc_gga_exch.cpp ../xc_lda_exch.cpp ../xc_hcth.cpp ../xc_pot.cpp + ../xc_functional_ncgga_sf.cpp + ../xc_ncgga_radial.cpp ../libxc_pot.cpp ../libxc_tools.cpp ../../../source_base/module_external/blas_connector_base.cpp ../../../source_base/module_external/blas_connector_vector.cpp ../../../source_base/module_external/blas_connector_matrix.cpp @@ -86,6 +88,35 @@ AddTest( ${FFT_SRC} ) +AddTest( + TARGET MODULE_HAMILT_XCTest_NCGGA_RADIAL + SOURCES test_xc_ncgga_radial.cpp ../xc_ncgga_radial.cpp +) + +AddTest( + TARGET MODULE_HAMILT_XCTest_NCGGA_DISCRETE_FD + LIBS xc_ planewave parameter MPI::MPI_CXX Libxc::xc psi device container base + SOURCES test_xc_functional_ncgga_sf.cpp +) +set_tests_properties( + MODULE_HAMILT_XCTest_NCGGA_DISCRETE_FD + PROPERTIES ENVIRONMENT "OMP_NUM_THREADS=1" +) +add_test( + NAME MODULE_HAMILT_XCTest_NCGGA_DISCRETE_FD_MPI2 + COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} 2 + ${MPIEXEC_PREFLAGS} + $ + ${MPIEXEC_POSTFLAGS} +) +set_tests_properties( + MODULE_HAMILT_XCTest_NCGGA_DISCRETE_FD_MPI2 + PROPERTIES + ENVIRONMENT "OMP_NUM_THREADS=1" + PROCESSORS 2 + WORKING_DIRECTORY $ +) + AddTest( TARGET MODULE_HAMILT_XCTest_SCANL_LAPL LIBS parameter MPI::MPI_CXX Libxc::xc @@ -126,3 +157,10 @@ AddTest( ../../../source_base/module_fft/fft_cpu.cpp ${FFT_SRC} ) + +# gradcorr dispatches the noncollinear stress path even in focused FFT tests. +foreach(xc_gradient_test MODULE_HAMILT_XCTest_GRADCORR MODULE_HAMILT_XCTest_LAPL) + target_sources(${xc_gradient_test} PRIVATE + ../xc_functional_ncgga_sf.cpp ../xc_ncgga_radial.cpp + ../libxc_pot.cpp ../libxc_tools.cpp ../../../source_base/timer.cpp) +endforeach() diff --git a/source/source_hamilt/module_xc/test/test_libxc_tools.cpp b/source/source_hamilt/module_xc/test/test_libxc_tools.cpp new file mode 100644 index 00000000000..0714c6a3d6d --- /dev/null +++ b/source/source_hamilt/module_xc/test/test_libxc_tools.cpp @@ -0,0 +1,228 @@ +#include +#include "../libxc_abacus.h" + +#include "gtest/gtest.h" +#include +#include +#include +#include + +#ifdef __LIBXC + +TEST(LibxcSanitizer, FullPolarizationUsesTheWeightedEnergyDerivative) +{ + const std::array functional_ids = {{XC_LDA_X, XC_LDA_C_PZ}}; + const double density = 0.45; + + for (std::size_t ifunc = 0; ifunc < functional_ids.size(); ++ifunc) + { + xc_func_type func; + ASSERT_EQ(xc_func_init(&func, functional_ids[ifunc], XC_POLARIZED), 0); + xc_func_set_dens_threshold(&func, 1.0e-6); + + const std::function evaluate + = [&func](const double rho_up, + const double rho_down, + XC_Functional_Libxc::LibxcWeightedDerivatives* const weighted) { + const std::vector rho = {rho_up, rho_down}; + const std::vector mask = {1.0, 1.0}; + std::vector exc(1, 0.0); + std::vector vrho(2, 0.0); + xc_lda_exc_vxc(&func, 1, rho.data(), exc.data(), vrho.data()); + if (weighted != nullptr) + { + *weighted = XC_Functional_Libxc::make_libxc_weighted_derivatives(func, + 2, + 1, + mask, + rho, + std::vector(), + exc, + vrho, + std::vector()); + } + return (rho_up + rho_down) * exc[0]; + }; + + XC_Functional_Libxc::LibxcWeightedDerivatives weighted; + const double energy = evaluate(density, 0.0, &weighted); + ASSERT_EQ(weighted.drho.size(), 2U); + EXPECT_DOUBLE_EQ(weighted.energy_sum, energy); + + const double steps[] = {1.0e-3, 5.0e-4, 2.5e-4, 1.25e-4}; + std::array errors = {{0.0, 0.0, 0.0, 0.0}}; + for (std::size_t ieps = 0; ieps < 4; ++ieps) + { + const double step = steps[ieps]; + const double finite_difference + = (evaluate(density + step, 0.0, nullptr) - evaluate(density - step, 0.0, nullptr)) / (2.0 * step); + errors[ieps] = std::abs(weighted.drho[0] - finite_difference); + EXPECT_LE(errors[ieps], 2.0e-7 * std::max(1.0, std::abs(weighted.drho[0]))) + << "functional_id=" << functional_ids[ifunc] << ", step=" << step; + } + EXPECT_LE(errors[1], 0.4 * errors[0] + 1.0e-12); + EXPECT_LE(errors[2], 0.4 * errors[1] + 1.0e-12); + + const double inactive_density = 0.5 * func.dens_threshold; + const double inactive_step = 0.2 * func.dens_threshold; + XC_Functional_Libxc::LibxcWeightedDerivatives inactive_weighted; + evaluate(density, inactive_density, &inactive_weighted); + const double inactive_finite_difference = (evaluate(density, inactive_density + inactive_step, nullptr) + - evaluate(density, inactive_density - inactive_step, nullptr)) + / (2.0 * inactive_step); + EXPECT_NEAR(inactive_finite_difference, + inactive_weighted.drho[1], + 2.0e-8 * std::max(1.0, std::abs(inactive_finite_difference))); + xc_func_end(&func); + } +} + +TEST(LibxcSanitizer, GgaSigmaReverseMatchesTheWeightedEnergy) +{ + xc_func_type func; + ASSERT_EQ(xc_func_init(&func, XC_GGA_C_PBE, XC_POLARIZED), 0); + xc_func_set_dens_threshold(&func, 1.0e-6); + xc_func_set_sigma_threshold(&func, 1.0e-2); + + const std::vector mask = {1.0, 0.0}; + const std::vector density = {0.40, 0.20}; + const double sigma_floor = func.sigma_threshold * func.sigma_threshold; + const std::array, 5> sigma_states = {{{{0.040, 0.010, 0.030}}, + {{0.5 * sigma_floor, 0.0, 0.030}}, + {{0.040, 0.200, 0.030}}, + {{0.040, -0.200, 0.030}}, + {{0.5 * sigma_floor, 0.200, 0.030}}}}; + + const std::function&, const std::vector&, XC_Functional_Libxc::LibxcWeightedDerivatives*)> evaluate + = [&func, &mask](const std::vector& rho, + const std::vector& sigma, + XC_Functional_Libxc::LibxcWeightedDerivatives* const weighted) { + std::vector exc(1, 0.0); + std::vector vrho(2, 0.0); + std::vector vsigma(3, 0.0); + xc_gga_exc_vxc(&func, 1, rho.data(), sigma.data(), exc.data(), vrho.data(), vsigma.data()); + if (weighted != nullptr) + { + *weighted + = XC_Functional_Libxc::make_libxc_weighted_derivatives(func, 2, 1, mask, rho, sigma, exc, vrho, vsigma); + } + return (mask[0] * rho[0] + mask[1] * rho[1]) * exc[0]; + }; + + for (std::size_t icase = 0; icase < sigma_states.size(); ++icase) + { + std::vector sigma(sigma_states[icase].begin(), sigma_states[icase].end()); + XC_Functional_Libxc::LibxcWeightedDerivatives weighted; + const double energy = evaluate(density, sigma, &weighted); + EXPECT_DOUBLE_EQ(weighted.energy_sum, energy); + ASSERT_EQ(weighted.drho.size(), 2U); + ASSERT_EQ(weighted.dsigma.size(), 3U); + + if (icase == 0) + { + for (int component = 0; component < 2; ++component) + { + std::vector perturbed_density = density; + const double step = 1.0e-6; + perturbed_density[component] += step; + const double energy_plus = evaluate(perturbed_density, sigma, nullptr); + perturbed_density[component] -= 2.0 * step; + const double energy_minus = evaluate(perturbed_density, sigma, nullptr); + const double finite_difference = (energy_plus - energy_minus) / (2.0 * step); + EXPECT_NEAR(finite_difference, + weighted.drho[component], + 2.0e-7 * std::max(1.0, std::abs(weighted.drho[component]))); + } + } + + if (icase == 1 || icase == 4) + { + EXPECT_DOUBLE_EQ(weighted.dsigma[0], 0.0); + } + if (icase >= 2) + { + EXPECT_DOUBLE_EQ(weighted.dsigma[1], 0.0); + } + + for (int component = 0; component < 3; ++component) + { + const double step = 1.0e-6; + sigma[component] += step; + const double energy_plus = evaluate(density, sigma, nullptr); + sigma[component] -= 2.0 * step; + const double energy_minus = evaluate(density, sigma, nullptr); + sigma[component] += step; + const double finite_difference = (energy_plus - energy_minus) / (2.0 * step); + EXPECT_NEAR(finite_difference, + weighted.dsigma[component], + 2.0e-7 * std::max(1.0, std::abs(weighted.dsigma[component]))) + << "case=" << icase << ", sigma component=" << component; + } + } + xc_func_end(&func); +} + +TEST(LibxcSanitizer, UnpolarizedSelfSigmaReverseMatchesTheWeightedEnergy) +{ + xc_func_type func; + ASSERT_EQ(xc_func_init(&func, XC_GGA_C_PBE, XC_UNPOLARIZED), 0); + xc_func_set_dens_threshold(&func, 1.0e-6); + xc_func_set_sigma_threshold(&func, 1.0e-2); + + const std::vector mask = {1.0}; + const std::vector density = {0.40}; + const double sigma_floor = func.sigma_threshold * func.sigma_threshold; + + const std::function evaluate + = [&func, &mask, &density](const double sigma_value, + XC_Functional_Libxc::LibxcWeightedDerivatives* const weighted) { + const std::vector sigma = {sigma_value}; + std::vector exc(1, 0.0); + std::vector vrho(1, 0.0); + std::vector vsigma(1, 0.0); + xc_gga_exc_vxc(&func, 1, density.data(), sigma.data(), exc.data(), vrho.data(), vsigma.data()); + if (weighted != nullptr) + { + *weighted = XC_Functional_Libxc::make_libxc_weighted_derivatives(func, + 1, + 1, + mask, + density, + sigma, + exc, + vrho, + vsigma); + } + return mask[0] * density[0] * exc[0]; + }; + + const double below_floor_sigma = 0.5 * sigma_floor; + const double below_floor_step = 0.2 * sigma_floor; + XC_Functional_Libxc::LibxcWeightedDerivatives below_floor_weighted; + const double below_floor_energy = evaluate(below_floor_sigma, &below_floor_weighted); + ASSERT_EQ(below_floor_weighted.dsigma.size(), 1U); + EXPECT_DOUBLE_EQ(below_floor_weighted.energy_sum, below_floor_energy); + EXPECT_DOUBLE_EQ(below_floor_weighted.dsigma[0], 0.0); + const double below_floor_finite_difference = (evaluate(below_floor_sigma + below_floor_step, nullptr) + - evaluate(below_floor_sigma - below_floor_step, nullptr)) + / (2.0 * below_floor_step); + EXPECT_NEAR(below_floor_finite_difference, below_floor_weighted.dsigma[0], 1.0e-12); + + const double above_floor_sigma = 0.040; + const double above_floor_step = 1.0e-6; + XC_Functional_Libxc::LibxcWeightedDerivatives above_floor_weighted; + const double above_floor_energy = evaluate(above_floor_sigma, &above_floor_weighted); + ASSERT_EQ(above_floor_weighted.dsigma.size(), 1U); + EXPECT_DOUBLE_EQ(above_floor_weighted.energy_sum, above_floor_energy); + EXPECT_GT(std::abs(above_floor_weighted.dsigma[0]), 1.0e-12); + const double above_floor_finite_difference = (evaluate(above_floor_sigma + above_floor_step, nullptr) + - evaluate(above_floor_sigma - above_floor_step, nullptr)) + / (2.0 * above_floor_step); + EXPECT_NEAR(above_floor_finite_difference, + above_floor_weighted.dsigma[0], + 2.0e-7 * std::max(1.0, std::abs(above_floor_weighted.dsigma[0]))); + + xc_func_end(&func); +} + +#endif diff --git a/source/source_hamilt/module_xc/test/test_xc3.cpp b/source/source_hamilt/module_xc/test/test_xc3.cpp index d472d972cc9..cd7d0e52e34 100644 --- a/source/source_hamilt/module_xc/test/test_xc3.cpp +++ b/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, hybrid_alpha, hse_omega); - XC_Functional::gradcorr(et1,vt1,v1,&chr,&rhopw,&ucell,stress1,true,nspin1,domag,domag_z, 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, hybrid_alpha, hse_omega); - XC_Functional::gradcorr(et2,vt2,v2,&chr,&rhopw,&ucell,stress2,true,nspin2,domag,domag_z, 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, 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); } }; @@ -236,16 +237,17 @@ class XCTest_GRADCORR_HF : public XCTest XC_Functional::set_xc_type("PBE"); XC_Functional::set_xc_type("HF"); + const int gga_grad = 0; const double hybrid_alpha = 1.0; const double hse_omega = 0.0; - XC_Functional::gradcorr(et1,vt1,v1,&chr,&rhopw,&ucell,stress1,false,nspin1,domag,domag_z, hybrid_alpha, hse_omega); - XC_Functional::gradcorr(et1,vt1,v1,&chr,&rhopw,&ucell,stress1,true, nspin1,domag,domag_z, 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, hybrid_alpha, hse_omega); - XC_Functional::gradcorr(et2,vt2,v2,&chr,&rhopw,&ucell,stress2,true, nspin2,domag,domag_z, 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, hybrid_alpha, hse_omega); - XC_Functional::gradcorr(et4,vt4,v4,&chr,&rhopw,&ucell,stress4,true, nspin4,domag_true,domag_z, 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); + XC_Functional::gradcorr(et4,vt4,v4,&chr,&rhopw,&ucell,stress4,true, nspin4,domag_true,domag_z,gga_grad,hybrid_alpha,hse_omega); } }; @@ -343,4 +345,4 @@ TEST_F(XCTest_GRADWFC, set_xc_type) EXPECT_NEAR(grad[i+j*5].imag(),0,1e-8); } } -} \ No newline at end of file +} diff --git a/source/source_hamilt/module_xc/test/test_xc5.cpp b/source/source_hamilt/module_xc/test/test_xc5.cpp index 01aa214feae..ce0f975f152 100644 --- a/source/source_hamilt/module_xc/test/test_xc5.cpp +++ b/source/source_hamilt/module_xc/test/test_xc5.cpp @@ -8,6 +8,13 @@ #include "source_cell/cal_ux.h" #include "../../../source_base/parallel_reduce.h" +#include +#include +#include +#include +#include +#include + /************************************************ * unit test of functionals ***********************************************/ @@ -22,7 +29,8 @@ class XCTest_VXC : public XCTest { protected: - double et1 = 0, vt1 = 0; + double et1 = 0; + double vt1 = 0; ModuleBase::matrix v1; double et2 = 0, vt2 = 0; @@ -82,13 +90,13 @@ class XCTest_VXC : public XCTest const double hybrid_alpha = XC_Functional::get_hybrid_alpha(); const double hse_omega = XC_Functional::get_hse_omega(); std::tuple etxc_vtxc_v - = XC_Functional::v_xc(rhopw.nrxx,&chr,&ucell,nspin1,domag,domag_z, hybrid_alpha, hse_omega); + = XC_Functional::v_xc(rhopw.nrxx,&chr,&ucell,nspin1,domag,domag_z,0, hybrid_alpha, hse_omega); et1 = std::get<0>(etxc_vtxc_v); vt1 = std::get<1>(etxc_vtxc_v); v1 = std::get<2>(etxc_vtxc_v); etxc_vtxc_v - = XC_Functional::v_xc(rhopw.nrxx,&chr,&ucell,nspin2,domag,domag_z, hybrid_alpha, hse_omega); + = XC_Functional::v_xc(rhopw.nrxx,&chr,&ucell,nspin2,domag,domag_z,0, hybrid_alpha, hse_omega); et2 = std::get<0>(etxc_vtxc_v); vt2 = std::get<1>(etxc_vtxc_v); v2 = std::get<2>(etxc_vtxc_v); @@ -126,7 +134,8 @@ class XCTest_VXC_Libxc : public XCTest { protected: - double et1 = 0, vt1 = 0; + double et1 = 0; + double vt1 = 0; ModuleBase::matrix v1; double et2 = 0, vt2 = 0; @@ -186,13 +195,13 @@ class XCTest_VXC_Libxc : public XCTest const double hybrid_alpha = XC_Functional::get_hybrid_alpha(); const double hse_omega = XC_Functional::get_hse_omega(); std::tuple etxc_vtxc_v - = XC_Functional::v_xc(rhopw.nrxx,&chr,&ucell,nspin1,domag,domag_z, hybrid_alpha, hse_omega); + = XC_Functional::v_xc(rhopw.nrxx,&chr,&ucell,nspin1,domag,domag_z,0, hybrid_alpha, hse_omega); et1 = std::get<0>(etxc_vtxc_v); vt1 = std::get<1>(etxc_vtxc_v); v1 = std::get<2>(etxc_vtxc_v); etxc_vtxc_v - = XC_Functional::v_xc(rhopw.nrxx,&chr,&ucell,nspin2,domag,domag_z, hybrid_alpha, hse_omega); + = XC_Functional::v_xc(rhopw.nrxx,&chr,&ucell,nspin2,domag,domag_z,0, hybrid_alpha, hse_omega); et2 = std::get<0>(etxc_vtxc_v); vt2 = std::get<1>(etxc_vtxc_v); v2 = std::get<2>(etxc_vtxc_v); @@ -230,7 +239,8 @@ class XCTest_VXC_meta : public XCTest { protected: - double et1 = 0, vt1 = 0; + double et1 = 0; + double vt1 = 0; ModuleBase::matrix v1,vtau1; double et2 = 0, vt2 = 0; @@ -353,6 +363,370 @@ TEST_F(XCTest_VXC_meta, set_xc_type) EXPECT_NEAR(vtau2(1,4),0.0311787189,1.0e-8); } +/************************************************ + * unit tests for the gga_grad keyword (nspin=4 + * noncollinear GGA gradient methods) + * + * Method 2 differentiates the complete discrete local-spin-map and FFT + * gradient graph. Its reverse is tested on a real PW grid in + * test_xc_functional_ncgga_sf.cpp. + ************************************************/ + +namespace +{ +constexpr int gga_grad_nrxx = 5; + +// build a mock 4-component charge on the mocked 5-point grid. +// pattern 0: m = (0,0,mz), mz>0, so m_hat = (0,0,1) everywhere +// pattern 1: m direction varies from point to point +// pattern 2: m . ux changes sign across the grid (ux = (0,1,2) in the mock) +struct Ns4Charge +{ + ModulePW::PW_Basis rhopw; + UnitCell ucell; + Charge chr; + + Ns4Charge(const int pattern) + { + rhopw.nrxx = gga_grad_nrxx; + rhopw.npw = gga_grad_nrxx; + rhopw.nmaxgr = gga_grad_nrxx; + rhopw.gcar = new ModuleBase::Vector3[gga_grad_nrxx]; + rhopw.nxyz = 1; + + ucell.tpiba = 1; + ucell.omega = 1; + ucell.magnet.lsign_ = true; + unitcell::cal_ux(ucell, 4); + + chr.rhopw = &(rhopw); + chr.rho = new double*[4]; + for (int is = 0; is < 4; ++is) + { + chr.rho[is] = new double[gga_grad_nrxx]; + } + chr.rhog = new std::complex*[2]; + chr.rhog[0] = new std::complex[gga_grad_nrxx]; + chr.rhog[1] = new std::complex[gga_grad_nrxx]; + chr.rho_core = new double[gga_grad_nrxx]; + chr.rhog_core = new std::complex[gga_grad_nrxx]; + + for (int i = 0; i < gga_grad_nrxx; ++i) + { + chr.rho[0][i] = 2.0 + i; + if (pattern == 1) + { + chr.rho[1][i] = 0.10 * (i + 1); + chr.rho[2][i] = 0.05 * (gga_grad_nrxx - i); + chr.rho[3][i] = 0.20 * (i + 1); + } + else if (pattern == 2) + { + chr.rho[1][i] = 0.0; + chr.rho[2][i] = (i % 2 == 0) ? 0.3 : -0.3; + chr.rho[3][i] = 0.05; + } + else + { + chr.rho[1][i] = 0.0; + chr.rho[2][i] = 0.0; + chr.rho[3][i] = 0.2 * (i + 1); + } + chr.rhog[0][i] = chr.rho[0][i]; + chr.rhog[1][i] = chr.rho[1][i]; + chr.rho_core[i] = 0; + chr.rhog_core[i] = 0; + rhopw.gcar[i] = 1; + } + } +}; + + + +struct Ns2LocalCharge +{ + ModulePW::PW_Basis rhopw; + Charge chr; + + Ns2LocalCharge() + { + rhopw.nrxx = 1; + rhopw.npw = 1; + rhopw.nmaxgr = 1; + rhopw.nxyz = 1; + rhopw.gcar = new ModuleBase::Vector3[1]; + rhopw.gcar[0] = 0.0; + + chr.rhopw = &rhopw; + chr.rho = new double*[2]; + chr.rhog = new std::complex*[2]; + for (int is = 0; is < 2; ++is) + { + chr.rho[is] = new double[1]; + chr.rhog[is] = new std::complex[1]; + chr.rhog[is][0] = 0.0; + } + chr.rho_core = new double[1]; + chr.rhog_core = new std::complex[1]; + chr.rho_core[0] = 0.0; + chr.rhog_core[0] = 0.0; + } +}; + +// run XC_Functional::v_xc for nspin=4 with noncollinear magnetism +std::tuple run_vxc_nspin4( + const std::string& functional, + const int pattern, + const int gga_grad) +{ + Ns4Charge mock(pattern); + XC_Functional::set_xc_type(functional); + return XC_Functional::v_xc(gga_grad_nrxx, + &mock.chr, + &mock.ucell, + 4, + true, + false, + gga_grad, + XC_Functional::get_hybrid_alpha(), + XC_Functional::get_hse_omega()); +} + +// compare two (etxc, vtxc, v) results +void expect_vxc_equal(const std::tuple& a, + const std::tuple& b, + const double tol) +{ + EXPECT_NEAR(std::get<0>(a), std::get<0>(b), tol); + EXPECT_NEAR(std::get<1>(a), std::get<1>(b), tol); + const ModuleBase::matrix& va = std::get<2>(a); + const ModuleBase::matrix& vb = std::get<2>(b); + ASSERT_EQ(va.nr, vb.nr); + ASSERT_EQ(va.nc, vb.nc); + for (int ir = 0; ir < va.nr; ++ir) + { + for (int ic = 0; ic < va.nc; ++ic) + { + EXPECT_NEAR(va(ir, ic), vb(ir, ic), tol); + } + } +} +} // namespace + +// m_hat = m/|m|, zero where |m| ~ 0 + + +// v_tot = 0.5*(v_up+v_dn), v_mu = 0.5*(v_up-v_dn)*m_hat_mu + + +// original conversion: has_mag=false leaves magnetic channels zero +TEST(GgaGradTools, ConvertVNspin4HasMag) +{ + Ns4Charge mock(0); + std::vector amag(gga_grad_nrxx); + for (int ir = 0; ir < gga_grad_nrxx; ++ir) + { + amag[ir] = mock.chr.rho[3][ir]; + } + ModuleBase::matrix v(2, gga_grad_nrxx); + for (int ir = 0; ir < gga_grad_nrxx; ++ir) + { + v(0, ir) = 1.0 + ir; + v(1, ir) = 0.5 * ir; + } + + const ModuleBase::matrix v_nomag + = XC_Functional_Libxc::convert_v_nspin4(gga_grad_nrxx, &mock.chr, amag, v, false); + for (int ir = 0; ir < gga_grad_nrxx; ++ir) + { + EXPECT_NEAR(v_nomag(0, ir), 0.5 * (v(0, ir) + v(1, ir)), 1e-14); + EXPECT_NEAR(v_nomag(1, ir), 0.0, 1e-14); + EXPECT_NEAR(v_nomag(2, ir), 0.0, 1e-14); + EXPECT_NEAR(v_nomag(3, ir), 0.0, 1e-14); + } + + const ModuleBase::matrix v_mag + = XC_Functional_Libxc::convert_v_nspin4(gga_grad_nrxx, &mock.chr, amag, v, true); + for (int ir = 0; ir < gga_grad_nrxx; ++ir) + { + const double vs = 0.5 * (v(0, ir) - v(1, ir)); + EXPECT_NEAR(v_mag(3, ir), vs * mock.chr.rho[3][ir] / amag[ir], 1e-14); + } +} + + + +// gga_grad=0 keeps the original built-in algorithm and must not crash +TEST(GgaGradVxc, BuiltinOriginalAlgorithmRuns) +{ + const std::tuple r0 = run_vxc_nspin4("PBE", 1, 0); + EXPECT_EQ(std::get<2>(r0).nr, 4); + EXPECT_TRUE(std::isfinite(std::get<0>(r0))); + EXPECT_TRUE(std::isfinite(std::get<1>(r0))); +} + +// noncolin_rho with lsign=true defines up/down w.r.t. the global axis ux +// through sign(m . ux); with lsign=false, up is always the local |m| +TEST(GgaGradTools, NoncolinRhoGlobalAxis) +{ + Ns4Charge mock(2); // pattern 2: m . ux changes sign across the grid + const double* ux = mock.ucell.magnet.ux_; // (0,1,2) in the mock + + std::vector rup(gga_grad_nrxx), rdn(gga_grad_nrxx), neg(gga_grad_nrxx); + XC_Functional::noncolin_rho( + rup.data(), rdn.data(), neg.data(), mock.chr.rho, gga_grad_nrxx, ux, true); + for (int ir = 0; ir < gga_grad_nrxx; ++ir) + { + const double mx = mock.chr.rho[1][ir]; + const double my = mock.chr.rho[2][ir]; + const double mz = mock.chr.rho[3][ir]; + const double amag = std::sqrt(mx * mx + my * my + mz * mz); + const double sign = (mx * ux[0] + my * ux[1] + mz * ux[2] > 0) ? 1.0 : -1.0; + EXPECT_NEAR(rup[ir], 0.5 * (mock.chr.rho[0][ir] + sign * amag), 1e-14); + EXPECT_NEAR(rdn[ir], 0.5 * (mock.chr.rho[0][ir] - sign * amag), 1e-14); + } + // the sign really flips on this grid, i.e. the global axis matters here + EXPECT_NEAR(neg[0], 1.0, 1e-14); + EXPECT_NEAR(neg[1], -1.0, 1e-14); + + XC_Functional::noncolin_rho( + rup.data(), rdn.data(), neg.data(), mock.chr.rho, gga_grad_nrxx, ux, false); + for (int ir = 0; ir < gga_grad_nrxx; ++ir) + { + const double mx = mock.chr.rho[1][ir]; + const double my = mock.chr.rho[2][ir]; + const double mz = mock.chr.rho[3][ir]; + const double amag = std::sqrt(mx * mx + my * my + mz * mz); + EXPECT_NEAR(rup[ir], 0.5 * (mock.chr.rho[0][ir] + amag), 1e-14); + EXPECT_NEAR(rdn[ir], 0.5 * (mock.chr.rho[0][ir] - amag), 1e-14); + } +} + +// gga_grad=1 must ignore the global magnetization direction: with lsign_=true +// it has to give the same gradcorr result as gga_grad=0 with lsign_=false +TEST(GgaGradVxc, BuiltinGgaGrad1IgnoresGlobalAxis) +{ + XC_Functional::set_xc_type("PBE"); + const double hybrid_alpha = XC_Functional::get_hybrid_alpha(); + const double hse_omega = XC_Functional::get_hse_omega(); + + Ns4Charge mock_a(2); // lsign_ = true + double et1 = 0; + double vt1 = 0; + ModuleBase::matrix v1(4, gga_grad_nrxx); + std::vector dum; + XC_Functional::gradcorr(et1, vt1, v1, &mock_a.chr, &mock_a.rhopw, &mock_a.ucell, + dum, false, 4, true, false, 1, hybrid_alpha, hse_omega); + + Ns4Charge mock_b(2); + mock_b.ucell.magnet.lsign_ = false; + double et0 = 0; + double vt0 = 0; + ModuleBase::matrix v0(4, gga_grad_nrxx); + XC_Functional::gradcorr(et0, vt0, v0, &mock_b.chr, &mock_b.rhopw, &mock_b.ucell, + dum, false, 4, true, false, 0, hybrid_alpha, hse_omega); + + EXPECT_NEAR(et0, et1, 1e-12); + EXPECT_NEAR(vt0, vt1, 1e-12); + for (int is = 0; is < 4; ++is) + { + for (int ir = 0; ir < gga_grad_nrxx; ++ir) + { + EXPECT_NEAR(v0(is, ir), v1(is, ir), 1e-12); + } + } +} + + + +TEST(GgaGradVxc, LibxcNspin4NearSaturationDifferentiatesTheWeightedEnergy) +{ + constexpr double density_threshold = 1.0e-6; + Ns4Charge mock(0); + for (int ir = 0; ir < gga_grad_nrxx; ++ir) + { + mock.chr.rho[0][ir] = 0.45; + mock.chr.rho[1][ir] = 0.0; + mock.chr.rho[2][ir] = 0.0; + mock.chr.rho[3][ir] = 0.45 - density_threshold; + } + const std::vector func_ids = {XC_LDA_X}; + const int gga_grad_modes[] = {2}; + for (const int gga_grad : gga_grad_modes) + { + const std::function()> evaluate = [&, gga_grad]() + { + return XC_Functional_Libxc::v_xc_libxc(func_ids, + gga_grad_nrxx, + mock.ucell.omega, + mock.ucell.tpiba, + &mock.chr, + 4, + true, + false, + gga_grad, + nullptr, + 0.0, + 0.0); + }; + + const std::tuple reference = evaluate(); + const int components[] = {0, 3}; + const double steps[] = {8.0e-8, 4.0e-8, 2.0e-8}; + for (const int component : components) + { + double analytic = std::get<2>(reference)(component, 0); + Parallel_Reduce::reduce_pool(analytic); + if (std::getenv("ABACUS_XC_FD_TRACE") != nullptr) + { + std::cout << std::setprecision(17) + << "XC_SANITIZER_REFERENCE case=nspin4_near_saturation" + << " gga_grad=" << gga_grad + << " component=" << component + << " energy=" << std::get<0>(reference) + << " vtxc=" << std::get<1>(reference) + << " analytic=" << analytic << std::endl; + } + const double original = mock.chr.rho[component][0]; + for (const double step : steps) + { + mock.chr.rho[component][0] = original + step; + const double energy_plus = std::get<0>(evaluate()); + mock.chr.rho[component][0] = original - step; + const double energy_minus = std::get<0>(evaluate()); + mock.chr.rho[component][0] = original; + + const double finite_difference = (energy_plus - energy_minus) / (2.0 * step); + if (std::getenv("ABACUS_XC_FD_TRACE") != nullptr) + { + std::cout << std::setprecision(17) + << "XC_SANITIZER_FD case=nspin4_near_saturation" + << " gga_grad=" << gga_grad + << " component=" << component + << " eps=" << step + << " analytic=" << analytic + << " finite_difference=" << finite_difference + << " absolute_error=" << std::abs(analytic - finite_difference) + << std::endl; + } + const double scale = std::max(1.0, std::max(std::abs(analytic), + std::abs(finite_difference))); + EXPECT_NEAR(analytic, finite_difference, 2.0e-8 * scale) + << "gga_grad=" << gga_grad + << ", component=" << component + << ", step=" << step; + } + } + } +} + +// for LIBXC, gga_grad=0 and 1 both keep the original collinear algorithm +TEST(GgaGradVxc, LibxcZeroEqualsOne) +{ + const std::tuple r0 = run_vxc_nspin4("GGA_X_PBE+GGA_C_PBE", 1, 0); + const std::tuple r1 = run_vxc_nspin4("GGA_X_PBE+GGA_C_PBE", 1, 1); + expect_vxc_equal(r0, r1, 1e-12); +} int main(int argc, char **argv) { @@ -361,4 +735,4 @@ int main(int argc, char **argv) int result = RUN_ALL_TESTS(); MPI_Finalize(); return result; -} \ No newline at end of file +} diff --git a/source/source_hamilt/module_xc/test/test_xc_functional_ncgga_sf.cpp b/source/source_hamilt/module_xc/test/test_xc_functional_ncgga_sf.cpp new file mode 100644 index 00000000000..3166ec625fa --- /dev/null +++ b/source/source_hamilt/module_xc/test/test_xc_functional_ncgga_sf.cpp @@ -0,0 +1,295 @@ +#include "../xc_functional.h" +#include "../libxc_abacus.h" +#include "source_base/constants.h" +#include "source_base/parallel_comm.h" +#include "source_base/parallel_global.h" +#include "source_base/parallel_reduce.h" +#include "source_basis/module_pw/pw_basis.h" +#include "source_cell/unitcell.h" +#include "source_estate/module_charge/charge.h" +#include "gtest/gtest.h" +#include +#include +#include +#include +#include +#include + +// The focused target links PW/XC objects; the fixture owns charge storage. +Charge::Charge() {} +Charge::~Charge() {} +UnitCell::UnitCell() {} +UnitCell::~UnitCell() {} +Magnetism::Magnetism() {} +Magnetism::~Magnetism() {} +SepPot::SepPot() {} +SepPot::~SepPot() {} +Sep_Cell::Sep_Cell() noexcept {} +Sep_Cell::~Sep_Cell() noexcept {} + +namespace +{ +int test_rank = 0; +int test_size = 1; +const int noncollinear_spin = 4; +const int variational_gga = 2; +const double hybrid_alpha = 0.0; +const double hse_omega = 0.11; +typedef std::tuple VxcResult; +enum Branch { Smooth, Negative, Saturated, Radial }; + +// Each branch runs through both the built-in and LibXC production dispatch. +class RealPwNcgga : public testing::TestWithParam> +{ + protected: + ModulePW::PW_Basis pw; + UnitCell cell; + Charge charge; + std::array, 4> density; + std::array, 4> direction; + std::array pointers; + std::vector core; + std::vector> core_g; + + double sum(double value) + { + Parallel_Reduce::reduce_pool(value); + return value; + } + + void SetUp() override + { + pw.initmpi(test_size, test_rank, MPI_COMM_WORLD); + // Unequal MPI2 slabs and a finite PW cutoff exercise the projection. + pw.initgrids(7.0, ModuleBase::Matrix3(1, 0, 0, 0, 1, 0, 0, 0, 1), 24, 10, 9); + pw.initparameters(false, 80.0, 2, false); + pw.setuptransform(); + pw.collect_local_pw(); + cell.omega = pw.omega; + cell.tpiba = pw.tpiba; + cell.magnet.lsign_ = false; + for (int c = 0; c < noncollinear_spin; ++c) + { + density[c].resize(pw.nrxx); + direction[c].resize(pw.nrxx); + pointers[c] = density[c].data(); + } + core.resize(pw.nrxx); + core_g.resize(pw.npw); + charge.rhopw = &pw; + charge.nrxx = pw.nrxx; + charge.nxyz = pw.nxyz; + charge.ngmc = pw.npw; + charge.nspin = noncollinear_spin; + charge.rho = pointers.data(); + charge.rho_core = core.data(); + charge.rhog_core = core_g.data(); + const int branch = std::get<1>(GetParam()); + const double totals[] = {2.2, -1.55, 0.43, 0.030}; + const double magnitudes[] = {0.62, 0.45, 0.95, 4.0e-4}; + for (int ir = 0; ir < pw.nrxx; ++ir) + { + const double x = ModuleBase::TWO_PI * (ir / (pw.ny * pw.nplane)) / pw.nx; + const double y = ModuleBase::TWO_PI * ((ir / pw.nplane) % pw.ny) / pw.ny; + const double z = ModuleBase::TWO_PI * (ir % pw.nplane + pw.startz_current) / pw.nz; + density[0][ir] = totals[branch] * (1.0 + 0.08 * std::sin(x) + 0.03 * std::cos(y + z)); + const double magnitude = magnitudes[branch] * (1.0 + 0.1 * std::cos(5 * x) + 0.03 * std::sin(y)); + const double theta = 0.7 + 0.32 * std::sin(3 * x) + 0.18 * std::cos(7 * x) + 0.1 * std::cos(y); + const double phi = 0.4 + 0.27 * std::cos(4 * x) - 0.16 * std::sin(6 * x) + 0.09 * std::sin(z); + density[1][ir] = magnitude * std::sin(theta) * std::cos(phi); + density[2][ir] = magnitude * std::sin(theta) * std::sin(phi); + density[3][ir] = magnitude * std::cos(theta); + core[ir] = 0.05 * totals[branch] * (1.0 + 0.1 * std::cos(3 * x + z)); + if (branch == Negative) ASSERT_LT(density[0][ir] + core[ir], 0.0); + if (branch == Saturated) ASSERT_GT(magnitude, density[0][ir] + core[ir]); + if (branch == Radial) ASSERT_LT(magnitude, 1e-3); + for (int c = 0; c < noncollinear_spin; ++c) + direction[c][ir] = 0.1 + 0.2 * std::cos((c + 2) * x + y) + 0.1 * std::sin(z + c); + } + pw.real2recip(core.data(), core_g.data()); + XC_Functional::set_xc_type(std::get<0>(GetParam()) ? "GGA_X_PBE+GGA_C_PBE" : "PBE"); + } + + VxcResult evaluate() + { + return XC_Functional::v_xc(pw.nrxx, &charge, &cell, noncollinear_spin, + true, false, variational_gga, hybrid_alpha, hse_omega); + } + + double inner(const ModuleBase::matrix& potential, const bool perturbation) + { + double value = 0.0; + for (int c = 0; c < noncollinear_spin; ++c) + for (int ir = 0; ir < pw.nrxx; ++ir) + value += potential(c, ir) * (perturbation ? direction[c][ir] : density[c][ir]); + return sum(value) * pw.omega / pw.nxyz; + } +}; + +TEST_P(RealPwNcgga, GradientAndDivergenceAreAdjoints) +{ + std::vector> reciprocal(pw.npw); + std::vector> field(pw.nrxx); + std::vector> gradient(pw.nrxx); + std::vector divergence(pw.nrxx); + for (int ir = 0; ir < pw.nrxx; ++ir) + field[ir] = ModuleBase::Vector3(direction[1][ir], direction[2][ir], direction[3][ir]); + pw.real2recip(direction[0].data(), reciprocal.data()); + XC_Functional::grad_rho(reciprocal.data(), gradient.data(), &pw, pw.tpiba); + XC_Functional::grad_dot(field.data(), divergence.data(), &pw, pw.tpiba); + double identity = 0.0; + double norm = 0.0; + for (int ir = 0; ir < pw.nrxx; ++ir) + { + const double left = field[ir] * gradient[ir]; + const double right = divergence[ir] * direction[0][ir]; + identity += left + right; + norm += std::abs(left) + std::abs(right); + } + const double global_norm = sum(norm); + EXPECT_GT(global_norm, 1e-4); + EXPECT_LE(std::abs(sum(identity)), 5e-11 * global_norm); +} + +TEST_P(RealPwNcgga, DensityAndCoreDerivatives) +{ + const VxcResult reference = evaluate(); + const ModuleBase::matrix& potential = std::get<2>(reference); + EXPECT_NEAR(std::get<1>(reference), inner(potential, false), 2e-12 * std::max(1.0, std::abs(std::get<1>(reference)))); + // Four independent spin channels plus a core-density perturbation. + const bool radial = std::get<1>(GetParam()) == Radial; + const double step = radial ? 2e-4 : 2e-3; + const int refinements = radial ? 10 : 3; + for (int channel = 0; channel <= noncollinear_spin; ++channel) + { + SCOPED_TRACE(channel); + const bool is_core = channel == noncollinear_spin; + const int c = is_core ? 0 : channel; + std::vector& values = is_core ? core : density[c]; + const std::vector original = values; + double analytic = 0.0; + for (int ir = 0; ir < pw.nrxx; ++ir) + analytic += potential(c, ir) * direction[c][ir]; + analytic = sum(analytic) * pw.omega / pw.nxyz; + double previous = 0.0; + double best = 1.0; + for (int level = 0; level < refinements; ++level) + { + const double epsilon = step / (1 << level); + double energies[2]; + for (int sign = 0; sign < 2; ++sign) + { + for (int ir = 0; ir < pw.nrxx; ++ir) + values[ir] = original[ir] + (sign == 0 ? epsilon : -epsilon) * direction[c][ir]; + if (is_core) pw.real2recip(core.data(), core_g.data()); + energies[sign] = std::get<0>(evaluate()); + } + std::copy(original.begin(), original.end(), values.begin()); + if (is_core) pw.real2recip(core.data(), core_g.data()); + const double error = std::abs((energies[0] - energies[1]) / (2 * epsilon) - analytic); + const double scale = std::max(1.0, std::abs(analytic)); + if (level > 0 && level <= 2) EXPECT_LE(error, 0.4 * previous + 5e-9 * scale); + best = std::min(best, error); + if (level == refinements - 1) EXPECT_LE(best, 3e-8 * scale); + previous = error; + } + } +} + +TEST_P(RealPwNcgga, SixComponentStressDerivative) +{ + double energy = 0.0; + double vtxc = 0.0; + ModuleBase::matrix potential; + std::vector stress; + XC_Functional::gradcorr(energy, vtxc, potential, &charge, &pw, &cell, stress, + true, noncollinear_spin, true, false, variational_gga, hybrid_alpha, hse_omega); + ASSERT_EQ(stress.size(), 9U); + const std::vector> original(pw.gcar, pw.gcar + pw.npw); + for (int row = 0; row < 3; ++row) + for (int column = 0; column <= row; ++column) + { + SCOPED_TRACE(row * 3 + column); + const double analytic = sum(stress[row * 3 + column]) / pw.nxyz; + double previous = 0.0; + double best = 1.0; + for (int level = 0; level < 4; ++level) + { + const double epsilon = 2e-3 / (1 << level); + double energies[2]; + for (int sign = 0; sign < 2; ++sign) + { + // Inverse of I + epsilon e_row e_column, holding density fixed. + const double strain = sign == 0 ? epsilon : -epsilon; + const double factor = row == column ? strain / (1 + strain) : strain; + for (int ig = 0; ig < pw.npw; ++ig) + { + pw.gcar[ig] = original[ig]; + pw.gcar[ig][column] -= factor * original[ig][row]; + } + energies[sign] = std::get<0>(evaluate()); + } + std::copy(original.begin(), original.end(), pw.gcar); + const double error = std::abs(-(energies[0] - energies[1]) / (2 * epsilon * pw.omega) - analytic); + const double tolerance = 3e-4 * std::max(1e-10, std::abs(analytic)) + 1e-11; + if (level > 0 && previous > 1e-10 && error > 1e-10) + EXPECT_LE(error, 0.4 * previous + tolerance); + best = std::min(best, error); + if (level == 3) EXPECT_LE(best, tolerance); + previous = error; + } + } +} + +TEST_P(RealPwNcgga, SpinRotationAndInversion) +{ + const VxcResult reference = evaluate(); + // Cyclic permutation is a proper global rotation; inversion is checked separately. + for (int inversion = 0; inversion < 2; ++inversion) + { + const std::array, 4> original = density; + for (int c = 1; c < noncollinear_spin; ++c) + for (int ir = 0; ir < pw.nrxx; ++ir) + density[c][ir] = (inversion == 0 ? 1 : -1) * original[c % 3 + 1][ir]; + const VxcResult rotated = evaluate(); + EXPECT_NEAR(std::get<0>(reference), std::get<0>(rotated), 1e-10); + for (int c = 0; c < noncollinear_spin; ++c) + for (int ir = 0; ir < pw.nrxx; ++ir) + EXPECT_NEAR(std::get<2>(rotated)(c, ir), + (c == 0 ? 1 : (inversion == 0 ? 1 : -1)) + * std::get<2>(reference)(c == 0 ? 0 : c % 3 + 1, ir), 1e-11); + for (int c = 0; c < noncollinear_spin; ++c) + std::copy(original[c].begin(), original[c].end(), density[c].begin()); + } +} + +TEST_P(RealPwNcgga, ZeroMagnetizationIsFinite) +{ + for (int c = 1; c < noncollinear_spin; ++c) + std::fill(density[c].begin(), density[c].end(), 0.0); + const VxcResult result = evaluate(); + EXPECT_TRUE(std::isfinite(std::get<0>(result))); + for (int c = 1; c < noncollinear_spin; ++c) + for (int ir = 0; ir < pw.nrxx; ++ir) + EXPECT_NEAR(std::get<2>(result)(c, ir), 0.0, 1e-13); +} + +INSTANTIATE_TEST_SUITE_P(ProductionBranches, RealPwNcgga, + testing::Combine(testing::Bool(), testing::Values(Smooth, Negative, Saturated, Radial))); +} // namespace + +int main(int argc, char** argv) +{ + int threads = 1; + Parallel_Global::read_pal_param(argc, argv, test_size, threads, test_rank); + POOL_WORLD = MPI_COMM_WORLD; + KP_WORLD = MPI_COMM_NULL; + INT_BGROUP = MPI_COMM_NULL; + BP_WORLD = MPI_COMM_NULL; + GRID_WORLD = MPI_COMM_NULL; + DIAG_WORLD = MPI_COMM_NULL; + testing::InitGoogleTest(&argc, argv); + const int result = RUN_ALL_TESTS(); + Parallel_Global::finalize_mpi(); + return result; +} diff --git a/source/source_hamilt/module_xc/test/test_xc_ncgga_radial.cpp b/source/source_hamilt/module_xc/test/test_xc_ncgga_radial.cpp new file mode 100644 index 00000000000..3ad9a9d19f0 --- /dev/null +++ b/source/source_hamilt/module_xc/test/test_xc_ncgga_radial.cpp @@ -0,0 +1,278 @@ +#include "../xc_ncgga_radial.h" + +#include "gtest/gtest.h" +#include +#include +#include +#include +#include +#include + +namespace +{ + +double dot(const std::array& left, const std::array& right) +{ + double result = 0.0; + for (int component = 0; component < 3; ++component) + { + result += left[component] * right[component]; + } + return result; +} + +std::array shifted(const std::array& point, + const std::array& direction, + const double step) +{ + std::array result = point; + for (int component = 0; component < 3; ++component) + { + result[component] += step * direction[component]; + } + return result; +} + +double directional_hessian(const ModuleXC::NcggaRadialPoint& point, const std::array& direction) +{ + double result = 0.0; + for (int row = 0; row < 3; ++row) + { + for (int column = 0; column < 3; ++column) + { + result += direction[row] * point.jacobian(row, column) * direction[column]; + } + } + return result; +} + +ModuleXC::NcggaSpinMapPoint make_spin_map(const std::array& state, const double eta) +{ + const std::array magnetization = {{state[1], state[2], state[3]}}; + return ModuleXC::make_ncgga_spin_map_point(state[0], ModuleXC::make_ncgga_radial_point(magnetization, eta)); +} + +TEST(NcggaRadial, RejectsNonPositiveEta) +{ + const std::array magnetization = {{0.1, 0.0, 0.0}}; + EXPECT_THROW(ModuleXC::make_ncgga_radial_point(magnetization, 0.0), std::invalid_argument); + EXPECT_THROW(ModuleXC::make_ncgga_radial_point(magnetization, -0.5), std::invalid_argument); +} + +TEST(NcggaRadial, MatchesZeroInteriorJoinAndExteriorValues) +{ + const double eta = 0.5; + const ModuleXC::NcggaRadialPoint zero = ModuleXC::make_ncgga_radial_point({{0.0, 0.0, 0.0}}, eta); + EXPECT_DOUBLE_EQ(zero.value, 0.0); + for (int row = 0; row < 3; ++row) + { + EXPECT_DOUBLE_EQ(zero.direction[row], 0.0); + EXPECT_DOUBLE_EQ(zero.gradient[row], 0.0); + for (int column = 0; column < 3; ++column) + { + EXPECT_DOUBLE_EQ(zero.jacobian(row, column), 0.0); + } + } + + const ModuleXC::NcggaRadialPoint interior = ModuleXC::make_ncgga_radial_point({{0.5 * eta, 0.0, 0.0}}, eta); + EXPECT_DOUBLE_EQ(interior.value, 11.0 * eta / 32.0); + EXPECT_DOUBLE_EQ(interior.gradient[0], 23.0 / 16.0); + EXPECT_DOUBLE_EQ(interior.jacobian(0, 0), 3.0 / (2.0 * eta)); + EXPECT_DOUBLE_EQ(interior.jacobian(1, 1), 23.0 / (8.0 * eta)); + EXPECT_DOUBLE_EQ(interior.jacobian(2, 2), 23.0 / (8.0 * eta)); + + const ModuleXC::NcggaRadialPoint join = ModuleXC::make_ncgga_radial_point({{eta, 0.0, 0.0}}, eta); + EXPECT_DOUBLE_EQ(join.value, eta); + EXPECT_DOUBLE_EQ(join.gradient[0], 1.0); + EXPECT_DOUBLE_EQ(join.jacobian(0, 0), 0.0); + EXPECT_DOUBLE_EQ(join.jacobian(1, 1), 1.0 / eta); + EXPECT_DOUBLE_EQ(join.jacobian(2, 2), 1.0 / eta); + + const ModuleXC::NcggaRadialPoint exterior = ModuleXC::make_ncgga_radial_point({{3.0, 4.0, 0.0}}, eta); + EXPECT_DOUBLE_EQ(exterior.value, 5.0); + EXPECT_DOUBLE_EQ(exterior.gradient[0], 0.6); + EXPECT_DOUBLE_EQ(exterior.gradient[1], 0.8); + EXPECT_NEAR(exterior.jacobian(0, 0), 0.128, 1.0e-15); + EXPECT_NEAR(exterior.jacobian(0, 1), -0.096, 1.0e-15); + EXPECT_NEAR(exterior.jacobian(1, 1), 0.072, 1.0e-15); + EXPECT_NEAR(exterior.jacobian(2, 2), 0.2, 1.0e-15); +} + +TEST(NcggaRadial, CompactBranchConvergesC2AtJoin) +{ + const double eta = 0.5; + const double coarse_delta = 1.0e-3; + const double fine_delta = 0.5 * coarse_delta; + const ModuleXC::NcggaRadialPoint coarse + = ModuleXC::make_ncgga_radial_point({{eta * (1.0 - coarse_delta), 0.0, 0.0}}, eta); + const ModuleXC::NcggaRadialPoint fine + = ModuleXC::make_ncgga_radial_point({{eta * (1.0 - fine_delta), 0.0, 0.0}}, eta); + + EXPECT_LE(std::abs(fine.value - eta), 0.51 * std::abs(coarse.value - eta)); + EXPECT_LE(std::abs(fine.gradient[0] - 1.0), 0.26 * std::abs(coarse.gradient[0] - 1.0)); + EXPECT_LE(std::abs(fine.jacobian(0, 0)), 0.51 * std::abs(coarse.jacobian(0, 0))); + EXPECT_LE(std::abs(fine.jacobian(1, 1) - 1.0 / eta), 0.51 * std::abs(coarse.jacobian(1, 1) - 1.0 / eta)); +} + +TEST(NcggaRadial, CompactBranchConvergesC2AtOrigin) +{ + const double eta = 0.5; + const double coarse_radius = 1.0e-3 * eta; + const double fine_radius = 0.5 * coarse_radius; + const ModuleXC::NcggaRadialPoint coarse = ModuleXC::make_ncgga_radial_point({{coarse_radius, 0.0, 0.0}}, eta); + const ModuleXC::NcggaRadialPoint fine = ModuleXC::make_ncgga_radial_point({{fine_radius, 0.0, 0.0}}, eta); + + EXPECT_LE(fine.value, 0.13 * coarse.value); + EXPECT_LE(std::abs(fine.gradient[0]), 0.26 * std::abs(coarse.gradient[0])); + EXPECT_LE(std::abs(fine.jacobian(0, 0)), 0.51 * std::abs(coarse.jacobian(0, 0))); + EXPECT_LE(std::abs(fine.jacobian(1, 1)), 0.51 * std::abs(coarse.jacobian(1, 1))); +} + +TEST(NcggaRadial, GradientMatchesCentralFiniteDifference) +{ + const double eta = 0.5; + const double step = 2.0e-6 * eta; + const std::array direction = {{0.31, -0.27, 0.19}}; + const std::array interior = {{0.12, -0.08, 0.05}}; + const std::array exterior = {{0.62, -0.41, 0.23}}; + const std::array join = {{eta, 0.0, 0.0}}; + const std::array, 3> points = {{interior, exterior, join}}; + + for (std::size_t sample = 0; sample < points.size(); ++sample) + { + const ModuleXC::NcggaRadialPoint point = ModuleXC::make_ncgga_radial_point(points[sample], eta); + const double finite_difference + = (ModuleXC::make_ncgga_radial_point(shifted(points[sample], direction, step), eta).value + - ModuleXC::make_ncgga_radial_point(shifted(points[sample], direction, -step), eta).value) + / (2.0 * step); + EXPECT_NEAR(finite_difference, dot(point.gradient, direction), 2.0e-10); + } +} + +TEST(NcggaRadial, JacobianMatchesGradientFiniteDifference) +{ + const double eta = 0.5; + const double step = 2.0e-6 * eta; + const std::array direction = {{0.31, -0.27, 0.19}}; + const std::array interior = {{0.12, -0.08, 0.05}}; + const std::array exterior = {{0.62, -0.41, 0.23}}; + const std::array, 2> points = {{interior, exterior}}; + + for (std::size_t sample = 0; sample < points.size(); ++sample) + { + const ModuleXC::NcggaRadialPoint point = ModuleXC::make_ncgga_radial_point(points[sample], eta); + const ModuleXC::NcggaRadialPoint plus + = ModuleXC::make_ncgga_radial_point(shifted(points[sample], direction, step), eta); + const ModuleXC::NcggaRadialPoint minus + = ModuleXC::make_ncgga_radial_point(shifted(points[sample], direction, -step), eta); + for (int row = 0; row < 3; ++row) + { + double analytic = 0.0; + for (int column = 0; column < 3; ++column) + { + analytic += point.jacobian(row, column) * direction[column]; + EXPECT_NEAR(point.jacobian(row, column), point.jacobian(column, row), 1.0e-15); + } + const double finite_difference = (plus.gradient[row] - minus.gradient[row]) / (2.0 * step); + EXPECT_NEAR(finite_difference, analytic, 2.0e-9); + } + } +} + +TEST(NcggaRadial, DirectionalSecondDerivativeMatchesValueFiniteDifference) +{ + const double eta = 0.5; + const double step = 2.0e-4 * eta; + const std::array direction = {{0.31, -0.27, 0.19}}; + const std::array interior = {{0.12, -0.08, 0.05}}; + const std::array exterior = {{0.62, -0.41, 0.23}}; + const std::array, 2> points = {{interior, exterior}}; + + for (std::size_t sample = 0; sample < points.size(); ++sample) + { + const ModuleXC::NcggaRadialPoint point = ModuleXC::make_ncgga_radial_point(points[sample], eta); + const double plus = ModuleXC::make_ncgga_radial_point(shifted(points[sample], direction, step), eta).value; + const double minus = ModuleXC::make_ncgga_radial_point(shifted(points[sample], direction, -step), eta).value; + const double finite_difference = (plus - 2.0 * point.value + minus) / (step * step); + EXPECT_NEAR(finite_difference, directional_hessian(point, direction), 2.0e-7); + } +} + +TEST(NcggaRadial, SpinDensityMapJacobianMatchesFiniteDifferenceInEachBranch) +{ + const double eta = 0.5; + const double step = 2.0e-7; + const std::array, 5> states = {{{{1.20, 0.30, -0.20, 0.10}}, + {{-1.20, 0.30, -0.20, 0.10}}, + {{0.20, 0.42, -0.31, 0.16}}, + {{-0.20, 0.42, -0.31, 0.16}}, + {{0.90, 0.12, -0.08, 0.05}}}}; + + for (std::size_t sample = 0; sample < states.size(); ++sample) + { + const ModuleXC::NcggaSpinMapPoint point = make_spin_map(states[sample], eta); + EXPECT_GE(point.spin_density[0], 0.0); + EXPECT_GE(point.spin_density[1], 0.0); + EXPECT_NEAR(point.spin_density[0] + point.spin_density[1], std::abs(states[sample][0]), 1.0e-15); + EXPECT_NEAR(point.spin_density[0] - point.spin_density[1], point.clipped_magnitude, 1.0e-15); + + double maximum_error = 0.0; + for (int channel = 0; channel < 4; ++channel) + { + std::array plus_state = states[sample]; + std::array minus_state = states[sample]; + plus_state[channel] += step; + minus_state[channel] -= step; + const ModuleXC::NcggaSpinMapPoint plus = make_spin_map(plus_state, eta); + const ModuleXC::NcggaSpinMapPoint minus = make_spin_map(minus_state, eta); + for (int spin = 0; spin < 2; ++spin) + { + const double finite_difference = (plus.spin_density[spin] - minus.spin_density[spin]) / (2.0 * step); + maximum_error = std::max(maximum_error, std::abs(finite_difference - point.jacobian(spin, channel))); + EXPECT_NEAR(finite_difference, point.jacobian(spin, channel), 3.0e-9) + << "sample=" << sample << " spin=" << spin << " channel=" << channel; + } + } + std::cout << std::setprecision(17) << "NCGGA_SPIN_MAP_FD sample=" << sample + << " total_density=" << states[sample][0] << " absolute_density=" << point.absolute_density + << " radial_value=" << point.radial.value << " clipped_magnitude=" << point.clipped_magnitude + << " saturated=" << point.saturated << " step=" << step << " maximum_error=" << maximum_error << '\n'; + } +} + +TEST(NcggaRadial, SpinDensityMapDefinesAbsAndSaturationKinkConventions) +{ + const double eta = 0.1; + const ModuleXC::NcggaRadialPoint nonzero_radial = ModuleXC::make_ncgga_radial_point({{0.30, 0.0, 0.0}}, eta); + const ModuleXC::NcggaSpinMapPoint zero_density = ModuleXC::make_ncgga_spin_map_point(0.0, nonzero_radial); + EXPECT_TRUE(zero_density.saturated); + EXPECT_DOUBLE_EQ(zero_density.spin_density[0], 0.0); + EXPECT_DOUBLE_EQ(zero_density.spin_density[1], 0.0); + for (int spin = 0; spin < 2; ++spin) + { + for (int channel = 0; channel < 4; ++channel) + { + EXPECT_DOUBLE_EQ(zero_density.jacobian(spin, channel), 0.0); + } + } + + const ModuleXC::NcggaRadialPoint equality_radial = ModuleXC::make_ncgga_radial_point({{0.50, 0.0, 0.0}}, eta); + const ModuleXC::NcggaSpinMapPoint positive = ModuleXC::make_ncgga_spin_map_point(0.50, equality_radial); + const ModuleXC::NcggaSpinMapPoint negative = ModuleXC::make_ncgga_spin_map_point(-0.50, equality_radial); + EXPECT_TRUE(positive.saturated); + EXPECT_TRUE(negative.saturated); + EXPECT_DOUBLE_EQ(positive.jacobian(0, 0), 1.0); + EXPECT_DOUBLE_EQ(positive.jacobian(1, 0), 0.0); + EXPECT_DOUBLE_EQ(negative.jacobian(0, 0), -1.0); + EXPECT_DOUBLE_EQ(negative.jacobian(1, 0), 0.0); + for (int spin = 0; spin < 2; ++spin) + { + for (int channel = 1; channel < 4; ++channel) + { + EXPECT_DOUBLE_EQ(positive.jacobian(spin, channel), 0.0); + EXPECT_DOUBLE_EQ(negative.jacobian(spin, channel), 0.0); + } + } +} + +} // namespace diff --git a/source/source_hamilt/module_xc/xc_functional.h b/source/source_hamilt/module_xc/xc_functional.h index 4dce5e3e99f..f72510850da 100644 --- a/source/source_hamilt/module_xc/xc_functional.h +++ b/source/source_hamilt/module_xc/xc_functional.h @@ -72,6 +72,7 @@ class XC_Functional const int nspin, const bool domag, const bool domag_z, + const int gga_grad, const double hybrid_alpha, const double hse_omega); @@ -262,6 +263,7 @@ class XC_Functional const int nspin, const bool domag, const bool domag_z, + const int gga_grad, const double hybrid_alpha, const double hse_omega); diff --git a/source/source_hamilt/module_xc/xc_functional_ncgga_sf.cpp b/source/source_hamilt/module_xc/xc_functional_ncgga_sf.cpp new file mode 100644 index 00000000000..92808791ca9 --- /dev/null +++ b/source/source_hamilt/module_xc/xc_functional_ncgga_sf.cpp @@ -0,0 +1,347 @@ +#include "xc_functional_ncgga_sf.h" + +#include "source_base/parallel_reduce.h" +#include "source_base/timer.h" +#include "source_base/vector3.h" +#include "source_basis/module_pw/pw_basis.h" +#include "source_estate/module_charge/charge.h" +#include "xc_functional.h" +#include "xc_ncgga_radial.h" + +#include +#include +#include + +namespace ModuleXC +{ +namespace NCGGA_SF_Builtin +{ + +std::tuple v_xc_ncgga_sf_builtin(const int& nrxx, + const double& omega, + const double tpiba, + const Charge* const chr) +{ + ModuleBase::TITLE("XC_Functional", "v_xc_ncgga_sf_builtin"); + ModuleBase::timer::start("XC_Functional", "v_xc_ncgga_sf_builtin"); + + // Regularized projected local-collinear energy: rho_s = N_s(n + rho_core, m), + // g_s = sum_A J_sA G_h x_A. Reverse the same discrete graph, including + // both -D_h(sum_s J_sA h_s) and the local Hessian response of N_s. + ModulePW::PW_Basis* rhopw = chr->rhopw; + const int npw = rhopw->npw; + const double e2 = ModuleBase::e2; + constexpr double vanishing = 1e-10; + constexpr double epsr = 1e-6; + const bool is_gga = (XC_Functional::get_func_type() == 2 || XC_Functional::get_func_type() == 4); + + std::vector rhotmp1(nrxx); + std::vector rhotmp2(nrxx); + std::vector spin_map(nrxx); + + for (int ir = 0; ir < nrxx; ++ir) + { + const double mx = chr->rho[1][ir]; + const double my = chr->rho[2][ir]; + const double mz = chr->rho[3][ir]; + + const std::array magnetization = {{mx, my, mz}}; + const NcggaRadialPoint radial = make_ncgga_radial_point(magnetization, ncgga_lca_radial_eta()); + spin_map[ir] = make_ncgga_spin_map_point(chr->rho[0][ir] + chr->rho_core[ir], radial); + rhotmp1[ir] = spin_map[ir].spin_density[0]; + rhotmp2[ir] = spin_map[ir].spin_density[1]; + } + + std::vector> rhogsum1(npw); + std::vector> tmp_recip(npw); + rhopw->real2recip(chr->rho[0], rhogsum1.data()); + for (int ig = 0; ig < npw; ++ig) + rhogsum1[ig] += chr->rhog_core[ig]; + + std::vector> gdr1(nrxx); + std::vector> gdr2(nrxx); + std::vector> grad_rho(nrxx); + std::array>, 3> grad_m; + for (int mu = 0; mu < 3; ++mu) + { + grad_m[mu].resize(nrxx); + } + std::vector> gdr_mag(nrxx); + XC_Functional::grad_rho(rhogsum1.data(), gdr1.data(), rhopw, tpiba); + + for (int ir = 0; ir < nrxx; ++ir) + { + grad_rho[ir] = gdr1[ir]; + + gdr1[ir] = spin_map[ir].jacobian(0, 0) * grad_rho[ir]; + gdr2[ir] = spin_map[ir].jacobian(1, 0) * grad_rho[ir]; + } + for (int is = 1; is <= 3; ++is) + { + rhopw->real2recip(chr->rho[is], tmp_recip.data()); + XC_Functional::grad_rho(tmp_recip.data(), gdr_mag.data(), rhopw, tpiba); + grad_m[is - 1] = gdr_mag; + for (int ir = 0; ir < nrxx; ++ir) + { + + gdr1[ir] += spin_map[ir].jacobian(0, is) * gdr_mag[ir]; + gdr2[ir] += spin_map[ir].jacobian(1, is) * gdr_mag[ir]; + } + } + + double etxc = 0; + double vtxc = 0; + ModuleBase::matrix v(4, nrxx); + + for (int ir = 0; ir < nrxx; ++ir) + { + const double arho = spin_map[ir].absolute_density; + if (arho <= vanishing) + continue; + + double zeta = spin_map[ir].clipped_magnitude / arho; + if (std::abs(zeta) > 1.0) + zeta = (zeta > 0) ? 1.0 : -1.0; + double exc = 0; + double vxc[2] = {0, 0}; + XC_Functional::xc_spin(arho, zeta, exc, vxc[0], vxc[1]); + + for (int channel = 0; channel < 4; ++channel) + { + v(channel, ir) + = e2 * (spin_map[ir].jacobian(0, channel) * vxc[0] + spin_map[ir].jacobian(1, channel) * vxc[1]); + } + + etxc += e2 * exc * arho; + } + + // Step 4: GGA contribution and variational divergence correction. + if (is_gga) + { + double etxcgc = 0; + std::vector vup_gga(nrxx, 0); + std::vector vdw_gga(nrxx, 0); + std::vector> h1(nrxx); + std::vector> h2(nrxx); + for (int ir = 0; ir < nrxx; ++ir) + { + double sx = 0; + double v1xup = 0; + double v1xdw = 0; + double v2xup = 0; + double v2xdw = 0; + double sc = 0; + double v1cup = 0; + double v1cdw = 0; + double v2c = 0; + double grho2a = gdr1[ir] * gdr1[ir]; + double grho2b = gdr2[ir] * gdr2[ir]; + + const double rh = rhotmp1[ir] + rhotmp2[ir]; + + XC_Functional::gcx_spin(rhotmp1[ir], rhotmp2[ir], grho2a, grho2b, sx, v1xup, v1xdw, v2xup, v2xdw); + + if (rh > epsr) + { + const double zeta_input = std::fabs((rhotmp1[ir] - rhotmp2[ir]) / rh); + double zeta = zeta_input; + const double grh2 = (gdr1[ir] + gdr2[ir]) * (gdr1[ir] + gdr2[ir]); + XC_Functional::gcc_spin(rh, zeta, grh2, sc, v1cup, v1cdw, v2c); + if (zeta_input > 1.0 - epsr) + { + // gcc_spin evaluates this branch at a fixed clipped zeta. + // Reverse that actual branch instead of differentiating + // through the discarded input polarization. + const double fixed_zeta_density_derivative = 0.5 * ((1.0 + zeta) * v1cup + (1.0 - zeta) * v1cdw); + v1cup = fixed_zeta_density_derivative; + v1cdw = fixed_zeta_density_derivative; + } + } + + vup_gga[ir] = e2 * (v1xup + v1cup); + vdw_gga[ir] = e2 * (v1xdw + v1cdw); + + const double v2cup = v2c; + const double v2cdw = v2c; + const double v2cud = v2c; + h1[ir] = e2 * ((v2xup + v2cup) * gdr1[ir] + v2cud * gdr2[ir]); + h2[ir] = e2 * ((v2xdw + v2cdw) * gdr2[ir] + v2cud * gdr1[ir]); + + etxcgc += e2 * (sx + sc); + } + + for (int ir = 0; ir < nrxx; ++ir) + { + + for (int channel = 0; channel < 4; ++channel) + { + v(channel, ir) += spin_map[ir].jacobian(0, channel) * vup_gga[ir] + + spin_map[ir].jacobian(1, channel) * vdw_gga[ir]; + } + } + + std::vector dh(nrxx); + std::vector> tmp_h(nrxx); + + // Exact reverse of g_s=sum_A J_sA G_h(x_A): + // v_B = -D_h(sum_s J_sB h_s) + // + sum_s,A dJ_sA/dx_B h_s.G_h(x_A). + for (int channel = 0; channel < 4; ++channel) + { + for (int ir = 0; ir < nrxx; ++ir) + { + tmp_h[ir] = spin_map[ir].jacobian(0, channel) * h1[ir] + spin_map[ir].jacobian(1, channel) * h2[ir]; + } + XC_Functional::grad_dot(tmp_h.data(), dh.data(), rhopw, tpiba); + for (int ir = 0; ir < nrxx; ++ir) + { + v(channel, ir) -= dh[ir]; + if (channel == 0 || spin_map[ir].saturated) + { + continue; + } + const ModuleBase::Vector3 spin_flux = 0.5 * (h1[ir] - h2[ir]); + double local_response = 0.0; + for (int nu = 0; nu < 3; ++nu) + { + local_response += spin_map[ir].radial.jacobian(nu, channel - 1) * (spin_flux * grad_m[nu][ir]); + } + v(channel, ir) += local_response; + } + } + + etxc += etxcgc; + } + + // vtxc uses the same completed four-component potential returned to the + // caller. This unifies the bookkeeping for both modes. + vtxc = 0.0; + for (int ir = 0; ir < nrxx; ++ir) + { + for (int is = 0; is < 4; ++is) + { + vtxc += v(is, ir) * chr->rho[is][ir]; + } + } + +#ifdef __MPI + Parallel_Reduce::reduce_pool(etxc); + Parallel_Reduce::reduce_pool(vtxc); +#endif + etxc *= omega / rhopw->nxyz; + vtxc *= omega / rhopw->nxyz; + + ModuleBase::timer::end("XC_Functional", "v_xc_ncgga_sf_builtin"); + return std::make_tuple(etxc, vtxc, std::move(v)); +} + +void gradcorr_ncgga_lca_builtin(const Charge* const chr, + ModulePW::PW_Basis* rhopw, + const double tpiba, + std::vector& stress_gga) +{ + stress_gga.assign(9, 0.0); + + const int nrxx = rhopw->nrxx; + const int npw = rhopw->npw; + const double e2 = ModuleBase::e2; + constexpr double epsr = 1.0e-6; + + // Rebuild the same complete local map used by the gga_grad=2 energy: + // rho_s = N_s(n + rho_core, m), + // g_s = sum_A dN_s/dx_A G_h x_A. + // Metric differentiation keeps the real-grid values x_A fixed, so the + // map Jacobian is unchanged and every G_h x_A transforms covariantly. + std::vector spin_map(nrxx); + std::array>, 4> field_gradient; + for (int channel = 0; channel < 4; ++channel) + { + field_gradient[channel].resize(nrxx); + } + std::array>, 2> spin_gradient; + for (int spin = 0; spin < 2; ++spin) + { + spin_gradient[spin].resize(nrxx); + } + + std::vector> reciprocal(npw); + rhopw->real2recip(chr->rho[0], reciprocal.data()); + for (int ig = 0; ig < npw; ++ig) + { + reciprocal[ig] += chr->rhog_core[ig]; + } + XC_Functional::grad_rho(reciprocal.data(), field_gradient[0].data(), rhopw, tpiba); + + for (int channel = 1; channel < 4; ++channel) + { + rhopw->real2recip(chr->rho[channel], reciprocal.data()); + XC_Functional::grad_rho(reciprocal.data(), field_gradient[channel].data(), rhopw, tpiba); + } + + for (int ir = 0; ir < nrxx; ++ir) + { + const std::array magnetization = {{chr->rho[1][ir], chr->rho[2][ir], chr->rho[3][ir]}}; + spin_map[ir] = make_ncgga_spin_map_point(chr->rho[0][ir] + chr->rho_core[ir], + make_ncgga_radial_point(magnetization, ncgga_lca_radial_eta())); + for (int spin = 0; spin < 2; ++spin) + { + for (int channel = 0; channel < 4; ++channel) + { + spin_gradient[spin][ir] += spin_map[ir].jacobian(spin, channel) * field_gradient[channel][ir]; + } + } + } + + for (int ir = 0; ir < nrxx; ++ir) + { + const double rho_up = spin_map[ir].spin_density[0]; + const double rho_down = spin_map[ir].spin_density[1]; + const ModuleBase::Vector3& grad_up = spin_gradient[0][ir]; + const ModuleBase::Vector3& grad_down = spin_gradient[1][ir]; + + double sx = 0.0; + double v1xup = 0.0; + double v1xdw = 0.0; + double v2xup = 0.0; + double v2xdw = 0.0; + XC_Functional::gcx_spin(rho_up, + rho_down, + grad_up * grad_up, + grad_down * grad_down, + sx, + v1xup, + v1xdw, + v2xup, + v2xdw); + + double sc = 0.0; + double v1cup = 0.0; + double v1cdw = 0.0; + double v2c = 0.0; + const double rho = rho_up + rho_down; + if (rho > epsr) + { + double zeta = std::fabs((rho_up - rho_down) / rho); + const ModuleBase::Vector3 grad_rho = grad_up + grad_down; + XC_Functional::gcc_spin(rho, zeta, grad_rho * grad_rho, sc, v1cup, v1cdw, v2c); + } + + const ModuleBase::Vector3 h_up = e2 * ((v2xup + v2c) * grad_up + v2c * grad_down); + const ModuleBase::Vector3 h_down = e2 * ((v2xdw + v2c) * grad_down + v2c * grad_up); + const double grad_up_component[3] = {grad_up.x, grad_up.y, grad_up.z}; + const double grad_down_component[3] = {grad_down.x, grad_down.y, grad_down.z}; + const double h_up_component[3] = {h_up.x, h_up.y, h_up.z}; + const double h_down_component[3] = {h_down.x, h_down.y, h_down.z}; + for (int row = 0; row < 3; ++row) + { + for (int column = 0; column <= row; ++column) + { + stress_gga[row * 3 + column] += h_up_component[row] * grad_up_component[column] + + h_down_component[row] * grad_down_component[column]; + } + } + } +} + +} // namespace NCGGA_SF_Builtin +} // namespace ModuleXC diff --git a/source/source_hamilt/module_xc/xc_functional_ncgga_sf.h b/source/source_hamilt/module_xc/xc_functional_ncgga_sf.h new file mode 100644 index 00000000000..c4d71a28ad7 --- /dev/null +++ b/source/source_hamilt/module_xc/xc_functional_ncgga_sf.h @@ -0,0 +1,37 @@ +#ifndef XC_FUNCTIONAL_NCGGA_SF_H +#define XC_FUNCTIONAL_NCGGA_SF_H + +#include "source_base/matrix.h" + +#include +#include + +class Charge; +namespace ModulePW +{ +class PW_Basis; +} + +namespace ModuleXC +{ +namespace NCGGA_SF_Builtin +{ + +// Exact discrete reverse of the regularized projected LCA graph (gga_grad=2). +std::tuple v_xc_ncgga_sf_builtin(const int& nrxx, + const double& omega, + const double tpiba, + const Charge* const chr); + +// Gradient-metric stress of the exact gga_grad=2 projected-LCA graph. The +// returned lower triangle is the unnormalised real-grid sum; Stress_Func +// applies the existing pool reduction and 1/nxyz normalisation. +void gradcorr_ncgga_lca_builtin(const Charge* const chr, + ModulePW::PW_Basis* rhopw, + const double tpiba, + std::vector& stress_gga); + +} // namespace NCGGA_SF_Builtin +} // namespace ModuleXC + +#endif diff --git a/source/source_hamilt/module_xc/xc_grad.cpp b/source/source_hamilt/module_xc/xc_grad.cpp index c7bf9a20136..2df043824a8 100644 --- a/source/source_hamilt/module_xc/xc_grad.cpp +++ b/source/source_hamilt/module_xc/xc_grad.cpp @@ -14,11 +14,13 @@ // noncolin_rho. #include "xc_functional.h" +#include "xc_functional_ncgga_sf.h" #include "xc_grad_internal.h" #include "source_base/timer.h" #ifdef __LIBXC #include +#include "libxc_abacus.h" #endif void XC_Functional::gradcorr( @@ -33,6 +35,7 @@ void XC_Functional::gradcorr( const int nspin, const bool domag, const bool domag_z, + const int gga_grad, const double hybrid_alpha_in, const double hse_omega_in) { @@ -63,6 +66,32 @@ void XC_Functional::gradcorr( return; } + if (is_stress && !use_libxc && nspin == 4 && (domag || domag_z) && gga_grad == 2) + { + ModuleXC::NCGGA_SF_Builtin::gradcorr_ncgga_lca_builtin( + chr, rhopw, ucell->tpiba, stress_gga); + return; + } + + + +#ifdef __LIBXC + if (is_stress && use_libxc && nspin == 4 && (domag || domag_z) + && gga_grad == 2) + { + XC_Functional_Libxc::gradcorr_ncgga_sf_libxc( + func_id, + rhopw->nrxx, + ucell->tpiba, + chr, + &scaling_factor_xc, + hybrid_alpha_in, + hse_omega_in, + stress_gga); + return; + } +#endif + bool igcc_is_lyp = false; // func_id may hold a single entry (e.g. PBE0 -> {XC_HYB_GGA_XC_PBEH}), so guard the index. if( func_id.size() > 1 && func_id[1] == XC_GGA_C_LYP) @@ -108,6 +137,7 @@ void XC_Functional::gradcorr( params.igcc_is_lyp = igcc_is_lyp; params.domag = domag; params.domag_z = domag_z; + params.gga_grad = gga_grad; params.hybrid_alpha = hybrid_alpha_in; params.hse_omega = hse_omega_in; params.use_libxc = use_libxc; diff --git a/source/source_hamilt/module_xc/xc_grad_internal.h b/source/source_hamilt/module_xc/xc_grad_internal.h index 8d0e54abedd..76d27832a89 100644 --- a/source/source_hamilt/module_xc/xc_grad_internal.h +++ b/source/source_hamilt/module_xc/xc_grad_internal.h @@ -41,6 +41,7 @@ struct GradCorrParams bool igcc_is_lyp; bool domag; bool domag_z; + int gga_grad; double hybrid_alpha; double hse_omega; bool use_libxc; diff --git a/source/source_hamilt/module_xc/xc_grad_prepare.cpp b/source/source_hamilt/module_xc/xc_grad_prepare.cpp index 756e7c9ff2e..b60d0d45207 100644 --- a/source/source_hamilt/module_xc/xc_grad_prepare.cpp +++ b/source/source_hamilt/module_xc/xc_grad_prepare.cpp @@ -150,7 +150,10 @@ void gradcorr_prepare_rho( } } } - XC_Functional::noncolin_rho(buf.rhotmp1.data(), buf.rhotmp2.data(), buf.neg.data(), chr->rho, rhopw->nrxx, ucell->magnet.ux_, ucell->magnet.lsign_); + // Mode 1 ignores the global quantization axis to remain continuous + // when the magnetic moments tilt away from a collinear state. + const bool use_global_axis = ucell->magnet.lsign_ && params.gga_grad != 1; + XC_Functional::noncolin_rho(buf.rhotmp1.data(), buf.rhotmp2.data(), buf.neg.data(), chr->rho, rhopw->nrxx, ucell->magnet.ux_, use_global_axis); rhopw->real2recip(buf.rhotmp1.data(), buf.rhogsum1.data()); rhopw->real2recip(buf.rhotmp2.data(), buf.rhogsum2.data()); #ifdef _OPENMP diff --git a/source/source_hamilt/module_xc/xc_ncgga_radial.cpp b/source/source_hamilt/module_xc/xc_ncgga_radial.cpp new file mode 100644 index 00000000000..d01f4159b09 --- /dev/null +++ b/source/source_hamilt/module_xc/xc_ncgga_radial.cpp @@ -0,0 +1,95 @@ +#include "xc_ncgga_radial.h" + +#include +#include +#include + +namespace ModuleXC +{ + +double NcggaRadialPoint::jacobian(const int row, const int column) const +{ + const double identity = (row == column) ? 1.0 : 0.0; + return transverse_hessian * identity + (radial_hessian - transverse_hessian) * direction[row] * direction[column]; +} + +NcggaRadialPoint make_ncgga_radial_point(const std::array& magnetization, const double eta) +{ + if (!(eta > 0.0)) + { + throw std::invalid_argument("noncollinear GGA radial eta must be positive"); + } + + NcggaRadialPoint point; + const double magnitude = std::sqrt(magnetization[0] * magnetization[0] + magnetization[1] * magnetization[1] + + magnetization[2] * magnetization[2]); + if (magnitude == 0.0) + { + return point; + } + + for (int component = 0; component < 3; ++component) + { + point.direction[component] = magnetization[component] / magnitude; + } + + if (magnitude < eta) + { + const double x = magnitude / eta; + const double x2 = x * x; + const double x3 = x2 * x; + point.value = eta * x3 * (3.0 * x2 - 8.0 * x + 6.0); + point.transverse_hessian = x * (15.0 * x2 - 32.0 * x + 18.0) / eta; + point.radial_hessian = x * (60.0 * x2 - 96.0 * x + 36.0) / eta; + } + else + { + point.value = magnitude; + point.transverse_hessian = 1.0 / magnitude; + point.radial_hessian = 0.0; + } + + for (int component = 0; component < 3; ++component) + { + point.gradient[component] = point.transverse_hessian * magnetization[component]; + } + return point; +} + +double ncgga_lca_radial_eta() +{ + return 1.0e-3; +} + +double NcggaSpinMapPoint::jacobian(const int spin, const int channel) const +{ + if (channel == 0) + { + if (saturated) + { + return spin == 0 ? density_sign : 0.0; + } + return 0.5 * density_sign; + } + if (saturated) + { + return 0.0; + } + const double spin_sign = spin == 0 ? 0.5 : -0.5; + return spin_sign * radial.gradient[channel - 1]; +} + +NcggaSpinMapPoint make_ncgga_spin_map_point(const double total_density, const NcggaRadialPoint& radial) +{ + NcggaSpinMapPoint point; + point.radial = radial; + point.absolute_density = std::abs(total_density); + point.clipped_magnitude = std::min(radial.value, point.absolute_density); + point.spin_density[0] = 0.5 * (point.absolute_density + point.clipped_magnitude); + point.spin_density[1] = 0.5 * (point.absolute_density - point.clipped_magnitude); + point.density_sign = total_density > 0.0 ? 1.0 : total_density < 0.0 ? -1.0 : 0.0; + point.saturated = !(radial.value < point.absolute_density); + return point; +} + +} // namespace ModuleXC diff --git a/source/source_hamilt/module_xc/xc_ncgga_radial.h b/source/source_hamilt/module_xc/xc_ncgga_radial.h new file mode 100644 index 00000000000..13da072ed22 --- /dev/null +++ b/source/source_hamilt/module_xc/xc_ncgga_radial.h @@ -0,0 +1,60 @@ +#ifndef XC_NCGGA_RADIAL_H +#define XC_NCGGA_RADIAL_H + +#include + +namespace ModuleXC +{ + +struct NcggaRadialPoint +{ + // Value, gradient, and Hessian of the same radial scalar map. At zero, + // direction is represented by the zero vector because the scalar map has + // a unique zero gradient and zero Hessian there. + double value = 0.0; + std::array direction = {{0.0, 0.0, 0.0}}; + std::array gradient = {{0.0, 0.0, 0.0}}; + double transverse_hessian = 0.0; + double radial_hessian = 0.0; + + double jacobian(const int row, const int column) const; +}; + +// For r = |magnetization| and x = r / eta, the returned scalar is +// eta * x^3 * (3 x^2 - 8 x + 6), r < eta, +// r, r >= eta. +// The splice is C2 at both r = 0 and r = eta. Eta is explicit so this +// mathematical primitive does not choose policy for any XC mode. +NcggaRadialPoint make_ncgga_radial_point(const std::array& magnetization, const double eta); + +// The C2 regularization scale is part of the gga_grad=2 LCA functional, not a +// divide-by-zero guard. Other noncollinear modes do not inherit this policy. +double ncgga_lca_radial_eta(); + +struct NcggaSpinMapPoint +{ + NcggaRadialPoint radial; + double absolute_density = 0.0; + double clipped_magnitude = 0.0; + std::array spin_density = {{0.0, 0.0}}; + double density_sign = 0.0; + bool saturated = true; + + // spin is 0/1 for up/down; channel is 0 for the raw total density and + // 1..3 for mx,my,mz. + double jacobian(const int spin, const int channel) const; +}; + +// Compose a raw total density t=n+rho_core with one radial magnetization map: +// a = |t|, +// c = min(radial.value, a), +// rho_up = (a+c)/2, +// rho_down= (a-c)/2. +// The Jacobian follows that exact graph. At the abs kink t=0 it selects zero; +// at the clipping kink radial.value=a it selects the saturated branch. Eta and +// the radial policy remain explicit in make_ncgga_radial_point. +NcggaSpinMapPoint make_ncgga_spin_map_point(const double total_density, const NcggaRadialPoint& radial); + +} // namespace ModuleXC + +#endif diff --git a/source/source_hamilt/module_xc/xc_pot.cpp b/source/source_hamilt/module_xc/xc_pot.cpp index dc6a13dabdf..f8b54bde2dd 100644 --- a/source/source_hamilt/module_xc/xc_pot.cpp +++ b/source/source_hamilt/module_xc/xc_pot.cpp @@ -8,6 +8,8 @@ #include "source_base/timer.h" #include "xc_functional.h" +#include "xc_functional_ncgga_sf.h" + #ifdef __LIBXC #include "libxc_abacus.h" #ifdef __EXX @@ -15,6 +17,7 @@ #endif #endif + // [etxc, vtxc, v] = XC_Functional::v_xc(...) std::tuple XC_Functional::v_xc( const int& nrxx, @@ -23,11 +26,13 @@ std::tuple XC_Functional::v_xc( const int nspin, const bool domag, const bool domag_z, + const int gga_grad, const double hybrid_alpha, const double hse_omega) { ModuleBase::TITLE("XC_Functional", "v_xc"); + if (use_libxc) { #ifdef __LIBXC @@ -39,6 +44,7 @@ std::tuple XC_Functional::v_xc( nspin, domag, domag_z, + gga_grad, &(scaling_factor_xc), hybrid_alpha, hse_omega); @@ -47,6 +53,12 @@ std::tuple XC_Functional::v_xc( #endif } + // Evaluate the regularized projected local-collinear graph. + if (nspin == 4 && (domag || domag_z) && gga_grad == 2) + { + return ModuleXC::NCGGA_SF_Builtin::v_xc_ncgga_sf_builtin(nrxx, ucell->omega, ucell->tpiba, chr); + } + ModuleBase::timer::start("XC_Functional", "v_xc"); //Exchange-Correlation potential Vxc(r) from n(r) @@ -183,7 +195,7 @@ std::tuple XC_Functional::v_xc( // the dummy variable dum contains gradient correction to stress // which is not used here std::vector dum; - gradcorr(etxc, vtxc, v, chr, chr->rhopw, ucell, dum, false, nspin, domag, domag_z, hybrid_alpha, hse_omega); + gradcorr(etxc, vtxc, v, chr, chr->rhopw, ucell, dum, false, nspin, domag, domag_z, gga_grad, hybrid_alpha, hse_omega); // parallel code : collect vtxc,etxc // mohan add 2008-06-01 diff --git a/source/source_io/module_hs/write_h_terms.cpp b/source/source_io/module_hs/write_h_terms.cpp index f85aa4b97a7..d3785b806df 100644 --- a/source/source_io/module_hs/write_h_terms.cpp +++ b/source/source_io/module_hs/write_h_terms.cpp @@ -350,7 +350,13 @@ void write_h_vxc(WriteHParams& params) #else const double hse_omega = 0.0; #endif - std::tie(etxc, vtxc, v_xc) = XC_Functional::v_xc(nrxx, chg, &ucell, PARAM.inp.nspin, PARAM.globalv.domag, PARAM.globalv.domag_z, hybrid_alpha, hse_omega); + std::tie(etxc, vtxc, v_xc) = XC_Functional::v_xc(nrxx, chg, &ucell, + PARAM.inp.nspin, + PARAM.globalv.domag, + PARAM.globalv.domag_z, + PARAM.inp.gga_grad, + hybrid_alpha, + hse_omega); for (int ispin = 0; ispin < nspin_out; ispin++) { diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h index c5afaaaae1e..1b1f490bafd 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -85,6 +85,7 @@ struct Input_para double nelec_delta = 0.0; ///< change in the number of total electrons double nupdown = 0.0; std::string dft_functional = "default"; ///< input DFT functional. + int gga_grad = 0; ///< Noncollinear GGA: 0 original, 1 local axis, 2 regularized projected LCA. double xc_temperature = 0.0; ///< only relevant if finite temperature functional is used double pseudo_rcut = 15.0; ///< cut-off radius for calculating msh bool pseudo_mesh = false; ///< 0: use msh to normalize radial wave functions; 1: diff --git a/source/source_io/module_parameter/read_inp_estruc.cpp b/source/source_io/module_parameter/read_inp_estruc.cpp index 6f96654199a..49c5fec46b3 100644 --- a/source/source_io/module_parameter/read_inp_estruc.cpp +++ b/source/source_io/module_parameter/read_inp_estruc.cpp @@ -488,6 +488,27 @@ The other way is only available when compiling with LIBXC, and it allows for sup }; this->add_item(item); } + { + Input_Item item("gga_grad"); + item.annotation = "Noncollinear GGA gradient method: 0 original, 1 local axis, 2 regularized projected LCA"; + item.category = "Electronic structure"; + item.type = "Integer"; + item.description = R"(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.)"; + item.default_value = "0"; + read_sync_int(input.gga_grad); + item.check_value = [](const Input_Item& item, const Parameter& para) { + if (para.input.gga_grad < 0 || para.input.gga_grad > 2) + { + ModuleBase::WARNING_QUIT("ReadInput", "gga_grad must be 0, 1, or 2."); + } + }; + this->add_item(item); + } { Input_Item item("smearing_method"); item.annotation = "type of smearing_method: gauss; fd; fixed; mp; mp2; mv"; diff --git a/source/source_io/test/read_input_ptest.cpp b/source/source_io/test/read_input_ptest.cpp index b6a8657fb48..e1dd3e230ce 100644 --- a/source/source_io/test/read_input_ptest.cpp +++ b/source/source_io/test/read_input_ptest.cpp @@ -95,6 +95,8 @@ TEST_F(InputParaTest, ParaRead) EXPECT_DOUBLE_EQ(param.inp.min_dist_coef, 0.2); EXPECT_EQ(param.inp.gint_precision, "double"); EXPECT_EQ(param.inp.dft_functional, "hse"); + EXPECT_EQ(Parameter().inp.gga_grad, 0); + EXPECT_EQ(param.inp.gga_grad, 2); EXPECT_DOUBLE_EQ(param.inp.xc_temperature, 0.0); EXPECT_EQ(param.inp.nspin, 1); EXPECT_DOUBLE_EQ(param.inp.nelec, 0.0); @@ -484,6 +486,37 @@ TEST_F(InputParaTest, ParaRead) EXPECT_DOUBLE_EQ(param.inp.rdmft_power_alpha, 0.656); } +TEST_F(InputParaTest, GgaGradAcceptedRange) +{ + ModuleIO::ReadInput readinput(0); + bool found = false; + for (const std::pair& entry: readinput.get_input_lists()) + { + if (entry.first != "gga_grad") + { + continue; + } + found = true; + ModuleIO::Input_Item item(entry.second); + for (int mode = 0; mode <= 2; ++mode) + { + Parameter param; + item.str_values = {std::to_string(mode)}; + item.read_value(item, param); + EXPECT_EQ(param.inp.gga_grad, mode); + item.check_value(item, param); + } + for (const int mode: {-1, 3}) + { + Parameter param; + item.str_values = {std::to_string(mode)}; + item.read_value(item, param); + EXPECT_EXIT(item.check_value(item, param), testing::ExitedWithCode(1), ""); + } + } + EXPECT_TRUE(found); +} + TEST_F(InputParaTest, TypedTDFieldLists) { ModuleIO::ReadInput readinput(GlobalV::MY_RANK); diff --git a/source/source_io/test/support/INPUT b/source/source_io/test/support/INPUT index df78fb591ee..d72742cd4ee 100644 --- a/source/source_io/test/support/INPUT +++ b/source/source_io/test/support/INPUT @@ -396,3 +396,5 @@ sccut 4 #Maximal step size for lambda in eV/uB #Parameters (23. Time-dependent orbital-free DFT) of_cd 0 #0: no CD potential; 1: add CD potential of_mCD_alpha 1.0 # parameter of modified CD potential + +gga_grad 2 diff --git a/source/source_lcao/force_stress_lcao.cpp b/source/source_lcao/force_stress_lcao.cpp index 24c0c10535c..211d737ed77 100644 --- a/source/source_lcao/force_stress_lcao.cpp +++ b/source/source_lcao/force_stress_lcao.cpp @@ -137,7 +137,7 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, // calculate basic terms in Force, same method with PW base this->calForcePwPart(ucell, parts.fvl_dvl, parts.fewalds, parts.fcc, parts.fscc, pelec->f_en.etxc, pelec->vnew, pelec->vnew_exist, pelec->charge, rhopw, - locpp, sf, cfg.device); + locpp, sf, cfg); } // total stress : ModuleBase::matrix scs @@ -161,7 +161,7 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, // calculate basic terms in Stress, similar method with PW base this->sc_pw.stress_pw_terms(ucell, sparts.sigmadvl, sparts.sigmahar, sparts.sigmaewa, sparts.sigmacc, sparts.sigmaxc, pelec->f_en.etxc, pelec->charge, - rhopw, locpp, sf); + rhopw, locpp, sf, cfg.nspin, cfg.domag, cfg.domag_z, cfg.gga_grad, cfg.gamma_only_pw); } // Calculate operator-based force/stress terms (kinetic, overlap, // nonlocal, rt-TDDFT hybrid gauge, local Pulay term and DeltaSpin). @@ -394,16 +394,17 @@ void Force_Stress_LCAO::calForcePwPart(UnitCell& ucell, ModulePW::PW_Basis* rhopw, const pseudopot_cell_vl& locpp, const Structure_Factor& sf, - const std::string& device) + const FSCalcConfig& cfg) { ModuleBase::TITLE("Force_Stress_LCAO", "calForcePwPart"); #ifdef __CUDA - if (device == "gpu") + if (cfg.device == "gpu") { Forces 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, + cfg.nspin, cfg.domag, cfg.domag_z, cfg.gga_grad); f_pw.cal_force_scc(fscc, rhopw, vnew, vnew_exist, locpp.numeric, ucell); } else @@ -412,7 +413,8 @@ void Force_Stress_LCAO::calForcePwPart(UnitCell& ucell, Forces 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, + cfg.nspin, cfg.domag, cfg.domag_z, cfg.gga_grad); f_pw.cal_force_scc(fscc, rhopw, vnew, vnew_exist, locpp.numeric, ucell); } diff --git a/source/source_lcao/force_stress_lcao.h b/source/source_lcao/force_stress_lcao.h index c075e876364..bb48c66900f 100644 --- a/source/source_lcao/force_stress_lcao.h +++ b/source/source_lcao/force_stress_lcao.h @@ -39,6 +39,10 @@ struct FSCalcConfig bool t_in_h; bool sc_mag_switch; std::string device; + bool domag; + bool domag_z; + int gga_grad; + bool gamma_only_pw; }; // Force/stress component matrices assembled by getForceStress. Grouping them @@ -168,7 +172,7 @@ class Force_Stress_LCAO ModulePW::PW_Basis* rhopw, const pseudopot_cell_vl& locpp, const Structure_Factor& sf, - const std::string& device); + const FSCalcConfig& cfg); static double force_invalid_threshold_ev; }; diff --git a/source/source_lcao/module_operator_lcao/veff_dh.cpp b/source/source_lcao/module_operator_lcao/veff_dh.cpp index 34eea550af1..9e69d7e4970 100644 --- a/source/source_lcao/module_operator_lcao/veff_dh.cpp +++ b/source/source_lcao/module_operator_lcao/veff_dh.cpp @@ -414,7 +414,8 @@ void Veff>::cal_dH(std::arraynrxx; // finite-difference XC: delta V^XC(r) = V^XC[rho0 + drho](r) - V^XC[rho0](r) - elecstate::PotXC_FDM dvxcr_fdm_op(rho_basis, chg, this->ucell); + elecstate::PotXC_FDM dvxcr_fdm_op(chg->nspin, PARAM.globalv.domag, PARAM.globalv.domag_z, + PARAM.inp.gga_grad, rho_basis, chg, this->ucell); std::vector chg_drho(3); for (int d = 0; d < 3; ++d) diff --git a/source/source_lcao/module_rdmft/rdmft_tools.cpp b/source/source_lcao/module_rdmft/rdmft_tools.cpp index 84998880611..45103b62058 100644 --- a/source/source_lcao/module_rdmft/rdmft_tools.cpp +++ b/source/source_lcao/module_rdmft/rdmft_tools.cpp @@ -281,7 +281,7 @@ void Veff_rdmft, double>::contributeHR() ModuleBase::matrix vofk = *vloc_; vofk.zero_out(); ModuleBase::matrix v_matrix_XC(this->nspin, charge_->nrxx); - elecstate::PotXC potXC(rho_basis_, etxc, vtxc, &vofk); + elecstate::PotXC potXC(PARAM.globalv.domag, PARAM.globalv.domag_z, PARAM.inp.gga_grad, rho_basis_, etxc, vtxc, &vofk); potXC.cal_v_eff(charge_, ucell, v_matrix_XC); // if need meta-GGA, go to study veff_lcao.cpp and modify the code @@ -363,7 +363,7 @@ void Veff_rdmft::contributeHR() ModuleBase::matrix vofk = *vloc_; vofk.zero_out(); ModuleBase::matrix v_matrix_XC(this->nspin, charge_->nrxx); - elecstate::PotXC potXC(rho_basis_, etxc, vtxc, &vofk); + elecstate::PotXC potXC(PARAM.globalv.domag, PARAM.globalv.domag_z, PARAM.inp.gga_grad, rho_basis_, etxc, vtxc, &vofk); potXC.cal_v_eff(charge_, ucell, v_matrix_XC); for(int is=0; isnspin; ++is) diff --git a/source/source_pw/module_ofdft/of_stress_pw.cpp b/source/source_pw/module_ofdft/of_stress_pw.cpp index 6b988baf786..20189bd6f62 100644 --- a/source/source_pw/module_ofdft/of_stress_pw.cpp +++ b/source/source_pw/module_ofdft/of_stress_pw.cpp @@ -7,7 +7,12 @@ // Since the kinetic stress of OFDFT is calculated by kinetic functionals in esolver_of.cpp, here we regard it as an // input variable. -void OF_Stress_PW::cal_stress(ModuleBase::matrix& sigmatot, +void OF_Stress_PW::cal_stress(const int nspin, + const bool domag, + const bool domag_z, + const int gga_grad, + const bool gamma_only_pw, + ModuleBase::matrix& sigmatot, ModuleBase::matrix& kinetic_stress, UnitCell& ucell, const vdw::VdwResult* vdw_result, @@ -74,13 +79,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, + nspin, domag, domag_z, 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, + nspin, domag, domag_z, gga_grad, gamma_only_pw); // vdW term prepared before SCF for this ionic configuration. if (vdw_result != nullptr) diff --git a/source/source_pw/module_ofdft/of_stress_pw.h b/source/source_pw/module_ofdft/of_stress_pw.h index 17fc50f4775..edb62c9e04c 100644 --- a/source/source_pw/module_ofdft/of_stress_pw.h +++ b/source/source_pw/module_ofdft/of_stress_pw.h @@ -17,7 +17,12 @@ class OF_Stress_PW : public Stress_Func : pelec(pelec_in), rhopw(rhopw_in){}; // calculate the stress in OFDFT - void cal_stress(ModuleBase::matrix& sigmatot, + void cal_stress(const int nspin, + const bool domag, + const bool domag_z, + const int gga_grad, + const bool gamma_only_pw, + ModuleBase::matrix& sigmatot, ModuleBase::matrix& kinetic_stress, UnitCell& ucell, const vdw::VdwResult* vdw_result, diff --git a/source/source_pw/module_pwdft/force_pw.cpp b/source/source_pw/module_pwdft/force_pw.cpp index 85d27eb0013..d37e2379954 100644 --- a/source/source_pw/module_pwdft/force_pw.cpp +++ b/source/source_pw/module_pwdft/force_pw.cpp @@ -27,7 +27,12 @@ template -void Forces::cal_force(UnitCell& ucell, +void Forces::cal_force(const int nspin, + const bool domag, + const bool domag_z, + const int gga_grad, + const bool use_onsite_projection, + UnitCell& ucell, ModuleBase::matrix& force, const vdw::VdwResult* vdw_result, const elecstate::ElecState& elec, @@ -77,14 +82,15 @@ void Forces::cal_force(UnitCell& ucell, // DFT+U and DeltaSpin // here maybe a bug when OFDFT calls +U, mohan add 20251107 - if(PARAM.inp.dft_plus_u || PARAM.inp.sc_mag_switch) + if(use_onsite_projection) { this->cal_force_onsite(forceonsite, wg, wfc_basis, ucell, *p_dftu, psi_in); } } // 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, + nspin, domag, domag_z, gga_grad); // force due to core charge this->cal_force_scc(forcescc, rho_basis, elec.vnew, elec.vnew_exist, locpp->numeric, ucell); @@ -135,7 +141,7 @@ void Forces::cal_force(UnitCell& ucell, if (PARAM.inp.imp_sol) { forcesol.create(this->nat, 3); - solvent.cal_force_sol(ucell, rho_basis, locpp->vloc, PARAM.inp.nspin, forcesol); + solvent.cal_force_sol(ucell, rho_basis, locpp->vloc, nspin, forcesol); if (PARAM.inp.test_force) { ModuleIO::print_force(GlobalV::ofs_running, ucell, "IMP_SOL FORCE (Ry/Bohr)", forcesol); @@ -174,7 +180,7 @@ void Forces::cal_force(UnitCell& ucell, force(iat, ipol) = force(iat, ipol) + forcesol(iat, ipol); } - if(PARAM.inp.dft_plus_u || PARAM.inp.sc_mag_switch) + if(use_onsite_projection) { force(iat, ipol) += forceonsite(iat, ipol); } @@ -261,7 +267,7 @@ void Forces::cal_force(UnitCell& ucell, forcesol, false); } - if (PARAM.inp.dft_plus_u || PARAM.inp.sc_mag_switch) + if (use_onsite_projection) { ModuleIO::print_force(GlobalV::ofs_running, ucell, diff --git a/source/source_pw/module_pwdft/force_pw.h b/source/source_pw/module_pwdft/force_pw.h index 3a2f3c5fa9f..6aeaff6069e 100644 --- a/source/source_pw/module_pwdft/force_pw.h +++ b/source/source_pw/module_pwdft/force_pw.h @@ -45,7 +45,12 @@ class Forces Forces(const int nat_in) : nat(nat_in){}; ~Forces(){}; - void cal_force(UnitCell& ucell, + void cal_force(const int nspin, + const bool domag, + const bool domag_z, + const int gga_grad, + const bool use_onsite_projection, + UnitCell& ucell, ModuleBase::matrix& force, const vdw::VdwResult* vdw_result, const elecstate::ElecState& elec, @@ -77,7 +82,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} [ diff --git a/source/source_pw/module_pwdft/force_pw_cc.cpp b/source/source_pw/module_pwdft/force_pw_cc.cpp index 56575ea0925..a0d0b944e4c 100644 --- a/source/source_pw/module_pwdft/force_pw_cc.cpp +++ b/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,7 +31,11 @@ void Forces::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) { ModuleBase::TITLE("Forces", "cal_force_cc"); // recalculate the exchange-correlation potential. @@ -53,7 +56,7 @@ void Forces::cal_force_cc(ModuleBase::matrix& forcecc, return; } - ModuleBase::matrix v(PARAM.inp.nspin, rho_basis->nrxx); + ModuleBase::matrix v(nspin, rho_basis->nrxx); const double hybrid_alpha = XC_Functional::get_hybrid_alpha(); #ifdef __EXX @@ -66,7 +69,7 @@ void Forces::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, - PARAM.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); @@ -77,11 +80,12 @@ void Forces::cal_force_cc(ModuleBase::matrix& forcecc, } else { - unitcell::cal_ux(ucell_in, PARAM.inp.nspin); + unitcell::cal_ux(ucell_in, nspin); const auto etxc_vtxc_v = XC_Functional::v_xc(rho_basis->nrxx, chr, &ucell_in, - PARAM.inp.nspin, - PARAM.globalv.domag, - PARAM.globalv.domag_z, + nspin, + domag, + domag_z, + gga_grad, hybrid_alpha, hse_omega); @@ -92,7 +96,7 @@ void Forces::cal_force_cc(ModuleBase::matrix& forcecc, const ModuleBase::matrix vxc = v; std::complex* psiv = new std::complex[rho_basis->nmaxgr]; - if (PARAM.inp.nspin == 1 || PARAM.inp.nspin == 4) + if (nspin == 1 || nspin == 4) { #ifdef _OPENMP #pragma omp parallel for schedule(static, 1024) @@ -368,4 +372,4 @@ void Forces::deriv_drhoc template class Forces; #if ((defined __CUDA) || (defined __ROCM)) template class Forces; -#endif \ No newline at end of file +#endif diff --git a/source/source_pw/module_pwdft/stress_cc.cpp b/source/source_pw/module_pwdft/stress_cc.cpp index ffe5dc85e09..a0b6cb3acb3 100644 --- a/source/source_pw/module_pwdft/stress_cc.cpp +++ b/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,18 +14,23 @@ template void Stress_Func::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) { ModuleBase::TITLE("Stress","stress_cc"); ModuleBase::timer::start("Stress","stress_cc"); FPTYPE fact=1.0; - if(is_pw&&PARAM.globalv.gamma_only_pw) + if(is_pw&&gamma_only_pw) { fact = 2.0; //is_pw:PW basis, gamma_only need to FPTYPE. } @@ -62,7 +66,7 @@ void Stress_Func::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, - PARAM.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); @@ -73,11 +77,12 @@ void Stress_Func::stress_cc(ModuleBase::matrix& sigma, } else { - unitcell::cal_ux(ucell, PARAM.inp.nspin); + unitcell::cal_ux(ucell, nspin); const auto etxc_vtxc_v = XC_Functional::v_xc(rho_basis->nrxx, chr, &ucell, - PARAM.inp.nspin, - PARAM.globalv.domag, - PARAM.globalv.domag_z, + nspin, + domag, + domag_z, + gga_grad, hybrid_alpha, hse_omega); // etxc = std::get<0>(etxc_vtxc_v); // may delete? @@ -87,7 +92,7 @@ void Stress_Func::stress_cc(ModuleBase::matrix& sigma, std::complex* psic = new std::complex[rho_basis->nmaxgr]; - if(PARAM.inp.nspin==1||PARAM.inp.nspin==4) + if(nspin==1||nspin==4) { #ifdef _OPENMP #pragma omp parallel for schedule(static, 1024) diff --git a/source/source_pw/module_pwdft/stress_func.h b/source/source_pw/module_pwdft/stress_func.h index d6020e4b6b5..4eb758999b8 100644 --- a/source/source_pw/module_pwdft/stress_func.h +++ b/source/source_pw/module_pwdft/stress_func.h @@ -119,7 +119,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, @@ -136,7 +141,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, @@ -164,7 +173,12 @@ class Stress_Func const Charge* const chr, ModulePW::PW_Basis* rhopw, const pseudopot_cell_vl& locpp, - const Structure_Factor& sf); + const Structure_Factor& sf, + const int nspin, + const bool domag, + const bool domag_z, + const int gga_grad, + const bool gamma_only_pw); // 7) the stress from the non-local pseudopotentials /** diff --git a/source/source_pw/module_pwdft/stress_gga.cpp b/source/source_pw/module_pwdft/stress_gga.cpp index bb175514803..1716bd093db 100644 --- a/source/source_pw/module_pwdft/stress_gga.cpp +++ b/source/source_pw/module_pwdft/stress_gga.cpp @@ -1,14 +1,17 @@ #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 void Stress_Func::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) { ModuleBase::TITLE("Stress","stress_gga"); ModuleBase::timer::start("Stress","stress_gga"); @@ -31,7 +34,7 @@ void Stress_Func::stress_gga(const UnitCell& ucell, XC_Functional::gradcorr( dum1, dum2, dum3, chr, rho_basis, &ucell, stress_gga, is_stress, - PARAM.inp.nspin, PARAM.globalv.domag, PARAM.globalv.domag_z, + nspin, domag, domag_z, gga_grad, hybrid_alpha, hse_omega); for(int l = 0;l< 3;l++) diff --git a/source/source_pw/module_pwdft/stress_loc.cpp b/source/source_pw/module_pwdft/stress_loc.cpp index 1f4c0a8b2ca..94f1746dbe0 100644 --- a/source/source_pw/module_pwdft/stress_loc.cpp +++ b/source/source_pw/module_pwdft/stress_loc.cpp @@ -314,7 +314,12 @@ void Stress_Func::stress_pw_terms(UnitCell& ucell, const Charge* const chr, ModulePW::PW_Basis* rhopw, const pseudopot_cell_vl& locpp, - const Structure_Factor& sf) + const Structure_Factor& sf, + const int nspin, + const bool domag, + const bool domag_z, + const int gga_grad, + const bool gamma_only_pw) { ModuleBase::TITLE("Stress", "stress_pw_terms"); @@ -328,7 +333,8 @@ void Stress_Func::stress_pw_terms(UnitCell& ucell, this->stress_ewa(ucell, sigmaewa, rhopw, 0); // remain problem // stress due to core correlation. - this->stress_cc(sigmacc, rhopw, ucell, &sf, 0, locpp.numeric, chr); + this->stress_cc(sigmacc, rhopw, ucell, &sf, 0, locpp.numeric, chr, + nspin, domag, domag_z, gga_grad, gamma_only_pw); // stress due to self-consistent charge. for (int i = 0; i < 3; i++) @@ -336,7 +342,7 @@ void Stress_Func::stress_pw_terms(UnitCell& ucell, sigmaxc(i, i) = -etxc / ucell.omega; } // Exchange-correlation for PBE - this->stress_gga(ucell, sigmaxc, rhopw, chr); + this->stress_gga(ucell, sigmaxc, rhopw, chr, nspin, domag, domag_z, gga_grad); return; } diff --git a/source/source_pw/module_pwdft/stress_pw.cpp b/source/source_pw/module_pwdft/stress_pw.cpp index bcf6aff8db2..b2cea68404f 100644 --- a/source/source_pw/module_pwdft/stress_pw.cpp +++ b/source/source_pw/module_pwdft/stress_pw.cpp @@ -9,7 +9,12 @@ #include "source_hamilt/module_xc/general_exx_info.h" // for General_Exx_Info type template -void Stress_PW::cal_stress(ModuleBase::matrix& sigmatot, +void Stress_PW::cal_stress(const int nspin, + const bool domag, + const bool domag_z, + const int gga_grad, + const bool gamma_only_pw, + ModuleBase::matrix& sigmatot, UnitCell& ucell, const vdw::VdwResult* vdw_result, Plus_U_Base& dftu, @@ -91,7 +96,8 @@ void Stress_PW::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, + nspin, domag, domag_z, gga_grad); if (XC_Functional::get_ked_flag()) { this->stress_mgga(ucell, @@ -108,7 +114,8 @@ void Stress_PW::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, + nspin, domag, domag_z, gga_grad, 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); diff --git a/source/source_pw/module_pwdft/stress_pw.h b/source/source_pw/module_pwdft/stress_pw.h index f203a642e71..19f4bcb87ea 100644 --- a/source/source_pw/module_pwdft/stress_pw.h +++ b/source/source_pw/module_pwdft/stress_pw.h @@ -21,8 +21,13 @@ class Stress_PW : public Stress_Func Stress_PW(const elecstate::ElecState* pelec_in) : pelec(pelec_in){}; // calculate the stress in PW basis - void cal_stress(ModuleBase::matrix& smearing_sigmatot, - UnitCell& ucell, + void cal_stress(const int nspin, + const bool domag, + const bool domag_z, + const int gga_grad, + const bool gamma_only_pw, + ModuleBase::matrix& smearing_sigmatot, + UnitCell& ucell, const vdw::VdwResult* vdw_result, Plus_U_Base& dftu, const pseudopot_cell_vl& locpp, diff --git a/source/source_pw/module_stodft/sto_forces.cpp b/source/source_pw/module_stodft/sto_forces.cpp index e092b8f9327..3f455dec4a5 100644 --- a/source/source_pw/module_stodft/sto_forces.cpp +++ b/source/source_pw/module_stodft/sto_forces.cpp @@ -16,7 +16,11 @@ #include "source_hamilt/module_xc/xc_functional.h" template -void Sto_Forces::cal_stoforce(ModuleBase::matrix& force, +void Sto_Forces::cal_stoforce(const int nspin, + const bool domag, + const bool domag_z, + const int gga_grad, + ModuleBase::matrix& force, const elecstate::ElecState& elec, ModulePW::PW_Basis* rho_basis, ModuleSymmetry::Symmetry* p_symm, @@ -44,7 +48,8 @@ void Sto_Forces::cal_stoforce(ModuleBase::matrix& force, this->cal_force_loc(ucell, forcelc, rho_basis, locpp.vloc, chr); this->cal_force_ew(ucell,forceion, rho_basis, p_sf); this->cal_sto_force_nl(forcenl, wg, pkv, wfc_basis, p_sf, nlpp, ucell, psi, stowf); - this->cal_force_cc(forcecc, rho_basis, chr, locpp.numeric, ucell); + this->cal_force_cc(forcecc, rho_basis, chr, locpp.numeric, ucell, + nspin, domag, domag_z, gga_grad); this->cal_force_scc(forcescc, rho_basis, elec.vnew, elec.vnew_exist, locpp.numeric, ucell); // impose total force = 0 diff --git a/source/source_pw/module_stodft/sto_forces.h b/source/source_pw/module_stodft/sto_forces.h index d8416b916bc..4b631f72a5f 100644 --- a/source/source_pw/module_stodft/sto_forces.h +++ b/source/source_pw/module_stodft/sto_forces.h @@ -21,7 +21,11 @@ class Sto_Forces : public Forces Sto_Forces(const int nat_in) : Forces(nat_in){}; ~Sto_Forces(){}; - void cal_stoforce(ModuleBase::matrix& force, + void cal_stoforce(const int nspin, + const bool domag, + const bool domag_z, + const int gga_grad, + ModuleBase::matrix& force, const elecstate::ElecState& elec, ModulePW::PW_Basis* rho_basis, ModuleSymmetry::Symmetry* p_symm, diff --git a/source/source_pw/module_stodft/sto_stress_pw.cpp b/source/source_pw/module_stodft/sto_stress_pw.cpp index b85ab58b503..68a42d67006 100644 --- a/source/source_pw/module_stodft/sto_stress_pw.cpp +++ b/source/source_pw/module_stodft/sto_stress_pw.cpp @@ -9,7 +9,12 @@ #include "source_io/module_parameter/parameter.h" template -void Sto_Stress_PW::cal_stress(ModuleBase::matrix& sigmatot, +void Sto_Stress_PW::cal_stress(const int nspin, + const bool domag, + const bool domag_z, + const int gga_grad, + const bool gamma_only_pw, + ModuleBase::matrix& sigmatot, const elecstate::ElecState& elec, ModulePW::PW_Basis* rho_basis, ModuleSymmetry::Symmetry* p_symm, @@ -50,13 +55,15 @@ void Sto_Stress_PW::cal_stress(ModuleBase::matrix& sigmatot, { sigmaxc(i, i) = -(elec.f_en.etxc - elec.f_en.vtxc) / this->ucell->omega; } - this->stress_gga(ucell_in, sigmaxc, rho_basis, chr); + this->stress_gga(ucell_in, sigmaxc, rho_basis, chr, + nspin, domag, domag_z, gga_grad); // local contribution this->stress_loc(ucell_in, sigmaloc, rho_basis, locpp->vloc, p_sf, true, chr); // nlcc - this->stress_cc(sigmaxcc, rho_basis, ucell_in, p_sf, true, locpp->numeric, chr); + this->stress_cc(sigmaxcc, rho_basis, ucell_in, p_sf, true, locpp->numeric, chr, + nspin, domag, domag_z, gga_grad, gamma_only_pw); // nonlocal this->sto_stress_nl(sigmanl, wg, p_sf, p_symm, p_kv, wfc_basis, *nlpp, ucell_in, psi_in, stowf); diff --git a/source/source_pw/module_stodft/sto_stress_pw.h b/source/source_pw/module_stodft/sto_stress_pw.h index c357cb1fb6f..cb5f41d2926 100644 --- a/source/source_pw/module_stodft/sto_stress_pw.h +++ b/source/source_pw/module_stodft/sto_stress_pw.h @@ -17,7 +17,12 @@ class Sto_Stress_PW : public Stress_Func ~Sto_Stress_PW(){}; // calculate the stress in PW basis - void cal_stress(ModuleBase::matrix& sigmatot, + void cal_stress(const int nspin, + const bool domag, + const bool domag_z, + const int gga_grad, + const bool gamma_only_pw, + ModuleBase::matrix& sigmatot, const elecstate::ElecState& elec, ModulePW::PW_Basis* rho_basis, ModuleSymmetry::Symmetry* p_symm,