diff --git a/source/source_hamilt/module_xc/libxc_abacus.h b/source/source_hamilt/module_xc/libxc_abacus.h index 0c912858bf..7a692b2382 100644 --- a/source/source_hamilt/module_xc/libxc_abacus.h +++ b/source/source_hamilt/module_xc/libxc_abacus.h @@ -105,11 +105,44 @@ namespace XC_Functional_Libxc const double tpiba, const Charge* const chr); + extern void cal_gdr_and_lapl( + const int nspin, + const std::size_t nrxx, + const std::vector &rho, + const double tpiba, + const Charge* const chr, + std::vector>> &gdr, + std::vector &lapl); + // converting grho (abacus=>libxc) extern std::vector convert_sigma( const std::vector>> &gdr); - // sgn for threshold mask + /// Calculate Laplacian of density using spectral method: ∇²ρ = IFFT(-|G|²·FFT(ρ)) + /// @see cal_lapl() implementation in libxc_tools.cpp for full documentation + extern std::vector cal_lapl( + const int nspin, + const std::size_t nrxx, + const std::vector &rho, + const double tpiba, + const Charge* const chr); + + /// Calculate Laplacian of density using finite-difference kernel (bounded at high G) + /// @see cal_lapl_fd() implementation in libxc_tools.cpp for full documentation + extern std::vector cal_lapl_fd( + const int nspin, + const std::size_t nrxx, + const std::vector &rho, + const Charge* const chr); + + /// Calculate density Hessian: H_ab = ∂²ρ/∂r_a∂r_b (6 independent components per spin) + /// @see cal_rho_hessian() implementation in libxc_tools.cpp for full documentation + extern std::vector> cal_rho_hessian( + const int nspin, + const std::size_t nrxx, + const std::vector &rho, + const Charge* const chr); + extern std::vector cal_sgn( const double rho_threshold, const double grho_threshold, @@ -209,25 +242,55 @@ namespace XC_Functional_Libxc // libxc_mgga_wrap.cpp //------------------- - // wrapper for the mGGA functionals + /// Wrapper for meta-GGA functionals (single spin channel). + /// Computes XC energy density and potentials including Laplacian-dependent terms. + /// + /// @param func_id LibXC functional IDs (exchange + correlation) + /// @param rho Electron density at grid point + /// @param grho |∇ρ|² at grid point + /// @param lapl_rho ∇²ρ at grid point (density Laplacian; pass 0.0 for SCAN) + /// @param atau Kinetic energy density τ at grid point + /// @param sxc [out] XC energy density + /// @param v1xc [out] ∂ε_xc/∂ρ + /// @param v2xc [out] 2·∂ε_xc/∂(|∇ρ|²) + /// @param v3xc [out] ∂ε_xc/∂τ + /// @param vlapl [out] ∂ε_xc/∂(∇²ρ) (Laplacian potential; 0.0 for SCAN) + /// @param hybrid_alpha Exact exchange mixing fraction (0.0 for pure GGA/meta-GGA) + /// @param hse_omega Range-separation parameter for HSE-type functionals + /// + /// @note For SCAN (non-Laplacian meta-GGA), pass lapl_rho=0.0 and vlapl=0.0. + /// For SCANL, pass the actual ∇²ρ value and vlapl will be populated. + /// @note The vlapl output is used by process_vlapl_potential() to compute the + /// FD Laplacian kernel contribution to the XC potential. extern void tau_xc( const std::vector &func_id, const double &rho, const double &grho, + const double &lapl_rho, const double &atau, double &sxc, - double &v1xc, - double &v2xc, - double &v3xc, - const double &hybrid_alpha, - const double &hse_omega); - + double &v1xc, + double &v2xc, + double &v3xc, + double &vlapl, + const double &hybrid_alpha, + const double &hse_omega); + + /// Wrapper for meta-GGA functionals (spin-polarized, two channels). + /// Same as tau_xc but handles spin-up/spin-down channels separately. + /// + /// @param laplup, lapldw [in] ∇²ρ_up, ∇²ρ_dn (Laplacian for each spin channel) + /// @param vlaplup, vlapldw [out] ∂ε_xc/∂(∇²ρ_up), ∂ε_xc/∂(∇²ρ_dn) + /// + /// @note See tau_xc() for other parameter descriptions. extern void tau_xc_spin( const std::vector &func_id, double rhoup, double rhodw, ModuleBase::Vector3 gdr1, ModuleBase::Vector3 gdr2, + double laplup, + double lapldw, double tauup, double taudw, double &sxc, @@ -236,10 +299,12 @@ namespace XC_Functional_Libxc double &v2xcup, double &v2xcdw, double &v2xcud, - double &v3xcup, - double &v3xcdw, - const double &hybrid_alpha, - const double &hse_omega); + double &v3xcup, + double &v3xcdw, + double &vlaplup, + double &vlapldw, + const double &hybrid_alpha, + const double &hse_omega); } // namespace XC_Functional_Libxc diff --git a/source/source_hamilt/module_xc/libxc_mgga_wrap.cpp b/source/source_hamilt/module_xc/libxc_mgga_wrap.cpp index 2c7d700f43..e079863cb4 100644 --- a/source/source_hamilt/module_xc/libxc_mgga_wrap.cpp +++ b/source/source_hamilt/module_xc/libxc_mgga_wrap.cpp @@ -17,11 +17,13 @@ void XC_Functional_Libxc::tau_xc( const std::vector& func_id, const double& rho, const double& grho, + const double& lapl_rho, const double& atau, double& sxc, double& v1xc, double& v2xc, double& v3xc, + double& vlapl, const double& hybrid_alpha, const double& hse_omega) { @@ -29,7 +31,6 @@ void XC_Functional_Libxc::tau_xc( double v1 = 0.0; double v2 = 0.0; double v3 = 0.0; - double lapl_rho = grho; double vlapl_rho = 0.0; std::vector funcs = XC_Functional_Libxc::init_func( /* func_id = */ func_id, @@ -41,6 +42,7 @@ void XC_Functional_Libxc::tau_xc( v1xc = 0.0; v2xc = 0.0; v3xc = 0.0; + vlapl = 0.0; for (xc_func_type& func : funcs) { @@ -52,12 +54,14 @@ void XC_Functional_Libxc::tau_xc( v1 *= (1.0 - hybrid_alpha); v2 *= (1.0 - hybrid_alpha); v3 *= (1.0 - hybrid_alpha); + vlapl_rho *= (1.0 - hybrid_alpha); } #endif sxc += s * rho; v2xc += v2 * 2.0; v1xc += v1; v3xc += v3; + vlapl += vlapl_rho; } XC_Functional_Libxc::finish_func(funcs); @@ -71,6 +75,8 @@ void XC_Functional_Libxc::tau_xc_spin( double rhodw, ModuleBase::Vector3 gdr1, ModuleBase::Vector3 gdr2, + double laplup, + double lapldw, double tauup, double taudw, double& sxc, @@ -81,6 +87,8 @@ void XC_Functional_Libxc::tau_xc_spin( double& v2xcud, double& v3xcup, double& v3xcdw, + double& vlaplup, + double& vlapldw, const double& hybrid_alpha, const double& hse_omega) { @@ -92,10 +100,13 @@ void XC_Functional_Libxc::tau_xc_spin( v2xcud = 0.0; v3xcup = 0.0; v3xcdw = 0.0; + vlaplup = 0.0; + vlapldw = 0.0; const std::array rho = {rhoup, rhodw}; const std::array grho = {gdr1.norm2(), gdr1 * gdr2, gdr2.norm2()}; const std::array tau = {tauup, taudw}; + const std::array lapl = {laplup, lapldw}; std::vector funcs = XC_Functional_Libxc::init_func( /* func_id = */ func_id, @@ -125,12 +136,10 @@ void XC_Functional_Libxc::tau_xc_spin( double s = 0.0; std::array v1xc = {0.0, 0.0}; std::array v3xc = {0.0, 0.0}; - std::array lapl = {0.0, 0.0}; std::array vlapl = {0.0, 0.0}; std::array v2xc = {0.0, 0.0, 0.0}; - // call Libxc function: xc_mgga_exc_vxc - xc_mgga_exc_vxc(&func, 1, rho.data(), grho.data(), lapl.data(), tau.data(), &s, - v1xc.data(), v2xc.data(), vlapl.data(), v3xc.data()); + xc_mgga_exc_vxc(&func, 1, rho.data(), grho.data(), lapl.data(), tau.data(), &s, + v1xc.data(), v2xc.data(), vlapl.data(), v3xc.data()); #ifdef __EXX if (func.info->number == XC_MGGA_X_SCAN && XC_Functional::get_func_type() == 5) @@ -143,6 +152,8 @@ void XC_Functional_Libxc::tau_xc_spin( v2xc[2] *= (1.0 - hybrid_alpha); v3xc[0] *= (1.0 - hybrid_alpha); v3xc[1] *= (1.0 - hybrid_alpha); + vlapl[0] *= (1.0 - hybrid_alpha); + vlapl[1] *= (1.0 - hybrid_alpha); } #endif @@ -154,6 +165,8 @@ void XC_Functional_Libxc::tau_xc_spin( v2xcdw += 2.0 * v2xc[2] * sgn[1]; v3xcup += v3xc[0] * sgn[0]; v3xcdw += v3xc[1] * sgn[1]; + vlaplup += vlapl[0] * sgn[0]; + vlapldw += vlapl[1] * sgn[1]; } } diff --git a/source/source_hamilt/module_xc/libxc_pot.cpp b/source/source_hamilt/module_xc/libxc_pot.cpp index 9d70a41274..9f4f873929 100644 --- a/source/source_hamilt/module_xc/libxc_pot.cpp +++ b/source/source_hamilt/module_xc/libxc_pot.cpp @@ -15,6 +15,144 @@ #include #include +#include +#include + +namespace { +/// Process the vlapl (Laplacian potential) contribution to the XC potential and energy +/// for meta-GGA functionals that depend on the density Laplacian (e.g., SCANL). +/// +/// The vlapl potential contribution to the XC potential is computed using a finite-difference +/// (FD) Laplacian kernel in reciprocal space: +/// +/// V_lapl(r) = e² · ∇²_FD(vlapl(r) · sgn(r)) +/// +/// where: +/// - e² = ModuleBase::e2 is the electron charge squared (conversion factor) +/// - vlapl(r) is the Laplacian potential from the XC functional +/// - sgn(r) is the spin sign (+1 for spin-up, -1 for spin-down, 0 if below threshold) +/// - ∇²_FD is the FD Laplacian operator (not the spectral -|G|²) +/// +/// The FD Laplacian kernel is used instead of the spectral Laplacian (-|G|²) to avoid +/// amplification of high-G noise that causes SCF divergence. The FD kernel matches the +/// spectral kernel at low G but remains bounded at high G. +/// +/// The FD Laplacian kernel in reciprocal space is: +/// gg_FD(G) = (1/(2π)²) Σ_{αβ} GGT[α][β] · FD_{αβ}(G) +/// +/// where: +/// - GGT is the metric tensor G·G^T +/// - FD_{αα} = 2·N_α²·(1 - cos(2π·m_α/N_α)) (diagonal terms) +/// - FD_{αβ} = N_α·N_β·sin(2π·m_α/N_α)·sin(2π·m_β/N_β) (off-diagonal terms) +/// - m_α are the Miller indices, N_α are the grid dimensions +/// +/// The contribution to vtxc (XC energy density integral) is: +/// vtxc += Σ_r V_lapl(r) · ρ(r) · sgn(r) +/// +/// @note This function is only called for meta-GGA functionals (func_type == 3 or 5). +/// @note If vlapl_max <= 1e-20, the function returns early (SCAN functional returns vlapl=0). +/// @note For hybrid meta-GGA (SCAN, func_type==5), the vlapl contribution is scaled by +/// (1 - hybrid_alpha) to account for the exact exchange mixing. +/// @warning The FD kernel is an approximation. For very high ecutwfc (>80 Ry), the FD kernel +/// may not accurately represent the Laplacian at high G, but it remains stable. +/// Consider using r2SCAN instead of SCANL for better numerical stability. +void process_vlapl_potential( + const int nspin, + const int nrxx, + const double tpiba, + const std::vector& vlapl, + const std::vector& sgn, + const std::vector& rho, + ModuleBase::matrix& v, + double& vtxc, + const Charge* chr, + const int func_type, + const double hybrid_alpha) +{ + double vlapl_max = 0.0; + for (int i = 0; i < nrxx * nspin; ++i) + vlapl_max = std::max(vlapl_max, std::abs(vlapl[i])); + if (vlapl_max <= 1e-20) + return; + + const int nx = chr->rhopw->nx; + const int ny = chr->rhopw->ny; + const int nz = chr->rhopw->nz; + const double tpiba2 = tpiba * tpiba; + const double twopi = 2.0 * M_PI; + + std::vector gg_fd(chr->rhopw->npw); +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 1024) +#endif + for (int ig = 0; ig < chr->rhopw->npw; ig++) + { + int m0 = chr->rhopw->gdirect[ig].x; + int m1 = chr->rhopw->gdirect[ig].y; + int m2 = chr->rhopw->gdirect[ig].z; + + double fd00 = 2.0 * nx * nx * (1.0 - std::cos(twopi * m0 / nx)); + double fd11 = 2.0 * ny * ny * (1.0 - std::cos(twopi * m1 / ny)); + double fd22 = 2.0 * nz * nz * (1.0 - std::cos(twopi * m2 / nz)); + double fd01 = nx * ny * std::sin(twopi * m0 / nx) * std::sin(twopi * m1 / ny); + double fd02 = nx * nz * std::sin(twopi * m0 / nx) * std::sin(twopi * m2 / nz); + double fd12 = ny * nz * std::sin(twopi * m1 / ny) * std::sin(twopi * m2 / nz); + + double ggt00 = chr->rhopw->GGT.e11; + double ggt01 = chr->rhopw->GGT.e12; + double ggt02 = chr->rhopw->GGT.e13; + double ggt11 = chr->rhopw->GGT.e22; + double ggt12 = chr->rhopw->GGT.e23; + double ggt22 = chr->rhopw->GGT.e33; + + gg_fd[ig] = (ggt00 * fd00 + ggt11 * fd11 + ggt22 * fd22 + + 2.0 * ggt01 * fd01 + 2.0 * ggt02 * fd02 + 2.0 * ggt12 * fd12) + / (twopi * twopi); + } + + for (int is = 0; is < nspin; ++is) + { + std::vector vlapl_sgn(nrxx); +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 1024) +#endif + for (int ir = 0; ir < nrxx; ++ir) + { + vlapl_sgn[ir] = vlapl[ir * nspin + is] * sgn[ir * nspin + is]; + } + + std::vector> vlapl_g(chr->rhopw->npw); + chr->rhopw->real2recip(vlapl_sgn.data(), vlapl_g.data()); + +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 1024) +#endif + for (int ig = 0; ig < chr->rhopw->npw; ig++) + { + vlapl_g[ig] *= -gg_fd[ig] * tpiba2; + } + + std::vector> vlapl_lapl(chr->rhopw->nmaxgr); + chr->rhopw->recip2real(vlapl_g.data(), vlapl_lapl.data()); + + double rvtxc = 0.0; +#ifdef _OPENMP +#pragma omp parallel for reduction(+:rvtxc) schedule(static, 256) +#endif + for (int ir = 0; ir < nrxx; ++ir) + { + double vlapl_val = ModuleBase::e2 * vlapl_lapl[ir].real(); + if (func_type == 5) + { + vlapl_val *= (1.0 - hybrid_alpha); + } + v(is, ir) += vlapl_val; + rvtxc += vlapl_val * rho[ir * nspin + is] * sgn[ir * nspin + is]; + } + vtxc += rvtxc; + } +} +} // anonymous namespace std::tuple XC_Functional_Libxc::v_xc_libxc( // Peize Lin update for nspin==4 at 2023.01.14 const std::vector &func_id, @@ -246,8 +384,9 @@ std::tuple XC_Functional_Li /* hse_omega = */ hse_omega); const std::vector rho = XC_Functional_Libxc::convert_rho(nspin, nrxx, chr); - const std::vector>> gdr - = XC_Functional_Libxc::cal_gdr(nspin, nrxx, rho, tpiba, chr); + std::vector>> gdr; + std::vector lapl; + XC_Functional_Libxc::cal_gdr_and_lapl(nspin, nrxx, rho, tpiba, chr, gdr, lapl); const std::vector sigma = XC_Functional_Libxc::convert_sigma(gdr); //converting kin_r @@ -328,7 +467,7 @@ std::tuple XC_Functional_Li nrxx_thread, rho.data() + ir_start * nspin, sigma.data() + ir_start * ((1==nspin)?1:3), - sigma.data() + ir_start * ((1==nspin)?1:3), + lapl.data() + ir_start * nspin, kin_r.data() + ir_start * nspin, exc.data() + ir_start, vrho.data() + ir_start * nspin, @@ -458,6 +597,23 @@ std::tuple XC_Functional_Li vofk(is,ir) += vtau[ir*nspin+is] * sgn[ir*nspin+is]; } } + + //process vlapl: compute ∇²(vlapl·sgn) using FD Laplacian kernel + //The FD kernel matches the spectral kernel at low G but is bounded at high G, + //preventing |G|² amplification that causes SCF divergence. + // + // For SCANL and other Laplacian-dependent meta-GGAs, the vlapl potential + // (∂ε_xc/∂(∇²ρ)) contributes to the XC potential via: + // V_xc += e² · ∇²_FD(vlapl · sgn) + // vtxc += Σ_r V_lapl(r) · ρ(r) · sgn(r) + // + // The FD Laplacian kernel is used instead of the spectral Laplacian (-|G|²) + // to avoid amplification of high-G noise that causes SCF divergence. + // See process_vlapl_potential() for full documentation of the FD kernel. + process_vlapl_potential( + nspin, nrxx, tpiba, vlapl, sgn, rho, v, vtxc, chr, + XC_Functional::get_func_type(), + XC_Functional::get_hybrid_alpha()); } //------------------------------------------------- diff --git a/source/source_hamilt/module_xc/libxc_setup.cpp b/source/source_hamilt/module_xc/libxc_setup.cpp index 84972a8063..daf01d9fe6 100644 --- a/source/source_hamilt/module_xc/libxc_setup.cpp +++ b/source/source_hamilt/module_xc/libxc_setup.cpp @@ -180,6 +180,22 @@ XC_Functional_Libxc::set_xc_type_libxc(const std::string& xc_func_in) ModuleBase::WARNING_QUIT("XC_Functional::set_xc_type_libxc", message); } + { + std::vector tmp_funcs = XC_Functional_Libxc::init_func(func_id, XC_UNPOLARIZED, 0.0, 0.0); + for (auto& f : tmp_funcs) + { + if (f.info->flags & XC_FLAGS_NEEDS_LAPLACIAN) + { + std::cout << " WARNING: XC functional \"" << f.info->name + << "\" requires Laplacian of density (nabla^2 rho)." + << " This may require a higher energy cutoff for numerical stability." + << std::endl; + break; + } + } + XC_Functional_Libxc::finish_func(tmp_funcs); + } + // return return std::make_pair(func_type, func_id); } diff --git a/source/source_hamilt/module_xc/libxc_tools.cpp b/source/source_hamilt/module_xc/libxc_tools.cpp index b66e1b1b26..182874d3e7 100644 --- a/source/source_hamilt/module_xc/libxc_tools.cpp +++ b/source/source_hamilt/module_xc/libxc_tools.cpp @@ -88,6 +88,45 @@ XC_Functional_Libxc::cal_gdr( return gdr; } +void XC_Functional_Libxc::cal_gdr_and_lapl( + const int nspin, + const std::size_t nrxx, + const std::vector &rho, + const double tpiba, + const Charge* const chr, + std::vector>> &gdr, + std::vector &lapl) +{ + gdr.resize(nspin); + lapl.assign(nrxx * nspin, 0.0); + for( int is=0; is!=nspin; ++is ) + { + std::vector rhor(nrxx); + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for(std::size_t ir=0; ir> rhog(chr->rhopw->npw); + chr->rhopw->real2recip(rhor.data(), rhog.data()); + + gdr[is].resize(nrxx); + XC_Functional::grad_rho(rhog.data(), gdr[is].data(), chr->rhopw, tpiba); + + std::vector lapl_spin(nrxx); + XC_Functional::laplacian_rho(rhog.data(), lapl_spin.data(), chr->rhopw, tpiba); + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for(std::size_t ir=0; irlibxc) std::vector XC_Functional_Libxc::convert_sigma( const std::vector>> &gdr) @@ -124,6 +163,231 @@ std::vector XC_Functional_Libxc::convert_sigma( return sigma; } +/// Calculate the Laplacian of the electron density for each spin channel. +/// Uses the spectral (reciprocal-space) Laplacian: ∇²ρ(r) = IFFT(-|G|² · FFT(ρ)) +/// +/// @param nspin Number of spin channels +/// @param nrxx Number of real-space grid points +/// @param rho Electron density (interleaved: [ρ_up(r0), ρ_dn(r0), ρ_up(r1), ...]) +/// @param tpiba 2π/a (reciprocal lattice scaling factor) +/// @param chr Charge density object (provides PW basis and FFT infrastructure) +/// @return Laplacian of density (same layout as rho) +/// +/// @note This computes the exact spectral Laplacian, which is used for the energy evaluation. +/// For the vlapl potential, use cal_lapl_fd() instead to avoid SCF divergence. +std::vector XC_Functional_Libxc::cal_lapl( + const int nspin, + const std::size_t nrxx, + const std::vector &rho, + const double tpiba, + const Charge* const chr) +{ + std::vector lapl(nrxx * nspin); + for (int is = 0; is != nspin; ++is) + { + std::vector rhor(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 + is]; + } + std::vector> rhog(chr->rhopw->npw); + chr->rhopw->real2recip(rhor.data(), rhog.data()); + + std::vector lapl_spin(nrxx); + XC_Functional::laplacian_rho(rhog.data(), lapl_spin.data(), chr->rhopw, tpiba); + + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + lapl[ir * nspin + is] = lapl_spin[ir]; + } + } + return lapl; +} + +/// Calculate the Laplacian of the electron density using a finite-difference (FD) kernel. +/// This is consistent with the FD Laplacian kernel used for vlapl potential processing. +/// +/// The FD Laplacian in reciprocal space is: +/// ∇²_FD ρ(G) = -gg_FD(G) · tpiba² · ρ(G) +/// +/// where gg_FD(G) is the FD kernel: +/// gg_FD(G) = (1/(2π)²) Σ_{αβ} GGT[α][β] · FD_{αβ}(G) +/// +/// with: +/// FD_{αα} = 2·N_α²·(1 - cos(2π·m_α/N_α)) (diagonal) +/// FD_{αβ} = N_α·N_β·sin(2π·m_α/N_α)·sin(2π·m_β/N_β) (off-diagonal) +/// +/// @note The FD kernel matches the spectral kernel at low G but remains bounded at high G, +/// preventing |G|² amplification that causes SCF divergence. +/// @note This function is currently not used in the main code path (the vlapl potential +/// is computed inline in process_vlapl_potential()). It is provided for consistency +/// and potential future use. +std::vector XC_Functional_Libxc::cal_lapl_fd( + const int nspin, + const std::size_t nrxx, + const std::vector &rho, + const Charge* const chr) +{ + const int nx = chr->rhopw->nx; + const int ny = chr->rhopw->ny; + const int nz = chr->rhopw->nz; + const double tpiba = chr->rhopw->tpiba; + const double tpiba2 = tpiba * tpiba; + const double twopi = 2.0 * M_PI; + + std::vector gg_fd(chr->rhopw->npw); + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for (int ig = 0; ig < chr->rhopw->npw; ++ig) + { + int m0 = chr->rhopw->gdirect[ig].x; + int m1 = chr->rhopw->gdirect[ig].y; + int m2 = chr->rhopw->gdirect[ig].z; + + double fd00 = 2.0 * nx * nx * (1.0 - std::cos(twopi * m0 / nx)); + double fd11 = 2.0 * ny * ny * (1.0 - std::cos(twopi * m1 / ny)); + double fd22 = 2.0 * nz * nz * (1.0 - std::cos(twopi * m2 / nz)); + double fd01 = nx * ny * std::sin(twopi * m0 / nx) * std::sin(twopi * m1 / ny); + double fd02 = nx * nz * std::sin(twopi * m0 / nx) * std::sin(twopi * m2 / nz); + double fd12 = ny * nz * std::sin(twopi * m1 / ny) * std::sin(twopi * m2 / nz); + + double ggt00 = chr->rhopw->GGT.e11; + double ggt01 = chr->rhopw->GGT.e12; + double ggt02 = chr->rhopw->GGT.e13; + double ggt11 = chr->rhopw->GGT.e22; + double ggt12 = chr->rhopw->GGT.e23; + double ggt22 = chr->rhopw->GGT.e33; + + gg_fd[ig] = (ggt00 * fd00 + ggt11 * fd11 + ggt22 * fd22 + + 2.0 * ggt01 * fd01 + 2.0 * ggt02 * fd02 + 2.0 * ggt12 * fd12) + / (twopi * twopi); + } + + std::vector lapl(nrxx * nspin); + for (int is = 0; is != nspin; ++is) + { + std::vector rhor(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 + is]; + } + std::vector> rhog(chr->rhopw->npw); + chr->rhopw->real2recip(rhor.data(), rhog.data()); + + std::vector> lapl_g(chr->rhopw->nmaxgr); + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for (int ig = 0; ig < chr->rhopw->npw; ++ig) + { + lapl_g[ig] = -rhog[ig] * gg_fd[ig] * tpiba2; + } + for (int ig = chr->rhopw->npw; ig < chr->rhopw->nmaxgr; ++ig) + { + lapl_g[ig] = 0.0; + } + + std::vector> lapl_r(chr->rhopw->nmaxgr); + chr->rhopw->recip2real(lapl_g.data(), lapl_r.data()); + + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + lapl[ir * nspin + is] = lapl_r[ir].real(); + } + } + return lapl; +} + +/// Calculate the Hessian of the electron density: H_ab(r) = ∂²ρ/∂r_a∂r_b +/// Returns 6 independent components per spin channel: xx, yy, zz, xy, yz, zx +/// +/// The Hessian is computed in reciprocal space and transformed to real space: +/// H_ab(G) = -G_a · G_b · ρ(G) +/// H_ab(r) = IFFT(H_ab(G)) · tpiba² +/// +/// where G_a, G_b are Cartesian components of the reciprocal lattice vector G. +/// +/// @param nspin Number of spin channels +/// @param nrxx Number of real-space grid points +/// @param rho Electron density (interleaved: [ρ_up(r0), ρ_dn(r0), ...]) +/// @param chr Charge density object (provides PW basis and FFT infrastructure) +/// @return Hessian components: hessian[is][ab * nrxx + ir] +/// where ab = 0:xx, 1:yy, 2:zz, 3:xy, 4:yz, 5:zx +/// +/// @note This is used by add_vlapl_stress_contribution() to compute the vlapl stress: +/// σ_ab^{vlapl} = (2/Ω) Σ_{s,ir} vlapl_s(r) · H_ab,s(r) · e² +/// @note The Hessian is computed using the spectral (exact) Laplacian, not the FD kernel. +/// This is appropriate for the stress calculation because the density is smooth +/// and well-resolved, and the spectral method is more accurate. +std::vector> XC_Functional_Libxc::cal_rho_hessian( + const int nspin, + const std::size_t nrxx, + const std::vector &rho, + const Charge* const chr) +{ + // 6 independent components: (0)xx (1)yy (2)zz (3)xy (4)yz (5)zx + std::vector> hessian(nspin, std::vector(6 * nrxx, 0.0)); + + for (int is = 0; is != nspin; ++is) + { + std::vector rhor(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 + is]; + } + std::vector> rhog(chr->rhopw->npw); + chr->rhopw->real2recip(rhor.data(), rhog.data()); + + // Hessian in G-space: H_ab(G) = -G_a * G_b * rho(G) + // We compute each of the 6 independent components + int ab_map[6][2] = {{0,0}, {1,1}, {2,2}, {0,1}, {1,2}, {0,2}}; + + for (int ab = 0; ab < 6; ++ab) + { + int a = ab_map[ab][0]; + int b = ab_map[ab][1]; + std::vector> hess_g(chr->rhopw->npw); + + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for (int ig = 0; ig < chr->rhopw->npw; ++ig) + { + hess_g[ig] = -rhog[ig] * chr->rhopw->gcar[ig][a] * chr->rhopw->gcar[ig][b]; + } + + std::vector> hess_r(chr->rhopw->nmaxgr); + chr->rhopw->recip2real(hess_g.data(), hess_r.data()); + + const double tpiba2 = chr->rhopw->tpiba * chr->rhopw->tpiba; + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + hessian[is][ab * nrxx + ir] = hess_r[ir].real() * tpiba2; + } + } + } + return hessian; +} + // sgn for threshold mask std::vector XC_Functional_Libxc::cal_sgn( const double rho_threshold, diff --git a/source/source_hamilt/module_xc/test/CMakeLists.txt b/source/source_hamilt/module_xc/test/CMakeLists.txt index 34634eae02..3b74511bad 100644 --- a/source/source_hamilt/module_xc/test/CMakeLists.txt +++ b/source/source_hamilt/module_xc/test/CMakeLists.txt @@ -34,6 +34,7 @@ AddTest( ../libxc_lda_wrap.cpp ../libxc_gga_wrap.cpp ../libxc_mgga_wrap.cpp + ../libxc_tools.cpp ../xc_gga_corr.cpp ../xc_lda_corr.cpp ../xc_gga_exch.cpp ../xc_lda_exch.cpp ../xc_hcth.cpp ../../../source_base/matrix.cpp @@ -48,15 +49,29 @@ AddTest( AddTest( TARGET MODULE_HAMILT_XCTest_SCAN - LIBS parameter MPI::MPI_CXX Libxc::xc + LIBS parameter MPI::MPI_CXX Libxc::xc ${math_libs} psi device container SOURCES test_xc4.cpp ../xc_functional.cpp ../xc_lda_wrap.cpp ../xc_gga_wrap.cpp ../libxc_setup.cpp ../libxc_lda_wrap.cpp ../libxc_gga_wrap.cpp ../libxc_mgga_wrap.cpp - ../xc_gga_corr.cpp ../xc_lda_corr.cpp - ../xc_gga_exch.cpp ../xc_lda_exch.cpp ../xc_hcth.cpp + ../libxc_tools.cpp + ../xc_grad.cpp + ../xc_pot.cpp + ../libxc_pot.cpp + ../xc_gga_corr.cpp ../xc_lda_corr.cpp + ../xc_gga_exch.cpp ../xc_lda_exch.cpp ../xc_hcth.cpp + ../../../source_base/matrix.cpp + ../../../source_base/memory_recorder.cpp + ../../../source_base/libm/branred.cpp + ../../../source_base/libm/sincos.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 + ../../../source_base/module_fft/fft_bundle.cpp + ../../../source_base/module_fft/fft_cpu.cpp + ${FFT_SRC} + ../../../source_base/timer.cpp + xc_reduce_mock.cpp ) AddTest( @@ -83,3 +98,44 @@ AddTest( ../../../source_base/module_fft/fft_cpu.cpp ${FFT_SRC} ) + +AddTest( + TARGET MODULE_HAMILT_XCTest_SCANL_LAPL + LIBS parameter MPI::MPI_CXX Libxc::xc + SOURCES test_xc6.cpp ../xc_functional.cpp ../xc_lda_wrap.cpp + ../xc_gga_wrap.cpp + ../libxc_setup.cpp + ../libxc_lda_wrap.cpp + ../libxc_gga_wrap.cpp + ../libxc_mgga_wrap.cpp + ../xc_gga_corr.cpp ../xc_lda_corr.cpp + ../xc_gga_exch.cpp ../xc_lda_exch.cpp ../xc_hcth.cpp +) + +if (USE_CUDA) +list(APPEND FFT_SRC ../../../source_base/module_fft/fft_cuda.cpp) +endif() +if (USE_ROCM) +list(APPEND FFT_SRC ../../../source_base/module_fft/fft_rocm.cpp) +endif() +AddTest( + TARGET MODULE_HAMILT_XCTest_LAPL + LIBS parameter MPI::MPI_CXX Libxc::xc psi device container + SOURCES test_xc7.cpp ../xc_grad.cpp ../xc_functional.cpp + ../xc_lda_wrap.cpp ../xc_gga_wrap.cpp + ../libxc_setup.cpp + ../libxc_lda_wrap.cpp + ../libxc_gga_wrap.cpp + ../libxc_mgga_wrap.cpp + ../libxc_tools.cpp + ../xc_gga_corr.cpp ../xc_lda_corr.cpp ../xc_gga_exch.cpp + ../xc_lda_exch.cpp ../xc_hcth.cpp + ../../../source_base/matrix.cpp + ../../../source_base/memory_recorder.cpp + ../../../source_base/libm/branred.cpp + ../../../source_base/libm/sincos.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 + ../../../source_base/module_fft/fft_bundle.cpp + ../../../source_base/module_fft/fft_cpu.cpp + ${FFT_SRC} +) diff --git a/source/source_hamilt/module_xc/test/test_xc4.cpp b/source/source_hamilt/module_xc/test/test_xc4.cpp index e1ef9a82c0..ba7afd9945 100644 --- a/source/source_hamilt/module_xc/test/test_xc4.cpp +++ b/source/source_hamilt/module_xc/test/test_xc4.cpp @@ -3,33 +3,8 @@ #include "gtest/gtest.h" #include "xctest.h" #include "../exx_info.h" +#include "xc3_mock.h" -/************************************************ -* unit test of functionals -***********************************************/ - -// For more information of the functions, check the comment of xc_functional.h -// mGGA Functional (SCAN) is tested under nspin = 1 - -namespace ModuleBase -{ - void WARNING_QUIT(const std::string &file,const std::string &description) {exit(1);} - void TITLE(const std::string &class_function_name,bool disable){}; - void TITLE(const std::string &class_name,const std::string &function_name,bool disable){}; -} - -namespace GlobalV -{ - std::string BASIS_TYPE = ""; - bool CAL_STRESS = false; - int CAL_FORCE = 0; - int NSPIN = 1; -} - -namespace GlobalC -{ - Exx_Info exx_info; -} class XCTest_SCAN : public XCTest { @@ -44,11 +19,12 @@ class XCTest_SCAN : public XCTest std::vector tau = {0,0.02403590412,0.01672229351,0.01340429824,0.01141731056}; for(int i=0;i<5;i++) - { - double e,v,v1,v2,v3; - double hybrid_alpha = 0.0; - double hse_omega = 0.0; - XC_Functional_Libxc::tau_xc(XC_Functional::get_func_id(), rho[i],grho[i],tau[i],e,v1,v2,v3,hybrid_alpha, hse_omega); + { + double e,v,v1,v2,v3,v4; + double hybrid_alpha = 0.0; + double hse_omega = 0.0; + double lapl_rho = 0.0; + XC_Functional_Libxc::tau_xc(XC_Functional::get_func_id(), rho[i],grho[i],lapl_rho,tau[i],e,v1,v2,v3,v4,hybrid_alpha, hse_omega); e_.push_back(e); v1_.push_back(v1); v2_.push_back(v2); diff --git a/source/source_hamilt/module_xc/test/test_xc6.cpp b/source/source_hamilt/module_xc/test/test_xc6.cpp new file mode 100644 index 0000000000..1168861a51 --- /dev/null +++ b/source/source_hamilt/module_xc/test/test_xc6.cpp @@ -0,0 +1,155 @@ +#include "../xc_functional.h" +#include "../libxc_abacus.h" +#include "gtest/gtest.h" +#include "xctest.h" +#include "../exx_info.h" +#include +#include +#include + +namespace ModuleBase +{ + void WARNING_QUIT(const std::string &file,const std::string &description) {exit(1);} + void TITLE(const std::string &class_function_name,bool disable){}; + void TITLE(const std::string &class_name,const std::string &function_name,bool disable){}; +} + +namespace GlobalV +{ + std::string BASIS_TYPE = ""; + bool CAL_STRESS = false; + int CAL_FORCE = 0; + int NSPIN = 1; +} + +namespace GlobalC +{ + Exx_Info exx_info; +} + +class XCTest_SCANL_Laplacian : public XCTest +{ + protected: + double e_base, v1_base, v2_base, v3_base, vlapl_base; + double e_modified, v1_modified, v2_modified, v3_modified, vlapl_modified; + double e_scaled, v1_scaled, v2_scaled, v3_scaled, vlapl_scaled; + + void SetUp() + { + XC_Functional::set_xc_type("MGGA_X_SCANL+MGGA_C_SCANL"); + + const double rho = 0.17E+01; + const double grho = 0.81E-11; + const double tau = 0.02403590412; + const double lapl_base = 0.15E+01; + double hybrid_alpha = 0.0; + double hse_omega = 0.0; + + XC_Functional_Libxc::tau_xc( + XC_Functional::get_func_id(), + rho, grho, lapl_base, tau, + e_base, v1_base, v2_base, v3_base, vlapl_base, hybrid_alpha, hse_omega + ); + + XC_Functional_Libxc::tau_xc( + XC_Functional::get_func_id(), + rho, grho, lapl_base + 1.0, tau, + e_modified, v1_modified, v2_modified, v3_modified, vlapl_modified, hybrid_alpha, hse_omega + ); + + XC_Functional_Libxc::tau_xc( + XC_Functional::get_func_id(), + rho, grho, 2.0 * lapl_base, tau, + e_scaled, v1_scaled, v2_scaled, v3_scaled, vlapl_scaled, hybrid_alpha, hse_omega + ); + } +}; + +TEST_F(XCTest_SCANL_Laplacian, laplacian_affects_energy) +{ + EXPECT_NE(e_base, e_modified); + EXPECT_NE(e_base, e_scaled); + + std::cout << std::scientific << std::setprecision(15); + std::cout << "\n=== SCAN-L Laplacian Sensitivity Test ===" << std::endl; + std::cout << "Base Laplacian: " << 0.15E+01 << std::endl; + std::cout << " E_xc = " << e_base << std::endl; + std::cout << " dE/dlapl ~= " << (e_modified - e_base) / 1.0 << std::endl; + std::cout << "Modified Laplacian:" << 0.15E+01 + 1.0 << std::endl; + std::cout << " E_xc = " << e_modified << std::endl; + std::cout << " Delta E = " << e_modified - e_base << std::endl; + std::cout << "Scaled Laplacian: " << 2.0 * 0.15E+01 << std::endl; + std::cout << " E_xc = " << e_scaled << std::endl; + std::cout << " Delta E = " << e_scaled - e_base << std::endl; + std::cout << "=========================================" << std::endl; +} + +TEST(XC_ScanL_Reference, MGGA_X_SCANL_direct_libxc) +{ + const double rho = 35.536521214608185; + const double sigma = 1.149382202334535e+05; + const double lapl = 8.411855859277239e+02; + const double tau = 1.389887953757970e+02; + + std::vector func_id = {XC_MGGA_X_SCANL}; + std::vector funcs = XC_Functional_Libxc::init_func(func_id, XC_UNPOLARIZED, 0.0, 0.0); + + double exc = 0.0, vrho = 0.0, vsigma = 0.0, vlapl = 0.0, vtau = 0.0; + xc_mgga_exc_vxc(&funcs[0], 1, &rho, &sigma, &lapl, &tau, + &exc, &vrho, &vsigma, &vlapl, &vtau); + + EXPECT_NEAR(exc, -2.569101943337681e+00, 5.0e-04); + EXPECT_NEAR(vrho, -2.912514021641506e+00, 1.0e-03); + EXPECT_NEAR(vsigma, -8.289754258579972e-05, 1.0e-12); + EXPECT_NEAR(vlapl, 3.971745124876924e-03, 1.0e-12); + EXPECT_EQ(vtau, 0.0); + + XC_Functional_Libxc::finish_func(funcs); +} + +TEST(XC_ScanL_Reference, MGGA_C_SCANL_direct_libxc) +{ + const double rho = 35.536521214608185; + const double sigma = 1.149382202334535e+05; + const double lapl = 8.411855859277239e+02; + const double tau = 1.389887953757970e+02; + + std::vector func_id = {XC_MGGA_C_SCANL}; + std::vector funcs = XC_Functional_Libxc::init_func(func_id, XC_UNPOLARIZED, 0.0, 0.0); + + double exc = 0.0, vrho = 0.0, vsigma = 0.0, vlapl = 0.0, vtau = 0.0; + xc_mgga_exc_vxc(&funcs[0], 1, &rho, &sigma, &lapl, &tau, + &exc, &vrho, &vsigma, &vlapl, &vtau); + + EXPECT_NEAR(exc, -5.250261832501450e-02, 1.0e-06); + EXPECT_NEAR(vrho, -1.323740339176898e-01, 1.0e-05); + EXPECT_NEAR(vsigma, 1.138021340988220e-05, 1.0e-10); + EXPECT_NEAR(vlapl, -3.867564402989276e-04, 1.0e-10); + EXPECT_EQ(vtau, 0.0); + + XC_Functional_Libxc::finish_func(funcs); +} + +TEST(XC_ScanL_Reference, tau_xc_wrapper_libxc) +{ + XC_Functional::set_xc_type("MGGA_X_SCANL+MGGA_C_SCANL"); + + const double rho = 35.536521214608185; + const double grho = 1.149382202334535e+05; + const double lapl = 8.411855859277239e+02; + const double tau = 1.389887953757970e+02; + + double sxc, v1xc, v2xc, v3xc, vlaplxc; + double hybrid_alpha = 0.0; + double hse_omega = 0.0; + XC_Functional_Libxc::tau_xc(XC_Functional::get_func_id(), rho, grho, lapl, tau, + sxc, v1xc, v2xc, v3xc, vlaplxc, hybrid_alpha, hse_omega); + + double exc_x_ref = -2.569101943337681e+00; + double exc_c_ref = -5.250261832501450e-02; + double sxc_ref = (exc_x_ref + exc_c_ref) * rho; + + EXPECT_NEAR(sxc, sxc_ref, 0.05); + EXPECT_NE(v1xc, 0.0); + EXPECT_NE(vlaplxc, 0.0); +} diff --git a/source/source_hamilt/module_xc/test/test_xc7.cpp b/source/source_hamilt/module_xc/test/test_xc7.cpp new file mode 100644 index 0000000000..300090e3e7 --- /dev/null +++ b/source/source_hamilt/module_xc/test/test_xc7.cpp @@ -0,0 +1,98 @@ +#include "gtest/gtest.h" +#include "xctest.h" +#include "../xc_functional.h" +#include "../exx_info.h" +#include "xc3_mock.h" +#include "source_base/matrix.h" + +class XCTest_LaplacianAnalytical : public XCTest +{ + protected: + ModulePW::PW_Basis rhopw; + const int npw = 5; + const int nrxx = 5; + const int nmaxgr = 5; + const double tpiba = 1.0; + + void SetUp() override + { + rhopw.nrxx = nrxx; + rhopw.npw = npw; + rhopw.nmaxgr = nmaxgr; + rhopw.gcar = new ModuleBase::Vector3[npw]; + for (int ig = 0; ig < npw; ig++) + rhopw.gcar[ig] = ModuleBase::Vector3(1.0, 1.0, 1.0); + } + + void TearDown() override + { + delete[] rhopw.gcar; + } +}; + +TEST_F(XCTest_LaplacianAnalytical, zero_input) +{ + std::vector> rhog(npw, 0.0); + std::vector lapl(nrxx, -1.0); + XC_Functional::laplacian_rho(rhog.data(), lapl.data(), &rhopw, tpiba); + for (int ir = 0; ir < nrxx; ir++) + EXPECT_DOUBLE_EQ(lapl[ir], 0.0); +} + +TEST_F(XCTest_LaplacianAnalytical, imaginary_rhog_constant) +{ + std::vector> rhog(npw); + for (int ig = 0; ig < npw; ig++) + rhog[ig] = std::complex(0.0, 1.0); + std::vector lapl(nrxx); + XC_Functional::laplacian_rho(rhog.data(), lapl.data(), &rhopw, tpiba); + double g2 = 3.0; + double expected = -g2 * tpiba * tpiba; + for (int ir = 0; ir < nrxx; ir++) + EXPECT_NEAR(lapl[ir], expected, 1e-14); +} + +TEST_F(XCTest_LaplacianAnalytical, linearity) +{ + std::vector> rhog_a(npw); + std::vector> rhog_b(npw); + std::vector> rhog_sum(npw); + for (int ig = 0; ig < npw; ig++) + { + rhog_a[ig] = std::complex(ig, 2.0 * ig); + rhog_b[ig] = std::complex(3.0 * ig, ig); + rhog_sum[ig] = rhog_a[ig] + rhog_b[ig]; + } + + std::vector lapl_a(nrxx), lapl_b(nrxx), lapl_sum(nrxx); + XC_Functional::laplacian_rho(rhog_a.data(), lapl_a.data(), &rhopw, tpiba); + XC_Functional::laplacian_rho(rhog_b.data(), lapl_b.data(), &rhopw, tpiba); + XC_Functional::laplacian_rho(rhog_sum.data(), lapl_sum.data(), &rhopw, tpiba); + + for (int ir = 0; ir < nrxx; ir++) + EXPECT_NEAR(lapl_sum[ir], lapl_a[ir] + lapl_b[ir], 1e-14); +} + +TEST_F(XCTest_LaplacianAnalytical, single_plane_wave) +{ + int n = 5; + rhopw.nrxx = n; + rhopw.npw = n; + rhopw.nmaxgr = n; + delete[] rhopw.gcar; + rhopw.gcar = new ModuleBase::Vector3[n]; + for (int ig = 0; ig < n; ig++) + rhopw.gcar[ig] = ModuleBase::Vector3(static_cast(ig), 0.0, 0.0); + + std::vector> rhog(n, 0.0); + rhog[0] = std::complex(0.0, 1.0); + std::vector lapl(n); + XC_Functional::laplacian_rho(rhog.data(), lapl.data(), &rhopw, tpiba); + + double g2 = 0.0; + for (int i = 0; i < 3; i++) + g2 += rhopw.gcar[0][i] * rhopw.gcar[0][i]; + double expected = -g2 * tpiba * tpiba; + for (int ir = 0; ir < n; ir++) + EXPECT_NEAR(lapl[ir], expected, 1e-14); +} diff --git a/source/source_hamilt/module_xc/test/xc3_mock.h b/source/source_hamilt/module_xc/test/xc3_mock.h index 57625c8c52..5e7d7ba748 100644 --- a/source/source_hamilt/module_xc/test/xc3_mock.h +++ b/source/source_hamilt/module_xc/test/xc3_mock.h @@ -229,9 +229,17 @@ namespace Parallel_Reduce #endif return; } + template<> + void Parallel_Reduce::reduce_pool(float& object) + { + #ifdef __MPI + float swap = object; + MPI_Allreduce(&swap , &object , 1, MPI_FLOAT , MPI_SUM , MPI_COMM_WORLD); + #endif + return; + } template void reduce_all(double& object); template void reduce_all(double* object, const int n); - template void reduce_pool(float& object); template void reduce_pool(float* object, const int n); template void reduce_pool(double* object, const int n); } diff --git a/source/source_hamilt/module_xc/test/xc_reduce_mock.cpp b/source/source_hamilt/module_xc/test/xc_reduce_mock.cpp new file mode 100644 index 0000000000..22902e3907 --- /dev/null +++ b/source/source_hamilt/module_xc/test/xc_reduce_mock.cpp @@ -0,0 +1,13 @@ +#include "../../../source_base/parallel_reduce.h" + +namespace Parallel_Reduce +{ +template<> +void reduce_all(double& object) +{ +#ifdef __MPI + double swap = object; + MPI_Allreduce(&swap, &object, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); +#endif +} +} diff --git a/source/source_hamilt/module_xc/xc_functional.cpp b/source/source_hamilt/module_xc/xc_functional.cpp index bef9d91b2b..4ee3a6b6f9 100644 --- a/source/source_hamilt/module_xc/xc_functional.cpp +++ b/source/source_hamilt/module_xc/xc_functional.cpp @@ -169,6 +169,31 @@ void XC_Functional::set_xc_type(const std::string xc_func_in) func_type = 5; use_libxc = true; } + // SCANL: SCAN + Laplacian (meta-GGA functional that depends on ∇²ρ) + // Unlike SCAN, SCANL depends on the density Laplacian, which introduces: + // 1. Additional potential term: vlapl = ∂ε_xc/∂(∇²ρ) + // 2. Additional stress contribution: σ_ab^{vlapl} = (2/Ω) Σ vlapl·H_ab·e² + // + // The vlapl potential is processed using a finite-difference (FD) Laplacian kernel + // in reciprocal space to avoid |G|² amplification that causes SCF divergence. + // See process_vlapl_potential() in libxc_pot.cpp for details. + // + // @warning SCANL is known to have convergence difficulties for metallic systems + // at high ecutwfc due to the Laplacian-dependent potential. + // Recommended settings: ecutwfc <= 80 Ry, mixing_tau = 1, Pulay mixing (beta 0.1-0.2). + // Consider using r2SCAN or SCAN for better numerical stability. + else if ( xc_func == "SCANL") + { + func_id.push_back(XC_MGGA_X_SCANL); + func_id.push_back(XC_MGGA_C_SCANL); + func_type = 3; + use_libxc = true; + std::cout << "\n WARNING: meta-GGA functional [scanl] depends on the density Laplacian." << std::endl; + std::cout << " The vlapl potential uses a finite-difference Laplacian kernel to avoid" << std::endl; + std::cout << " |G|^2 amplification in reciprocal space. For maximum numerical stability," << std::endl; + std::cout << " consider using r2SCAN or SCAN instead." << std::endl; + std::cout << " Recommended: ecutwfc <= 80 Ry, mixing_tau = 1, Pulay mixing (beta 0.1-0.2)\n" << std::endl; + } else if( xc_func == "LC_PBE") { func_id.push_back(XC_HYB_GGA_XC_LC_PBEOP); diff --git a/source/source_hamilt/module_xc/xc_functional.h b/source/source_hamilt/module_xc/xc_functional.h index ceb8135399..fcd06e0acd 100644 --- a/source/source_hamilt/module_xc/xc_functional.h +++ b/source/source_hamilt/module_xc/xc_functional.h @@ -250,6 +250,23 @@ class XC_Functional const ModulePW::PW_Basis* rho_basis, const double tpiba); + /// Compute the Laplacian of the electron density in real space. + /// Uses the spectral (reciprocal-space) Laplacian: ∇²ρ(r) = IFFT(-|G|² · FFT(ρ)) + /// + /// @param rhog Density in reciprocal space (input) + /// @param lapl [out] ∇²ρ in real space + /// @param rho_basis PW basis for the density + /// @param tpiba 2π/a (reciprocal lattice scaling factor) + /// + /// @note This computes the exact spectral Laplacian, which is accurate for smooth densities + /// but amplifies high-G noise. For the vlapl potential processing, use the FD Laplacian + /// kernel (see process_vlapl_potential() in libxc_pot.cpp) instead. + static void laplacian_rho( + const std::complex* rhog, + double* lapl, + const ModulePW::PW_Basis* rho_basis, + const double tpiba); + static void noncolin_rho( double* rhoout1, double* rhoout2, diff --git a/source/source_hamilt/module_xc/xc_grad.cpp b/source/source_hamilt/module_xc/xc_grad.cpp index c16dcc718b..ed75921446 100644 --- a/source/source_hamilt/module_xc/xc_grad.cpp +++ b/source/source_hamilt/module_xc/xc_grad.cpp @@ -22,6 +22,152 @@ #ifdef __EXX #include "source_hamilt/module_xc/exx_info.h" #endif + +namespace { +/// Compute and add the vlapl (Laplacian potential) contribution to the GGA stress tensor +/// for meta-GGA functionals that depend on the density Laplacian (e.g., SCANL). +/// +/// The vlapl stress contribution arises from the strain derivative of the Laplacian-dependent +/// exchange-correlation energy. When the XC functional depends on ∇²ρ, the stress tensor +/// acquires an additional term: +/// +/// σ_ab^{vlapl} = (2/Ω) Σ_{s,ir} vlapl_s(r) · H_ab,s(r) · e² +/// +/// where: +/// - Ω is the cell volume +/// - vlapl_s(r) is the Laplacian potential for spin s at grid point r +/// - H_ab,s(r) = ∂²ρ_s/∂r_a∂r_b is the density Hessian (second spatial derivatives) +/// - e² = ModuleBase::e2 is the electron charge squared (conversion factor) +/// - The sum runs over spin channels s and real-space grid points r +/// +/// The density Hessian H_ab is computed in G-space as: +/// H_ab(G) = -G_a · G_b · ρ(G) +/// and transformed to real space via inverse FFT. +/// +/// @note This function is only called for meta-GGA functionals (func_type == 3 for pure meta-GGA, +/// or func_type == 5 for hybrid meta-GGA) that depend on the Laplacian. +/// @note The vlapl values are recomputed here by calling tau_xc/tau_xc_spin at each grid point. +/// This is necessary because the main gradcorr() loop does not store vlapl values. +/// If performance becomes an issue, consider caching vlapl values during the main loop. +void add_vlapl_stress_contribution( + std::vector& stress_gga, + const int nspin0, + const double* rhotmp1, + const double* rhotmp2, + const ModuleBase::Vector3* gdr1, + const ModuleBase::Vector3* gdr2, + const std::vector& lapl1, + const std::vector& lapl2, + const Charge* chr, + const ModulePW::PW_Basis* rhopw, + const UnitCell* ucell, + const double epsr, + const double hybrid_alpha, + const double hse_omega) +{ + const std::size_t nrxx_s = rhopw->nrxx; + std::vector rho_interleaved(nrxx_s * nspin0); + if (nspin0 == 1) + { +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 1024) +#endif + for (std::size_t ir = 0; ir < nrxx_s; ++ir) + { + rho_interleaved[ir] = std::abs(rhotmp1[ir]); + } + } + else + { +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 1024) +#endif + for (std::size_t ir = 0; ir < nrxx_s; ++ir) + { + rho_interleaved[ir * 2] = std::abs(rhotmp1[ir]); + rho_interleaved[ir * 2 + 1] = std::abs(rhotmp2[ir]); + } + } + auto hessian = XC_Functional_Libxc::cal_rho_hessian(nspin0, nrxx_s, rho_interleaved, chr); + + std::vector vlapl_vals(nrxx_s * nspin0, 0.0); + if (nspin0 == 1) + { +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 256) +#endif + for (std::size_t ir = 0; ir < nrxx_s; ++ir) + { + double arho = std::abs(rhotmp1[ir]); + double grho2 = gdr1[ir].norm2(); + double atau = chr->kin_r[0][ir] / 2.0; + double lapl_val = lapl1.empty() ? 0.0 : lapl1[ir]; + double sxc_l = 0.0, v1xc_l = 0.0, v2xc_l = 0.0, v3xc_l = 0.0, vlapl_l = 0.0; + double ha = hybrid_alpha; + double hse_om = hse_omega; + if (arho > epsr) + XC_Functional_Libxc::tau_xc(XC_Functional::get_func_id(), arho, grho2, lapl_val, atau, + sxc_l, v1xc_l, v2xc_l, v3xc_l, vlapl_l, ha, hse_om); + vlapl_vals[ir] = vlapl_l; + } + } + else + { +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 128) +#endif + for (std::size_t ir = 0; ir < nrxx_s; ++ir) + { + double arho1 = std::abs(rhotmp1[ir]); + double arho2 = std::abs(rhotmp2[ir]); + double atau1 = chr->kin_r[0][ir] / 2.0; + double atau2 = chr->kin_r[1][ir] / 2.0; + double laplup_val = lapl1.empty() ? 0.0 : lapl1[ir]; + double lapldw_val = lapl2.empty() ? 0.0 : lapl2[ir]; + double sxc_l = 0.0; + double v1xcup = 0.0, v1xcdw = 0.0; + double v2xcup = 0.0, v2xcdw = 0.0, v2xcud = 0.0; + double v3xcup = 0.0, v3xcdw = 0.0; + double vlaplup = 0.0, vlapldw = 0.0; + double ha = hybrid_alpha; + double hse_om = hse_omega; + XC_Functional_Libxc::tau_xc_spin(XC_Functional::get_func_id(), arho1, arho2, gdr1[ir], gdr2[ir], + laplup_val, lapldw_val, atau1, atau2, + sxc_l, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud, + v3xcup, v3xcdw, vlaplup, vlapldw, ha, hse_om); + vlapl_vals[ir * 2] = vlaplup; + vlapl_vals[ir * 2 + 1] = vlapldw; + } + } + + double vlapl_stress[6] = {0.0}; + int ab_map[6][2] = {{0,0}, {1,1}, {2,2}, {0,1}, {1,2}, {0,2}}; + for (int is = 0; is < nspin0; ++is) + { + for (int ab = 0; ab < 6; ++ab) + { + double sum = 0.0; +#ifdef _OPENMP +#pragma omp parallel for reduction(+:sum) schedule(static, 256) +#endif + for (std::size_t ir = 0; ir < nrxx_s; ++ir) + { + sum += vlapl_vals[ir * nspin0 + is] * hessian[is][ab * nrxx_s + ir]; + } + vlapl_stress[ab] += sum; + } + } + + double inv_omega = 1.0 / ucell->omega; + for (int ab = 0; ab < 6; ++ab) + { + int l = ab_map[ab][0]; + int m = ab_map[ab][1]; + int ind = l * 3 + m; + stress_gga[ind] += 2.0 * vlapl_stress[ab] * ModuleBase::e2 * inv_omega; + } +} +} // anonymous namespace #endif void XC_Functional::gradcorr( @@ -70,6 +216,23 @@ void XC_Functional::gradcorr( assert(nspin0>0); const double fac = 1.0/ nspin0; + bool need_laplacian = false; +#ifdef USE_LIBXC + if (use_libxc && (func_type == 3 || func_type == 5)) + { + std::vector check_funcs = XC_Functional_Libxc::init_func(func_id, XC_UNPOLARIZED, 0.0, 0.0); + for (auto& f : check_funcs) + { + if (f.info->flags & XC_FLAGS_NEEDS_LAPLACIAN) + { + need_laplacian = true; + break; + } + } + XC_Functional_Libxc::finish_func(check_funcs); + } +#endif + if(is_stress) { stress_gga.resize(9); @@ -100,6 +263,7 @@ void XC_Functional::gradcorr( double* neg = nullptr; double** vsave = nullptr; double** vgg = nullptr; + std::vector lapl1, lapl2; // for spin unpolarized case, // calculate the gradient of (rho_core+rho) in reciprocal space. @@ -128,6 +292,12 @@ void XC_Functional::gradcorr( XC_Functional::grad_rho( rhogsum1 , gdr1, rhopw, ucell->tpiba); + if(need_laplacian) + { + lapl1.resize(rhopw->nrxx); + XC_Functional::laplacian_rho(rhogsum1, lapl1.data(), rhopw, ucell->tpiba); + } + // for spin polarized case; // calculate the gradient of (rho_core+rho) in reciprocal space. if(nspin==2) @@ -156,6 +326,12 @@ void XC_Functional::gradcorr( } XC_Functional::grad_rho( rhogsum2 , gdr2, rhopw, ucell->tpiba); + + if(need_laplacian) + { + lapl2.resize(rhopw->nrxx); + XC_Functional::laplacian_rho(rhogsum2, lapl2.data(), rhopw, ucell->tpiba); + } } if(nspin == 4&&(domag||domag_z)) @@ -302,8 +478,10 @@ void XC_Functional::gradcorr( if(func_type == 3 || func_type == 5) { double v3xc = 0.0; + double vlapl = 0.0; double atau = chr->kin_r[0][ir]/2.0; - XC_Functional_Libxc::tau_xc( func_id, arho, grho2a, atau, sxc, v1xc, v2xc, v3xc, hybrid_alpha_in, hse_omega_in); + double lapl_val = lapl1.empty() ? 0.0 : lapl1[ir]; + XC_Functional_Libxc::tau_xc( func_id, arho, grho2a, lapl_val, atau, sxc, v1xc, v2xc, v3xc, vlapl, hybrid_alpha_in, hse_omega_in); } else { @@ -366,12 +544,16 @@ void XC_Functional::gradcorr( { double v3xcup = 0.0; double v3xcdw = 0.0; + double vlaplup = 0.0; + double vlapldw = 0.0; double atau1 = chr->kin_r[0][ir]/2.0; double atau2 = chr->kin_r[1][ir]/2.0; + double laplup_val = lapl1.empty() ? 0.0 : lapl1[ir]; + double lapldw_val = lapl2.empty() ? 0.0 : lapl2[ir]; XC_Functional_Libxc::tau_xc_spin( func_id, rhotmp1[ir], rhotmp2[ir], gdr1[ir], gdr2[ir], - atau1, atau2, sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud, v3xcup, v3xcdw, hybrid_alpha_in, hse_omega_in); + laplup_val, lapldw_val, atau1, atau2, sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud, v3xcup, v3xcdw, vlaplup, vlapldw, hybrid_alpha_in, hse_omega_in); } else { @@ -526,6 +708,17 @@ void XC_Functional::gradcorr( stress_gga [ind] += local_stress_gga [ind]; } } + + // Add vlapl stress contribution for mGGA functionals + if (need_laplacian && use_libxc) + { +#ifdef USE_LIBXC + add_vlapl_stress_contribution( + stress_gga, nspin0, rhotmp1, rhotmp2, gdr1, gdr2, + lapl1, lapl2, chr, rhopw, ucell, epsr, + hybrid_alpha_in, hse_omega_in); +#endif + } } else { @@ -824,6 +1017,42 @@ void XC_Functional::grad_dot( return; } + +void XC_Functional::laplacian_rho( + const std::complex* rhog, + double* lapl, + const ModulePW::PW_Basis* rho_basis, + const double tpiba) +{ + std::complex* lapl_g = new std::complex[rho_basis->nmaxgr]; + +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 1024) +#endif + for (int ig = 0; ig < rho_basis->npw; ig++) + { + double gg = rho_basis->gcar[ig][0] * rho_basis->gcar[ig][0] + + rho_basis->gcar[ig][1] * rho_basis->gcar[ig][1] + + rho_basis->gcar[ig][2] * rho_basis->gcar[ig][2]; + lapl_g[ig] = -rhog[ig] * gg; + } + + rho_basis->recip2real(lapl_g, lapl_g); + + const double tpiba2 = tpiba * tpiba; +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 1024) +#endif + for (int ir = 0; ir < rho_basis->nrxx; ir++) + { + lapl[ir] = lapl_g[ir].real() * tpiba2; + } + + delete[] lapl_g; + return; +} + + void XC_Functional::noncolin_rho( double *rhoout1, double *rhoout2,