diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index f04f903448..ffbd65bc85 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -93,6 +93,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) @@ -1268,6 +1269,17 @@ - 4: For the case of noncollinear polarized, nspin will be automatically set to 4 without being specified by the user. - **Default**: 1 +### gga_grad + +- **Type**: Integer +- **Description**: Method used to evaluate the density gradient entering GGA exchange-correlation terms in noncollinear-spin (nspin=4) calculations. With rho_up/dn = (rho +/- |m|)/2 and m_hat = m/|m|: + - 0: original algorithm; the default. When the initial magnetic moments in STRU are all collinear, their common direction u is used as a global quantization axis and the up/down densities are defined w.r.t. this axis via the sign of m . u. This reproduces the nspin=2 collinear result exactly for collinear configurations, but makes the magnetic potential energy surface discontinuous when the initial moments are tilted slightly away from collinearity. + - 1: collinear approximation without the global direction: up/down are always defined w.r.t. the local |m|, so the magnetic potential energy surface stays continuous, at the price of giving up exact consistency with nspin=2 for collinear configurations. (For LIBXC functionals the global axis is never used, so 0 and 1 are equivalent.) + - 2: projected method. The gradients use the full chain rule, grad(rho_up/dn) = (grad(rho) +/- m_hat . grad(m))/2, but the divergence of h = df/d(grad rho) in the potential is projected onto m_hat: v_mu -= m_hat_mu * div((h_up - h_dn)/2), dropping the (h_up - h_dn) . grad(m_hat_mu) cross terms. + - 3: Scalmani-Frisch transformation (G. Scalmani and M. J. Frisch, J. Chem. Theory Comput. 8, 2193 (2012)). Same gradients as 2, but the full divergence is kept: v_mu -= div((h_up - h_dn)/2 * m_hat_mu), retaining all cross terms; the most accurate. + This parameter only takes effect for nspin=4 with GGA functionals (and magnetic calculation). +- **Default**: 0 + ### smearing_method - **Type**: String diff --git a/docs/parameters.yaml b/docs/parameters.yaml index 005fad3076..c5cbbe663c 100644 --- a/docs/parameters.yaml +++ b/docs/parameters.yaml @@ -665,6 +665,19 @@ parameters: default_value: "1" unit: "" availability: "" + - name: gga_grad + category: Electronic structure + type: Integer + description: | + Method used to evaluate the density gradient entering GGA exchange-correlation terms in noncollinear-spin (nspin=4) calculations. With rho_up/dn = (rho +/- |m|)/2 and m_hat = m/|m|: + * 0: original algorithm; the default. When the initial magnetic moments in STRU are all collinear, their common direction u is used as a global quantization axis and the up/down densities are defined w.r.t. this axis via the sign of m . u. This reproduces the nspin=2 collinear result exactly for collinear configurations, but makes the magnetic potential energy surface discontinuous when the initial moments are tilted slightly away from collinearity. + * 1: collinear approximation without the global direction: up/down are always defined w.r.t. the local |m|, so the magnetic potential energy surface stays continuous, at the price of giving up exact consistency with nspin=2 for collinear configurations. (For LIBXC functionals the global axis is never used, so 0 and 1 are equivalent.) + * 2: projected method. The gradients use the full chain rule, grad(rho_up/dn) = (grad(rho) +/- m_hat . grad(m))/2, but the divergence of h = df/d(grad rho) in the potential is projected onto m_hat: v_mu -= m_hat_mu * div((h_up - h_dn)/2), dropping the (h_up - h_dn) . grad(m_hat_mu) cross terms. + * 3: Scalmani-Frisch transformation (G. Scalmani and M. J. Frisch, J. Chem. Theory Comput. 8, 2193 (2012)). Same gradients as 2, but the full divergence is kept: v_mu -= div((h_up - h_dn)/2 * m_hat_mu), retaining all cross terms; the most accurate. + This parameter only takes effect for nspin=4 with GGA functionals (and magnetic calculation). + default_value: "0" + unit: "" + availability: "" - name: smearing_method category: Electronic structure type: String diff --git a/source/source_estate/module_pot/pot_xc.cpp b/source/source_estate/module_pot/pot_xc.cpp index 19815f843a..de26641a38 100644 --- a/source/source_estate/module_pot/pot_xc.cpp +++ b/source/source_estate/module_pot/pot_xc.cpp @@ -16,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 = PARAM.inp.nspin; //---------------------------------------------------------- // calculate the exchange-correlation potential @@ -32,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); @@ -51,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, + nspin, PARAM.globalv.domag, PARAM.globalv.domag_z, + PARAM.inp.gga_grad, hybrid_alpha, hse_omega); *(this->etxc_) = std::get<0>(etxc_vtxc_v); diff --git a/source/source_estate/module_pot/pot_xc_fdm.cpp b/source/source_estate/module_pot/pot_xc_fdm.cpp index 03349fe4e7..1c69722c98 100644 --- a/source/source_estate/module_pot/pot_xc_fdm.cpp +++ b/source/source_estate/module_pot/pot_xc_fdm.cpp @@ -31,6 +31,7 @@ PotXC_FDM::PotXC_FDM( PARAM.inp.nspin, PARAM.globalv.domag, PARAM.globalv.domag_z, + PARAM.inp.gga_grad, hybrid_alpha, hse_omega); this->v_xc_0 = std::get<2>(etxc_vtxc_v_0); @@ -69,6 +70,7 @@ void PotXC_FDM::cal_v_eff( PARAM.inp.nspin, PARAM.globalv.domag, PARAM.globalv.domag_z, + PARAM.inp.gga_grad, hybrid_alpha, hse_omega); const ModuleBase::matrix &v_xc_01 = std::get<2>(etxc_vtxc_v_01); diff --git a/source/source_hamilt/module_xc/CMakeLists.txt b/source/source_hamilt/module_xc/CMakeLists.txt index 5f310d09d5..abc41012ac 100644 --- a/source/source_hamilt/module_xc/CMakeLists.txt +++ b/source/source_hamilt/module_xc/CMakeLists.txt @@ -2,6 +2,7 @@ add_library( xc_ OBJECT xc_functional.cpp + xc_functional_ncgga_sf.cpp xc_pot.cpp xc_grad.cpp xc_lda_wrap.cpp diff --git a/source/source_hamilt/module_xc/libxc_abacus.h b/source/source_hamilt/module_xc/libxc_abacus.h index bc584b7518..dccad5374a 100644 --- a/source/source_hamilt/module_xc/libxc_abacus.h +++ b/source/source_hamilt/module_xc/libxc_abacus.h @@ -24,23 +24,23 @@ namespace XC_Functional_Libxc //------------------- // sets functional type, which allows combination of LIBXC keyword connected by "+" - // for example, "XC_LDA_X+XC_LDA_C_PZ" + // for example: "XC_LDA_X+XC_LDA_C_PZ" extern std::pair> 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,6 +65,7 @@ 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); @@ -128,6 +129,9 @@ namespace XC_Functional_Libxc std::vector exc); // converting vtxc and v from vrho and vsigma (libxc=>abacus) + // use_sf: for nspin=4 magnetic GGA, apply the Scalmani-Frisch + // gradient correction instead of the collinear one; gga_grad then + // selects the projected (2) or full (3) divergence of h extern std::pair convert_vtxc_v( const xc_func_type &func, const int nspin, @@ -138,7 +142,9 @@ namespace XC_Functional_Libxc const std::vector &vrho, const std::vector &vsigma, const double tpiba, - const Charge* const chr); + const Charge* const chr, + const bool use_sf, + const int gga_grad); // dh for gga v extern std::vector> cal_dh( @@ -151,10 +157,42 @@ namespace XC_Functional_Libxc const Charge* const chr); // convert v for NSPIN=4 + // has_mag: whether the calculation has (noncollinear) magnetization, + // i.e. domag || domag_z extern ModuleBase::matrix convert_v_nspin4( const std::size_t nrxx, const Charge* const chr, const std::vector &amag, + const ModuleBase::matrix &v, + const bool has_mag); + + extern std::vector compute_mag_part_nspin4( + const std::size_t nrxx, + const Charge* const chr); + + extern std::vector>> cal_gdr_sf( + const int nspin, + const std::size_t nrxx, + const std::vector &rho, + const std::vector &mag_part, + const double tpiba, + const Charge* const chr); + + extern std::vector> cal_dh_sf( + const int nspin, + const std::size_t nrxx, + const std::vector &sgn, + const std::vector>> &gdr, + const std::vector &vsigma, + const std::vector &mag_part, + const int gga_grad, + const double tpiba, + const Charge* const chr); + + extern ModuleBase::matrix convert_v_nspin4_sf( + const std::size_t nrxx, + const Charge* const chr, + const std::vector &mag_part, const ModuleBase::matrix &v); diff --git a/source/source_hamilt/module_xc/libxc_pot.cpp b/source/source_hamilt/module_xc/libxc_pot.cpp index 6f605bbd48..41002b1974 100644 --- a/source/source_hamilt/module_xc/libxc_pot.cpp +++ b/source/source_hamilt/module_xc/libxc_pot.cpp @@ -4,7 +4,6 @@ #include "libxc_abacus.h" #include "source_estate/module_charge/charge.h" #include "source_base/global_variable.h" -#include "source_io/module_parameter/parameter.h" #include "source_base/parallel_reduce.h" #include "source_base/timer.h" #include "source_base/tool_title.h" @@ -25,6 +24,7 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( / 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) @@ -36,6 +36,12 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( / (nspin_in == 1 || ( nspin_in ==4 && !domag && !domag_z)) ? 1 : 2; + // For nspin=4 with noncollinear magnetism, gga_grad=2/3 selects the + // Scalmani-Frisch (SF) gradient decomposition; gga_grad=0/1 keeps the + // original collinear algorithm. + const bool has_mag = domag || domag_z; + const bool use_sf = (nspin_in == 4) && has_mag && (gga_grad == 2 || gga_grad == 3); + //---------------------------------------------------------- // xc_func_type is defined in Libxc package // to understand the usage of xc_func_type, @@ -64,8 +70,13 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( / }(); // converting rho + // For nspin=4, the charge density has 4 components: + // rho[0] = total charge, rho[1..3] = magnetization (mx, my, mz) + // libxc works with spin-up/spin-down densities: + // rho_up = 0.5*(rho[0] + |m|), rho_dn = 0.5*(rho[0] - |m|) std::vector rho; std::vector amag; + std::vector mag_part; if(1==nspin || 2==nspin_in) { rho = XC_Functional_Libxc::convert_rho(nspin, nrxx, chr); @@ -75,13 +86,32 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( / 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)); + // gga_grad=2/3: compute magnetization unit vector m_hat = m/|m| + // needed for the Scalmani-Frisch (SF) gradient decomposition: + // grad(rho_up) = 0.5 grad(rho[0]) + 0.5 m_hat_mu * grad(m_mu) + // grad(rho_dn) = 0.5 grad(rho[0]) - 0.5 m_hat_mu * grad(m_mu) + // gga_grad=0/1 (else): collinear approximation, uses grad(|m|) only + if(use_sf) + { + mag_part = XC_Functional_Libxc::compute_mag_part_nspin4(nrxx, chr); + } } std::vector>> gdr; std::vector sigma; if(is_gga) { - gdr = XC_Functional_Libxc::cal_gdr(nspin, nrxx, rho, tpiba, chr); + // gga_grad=2/3: use SF method to compute spin-up/spin-down gradients + // via the chain-rule decomposition. This is more accurate than the + // collinear approximation (gga_grad=0/1) for noncollinear magnetism. + if(use_sf) + { + gdr = XC_Functional_Libxc::cal_gdr_sf(nspin, nrxx, rho, mag_part, tpiba, chr); + } + else + { + gdr = XC_Functional_Libxc::cal_gdr(nspin, nrxx, rho, tpiba, chr); + } sigma = XC_Functional_Libxc::convert_sigma(gdr); } @@ -172,14 +202,24 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( / func, nspin, nrxx, sgn, rho, gdr, vrho, vsigma, - tpiba, chr); + tpiba, chr, use_sf, gga_grad); vtxc += std::get<0>(vtxc_v) * factor; v += std::get<1>(vtxc_v) * factor; } // end for( xc_func_type &func : funcs ) if(4==nspin_in) { - v = XC_Functional_Libxc::convert_v_nspin4(nrxx, chr, amag, v); + // gga_grad=2/3: convert libxc spin-up/spin-down potential back to + // nspin=4 representation via SF method (v_total, v_mag_hat) + // gga_grad=0/1: standard conversion using |m| decomposition only + if(use_sf) + { + v = XC_Functional_Libxc::convert_v_nspin4_sf(nrxx, chr, mag_part, v); + } + else + { + v = XC_Functional_Libxc::convert_v_nspin4(nrxx, chr, amag, v, has_mag); + } } //------------------------------------------------- diff --git a/source/source_hamilt/module_xc/libxc_tools.cpp b/source/source_hamilt/module_xc/libxc_tools.cpp index 916732fa8d..9ce2cf57dd 100644 --- a/source/source_hamilt/module_xc/libxc_tools.cpp +++ b/source/source_hamilt/module_xc/libxc_tools.cpp @@ -3,7 +3,6 @@ #include "libxc_abacus.h" #include "xc_functional.h" #include "source_estate/module_charge/charge.h" -#include "source_io/module_parameter/parameter.h" // converting rho (abacus=>libxc) std::vector XC_Functional_Libxc::convert_rho( @@ -32,7 +31,7 @@ XC_Functional_Libxc::convert_rho_amag_nspin4( const std::size_t nrxx, const Charge* const chr) { - assert(PARAM.inp.nspin==4); + 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 @@ -51,6 +50,28 @@ XC_Functional_Libxc::convert_rho_amag_nspin4( return std::make_tuple(std::move(rho), std::move(amag)); } +std::vector XC_Functional_Libxc::compute_mag_part_nspin4( + const std::size_t nrxx, + const Charge* const chr) +{ + std::vector mag_part(3 * nrxx, 0.0); + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + double mx = chr->rho[1][ir], my = chr->rho[2][ir], mz = chr->rho[3][ir]; + double amag = std::sqrt(mx * mx + my * my + mz * mz); + if (amag > 1e-12) + { + mag_part[ir] = mx / amag; + mag_part[ir + nrxx] = my / amag; + mag_part[ir + 2 * nrxx] = mz / amag; + } + } + return mag_part; +} + // calculating grho std::vector>> XC_Functional_Libxc::cal_gdr( @@ -191,7 +212,9 @@ std::pair XC_Functional_Libxc::convert_vtxc_v( const std::vector &vrho, const std::vector &vsigma, const double tpiba, - const Charge* const chr) + const Charge* const chr, + const bool use_sf, + const int gga_grad) { // assert(nrxx>0); // will cause error double vtxc = 0.0; @@ -213,22 +236,68 @@ std::pair XC_Functional_Libxc::convert_vtxc_v( if(func.info->family == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA) { - const std::vector> dh = XC_Functional_Libxc::cal_dh(nspin, nrxx, sgn, gdr, vsigma, tpiba, chr); - - double rvtxc = 0.0; - #ifdef _OPENMP - #pragma omp parallel for collapse(2) reduction(+:rvtxc) schedule(static, 256) - #endif - for( int is=0; is mag_part_tmp = XC_Functional_Libxc::compute_mag_part_nspin4(nrxx, chr); + const std::vector> dh = XC_Functional_Libxc::cal_dh_sf(nspin, nrxx, sgn, gdr, vsigma, mag_part_tmp, gga_grad, tpiba, chr); + + // vtxc contribution: sum dh[is] * rho[is] over the 4 nspin=4 channels + constexpr int nspin4 = 4; + double rvtxc = 0.0; + #ifdef _OPENMP + #pragma omp parallel for collapse(2) reduction(+:rvtxc) schedule(static, 256) + #endif + for (int is = 0; is < nspin4; ++is) + { + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + rvtxc += dh[is][ir] * chr->rho[is][ir]; + } + } + vtxc -= rvtxc; + + // v is kept in the spin-up/down representation here, so map the + // 4-component dh back with the inverse of the SF decomposition: + // v_up -= dh0 + sum_mu dh_mu * m_hat_mu + // v_dn -= dh0 - sum_mu dh_mu * m_hat_mu + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for (std::size_t ir = 0; ir < nrxx; ++ir) { - rvtxc += dh[is][ir] * rho[ir*nspin+is]; - v(is,ir) -= dh[is][ir]; + double dh_mag = 0.0; + for (int mu = 1; mu < nspin4; ++mu) + { + dh_mag += dh[mu][ir] * mag_part_tmp[ir + (mu - 1) * nrxx]; + } + v(0, ir) -= dh[0][ir] + dh_mag; + v(1, ir) -= dh[0][ir] - dh_mag; } } + else + { + const std::vector> dh = XC_Functional_Libxc::cal_dh(nspin, nrxx, sgn, gdr, vsigma, tpiba, chr); - vtxc -= rvtxc; + double rvtxc = 0.0; + #ifdef _OPENMP + #pragma omp parallel for collapse(2) reduction(+:rvtxc) schedule(static, 256) + #endif + for( int is=0; isfamily == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA)) return std::make_pair(vtxc, std::move(v)); @@ -290,24 +359,25 @@ 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 ModuleBase::matrix &v, + const bool has_mag) { //assert(nrxx>0); - assert(PARAM.inp.nspin==4); + constexpr int nspin4 = 4; constexpr double vanishing_charge = 1.0e-10; - ModuleBase::matrix v_nspin4(PARAM.inp.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]; } @@ -317,4 +387,203 @@ ModuleBase::matrix XC_Functional_Libxc::convert_v_nspin4( return v_nspin4; } +std::vector>> XC_Functional_Libxc::cal_gdr_sf( + const int nspin, + const std::size_t nrxx, + const std::vector &rho, + const std::vector &mag_part, + const double tpiba, + const Charge* const chr) +{ + std::vector>> gdr(nspin); + std::vector rhor(nrxx); + std::vector> rhog(chr->rhopw->npw); + std::vector> gdr_tmp(nrxx); + + // gradient of the TOTAL charge density (including the core charge), + // matching the up/down decomposition in convert_rho_amag_nspin4: + // grad(rho_up) = 0.5 grad(rho + rho_core) + 0.5 m_hat_mu * grad(m_mu) + // NOTE: rho[ir*nspin+0] is rho_up, not the total density; use chr->rho[0]. + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + rhor[ir] = chr->rho[0][ir] + chr->rho_core[ir]; + } + chr->rhopw->real2recip(rhor.data(), rhog.data()); + gdr[0].resize(nrxx); + XC_Functional::grad_rho(rhog.data(), gdr[0].data(), chr->rhopw, tpiba); + + gdr[1].resize(nrxx); + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + gdr_tmp[ir] = gdr[0][ir]; + gdr[0][ir] = 0.5 * gdr_tmp[ir]; + gdr[1][ir] = 0.5 * gdr_tmp[ir]; + } + + for (int is = 1; is <= 3; ++is) + { + chr->rhopw->real2recip(chr->rho[is], rhog.data()); + XC_Functional::grad_rho(rhog.data(), gdr_tmp.data(), chr->rhopw, tpiba); + const double* mp = mag_part.data() + (is - 1) * nrxx; + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + const ModuleBase::Vector3 g = 0.5 * gdr_tmp[ir] * mp[ir]; + gdr[0][ir] += g; + gdr[1][ir] -= g; + } + } + + return gdr; +} + +std::vector> XC_Functional_Libxc::cal_dh_sf( + const int nspin, + const std::size_t nrxx, + const std::vector &sgn, + const std::vector>> &gdr, + const std::vector &vsigma, + const std::vector &mag_part, + const int gga_grad, + const double tpiba, + const Charge* const chr) +{ + // h1/h2: functional derivative of the GGA energy density w.r.t. the + // spin-up/spin-down density gradients, h_s = df/d(grad rho_s): + // h_up = 2*(2*vsgm(0)*sgn_up^2*gdr_up + vsgm(1)*sgn_up*sgn_dn*gdr_dn) + // h_dn = 2*(2*vsgm(2)*sgn_dn^2*gdr_dn + vsgm(1)*sgn_up*sgn_dn*gdr_up) + std::vector> h1(nrxx), h2(nrxx); + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + h1[ir] = 2.0 * (gdr[0][ir] * vsigma[ir*3 ] * sgn[ir*2 ] * 2.0 + + gdr[1][ir] * vsigma[ir*3+1] * sgn[ir*2] * sgn[ir*2+1]); + h2[ir] = 2.0 * (gdr[1][ir] * vsigma[ir*3+2] * sgn[ir*2+1] * 2.0 + + gdr[0][ir] * vsigma[ir*3+1] * sgn[ir*2] * sgn[ir*2+1]); + } + + std::vector> tmp_h(nrxx); + std::vector dh0(nrxx), dh_mu(nrxx); + + // total-density channel (same for gga_grad=2 and 3): + // dh_0 = div( (h_up + h_dn)/2 ) + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + tmp_h[ir] = 0.5 * (h1[ir] + h2[ir]); + } + XC_Functional::grad_dot(tmp_h.data(), dh0.data(), chr->rhopw, tpiba); + + std::vector> dh_total(4, std::vector(nrxx, 0.0)); +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 1024) +#endif + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + dh_total[0][ir] = dh0[ir]; + } + + if (gga_grad == 2) + { + // gga_grad=2 (projected method): the divergence of the magnetic part + // of h is projected onto m_hat BEFORE differentiating, i.e. the + // cross terms (h_up - h_dn) . grad(m_hat_mu) are dropped: + // dh_mu = m_hat_mu * div( (h_up - h_dn)/2 ) + // Only one divergence (of a single vector field) is needed. + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + tmp_h[ir] = 0.5 * (h1[ir] - h2[ir]); + } + XC_Functional::grad_dot(tmp_h.data(), dh_mu.data(), chr->rhopw, tpiba); + for (int mu = 1; mu < 4; ++mu) + { + const double* mp_mu = mag_part.data() + (mu - 1) * nrxx; + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + dh_total[mu][ir] = mp_mu[ir] * dh_mu[ir]; + } + } + } + else // gga_grad == 3 + { + // gga_grad=3 (full Scalmani-Frisch): the divergence is taken AFTER + // multiplying by m_hat_mu, retaining the cross terms + // (h_up - h_dn) . grad(m_hat_mu): + // dh_mu = div( (h_up - h_dn)/2 * m_hat_mu ) + // = m_hat_mu * div((h_up-h_dn)/2) + (h_up-h_dn)/2 . grad(m_hat_mu) + for (int mu = 1; mu < 4; ++mu) + { + const double* mp_mu = mag_part.data() + (mu - 1) * nrxx; + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + tmp_h[ir] = 0.5 * (h1[ir] - h2[ir]) * mp_mu[ir]; + } + XC_Functional::grad_dot(tmp_h.data(), dh_mu.data(), chr->rhopw, tpiba); + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + dh_total[mu][ir] = dh_mu[ir]; + } + } + } + + return dh_total; +} + +ModuleBase::matrix XC_Functional_Libxc::convert_v_nspin4_sf( + const std::size_t nrxx, + const Charge* const chr, + const std::vector &mag_part, + const ModuleBase::matrix &v) +{ + // Only called for nspin=4 with noncollinear magnetism (use_sf in v_xc_libxc) + constexpr int nspin4 = 4; + constexpr double vanishing_charge = 1.0e-10; + ModuleBase::matrix v_nspin4(nspin4, nrxx); + + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + v_nspin4(0, ir) = 0.5 * (v(0, ir) + v(1, ir)); + } + + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + double amag = std::sqrt(std::pow(chr->rho[1][ir], 2) + + std::pow(chr->rho[2][ir], 2) + + std::pow(chr->rho[3][ir], 2)); + if (amag > vanishing_charge) + { + const double vs = 0.5 * (v(0, ir) - v(1, ir)); + for (int ipol = 1; ipol < nspin4; ++ipol) + { + v_nspin4(ipol, ir) = vs * mag_part[ir + (ipol - 1) * nrxx]; + } + } + } + return v_nspin4; +} + #endif diff --git a/source/source_hamilt/module_xc/test/CMakeLists.txt b/source/source_hamilt/module_xc/test/CMakeLists.txt index 34634eae02..ffcee4e3dd 100644 --- a/source/source_hamilt/module_xc/test/CMakeLists.txt +++ b/source/source_hamilt/module_xc/test/CMakeLists.txt @@ -71,6 +71,7 @@ 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 ../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 diff --git a/source/source_hamilt/module_xc/test/test_xc3.cpp b/source/source_hamilt/module_xc/test/test_xc3.cpp index 3d141accbc..b9eb8367f6 100644 --- a/source/source_hamilt/module_xc/test/test_xc3.cpp +++ b/source/source_hamilt/module_xc/test/test_xc3.cpp @@ -92,13 +92,13 @@ class XCTest_GRADCORR : public XCTest 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,0, hybrid_alpha, hse_omega); + XC_Functional::gradcorr(et1,vt1,v1,&chr,&rhopw,&ucell,stress1,true,nspin1,domag,domag_z,0, hybrid_alpha, hse_omega); - XC_Functional::gradcorr(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,0, hybrid_alpha, hse_omega); + XC_Functional::gradcorr(et2,vt2,v2,&chr,&rhopw,&ucell,stress2,true,nspin2,domag,domag_z,0, hybrid_alpha, hse_omega); - XC_Functional::gradcorr(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,0, hybrid_alpha, hse_omega); } }; diff --git a/source/source_hamilt/module_xc/test/test_xc5.cpp b/source/source_hamilt/module_xc/test/test_xc5.cpp index 01aa214fea..94131c9ea1 100644 --- a/source/source_hamilt/module_xc/test/test_xc5.cpp +++ b/source/source_hamilt/module_xc/test/test_xc5.cpp @@ -82,13 +82,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); @@ -186,13 +186,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); @@ -353,6 +353,432 @@ 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) + * + * Methods 2 and 3 share the same chain-rule spin-up/down gradients + * grad(rho_up/dn) = (grad(rho) +/- m_hat . grad(m))/2, m_hat = m/|m| + * but differ in the divergence of h = df/d(grad rho): + * gga_grad=2 (projected): v_mu -= m_hat_mu * div((h_up - h_dn)/2) + * gga_grad=3 (full SF): v_mu -= div((h_up - h_dn)/2 * m_hat_mu) + * The two are identical when grad(m_hat) = 0 (magnetization direction + * uniform in space) and differ otherwise. + ***********************************************/ + +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; + } + } +}; + +// 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 +TEST(GgaGradTools, ComputeMagPartNspin4) +{ + Charge chr; + chr.rho = new double*[4]; + for (int is = 0; is < 4; ++is) + { + chr.rho[is] = new double[3]; + } + chr.rho[1][0] = 3.0; chr.rho[2][0] = 0.0; chr.rho[3][0] = 4.0; // |m|=5 + chr.rho[1][1] = 0.0; chr.rho[2][1] = 0.0; chr.rho[3][1] = 0.0; // |m|=0 + chr.rho[1][2] = 1.0; chr.rho[2][2] = 2.0; chr.rho[3][2] = 2.0; // |m|=3 + + const std::vector mp = XC_Functional_Libxc::compute_mag_part_nspin4(3, &chr); + + EXPECT_NEAR(mp[0], 3.0 / 5.0, 1e-15); + EXPECT_NEAR(mp[3], 0.0, 1e-15); + EXPECT_NEAR(mp[6], 4.0 / 5.0, 1e-15); + EXPECT_NEAR(mp[1], 0.0, 1e-15); + EXPECT_NEAR(mp[4], 0.0, 1e-15); + EXPECT_NEAR(mp[7], 0.0, 1e-15); + EXPECT_NEAR(mp[2], 1.0 / 3.0, 1e-15); + EXPECT_NEAR(mp[5], 2.0 / 3.0, 1e-15); + EXPECT_NEAR(mp[8], 2.0 / 3.0, 1e-15); +} + +// v_tot = 0.5*(v_up+v_dn), v_mu = 0.5*(v_up-v_dn)*m_hat_mu +TEST(GgaGradTools, ConvertVNspin4Sf) +{ + Ns4Charge mock(0); // m_hat = (0,0,1) + const std::vector mag_part + = XC_Functional_Libxc::compute_mag_part_nspin4(gga_grad_nrxx, &mock.chr); + + 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 v4 + = XC_Functional_Libxc::convert_v_nspin4_sf(gga_grad_nrxx, &mock.chr, mag_part, v); + + for (int ir = 0; ir < gga_grad_nrxx; ++ir) + { + EXPECT_NEAR(v4(0, ir), 0.5 * (v(0, ir) + v(1, ir)), 1e-14); + EXPECT_NEAR(v4(1, ir), 0.0, 1e-14); + EXPECT_NEAR(v4(2, ir), 0.0, 1e-14); + EXPECT_NEAR(v4(3, ir), 0.5 * (v(0, ir) - v(1, ir)), 1e-14); + } +} + +// 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); + } +} + +class GgaGradDh : public testing::Test +{ + protected: + // dh from cal_dh_sf for gga_grad=2 and 3 on a given magnetization pattern + void run(const int pattern, + std::vector>& dh2, + std::vector>& dh3, + std::vector& mag_part) + { + Ns4Charge mock(pattern); + mag_part = XC_Functional_Libxc::compute_mag_part_nspin4(gga_grad_nrxx, &mock.chr); + + const std::tuple, std::vector> rho_amag + = XC_Functional_Libxc::convert_rho_amag_nspin4(2, gga_grad_nrxx, &mock.chr); + const std::vector& rho = std::get<0>(rho_amag); + const std::vector>> gdr + = XC_Functional_Libxc::cal_gdr_sf(2, gga_grad_nrxx, rho, mag_part, mock.ucell.tpiba, &mock.chr); + + std::vector sgn(gga_grad_nrxx * 2, 1.0); + std::vector vsigma(gga_grad_nrxx * 3); + for (int ir = 0; ir < gga_grad_nrxx; ++ir) + { + for (int j = 0; j < 3; ++j) + { + vsigma[ir * 3 + j] = 0.2 + 0.1 * ir + 0.05 * j; + } + } + + dh2 = XC_Functional_Libxc::cal_dh_sf( + 2, gga_grad_nrxx, sgn, gdr, vsigma, mag_part, 2, mock.ucell.tpiba, &mock.chr); + dh3 = XC_Functional_Libxc::cal_dh_sf( + 2, gga_grad_nrxx, sgn, gdr, vsigma, mag_part, 3, mock.ucell.tpiba, &mock.chr); + } +}; + +// uniform m_hat => grad(m_hat) = 0 => projected and full SF divergences agree +TEST_F(GgaGradDh, UniformDirectionMethodsAgree) +{ + std::vector> dh2, dh3; + std::vector mag_part; + run(0, dh2, dh3, mag_part); + + for (int is = 0; is < 4; ++is) + { + for (int ir = 0; ir < gga_grad_nrxx; ++ir) + { + EXPECT_NEAR(dh2[is][ir], dh3[is][ir], 1e-12); + } + } +} + +// for gga_grad=2, dh_mu = m_hat_mu * div((h_up-h_dn)/2), so the magnetic +// channels satisfy dh_mu = m_hat_mu * (m_hat . dh) exactly +TEST_F(GgaGradDh, ProjectedDivergenceIsProjection) +{ + std::vector> dh2, dh3; + std::vector mag_part; + run(1, dh2, dh3, mag_part); + + for (int ir = 0; ir < gga_grad_nrxx; ++ir) + { + double proj = 0.0; + for (int mu = 1; mu < 4; ++mu) + { + proj += mag_part[ir + (mu - 1) * gga_grad_nrxx] * dh2[mu][ir]; + } + for (int mu = 1; mu < 4; ++mu) + { + EXPECT_NEAR(dh2[mu][ir], mag_part[ir + (mu - 1) * gga_grad_nrxx] * proj, 1e-12); + } + } +} + +// NOTE on the mocked FFT: the mock derivative is pointwise +// (grad(f)[ir] ~ f[ir], div(h)[ir] ~ h[ir]), so multiplication by a +// scalar field commutes with the divergence and the projected (2) and +// full (3) SF divergences coincide under this mock. The two methods can +// only be distinguished with a real (nonlocal) FFT, e.g. in integration +// tests. The tests below therefore anchor the wiring (exact values, +// projection identities, crash-free dispatch) rather than the 2-vs-3 +// numerical difference. + +// built-in functionals: v_xc dispatches nspin=4 + gga_grad=2/3 to the SF builtin +TEST(GgaGradVxc, BuiltinUniformDirectionMethodsAgree) +{ + const auto r2 = run_vxc_nspin4("PBE", 0, 2); + const auto r3 = run_vxc_nspin4("PBE", 0, 3); + EXPECT_EQ(std::get<2>(r2).nr, 4); + expect_vxc_equal(r2, r3, 1e-10); +} + +// gga_grad=0 keeps the original built-in algorithm and must not crash +TEST(GgaGradVxc, BuiltinOriginalAlgorithmRuns) +{ + const auto 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], my = mock.chr.rho[2][ir], 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], my = mock.chr.rho[2][ir], 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, 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, 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); + } + } +} + +// regression anchors for the built-in SF path (gga_grad=3) +TEST(GgaGradVxc, BuiltinSfAnchoredValues) +{ + const auto r = run_vxc_nspin4("PBE", 1, 3); + const ModuleBase::matrix& v = std::get<2>(r); + EXPECT_NEAR(std::get<0>(r), -5.1906253324e+01, 1.0e-8); + EXPECT_NEAR(std::get<1>(r), -6.8184946486e+01, 1.0e-8); + EXPECT_NEAR(v(0, 0), -2.6284212276e+00, 1.0e-8); + EXPECT_NEAR(v(0, 4), -3.7476729308e+00, 1.0e-8); + EXPECT_NEAR(v(1, 0), -3.7774725540e-02, 1.0e-8); + EXPECT_NEAR(v(1, 4), -9.2204398939e-02, 1.0e-8); + EXPECT_NEAR(v(2, 0), -9.4436813851e-02, 1.0e-8); + EXPECT_NEAR(v(2, 4), -9.2204398939e-03, 1.0e-8); + EXPECT_NEAR(v(3, 0), -7.5549451081e-02, 1.0e-8); + EXPECT_NEAR(v(3, 4), -1.8440879788e-01, 1.0e-8); +} + +// LIBXC functionals: gga_grad=2/3 select the SF path in v_xc_libxc +TEST(GgaGradVxc, LibxcUniformDirectionMethodsAgree) +{ + const auto r2 = run_vxc_nspin4("GGA_X_PBE+GGA_C_PBE", 0, 2); + const auto r3 = run_vxc_nspin4("GGA_X_PBE+GGA_C_PBE", 0, 3); + EXPECT_EQ(std::get<2>(r2).nr, 4); + expect_vxc_equal(r2, r3, 1e-10); +} + +// for LIBXC, gga_grad=0 and 1 both keep the original collinear algorithm +TEST(GgaGradVxc, LibxcZeroEqualsOne) +{ + const auto r0 = run_vxc_nspin4("GGA_X_PBE+GGA_C_PBE", 1, 0); + const auto r1 = run_vxc_nspin4("GGA_X_PBE+GGA_C_PBE", 1, 1); + expect_vxc_equal(r0, r1, 1e-12); +} + +// regression anchors for the LIBXC SF path (gga_grad=3); this path goes +// through the SF branch of convert_vtxc_v +TEST(GgaGradVxc, LibxcSfAnchoredValues) +{ + const auto r = run_vxc_nspin4("GGA_X_PBE+GGA_C_PBE", 1, 3); + const ModuleBase::matrix& v = std::get<2>(r); + EXPECT_NEAR(std::get<0>(r), -5.1906239921e+01, 1.0e-8); + EXPECT_NEAR(std::get<1>(r), -6.8184928281e+01, 1.0e-8); + EXPECT_NEAR(v(0, 0), -2.6284203722e+00, 1.0e-8); + EXPECT_NEAR(v(0, 4), -3.7476719779e+00, 1.0e-8); + EXPECT_NEAR(v(1, 0), -3.7774739650e-02, 1.0e-8); + EXPECT_NEAR(v(1, 4), -9.2204427285e-02, 1.0e-8); + EXPECT_NEAR(v(2, 0), -9.4436849124e-02, 1.0e-8); + EXPECT_NEAR(v(2, 4), -9.2204427285e-03, 1.0e-8); + EXPECT_NEAR(v(3, 0), -7.5549479300e-02, 1.0e-8); + EXPECT_NEAR(v(3, 4), -1.8440885457e-01, 1.0e-8); +} + int main(int argc, char **argv) { diff --git a/source/source_hamilt/module_xc/xc_functional.h b/source/source_hamilt/module_xc/xc_functional.h index 540469b81e..9e1248a2e7 100644 --- a/source/source_hamilt/module_xc/xc_functional.h +++ b/source/source_hamilt/module_xc/xc_functional.h @@ -51,6 +51,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); @@ -225,6 +226,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); @@ -249,15 +251,13 @@ class XC_Functional double* dh, const ModulePW::PW_Basis* rho_basis, const double tpiba); - - static void noncolin_rho( - double* rhoout1, - double* rhoout2, - double* seg, - const double* const* const rho, - const int nrxx, - const double* ux_, - const bool lsign_); + static void noncolin_rho(double* rhoout1, + double* rhoout2, + double* neg, + const double* const* const rho, + const int nrxx, + const double* ux_, + const bool lsign_); //------------------- // xc_lda_exch.cpp 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 0000000000..5c54e8b7db --- /dev/null +++ b/source/source_hamilt/module_xc/xc_functional_ncgga_sf.cpp @@ -0,0 +1,376 @@ +#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_cell/unitcell.h" +#include "source_estate/module_charge/charge.h" +#include "xc_functional.h" + +#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, + const int gga_grad) +{ + ModuleBase::TITLE("XC_Functional", "v_xc_ncgga_sf_builtin"); + ModuleBase::timer::start("XC_Functional", "v_xc_ncgga_sf_builtin"); + + // Caller (XC_Functional::v_xc) guarantees nspin==4 with noncollinear magnetism + // and gga_grad==2 or 3. + + // ====================================================================== + // Scalmani-Frisch (SF) builtin for noncollinear GGA (gga_grad=2/3) + // + // Reference: Scalmani & Frisch, JCTC 8, 1069 (2012) + // + // Key formulas: + // rho_up = 0.5*(rho + |m|), rho_dn = 0.5*(rho - |m|) + // m_hat = m / |m| (magnetization unit vector) + // + // Gradient decomposition (chain rule via m_hat), identical for 2 and 3: + // grad(rho_up) = 0.5 grad(rho) + 0.5 m_hat_mu * grad(m_mu) + // grad(rho_dn) = 0.5 grad(rho) - 0.5 m_hat_mu * grad(m_mu) + // + // For GGA, the sigma = |grad(rho_up)|^2 etc. are computed from gdr1,gdr2, + // then passed to spin-polarized XC functionals (xc_spin, gcx_spin, gcc_spin). + // The resulting potential is converted back to nspin=4 representation: + // v_tot = 0.5*(v_up + v_dn) + // v_mag = 0.5*(v_up - v_dn) * m_hat + // + // The two methods differ in the divergence of h_s = df/d(grad rho_s): + // gga_grad=2 (projected): v_mu -= m_hat_mu * div( (h_up - h_dn)/2 ) + // the cross terms (h_up - h_dn) . grad(m_hat_mu) are dropped; + // only one divergence of a single vector field is computed. + // gga_grad=3 (full SF): v_mu -= div( (h_up - h_dn)/2 * m_hat_mu ) + // the divergence is taken after multiplying by m_hat_mu, retaining + // all cross terms; this is the most accurate method. + // ====================================================================== + + 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 double fac = 0.5; + const bool is_gga = (XC_Functional::get_func_type() == 2 || XC_Functional::get_func_type() == 4); + + // Step 1: decompose charge density into spin-up/spin-down and compute m_hat + std::vector rhotmp1(nrxx), rhotmp2(nrxx), amag(nrxx); + std::vector mag_part(3 * nrxx, 0.0); + + for (int ir = 0; ir < nrxx; ++ir) + { + const double mx = chr->rho[1][ir], my = chr->rho[2][ir], mz = chr->rho[3][ir]; + amag[ir] = std::sqrt(mx * mx + my * my + mz * mz); + rhotmp1[ir] = 0.5 * (chr->rho[0][ir] + amag[ir]); + rhotmp2[ir] = 0.5 * (chr->rho[0][ir] - amag[ir]); + if (amag[ir] > 1e-12) + { + mag_part[ir] = mx / amag[ir]; + mag_part[ir + nrxx] = my / amag[ir]; + mag_part[ir + 2 * nrxx] = mz / amag[ir]; + } + } + for (int ir = 0; ir < nrxx; ++ir) + { + rhotmp1[ir] += fac * chr->rho_core[ir]; + rhotmp2[ir] += fac * chr->rho_core[ir]; + } + + std::vector> rhogsum1(npw), 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), gdr2(nrxx); + std::vector> gdr_mag(nrxx); + XC_Functional::grad_rho(rhogsum1.data(), gdr1.data(), rhopw, tpiba); + + for (int ir = 0; ir < nrxx; ++ir) + { + gdr_mag[ir] = gdr1[ir]; + gdr1[ir] = 0.5 * gdr_mag[ir]; + gdr2[ir] = 0.5 * gdr_mag[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); + const double* mp = mag_part.data() + (is - 1) * nrxx; + for (int ir = 0; ir < nrxx; ++ir) + { + const ModuleBase::Vector3 g = 0.5 * gdr_mag[ir] * mp[ir]; + gdr1[ir] += g; + gdr2[ir] -= g; + } + } + + // Step 3: LDA contribution (xc_spin) and convert to nspin=4 potential + // v(0,:) = 0.5*(vxc[0] + vxc[1]) -- total density channel + // v(1..3,:) = 0.5*(vxc[0] - vxc[1]) * m_mu / |m| -- magnetization channels + double etxc = 0, vtxc = 0; + ModuleBase::matrix v(4, nrxx); + + for (int ir = 0; ir < nrxx; ++ir) + { + const double arho = std::abs(chr->rho[0][ir] + chr->rho_core[ir]); + if (arho <= vanishing) + continue; + + double zeta = amag[ir] / arho; + if (std::abs(zeta) > 1.0) + zeta = (zeta > 0) ? 1.0 : -1.0; + double exc = 0, vxc[2] = {0, 0}; + XC_Functional::xc_spin(arho, zeta, exc, vxc[0], vxc[1]); + + v(0, ir) = e2 * 0.5 * (vxc[0] + vxc[1]); + vtxc += v(0, ir) * chr->rho[0][ir]; + + if (amag[ir] > vanishing) + { + const double vs = e2 * 0.5 * (vxc[0] - vxc[1]); + const double inv_a = 1.0 / amag[ir]; + for (int mu = 1; mu < 4; ++mu) + { + v(mu, ir) = vs * chr->rho[mu][ir] * inv_a; + vtxc += v(mu, ir) * chr->rho[mu][ir]; + } + } + etxc += e2 * exc * arho; + } + + // Step 4: GGA contribution (gcx_spin + gcc_spin) and divergence correction (grad_dot) + // The gradient correction uses the same SF-decomposed gdr1, gdr2 and converts + // back to nspin=4 via m_hat, matching the LDA conversion above. + if (is_gga) + { + double etxcgc = 0, vtxcgc = 0; + std::vector vup_gga(nrxx, 0), vdw_gga(nrxx, 0); + std::vector> h1(nrxx), h2(nrxx); + + for (int ir = 0; ir < nrxx; ++ir) + { + double sx = 0, v1xup = 0, v1xdw = 0, v2xup = 0, v2xdw = 0; + double sc = 0, v1cup = 0, v1cdw = 0, v2c = 0; + const double grho2a = gdr1[ir] * gdr1[ir]; + const 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) + { + double zeta = (rhotmp1[ir] - rhotmp2[ir]) / rh; + zeta = std::fabs(zeta); + const double grh2 = (gdr1[ir] + gdr2[ir]) * (gdr1[ir] + gdr2[ir]); + XC_Functional::gcc_spin(rh, zeta, grh2, sc, v1cup, v1cdw, v2c); + } + + vup_gga[ir] = e2 * (v1xup + v1cup); + vdw_gga[ir] = e2 * (v1xdw + v1cdw); + + const double v2cup = v2c, v2cdw = v2c, v2cud = v2c; + h1[ir] = e2 * ((v2xup + v2cup) * gdr1[ir] + v2cud * gdr2[ir]); + h2[ir] = e2 * ((v2xdw + v2cdw) * gdr2[ir] + v2cud * gdr1[ir]); + + vtxcgc += vup_gga[ir] * (rhotmp1[ir] - chr->rho_core[ir] * fac); + vtxcgc += vdw_gga[ir] * (rhotmp2[ir] - chr->rho_core[ir] * fac); + etxcgc += e2 * (sx + sc); + } + + for (int ir = 0; ir < nrxx; ++ir) + { + v(0, ir) += 0.5 * (vup_gga[ir] + vdw_gga[ir]); + const double vdiff = 0.5 * (vup_gga[ir] - vdw_gga[ir]); + for (int mu = 1; mu < 4; ++mu) + { + v(mu, ir) += vdiff * mag_part[ir + (mu - 1) * nrxx]; + } + } + + std::vector dh(nrxx); + std::vector> tmp_h(nrxx); + + // total-density channel (same for gga_grad=2 and 3): + // v(0,:) -= div( (h_up + h_dn)/2 ) + for (int ir = 0; ir < nrxx; ++ir) + tmp_h[ir] = 0.5 * (h1[ir] + h2[ir]); + XC_Functional::grad_dot(tmp_h.data(), dh.data(), rhopw, tpiba); + for (int ir = 0; ir < nrxx; ++ir) + v(0, ir) -= dh[ir]; + double sum = 0; + for (int ir = 0; ir < nrxx; ++ir) + sum += dh[ir] * chr->rho[0][ir]; + vtxcgc -= sum; + + if (gga_grad == 2) + { + // gga_grad=2 (projected method): + // v(mu,:) -= m_hat_mu * div( (h_up - h_dn)/2 ) + // The divergence of (h_up - h_dn)/2 is computed once and then + // projected onto m_hat; the cross terms + // (h_up - h_dn) . grad(m_hat_mu) are dropped. + for (int ir = 0; ir < nrxx; ++ir) + tmp_h[ir] = 0.5 * (h1[ir] - h2[ir]); + XC_Functional::grad_dot(tmp_h.data(), dh.data(), rhopw, tpiba); + for (int mu = 1; mu < 4; ++mu) + { + const double* mp = mag_part.data() + (mu - 1) * nrxx; + double sum_mu = 0; + for (int ir = 0; ir < nrxx; ++ir) + { + const double dh_mu = mp[ir] * dh[ir]; + v(mu, ir) -= dh_mu; + sum_mu += dh_mu * chr->rho[mu][ir]; + } + vtxcgc -= sum_mu; + } + } + else // gga_grad == 3 + { + // gga_grad=3 (full Scalmani-Frisch): + // v(mu,:) -= div( (h_up - h_dn)/2 * m_hat_mu ) + // The divergence is taken after multiplying by m_hat_mu, + // retaining the cross terms (h_up - h_dn) . grad(m_hat_mu). + for (int mu = 1; mu < 4; ++mu) + { + const double* mp = mag_part.data() + (mu - 1) * nrxx; + for (int ir = 0; ir < nrxx; ++ir) + tmp_h[ir] = 0.5 * (h1[ir] - h2[ir]) * mp[ir]; + XC_Functional::grad_dot(tmp_h.data(), dh.data(), rhopw, tpiba); + for (int ir = 0; ir < nrxx; ++ir) + v(mu, ir) -= dh[ir]; + double sum_mu = 0; + for (int ir = 0; ir < nrxx; ++ir) + sum_mu += dh[ir] * chr->rho[mu][ir]; + vtxcgc -= sum_mu; + } + } + + etxc += etxcgc; + vtxc += vtxcgc; + } + +#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_sf_builtin(const Charge* const chr, ModulePW::PW_Basis* rhopw, + const UnitCell* ucell, 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; + constexpr double small = 1.0e-10; + const double fac = 0.5; + + std::vector rhotmp1(nrxx), rhotmp2(nrxx), amag_arr(nrxx); + std::vector mag_part(3 * nrxx, 0.0); + + for (int ir = 0; ir < nrxx; ++ir) + { + const double mx = chr->rho[1][ir], my = chr->rho[2][ir], mz = chr->rho[3][ir]; + double amag = std::sqrt(mx * mx + my * my + mz * mz); + amag_arr[ir] = amag; + rhotmp1[ir] = 0.5 * (chr->rho[0][ir] + amag); + rhotmp2[ir] = 0.5 * (chr->rho[0][ir] - amag); + if (amag > 1e-12) + { + mag_part[ir] = mx / amag; + mag_part[ir + nrxx] = my / amag; + mag_part[ir + 2 * nrxx] = mz / amag; + } + } + for (int ir = 0; ir < nrxx; ++ir) + { + rhotmp1[ir] += fac * chr->rho_core[ir]; + rhotmp2[ir] += fac * chr->rho_core[ir]; + } + + std::vector> rhogsum1(npw), 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), gdr2(nrxx); + std::vector> gdr_mag(nrxx); + XC_Functional::grad_rho(rhogsum1.data(), gdr1.data(), rhopw, ucell->tpiba); + + for (int ir = 0; ir < nrxx; ++ir) + { + gdr_mag[ir] = gdr1[ir]; + gdr1[ir] = 0.5 * gdr_mag[ir]; + gdr2[ir] = 0.5 * gdr_mag[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, ucell->tpiba); + const double* mp = mag_part.data() + (is - 1) * nrxx; + for (int ir = 0; ir < nrxx; ++ir) + { + const ModuleBase::Vector3 g = 0.5 * gdr_mag[ir] * mp[ir]; + gdr1[ir] += g; + gdr2[ir] -= g; + } + } + + for (int ir = 0; ir < nrxx; ++ir) + { + double sx = 0, v1xup = 0, v1xdw = 0, v2xup = 0, v2xdw = 0; + double sc = 0, v1cup = 0, v1cdw = 0, v2c = 0; + double v2cup = 0, v2cdw = 0, v2cud = 0; + const double grho2a = gdr1[ir] * gdr1[ir]; + const 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) + { + double zeta = (rhotmp1[ir] - rhotmp2[ir]) / rh; + zeta = std::fabs(zeta); + const double grh2 = (gdr1[ir] + gdr2[ir]) * (gdr1[ir] + gdr2[ir]); + XC_Functional::gcc_spin(rh, zeta, grh2, sc, v1cup, v1cdw, v2c); + v2cup = v2c; + v2cdw = v2c; + v2cud = v2c; + } + + double tt1[3] = {gdr1[ir].x, gdr1[ir].y, gdr1[ir].z}; + double tt2[3] = {gdr2[ir].x, gdr2[ir].y, gdr2[ir].z}; + for (int l = 0; l < 3; l++) + { + for (int m = 0; m < l + 1; m++) + { + int ind = l * 3 + m; + stress_gga[ind] += e2 * (tt1[l] * tt1[m] * v2xup + tt2[l] * tt2[m] * v2xdw); + stress_gga[ind] += e2 * (tt1[l] * tt1[m] * v2cup + tt2[l] * tt2[m] * v2cdw + + (tt1[l] * tt2[m] + tt2[l] * tt1[m]) * v2cud); + } + } + } +} + +} // 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 0000000000..34cc8276cd --- /dev/null +++ b/source/source_hamilt/module_xc/xc_functional_ncgga_sf.h @@ -0,0 +1,32 @@ +#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; +} +struct UnitCell; + +namespace ModuleXC +{ +namespace NCGGA_SF_Builtin +{ + +// gga_grad: 2 = projected divergence of h, 3 = full Scalmani-Frisch divergence +std::tuple v_xc_ncgga_sf_builtin( + const int& nrxx, const double& omega, const double tpiba, const Charge* const chr, + const int gga_grad); + +void gradcorr_ncgga_sf_builtin(const Charge* const chr, ModulePW::PW_Basis* rhopw, + const UnitCell* ucell, 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 12788b60a6..1a962a3f17 100644 --- a/source/source_hamilt/module_xc/xc_grad.cpp +++ b/source/source_hamilt/module_xc/xc_grad.cpp @@ -36,6 +36,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) { @@ -204,7 +205,18 @@ void XC_Functional::gradcorr( vgg[is] = new double[rhopw->nrxx]; } } - noncolin_rho(rhotmp1, rhotmp2, neg, chr->rho, rhopw->nrxx, ucell->magnet.ux_, ucell->magnet.lsign_); + // gga_grad=0 (original): when the initial magnetic moments in STRU are + // all collinear, lsign_ is set and their common direction ux_ is used + // as a global quantization axis: up/down are defined w.r.t. this axis + // through the sign of m . ux. This reproduces the nspin=2 collinear + // result exactly for collinear configurations, but makes the magnetic + // potential energy surface discontinuous when the moments are tilted + // slightly away from collinearity. + // gga_grad=1: ignore the global direction (neg = 1 everywhere), so the + // magnetic potential energy surface stays continuous, at the price of + // giving up exact consistency with nspin=2 for collinear cases. + const bool use_global_axis = ucell->magnet.lsign_ && (gga_grad != 1); + noncolin_rho(rhotmp1, rhotmp2, neg, chr->rho, rhopw->nrxx, ucell->magnet.ux_, use_global_axis); rhopw->real2recip(rhotmp1, rhogsum1); rhopw->real2recip(rhotmp2, rhogsum2); #ifdef _OPENMP diff --git a/source/source_hamilt/module_xc/xc_pot.cpp b/source/source_hamilt/module_xc/xc_pot.cpp index 1f8a3dcd32..a2e9a07fa2 100644 --- a/source/source_hamilt/module_xc/xc_pot.cpp +++ b/source/source_hamilt/module_xc/xc_pot.cpp @@ -6,9 +6,10 @@ #include "source_base/parallel_reduce.h" #include "source_base/timer.h" -#include "source_io/module_parameter/parameter.h" #include "xc_functional.h" +#include "xc_functional_ncgga_sf.h" + #ifdef __LIBXC #include "libxc_abacus.h" #ifdef __EXX @@ -16,6 +17,8 @@ #endif #endif + + // [etxc, vtxc, v] = XC_Functional::v_xc(...) std::tuple XC_Functional::v_xc( const int& nrxx, @@ -24,6 +27,7 @@ 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) { @@ -40,6 +44,7 @@ std::tuple XC_Functional::v_xc( nspin, domag, domag_z, + gga_grad, &(scaling_factor_xc), hybrid_alpha, hse_omega); @@ -48,6 +53,16 @@ std::tuple XC_Functional::v_xc( #endif } + // For non-libxc builds: gga_grad=2/3 uses the SF builtin that computes the + // full nspin=4 GGA potential via the Scalmani-Frisch transformation + // (gga_grad=2: projected divergence, gga_grad=3: full divergence). + // This path handles spin-up/spin-down decomposition and gradient + // corrections internally, returning (etxc, vtxc, v) directly. + if (nspin == 4 && (domag || domag_z) && (gga_grad == 2 || gga_grad == 3)) + { + return ModuleXC::NCGGA_SF_Builtin::v_xc_ncgga_sf_builtin(nrxx, ucell->omega, ucell->tpiba, chr, gga_grad); + } + ModuleBase::timer::start("XC_Functional", "v_xc"); //Exchange-Correlation potential Vxc(r) from n(r) @@ -184,7 +199,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 481f49cbed..404df69bf4 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 127e459548..8885d7e854 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -87,6 +87,7 @@ struct Input_para bool pseudo_mesh = false; ///< 0: use msh to normalize radial wave functions; 1: ///< use mesh, which is used in QE. int nspin = 1; ///< LDA ; LSDA ; non-linear spin + int gga_grad = 0; ///< GGA gradient method for noncollinear spin (nspin=4): 0=original algorithm (default), 1=collinear approx, 2=projected div(h), 3=Scalmani-Frisch transform int pw_diag_nmax = 50; double pw_diag_thr = 0.01; ///< used in cg method bool diago_smooth_ethr = false; ///< smooth ethr for iter methods diff --git a/source/source_io/module_parameter/read_input_item_elec_stru.cpp b/source/source_io/module_parameter/read_input_item_elec_stru.cpp index bc5ff4f4de..c370083a2c 100644 --- a/source/source_io/module_parameter/read_input_item_elec_stru.cpp +++ b/source/source_io/module_parameter/read_input_item_elec_stru.cpp @@ -501,6 +501,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 = "GGA gradient method for nspin=4: 0 original algorithm, 1 collinear approx, 2 projected div(h), 3 Scalmani-Frisch"; + item.category = "Electronic structure"; + item.type = "Integer"; + item.description = R"(Method used to evaluate the density gradient entering GGA exchange-correlation terms in noncollinear-spin (nspin=4) calculations. With rho_up/dn = (rho +/- |m|)/2 and m_hat = m/|m|: +* 0: original algorithm; the default. When the initial magnetic moments in STRU are all collinear, their common direction u is used as a global quantization axis and the up/down densities are defined w.r.t. this axis via the sign of m . u. This reproduces the nspin=2 collinear result exactly for collinear configurations, but makes the magnetic potential energy surface discontinuous when the initial moments are tilted slightly away from collinearity. +* 1: collinear approximation without the global direction: up/down are always defined w.r.t. the local |m|, so the magnetic potential energy surface stays continuous, at the price of giving up exact consistency with nspin=2 for collinear configurations. (For LIBXC functionals the global axis is never used, so 0 and 1 are equivalent.) +* 2: projected method. The gradients use the full chain rule, grad(rho_up/dn) = (grad(rho) +/- m_hat . grad(m))/2, but the divergence of h = df/d(grad rho) in the potential is projected onto m_hat: v_mu -= m_hat_mu * div((h_up - h_dn)/2), dropping the (h_up - h_dn) . grad(m_hat_mu) cross terms. +* 3: Scalmani-Frisch transformation (G. Scalmani and M. J. Frisch, J. Chem. Theory Comput. 8, 2193 (2012)). Same gradients as 2, but the full divergence is kept: v_mu -= div((h_up - h_dn)/2 * m_hat_mu), retaining all cross terms; the most accurate. +This parameter only takes effect for nspin=4 with GGA functionals (and magnetic calculation).)"; + 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 > 3) + { + ModuleBase::WARNING_QUIT("ReadInput", "gga_grad must be 0, 1, 2, or 3."); + } + }; + 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_pw/module_pwdft/forces_cc.cpp b/source/source_pw/module_pwdft/forces_cc.cpp index 2daf164bcf..b3fced1525 100644 --- a/source/source_pw/module_pwdft/forces_cc.cpp +++ b/source/source_pw/module_pwdft/forces_cc.cpp @@ -82,6 +82,7 @@ void Forces::cal_force_cc(ModuleBase::matrix& forcecc, PARAM.inp.nspin, PARAM.globalv.domag, PARAM.globalv.domag_z, + PARAM.inp.gga_grad, hybrid_alpha, hse_omega); diff --git a/source/source_pw/module_pwdft/stress_cc.cpp b/source/source_pw/module_pwdft/stress_cc.cpp index ffe5dc85e0..bf1a3ae6c9 100644 --- a/source/source_pw/module_pwdft/stress_cc.cpp +++ b/source/source_pw/module_pwdft/stress_cc.cpp @@ -78,6 +78,7 @@ void Stress_Func::stress_cc(ModuleBase::matrix& sigma, PARAM.inp.nspin, PARAM.globalv.domag, PARAM.globalv.domag_z, + PARAM.inp.gga_grad, hybrid_alpha, hse_omega); // etxc = std::get<0>(etxc_vtxc_v); // may delete? diff --git a/source/source_pw/module_pwdft/stress_gga.cpp b/source/source_pw/module_pwdft/stress_gga.cpp index bb17551480..536d97eeaf 100644 --- a/source/source_pw/module_pwdft/stress_gga.cpp +++ b/source/source_pw/module_pwdft/stress_gga.cpp @@ -31,7 +31,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, + PARAM.inp.nspin, PARAM.globalv.domag, PARAM.globalv.domag_z, PARAM.inp.gga_grad, hybrid_alpha, hse_omega); for(int l = 0;l< 3;l++)