diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index f04f903448..0ed36a4690 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -1268,6 +1268,18 @@ - 4: For the case of noncollinear polarized, nspin will be automatically set to 4 without being specified by the user. - **Default**: 1 +### gga_grad + +- **Type**: Integer +- **Description**: Method used to evaluate the density gradient entering GGA exchange-correlation terms in noncollinear-spin (nspin=4) calculations. With rho_up/dn = (rho +/- |m|)/2 and m_hat = m/|m|: + - 0: original algorithm, identical to the behavior before this keyword was introduced; the default. + - 1: collinear approximation, grad(rho_up/dn) = (grad(rho) +/- grad(|m|))/2; only the gradient of the magnetization magnitude is used (equivalent to 0). + - 2: projected method. The gradients use the full chain rule, grad(rho_up/dn) = (grad(rho) +/- m_hat . grad(m))/2, but the divergence of h = df/d(grad rho) in the potential is projected onto m_hat: v_mu -= m_hat_mu * div((h_up - h_dn)/2), dropping the (h_up - h_dn) . grad(m_hat_mu) cross terms. + - 3: Scalmani-Frisch transformation (G. Scalmani and M. J. Frisch, J. Chem. Theory Comput. 8, 2193 (2012)). Same gradients as 2, but the full divergence is kept: v_mu -= div((h_up - h_dn)/2 * m_hat_mu), retaining all cross terms; the most accurate. + + This parameter only takes effect for nspin=4 with GGA functionals (and magnetic calculation). +- **Default**: 0 + ### smearing_method - **Type**: String diff --git a/docs/parameters.yaml b/docs/parameters.yaml index 005fad3076..ebe352f69e 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, identical to the behavior before this keyword was introduced; the default. + * 1: collinear approximation, grad(rho_up/dn) = (grad(rho) +/- grad(|m|))/2; only the gradient of the magnetization magnitude is used (equivalent to 0). + * 2: projected method. The gradients use the full chain rule, grad(rho_up/dn) = (grad(rho) +/- m_hat . grad(m))/2, but the divergence of h = df/d(grad rho) in the potential is projected onto m_hat: v_mu -= m_hat_mu * div((h_up - h_dn)/2), dropping the (h_up - h_dn) . grad(m_hat_mu) cross terms. + * 3: Scalmani-Frisch transformation (G. Scalmani and M. J. Frisch, J. Chem. Theory Comput. 8, 2193 (2012)). Same gradients as 2, but the full divergence is kept: v_mu -= div((h_up - h_dn)/2 * m_hat_mu), retaining all cross terms; the most accurate. + This parameter only takes effect for nspin=4 with GGA functionals (and magnetic calculation). + default_value: "0" + unit: "" + availability: "" - name: smearing_method category: Electronic structure type: String diff --git a/source/source_estate/module_pot/pot_xc.cpp b/source/source_estate/module_pot/pot_xc.cpp index 19815f843a..a0f9770c54 100644 --- a/source/source_estate/module_pot/pot_xc.cpp +++ b/source/source_estate/module_pot/pot_xc.cpp @@ -54,6 +54,7 @@ void PotXC::cal_v_eff(const Charge*const chg, const UnitCell*const ucell, Module PARAM.inp.nspin, PARAM.globalv.domag, PARAM.globalv.domag_z, + PARAM.inp.gga_grad, hybrid_alpha, hse_omega); *(this->etxc_) = std::get<0>(etxc_vtxc_v); 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..4a016ed04d 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_gga_noncol_sf_builtin.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..e0cb52d467 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,199 @@ 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); + + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + rhor[ir] = rho[ir * nspin + 0]; + } + 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..e014d905cc 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_gga_noncol_sf_builtin.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..93e9dfc3c0 100644 --- a/source/source_hamilt/module_xc/test/test_xc3.cpp +++ b/source/source_hamilt/module_xc/test/test_xc3.cpp @@ -4,7 +4,6 @@ #include "../exx_info.h" #include "xc3_mock.h" #include "source_base/matrix.h" -#include "source_cell/cal_ux.h" /************************************************ * unit test of functionals @@ -54,7 +53,7 @@ class XCTest_GRADCORR : public XCTest ucell.tpiba = 1; ucell.magnet.lsign_ = true; - unitcell::cal_ux(ucell, 4); + elecstate::cal_ux(ucell, 4); chr.rho = new double*[4]; chr.rho[0] = new double[5]; diff --git a/source/source_hamilt/module_xc/test/test_xc5.cpp b/source/source_hamilt/module_xc/test/test_xc5.cpp index 01aa214fea..01a7826aca 100644 --- a/source/source_hamilt/module_xc/test/test_xc5.cpp +++ b/source/source_hamilt/module_xc/test/test_xc5.cpp @@ -5,7 +5,6 @@ #include "../exx_info.h" #include "xc3_mock.h" #include "source_base/matrix.h" -#include "source_cell/cal_ux.h" #include "../../../source_base/parallel_reduce.h" /************************************************ @@ -48,7 +47,7 @@ class XCTest_VXC : public XCTest ucell.tpiba = 1; ucell.magnet.lsign_ = true; - unitcell::cal_ux(ucell, 4); + elecstate::cal_ux(ucell, 4); ucell.omega = 1; chr.rhopw = &(rhopw); @@ -82,13 +81,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); @@ -152,7 +151,7 @@ class XCTest_VXC_Libxc : public XCTest ucell.tpiba = 1; ucell.magnet.lsign_ = true; - unitcell::cal_ux(ucell, 4); + elecstate::cal_ux(ucell, 4); ucell.omega = 1; chr.rhopw = &(rhopw); @@ -186,13 +185,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); @@ -254,7 +253,7 @@ class XCTest_VXC_meta : public XCTest ucell.tpiba = 1; ucell.magnet.lsign_ = true; - unitcell::cal_ux(ucell, 4); + elecstate::cal_ux(ucell, 4); ucell.omega = 1; chr.rhopw = &(rhopw); diff --git a/source/source_hamilt/module_xc/test/xc3_mock.h b/source/source_hamilt/module_xc/test/xc3_mock.h index afa1017ee5..57625c8c52 100644 --- a/source/source_hamilt/module_xc/test/xc3_mock.h +++ b/source/source_hamilt/module_xc/test/xc3_mock.h @@ -192,7 +192,7 @@ SepPot::~SepPot(){} Sep_Cell::Sep_Cell() noexcept {} Sep_Cell::~Sep_Cell() noexcept {} -namespace unitcell +namespace elecstate { void cal_ux(UnitCell& ucell, const int nspin) { diff --git a/source/source_hamilt/module_xc/xc_functional.h b/source/source_hamilt/module_xc/xc_functional.h index 540469b81e..bbf3a4094b 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); @@ -249,15 +250,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_gga_noncol_sf_builtin.cpp b/source/source_hamilt/module_xc/xc_functional_gga_noncol_sf_builtin.cpp new file mode 100644 index 0000000000..39c3109575 --- /dev/null +++ b/source/source_hamilt/module_xc/xc_functional_gga_noncol_sf_builtin.cpp @@ -0,0 +1,376 @@ +#include "xc_functional_gga_noncol_sf_builtin.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_gga_noncol_sf_builtin.h b/source/source_hamilt/module_xc/xc_functional_gga_noncol_sf_builtin.h new file mode 100644 index 0000000000..b763a6a47a --- /dev/null +++ b/source/source_hamilt/module_xc/xc_functional_gga_noncol_sf_builtin.h @@ -0,0 +1,32 @@ +#ifndef XC_FUNCTIONAL_GGA_NONCOL_SF_BUILTIN_H +#define XC_FUNCTIONAL_GGA_NONCOL_SF_BUILTIN_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_pot.cpp b/source/source_hamilt/module_xc/xc_pot.cpp index 1f8a3dcd32..da9eae7e77 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_gga_noncol_sf_builtin.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) 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..a31487eafe 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, identical to the behavior before this keyword was introduced; the default. +* 1: collinear approximation, grad(rho_up/dn) = (grad(rho) +/- grad(|m|))/2; only the gradient of the magnetization magnitude is used (equivalent to 0). +* 2: projected method. The gradients use the full chain rule, grad(rho_up/dn) = (grad(rho) +/- m_hat . grad(m))/2, but the divergence of h = df/d(grad rho) in the potential is projected onto m_hat: v_mu -= m_hat_mu * div((h_up - h_dn)/2), dropping the (h_up - h_dn) . grad(m_hat_mu) cross terms. +* 3: Scalmani-Frisch transformation (G. Scalmani and M. J. Frisch, J. Chem. Theory Comput. 8, 2193 (2012)). Same gradients as 2, but the full divergence is kept: v_mu -= div((h_up - h_dn)/2 * m_hat_mu), retaining all cross terms; the most accurate. +This parameter only takes effect for nspin=4 with GGA functionals (and magnetic calculation).)"; + 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?