From 812061e913c106ad7076bd02085a4afa49152bf2 Mon Sep 17 00:00:00 2001 From: linpeize Date: Mon, 26 Jan 2026 17:50:50 +0800 Subject: [PATCH 1/5] Fix(module_xc): compute real density laplacian and vlapl stress for meta-GGA functionals MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit SCANL, R2SCANL, and REVSCAN require the density laplacian (∇²ρ) but v_xc_meta was passing |∇ρ|² instead. This caused ~150 kbar stress error for Na with SCANL vs SCAN. Additionally, the vlapl stress contribution was completely missing. Laplacian computation method: ∇²ρ(r) = -∑_G |G|² ρ(G) e^{iGr} - ρ(r) is Fourier-transformed to ρ(G) via real2recip() - Each G-component is multiplied by -|G|² (i.e. -gg[ig] * tpiba²) - Inverse FFT via recip2real() yields ∇²ρ(r) on the real-space grid - Only the real part is retained (imaginary part should be ~0 for real ρ) This is exact within the plane-wave basis, not numerical differentiation. Spin handling: - cal_lapl(nspin, nrxx, rho, tpiba2, chr) computes ∇²ρ for each spin channel: rho is interleaved [rho_up[0], rho_dw[0], rho_up[1], ...] and output lapl is [lapl_up[0], lapl_dw[0], lapl_up[1], ...] - nspin=1: single channel (spin-unpolarized), loop runs once for is=0 - nspin=2: separate laplacian for spin-up and spin-down channels, each computed independently from its own real-space density - nspin=4 (noncollinear): blocked by set_xc_type for meta-GGA Hessian computation for vlapl stress: - cal_rho_hessian(nspin, nrxx, rho, chr) computes ∂²ρ/∂r_α∂r_β via reciprocal space: FFT[ -G_α G_β ρ(G) ] - Returns 6 independent components per spin (xx, yy, zz, xy, yz, zx) - Stress contribution: σ_αβ^vlapl = -2 ∫ vlapl(r) ∂²ρ/∂r_α∂r_β dr - Applied in gradcorr is_stress path for both nspin=1 and nspin=2 Changes: - Add cal_lapl() and cal_rho_hessian() as reusable functions - Fix v_xc_meta to call cal_lapl() and pass real laplacian to LibXC; handle vlapl output by computing ∇²(vlapl) via same G-space method - Fix tau_xc/tau_xc_spin wrappers to return vlapl derivative - Fix gradcorr stress/force paths: replace duplicated inline Laplacian code with cal_lapl() calls; add Hessian computation and vlapl stress contribution for meta-GGA functionals - Update test_xc4.cpp and test CMakeLists for new tau_xc signature --- source/module_esolver/esolver_ks.cpp | 1 + .../module_xc/test/CMakeLists.txt | 3 + .../module_xc/test/test_xc4.cpp | 5 +- .../module_xc/xc_functional.cpp | 67 +++++++- .../module_xc/xc_functional.h | 1 + .../module_xc/xc_functional_gradcorr.cpp | 156 +++++++++++++++++- .../module_xc/xc_functional_libxc.h | 23 ++- .../module_xc/xc_functional_libxc_tools.cpp | 92 +++++++++++ .../module_xc/xc_functional_libxc_vxc.cpp | 48 +++++- .../xc_functional_libxc_wrapper_tauxc.cpp | 25 +-- source/module_ri/Exx_LRI_interface.hpp | 15 +- 11 files changed, 394 insertions(+), 42 deletions(-) diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp index 5a8ec9a4d89..8afefa70dfc 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -185,6 +185,7 @@ void ESolver_KS::before_all_runners(UnitCell& ucell, const Input_para { XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func); } + GlobalV::ofs_running< funcs = XC_Functional_Libxc::init_func(func_id, XC_UNPOLARIZED); + for(const auto &func : funcs) + { + const xc_func_info_type *info = xc_func_get_info(&func); + ss<<" XC: "< get_func_id() { return func_id; } + static std::string output_info(); //------------------- // xc_functional_wrapper_xc.cpp diff --git a/source/module_hamilt_general/module_xc/xc_functional_gradcorr.cpp b/source/module_hamilt_general/module_xc/xc_functional_gradcorr.cpp index 42501465108..edfedf326e1 100644 --- a/source/module_hamilt_general/module_xc/xc_functional_gradcorr.cpp +++ b/source/module_hamilt_general/module_xc/xc_functional_gradcorr.cpp @@ -67,6 +67,10 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v, ModuleBase::Vector3* gdr2 = nullptr; ModuleBase::Vector3* h1 = nullptr; ModuleBase::Vector3* h2 = nullptr; + double* lapl1 = nullptr; + double* lapl2 = nullptr; + double* hess1 = nullptr; // nspin=1 Hessian (6 components) + double* hess2 = nullptr; // nspin=1 Hessian for spin-down channel double* neg = nullptr; double** vsave = nullptr; double** vgg = nullptr; @@ -95,6 +99,46 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v, XC_Functional::grad_rho( rhogsum1 , gdr1, rhopw, ucell->tpiba); + XC_Functional::grad_rho( rhogsum1 , gdr1, rhopw, ucell->tpiba); + + // compute laplacian of total density (rho[0] + rho_core) for nspin=1,2 + // and as the "up" channel equivalent for nspin=4 noncollinear case + { + std::vector rho_total(rhopw->nrxx); +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 1024) +#endif + for(int ir=0; irnrxx; ++ir) + rho_total[ir] = rhotmp1[ir]; + std::vector lapl_all = XC_Functional_Libxc::cal_lapl( + 1, rhopw->nrxx, rho_total, ucell->tpiba2, chr); + lapl1 = new double[rhopw->nrxx]; +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 1024) +#endif + for(int ir=0; irnrxx; ++ir) + lapl1[ir] = lapl_all[ir]; + } + + // compute Hessian of total density for meta-GGA stress + if(is_stress && (func_type == 3 || func_type == 5)) + { + std::vector rho_total(rhopw->nrxx); +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 1024) +#endif + for(int ir=0; irnrxx; ++ir) + rho_total[ir] = rhotmp1[ir]; + std::vector hess_all = XC_Functional_Libxc::cal_rho_hessian( + 1, rhopw->nrxx, rho_total, chr); + hess1 = new double[rhopw->nrxx * 6]; +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 1024) +#endif + for(int i=0; inrxx * 6; ++i) + hess1[i] = hess_all[i]; + } + // for spin polarized case; // calculate the gradient of (rho_core+rho) in reciprocal space. if(PARAM.inp.nspin==2) @@ -120,6 +164,43 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v, if(!is_stress) { h2 = new ModuleBase::Vector3[rhopw->nrxx]; } XC_Functional::grad_rho( rhogsum2 , gdr2, rhopw, ucell->tpiba); + + // compute laplacian of spin-down density (rho[1] + rho_core) + { + std::vector rho_dw(rhopw->nrxx); +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 1024) +#endif + for(int ir=0; irnrxx; ++ir) + rho_dw[ir] = rhotmp2[ir]; + std::vector lapl_all = XC_Functional_Libxc::cal_lapl( + 1, rhopw->nrxx, rho_dw, ucell->tpiba2, chr); + lapl2 = new double[rhopw->nrxx]; +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 1024) +#endif + for(int ir=0; irnrxx; ++ir) + lapl2[ir] = lapl_all[ir]; + } + + // compute Hessian of spin-down density for meta-GGA stress + if(is_stress && (func_type == 3 || func_type == 5)) + { + std::vector rho_dw(rhopw->nrxx); +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 1024) +#endif + for(int ir=0; irnrxx; ++ir) + rho_dw[ir] = rhotmp2[ir]; + std::vector hess_all = XC_Functional_Libxc::cal_rho_hessian( + 1, rhopw->nrxx, rho_dw, chr); + hess2 = new double[rhopw->nrxx * 6]; +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 1024) +#endif + for(int i=0; inrxx * 6; ++i) + hess2[i] = hess_all[i]; + } } if(PARAM.inp.nspin == 4&&(PARAM.globalv.domag||PARAM.globalv.domag_z)) @@ -189,6 +270,40 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v, XC_Functional::grad_rho( rhogsum1 , gdr1, rhopw, ucell->tpiba); XC_Functional::grad_rho( rhogsum2 , gdr2, rhopw, ucell->tpiba); + // re-compute laplacians after noncollinear rho update + { + std::vector rho_up(rhopw->nrxx); +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 1024) +#endif + for(int ir=0; irnrxx; ++ir) + rho_up[ir] = rhotmp1[ir]; + std::vector lapl_up = XC_Functional_Libxc::cal_lapl( + 1, rhopw->nrxx, rho_up, ucell->tpiba2, chr); + lapl1 = new double[rhopw->nrxx]; +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 1024) +#endif + for(int ir=0; irnrxx; ++ir) + lapl1[ir] = lapl_up[ir]; + } + { + std::vector rho_dw(rhopw->nrxx); +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 1024) +#endif + for(int ir=0; irnrxx; ++ir) + rho_dw[ir] = rhotmp2[ir]; + std::vector lapl_dw = XC_Functional_Libxc::cal_lapl( + 1, rhopw->nrxx, rho_dw, ucell->tpiba2, chr); + lapl2 = new double[rhopw->nrxx]; +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 1024) +#endif + for(int ir=0; irnrxx; ++ir) + lapl2[ir] = lapl_dw[ir]; + } + } const double epsr = 1.0e-6; @@ -249,9 +364,20 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v, #ifdef USE_LIBXC if(func_type == 3 || func_type == 5) //the gradcorr part to stress of mGGA { - double v3xc; + double v3xc, vlaplc; double atau = chr->kin_r[0][ir]/2.0; - XC_Functional_Libxc::tau_xc( func_id, arho, grho2a, atau, sxc, v1xc, v2xc, v3xc); + XC_Functional_Libxc::tau_xc( func_id, arho, grho2a, lapl1[ir], atau, sxc, v1xc, v2xc, v3xc, vlaplc); + // add vlapl stress contribution: σ_αβ^vlapl = -2 ∫ vlapl ∂²ρ/∂r_α∂r_β dr + // Hessian components order: xx, yy, zz, xy, yz, zx + for(int l = 0; l < 3; l++) + { + for(int m = 0; m < l+1; m++) + { + int ind = l*3 + m; + int ic = (l==0&&m==0) ? 0 : (l==1&&m==1) ? 1 : (l==2&&m==2) ? 2 : (l==0&&m==1) ? 3 : (l==1&&m==2) ? 4 : 5; + local_stress_gga[ind] -= 2.0 * vlaplc * hess1[ic * rhopw->nrxx + ir] * ModuleBase::e2; + } + } } else { @@ -308,13 +434,27 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v, double sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud; if(func_type == 3 || func_type == 5) //the gradcorr part to stress of mGGA { - double v3xcup, v3xcdw; + double v3xcup, v3xcdw, vlaplup, vlapldw; double atau1 = chr->kin_r[0][ir]/2.0; double atau2 = chr->kin_r[1][ir]/2.0; 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); + rhotmp1[ir], rhotmp2[ir], gdr1[ir], gdr2[ir], + lapl1[ir], lapl2[ir], + atau1, atau2, sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud, v3xcup, v3xcdw, vlaplup, vlapldw); + // add vlapl stress contribution for both spin channels + if(is_stress) + { + for(int l = 0; l < 3; l++) + { + for(int m = 0; m < l+1; m++) + { + int ind = l*3 + m; + int ic = (l==0&&m==0) ? 0 : (l==1&&m==1) ? 1 : (l==2&&m==2) ? 2 : (l==0&&m==1) ? 3 : (l==1&&m==2) ? 4 : 5; + local_stress_gga[ind] -= 2.0 * (vlaplup * hess1[ic * rhopw->nrxx + ir] + vlapldw * hess2[ic * rhopw->nrxx + ir]) * ModuleBase::e2; + } + } + } } else { @@ -582,6 +722,8 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v, delete[] rhotmp1; delete[] rhogsum1; delete[] gdr1; + delete[] lapl1; + delete[] hess1; if(!is_stress) { delete[] h1; } @@ -590,6 +732,8 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v, delete[] rhotmp2; delete[] rhogsum2; delete[] gdr2; + delete[] lapl2; + delete[] hess2; if(!is_stress) { delete[] h2; } } @@ -609,6 +753,8 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v, delete[] rhotmp2; delete[] rhogsum2; delete[] gdr2; + delete[] lapl2; + delete[] hess2; } return; diff --git a/source/module_hamilt_general/module_xc/xc_functional_libxc.h b/source/module_hamilt_general/module_xc/xc_functional_libxc.h index 12cf1fdd8ed..acadcd866f9 100644 --- a/source/module_hamilt_general/module_xc/xc_functional_libxc.h +++ b/source/module_hamilt_general/module_xc/xc_functional_libxc.h @@ -82,6 +82,22 @@ namespace XC_Functional_Libxc extern std::vector convert_sigma( const std::vector>> &gdr); + // calculating laplacian of density: ∇²ρ(r) + extern std::vector cal_lapl( + const int nspin, + const std::size_t nrxx, + const std::vector &rho, + const double tpiba2, + const Charge* const chr); + + // calculating Hessian of density: ∂²ρ/∂r_α∂r_β + // returns nspin * nrxx * 6 array (xx, yy, zz, xy, yz, zx) + extern std::vector cal_rho_hessian( + const int nspin, + const std::size_t nrxx, + const std::vector &rho, + const Charge* const chr); + // sgn for threshold mask extern std::vector cal_sgn( const double rho_threshold, @@ -166,16 +182,17 @@ namespace XC_Functional_Libxc // wrapper for the mGGA functionals extern void tau_xc( const std::vector &func_id, - const double &rho, const double &grho, const double &atau, double &sxc, - double &v1xc, double &v2xc, double &v3xc); + const double &rho, const double &grho, const double &lapl_rho, const double &atau, double &sxc, + double &v1xc, double &v2xc, double &v3xc, double &vlaplc); 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, double &v1xcup, double &v1xcdw, double &v2xcup, double &v2xcdw, double &v2xcud, - double &v3xcup, double &v3xcdw); + double &v3xcup, double &v3xcdw, double &vlaplup, double &vlapldw); } // namespace XC_Functional_Libxc #endif // USE_LIBXC diff --git a/source/module_hamilt_general/module_xc/xc_functional_libxc_tools.cpp b/source/module_hamilt_general/module_xc/xc_functional_libxc_tools.cpp index 567d04e62a8..f4297635cfd 100644 --- a/source/module_hamilt_general/module_xc/xc_functional_libxc_tools.cpp +++ b/source/module_hamilt_general/module_xc/xc_functional_libxc_tools.cpp @@ -82,6 +82,98 @@ XC_Functional_Libxc::cal_gdr( return gdr; } +// calculating laplacian of density: ∇²ρ(r) +std::vector XC_Functional_Libxc::cal_lapl( + const int nspin, + const std::size_t nrxx, + const std::vector &rho, + const double tpiba2, + const Charge* const chr) +{ + std::vector lapl(nspin * 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> rhog(chr->rhopw->npw); + chr->rhopw->real2recip(rhor.data(), rhog.data()); + + // ∇²ρ(r) = -∑_G |G|² ρ(G) e^{iGr} + std::vector> rhog_lapl(chr->rhopw->npw); + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for(int ig=0; igrhopw->npw; ++ig) + rhog_lapl[ig] = -rhog[ig] * chr->rhopw->gg[ig] * tpiba2; + + std::vector> aux(chr->rhopw->nmaxgr); + chr->rhopw->recip2real(rhog_lapl.data(), aux.data()); + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for(std::size_t ir=0; ir XC_Functional_Libxc::cal_rho_hessian( + const int nspin, + const std::size_t nrxx, + const std::vector &rho, + const Charge* const chr) +{ + std::vector hess(nspin * nrxx * 6, 0.0); + // ipol2xy maps 6 components to (x,y) pairs: 0:xx, 1:yy, 2:zz, 3:xy, 4:yz, 5:zx + const int ipol2xy[6][2] = {{0,0}, {1,1}, {2,2}, {0,1}, {1,2}, {2,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()); + + // compute all 6 Hessian components in one G-space pass + for(int ic=0; ic<6; ++ic) + { + const int ax = ipol2xy[ic][0]; + const int ay = ipol2xy[ic][1]; + std::vector> rhog_hess(chr->rhopw->npw); + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for(int ig=0; igrhopw->npw; ++ig) + { + const double gx = chr->rhopw->gcar[ig][ax]; + const double gy = chr->rhopw->gcar[ig][ay]; + rhog_hess[ig] = -rhog[ig] * gx * gy; + } + + std::vector> aux(chr->rhopw->nmaxgr); + chr->rhopw->recip2real(rhog_hess.data(), aux.data()); + #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) diff --git a/source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp b/source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp index 2456abe5ba9..5eaaa347ef5 100644 --- a/source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp +++ b/source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp @@ -199,6 +199,10 @@ std::tuple XC_Functional_Li = XC_Functional_Libxc::cal_gdr(nspin, nrxx, rho, tpiba, chr); const std::vector sigma = XC_Functional_Libxc::convert_sigma(gdr); + // compute laplacian of density: ∇²ρ(r) for each spin channel + const double tpiba2 = tpiba * tpiba; + const std::vector lapl = XC_Functional_Libxc::cal_lapl(nspin, nrxx, rho, tpiba2, chr); + //converting kin_r std::vector kin_r; kin_r.resize(nrxx*nspin); @@ -264,7 +268,7 @@ std::tuple XC_Functional_Li for ( xc_func_type &func : funcs ) { assert(func.info->family == XC_FAMILY_MGGA); - xc_mgga_exc_vxc(&func, nrxx, rho.data(), sigma.data(), sigma.data(), + xc_mgga_exc_vxc(&func, nrxx, rho.data(), sigma.data(), lapl.data(), kin_r.data(), exc.data(), vrho.data(), vsigma.data(), vlapl.data(), vtau.data()); //process etxc @@ -387,6 +391,48 @@ std::tuple XC_Functional_Li vofk(is,ir) += vtau[ir*nspin+is] * sgn[ir*nspin+is]; } } + + //process vlapl: ∇²(vlapl) contribution to potential + // vlapl = ∂ε/∂(∇²ρ), potential contribution = +∇²(vlapl) + std::vector vlapl_real(nrxx * nspin); + for( int is=0; is> vlapl_g(chr->rhopw->npw); + chr->rhopw->real2recip(&vlapl_real[is * nrxx], vlapl_g.data()); + + std::vector> vlapl_lapl_g(chr->rhopw->npw); +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 1024) +#endif + for(int ig=0; igrhopw->npw; ++ig) + vlapl_lapl_g[ig] = -vlapl_g[ig] * chr->rhopw->gg[ig] * tpiba2; + + std::vector> aux(chr->rhopw->nmaxgr); + chr->rhopw->recip2real(vlapl_lapl_g.data(), aux.data()); + + for( int ir=0; ir< nrxx; ++ir ) + { +#ifdef __EXX + double vlapl_contrib = aux[ir].real(); + if (func.info->number == XC_MGGA_X_SCAN && XC_Functional::get_func_type() == 5) + { + vlapl_contrib *= (1.0 - XC_Functional::get_hybrid_alpha()); + } + v(is,ir) += vlapl_contrib; +#else + v(is,ir) += aux[ir].real(); +#endif + } + } } //------------------------------------------------- diff --git a/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_tauxc.cpp b/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_tauxc.cpp index b265af9bb15..ee66590b801 100644 --- a/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_tauxc.cpp +++ b/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_tauxc.cpp @@ -12,19 +12,17 @@ //XC_POLARIZED, XC_UNPOLARIZED: internal flags used in LIBXC, denote the polarized(nspin=1) or unpolarized(nspin=2) calculations, definition can be found in xc.h from LIBXC void XC_Functional_Libxc::tau_xc( const std::vector &func_id, - const double &rho, const double &grho, const double &atau, double &sxc, - double &v1xc, double &v2xc, double &v3xc) + const double &rho, const double &grho, const double &lapl_rho, const double &atau, double &sxc, + double &v1xc, double &v2xc, double &v3xc, double &vlaplc) { - double s, v1, v2, v3; - double lapl_rho, vlapl_rho; - lapl_rho = grho; + double s, v1, v2, v3, vlapl; std::vector funcs = XC_Functional_Libxc::init_func(func_id, XC_UNPOLARIZED); - sxc = 0.0; v1xc = 0.0; v2xc = 0.0; v3xc = 0.0; + sxc = 0.0; v1xc = 0.0; v2xc = 0.0; v3xc = 0.0; vlaplc = 0.0; for(xc_func_type &func : funcs) { - xc_mgga_exc_vxc(&func,1,&rho,&grho,&lapl_rho,&atau,&s,&v1,&v2,&vlapl_rho,&v3); + xc_mgga_exc_vxc(&func,1,&rho,&grho,&lapl_rho,&atau,&s,&v1,&v2,&vlapl,&v3); #ifdef __EXX if (func.info->number == XC_MGGA_X_SCAN && XC_Functional::get_func_type() == 5) { @@ -32,12 +30,14 @@ void XC_Functional_Libxc::tau_xc( v1 *= (1.0 - GlobalC::exx_info.info_global.hybrid_alpha); v2 *= (1.0 - GlobalC::exx_info.info_global.hybrid_alpha); v3 *= (1.0 - GlobalC::exx_info.info_global.hybrid_alpha); + vlapl *= (1.0 - GlobalC::exx_info.info_global.hybrid_alpha); } #endif sxc += s * rho; v2xc += v2 * 2.0; v1xc += v1; v3xc += v3; + vlaplc += vlapl; } XC_Functional_Libxc::finish_func(funcs); @@ -49,16 +49,19 @@ void XC_Functional_Libxc::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, double &v1xcup, double &v1xcdw, double &v2xcup, double &v2xcdw, double &v2xcud, - double &v3xcup, double &v3xcdw) + double &v3xcup, double &v3xcdw, double &vlaplup, double &vlapldw) { sxc = v1xcup = v1xcdw = 0.0; v2xcup = v2xcdw = v2xcud = 0.0; v3xcup = v3xcdw = 0.0; + vlaplup = vlapldw = 0.0; const std::array rho = {rhoup, rhodw}; const std::array grho = {gdr1.norm2(), gdr1 * gdr2, gdr2.norm2()}; + const std::array lapl = {laplup, lapldw}; const std::array tau = {tauup, taudw}; std::vector funcs = XC_Functional_Libxc::init_func(func_id, XC_POLARIZED); @@ -79,9 +82,9 @@ void XC_Functional_Libxc::tau_xc_spin( } double s = 0.0; - std::array v1xc, v3xc, lapl, vlapl; + std::array v1xc, v3xc, vlapl; std::array v2xc; - // call Libxc function: xc_mgga_exc_vxc + // call Libxc function: xc_mgga_exc_vxc with real laplacian values xc_mgga_exc_vxc( &func, 1, rho.data(), grho.data(), lapl.data(), tau.data(), &s, v1xc.data(), v2xc.data(), vlapl.data(), v3xc.data()); sxc += s * (rho[0] * sgn[0] + rho[1] * sgn[1]); @@ -92,6 +95,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/module_ri/Exx_LRI_interface.hpp b/source/module_ri/Exx_LRI_interface.hpp index 724835d9974..0605c7e687d 100644 --- a/source/module_ri/Exx_LRI_interface.hpp +++ b/source/module_ri/Exx_LRI_interface.hpp @@ -147,20 +147,7 @@ void Exx_LRI_Interface::exx_beforescf(const int istep, } else { - if (ucell.atoms[0].ncpp.xc_func == "HF" || ucell.atoms[0].ncpp.xc_func == "PBE0" || ucell.atoms[0].ncpp.xc_func == "HSE") - { - XC_Functional::set_xc_type("pbe"); - } - else if (ucell.atoms[0].ncpp.xc_func == "SCAN0") - { - XC_Functional::set_xc_type("scan"); - } - // added by jghan, 2024-07-07 - else if ( ucell.atoms[0].ncpp.xc_func == "MULLER" || ucell.atoms[0].ncpp.xc_func == "POWER" - || ucell.atoms[0].ncpp.xc_func == "WP22" || ucell.atoms[0].ncpp.xc_func == "CWP22" ) - { - XC_Functional::set_xc_type("pbe"); - } + XC_Functional::set_xc_first_loop(ucell); } this->cal_exx_ions(ucell,PARAM.inp.out_ri_cv); From b096ead773fee79ac11d30c21f3b55743c5c850f Mon Sep 17 00:00:00 2001 From: dyzheng Date: Thu, 30 Apr 2026 17:30:57 +0800 Subject: [PATCH 2/5] Fix(module_xc): omit vlapl potential for SCF stability, add runtime warning MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit For meta-GGA functionals that depend on the Laplacian of the density (e.g., SCANL), the vlapl contribution to V_xc (e2*∇²(vlapl)) is omitted because it causes |G|² amplification in reciprocal space leading to SCF divergence. vtxc also omits vlapl for consistency with V. The vlapl energy is still included in etxc via exc, and the vlapl stress term (2*vlapl*Hess*e2) is computed in gradcorr. A runtime warning is printed when a Laplacian-dependent functional is detected, recommending r2SCAN or SCAN instead. See docs/scanl_laplacian_implementation_report.md for the full analysis including FD stress validation results. --- docs/scanl_laplacian_implementation_report.md | 297 ++++++++++++++++++ .../module_xc/xc_functional_gradcorr.cpp | 11 +- .../module_xc/xc_functional_libxc_vxc.cpp | 57 ++-- 3 files changed, 325 insertions(+), 40 deletions(-) create mode 100644 docs/scanl_laplacian_implementation_report.md diff --git a/docs/scanl_laplacian_implementation_report.md b/docs/scanl_laplacian_implementation_report.md new file mode 100644 index 00000000000..dbb48fb9d4a --- /dev/null +++ b/docs/scanl_laplacian_implementation_report.md @@ -0,0 +1,297 @@ +# SCANL Laplacian 依赖项实现与验证报告 + +## 1. 背景 + +SCANL（SCAN-L）是 SCAN 的 deorbitalized 版本，将 SCAN 中对动能密度 τ 的依赖替换为对密度 Laplacian ∇²ρ 的依赖 [Mejia-Rodriguez & Trickey, Phys. Rev. B **98**, 115161 (2018)]。这导致其 Kohn-Sham 势和应力公式中出现了 ∇²(∂ε/∂(∇²ρ)) 项，在平面波基组下引发严重的数值稳定性问题。 + +本文档记录了 ABACUS 中 SCANL 泛函 Laplacian 相关项的实现方案、有限差分（FD）应力验证结果，以及最终建议。 + +--- + +## 2. 理论公式 + +### 2.1 meta-GGA 能量泛函 + +meta-GGA 交换关联能量密度 ε_xc 依赖于四个变量： + +$$E_{xc} = \int \rho(\mathbf{r})\,\varepsilon_{xc}(\rho, |\nabla\rho|^2, \nabla^2\rho, \tau)\, d\mathbf{r}$$ + +对于 SCANL，∂ε/∂(∇²ρ) ≠ 0（与 SCAN 不同，SCAN 的此项为零）。 + +### 2.2 Kohn-Sham 势 + +由变分原理，XC 势为： + +$$V_{xc} = e^2 \frac{\delta E_{xc}}{\delta \rho} = e^2 \left[ \frac{\partial(\rho\varepsilon)}{\partial\rho} - \nabla \cdot \frac{\partial(\rho\varepsilon)}{\partial(\nabla\rho)} + \nabla^2 \frac{\partial(\rho\varepsilon)}{\partial(\nabla^2\rho)} + \frac{\partial(\rho\varepsilon)}{\partial\tau} \hat{\tau} \right]$$ + +其中 **vlapl 势项**为： + +$$V_{xc}^{\nabla^2\rho} = e^2\,\nabla^2\!\left(\frac{\partial(\rho\varepsilon)}{\partial(\nabla^2\rho)}\right)$$ + +在倒易空间中，∇² → -|G|²，因此： + +$$\tilde{V}_{xc}^{\nabla^2\rho}(\mathbf{G}) = -e^2\,|\mathbf{G}|^2\,\widetilde{\text{vlapl}}(\mathbf{G})$$ + +**关键问题**：|G|² 是高通滤波器，放大高频模式。vlapl(G) 中的数值噪声被 |G|² 放大后反馈到密度中，形成正反馈循环导致 SCF 发散。ecutwfc 越大，|G|_max 越大，放大效应越强，SCF 越容易发散。这并非 mixing 策略的问题，而是 Laplace 逆变换的数学病态性。 + +### 2.3 vtxc 与应力公式 + +vtxc（势对密度的积分）在应力计算中至关重要： + +$$\text{vtxc} = \int V_{xc}(\mathbf{r})\,\rho(\mathbf{r})\, d\mathbf{r}$$ + +应力中 XC 对角贡献为： + +$$\sigma_{ii}^{xc} = -\frac{E_{txc} - V_{txc}}{\Omega}$$ + +其中 $E_{txc} = \int e^2\,\rho\,\varepsilon_{xc}\, d\mathbf{r}$。 + +**vtxc 必须与实际使用的 V_xc 一致**。若 V 中不含 vlapl 而 vtxc 中含 vlapl，则 $-(E_{txc}-V_{txc})/\Omega$ 不正确。 + +### 2.4 vlapl 应力项 + +由能量泛函对应变的变分推导，vlapl 对应力的贡献为： + +$$\sigma_{\alpha\beta}^{\nabla^2\rho} = \frac{2}{\Omega} \sum_{\mathbf{r}} \frac{\partial(\rho\varepsilon)}{\partial(\nabla^2\rho)}\, H_{\alpha\beta}(\mathbf{r})\, e^2$$ + +其中 $H_{\alpha\beta} = \partial^2\rho / \partial r_\alpha \partial r_\beta$ 是密度的 Hessian 矩阵。 + +**注意**：LibXC 返回的 vlapl = ∂(ρε)/∂(∇²ρ) 已包含 ρ 因子，无需再乘 ρ。 + +--- + +## 3. 实现方案与测试结果 + +### 3.1 验证方法 + +采用有限差分（FD）应力验证：对 Si2 FCC 原胞（Γ 点）在 a₀±δ 下计算总能量，由中心差分公式得到应力： + +$$\sigma_{FD} = -\frac{E(a_0+\delta) - E(a_0-\delta)}{V(a_0+\delta) - V(a_0-\delta)}$$ + +与解析应力 σ_AB 比较，FD/AB 越接近 1.0 越好。 + +**单位转换**：ABACUS 总能量单位为 eV，应力原子单位为 Ry/Bohr³。FD 应力转换公式： + +$$\sigma_{FD}(\text{kbar}) = \sigma_{FD}(\text{eV/Bohr}^3) \times \frac{2\,\text{Ry}}{27.211386\,\text{eV}} \times 147105\,\frac{\text{kbar}}{\text{Ry/Bohr}^3}$$ + +### 3.2 测试参数 + +| 参数 | 值 | +|------|-----| +| 体系 | Si₂ FCC 原胞 | +| a₀ | 10.2 Bohr | +| δ | 0.001 Bohr | +| ecutwfc | 60 / 80 / 100 / 400 Ry | +| smearing_sigma | 0.02 | +| ks_solver | dav_subspace | +| mixing_beta | 0.1 | +| mixing_tau | 1 | +| scf_thr | 1e-7 | +| pw_seed | 1 | + +### 3.3 SCAN 基线验证 + +SCAN 不依赖 ∇²ρ（LibXC 返回 vlapl=0），用于验证 FD 方法和代码的正确性： + +| ecutwfc (Ry) | ecutwfc (eV) | FD (kbar) | AB (kbar) | FD/AB | 误差 | +|-------------|-------------|-----------|-----------|-------|------| +| 60 | 816 | 265.69 | 265.66 | 1.0001 | 0.01% | +| 80 | 1088 | 266.11 | 266.11 | 1.0000 | 0.00% | +| 100 | 1361 | 266.11 | 266.12 | 0.9999 | -0.01% | +| 400 | 5442 | 266.11 | 266.09 | 1.0001 | 0.01% | + +**结论**：SCAN 在所有 ecutwfc 下 FD/AB ≈ 1.000，验证了 FD 方法和应力代码的正确性。 + +### 3.4 各方案详情与结果 + +--- + +#### 方案 1：完整 e2·∇²(vlapl) 加入势和 vtxc（严格实现） + +**实现**： +- V_xc：含 e2·∇²(vlapl·sgn) +- vtxc：含 e2·∫∇²(vlapl·sgn)·ρ·sgn +- stress：含 2·vlapl·Hess·e2 + +**结果**：**SCF 发散**。无论 mixing_beta（0.01–0.7）、Kerker mixing、绝热演化（50 步 ramp）均无法收敛。 + +**原因**：V(G) = -e2·|G|²·vlapl(G)，|G|² 高通滤波放大高频噪声，正反馈增益超过任何混合策略的阻尼能力。ecutwfc 越大，|G|_max 越大，发散越严重。 + +--- + +#### 方案 2：sigma 伪装 lapl（原始代码方案） + +**实现**： +- `v_xc_libxc`：传 sigma（= |∇ρ|²）给 LibXC 的 Laplacian 参数位置 +- `tau_xc`/`tau_xc_spin`：传 grho（= |∇ρ|²）给 LibXC 的 Laplacian 参数位置 +- gradcorr stress：无 vlapl 应力项 + +**问题**：LibXC 在错误的 Laplacian 输入下计算 exc 和 vlapl，导致 etxc 和应力均不正确。 + +**FD 测试结果**： + +| ecutwfc (Ry) | FD (kbar) | AB (kbar) | FD/AB | 误差 | +|-------------|-----------|-----------|-------|------| +| 60 | 249.90 | 240.91 | 1.037 | 3.73% | +| 80 | 252.32 | 241.11 | 1.047 | 4.65% | +| 100 | 252.39 | 241.12 | 1.047 | 4.68% | +| 400 | 252.39 | 241.24 | 1.046 | 4.62% | + +FD/AB ≈ 1.04–1.05，存在系统性 4–5% 误差。SCANL 误差随 ecutwfc 趋于稳定（~4.6%），说明这是传假 Laplacian 导致的系统偏差，而非 ecutwfc 收敛问题。 + +--- + +#### 方案 3：V 和 vtxc 均不含 vlapl，stress 保留 vlapl 项 + +**实现**： +- `v_xc_libxc`：传真实 ∇²ρ 给 LibXC，获取正确的 exc、vrho、vsigma、vtau +- V_xc：不含 e2·∇²(vlapl) +- vtxc：不含 vlapl 贡献（与 V 保持一致） +- gradcorr stress：含 2·vlapl·Hess·e2 + +**物理逻辑**：vlapl 势无法加入 V（会发散），因此 vtxc 也不含 vlapl 以保持一致。etxc 中包含完整的 exc（含 Laplacian 贡献）。gradcorr 中的 vlapl 应力项从完整能量泛函严格推导。 + +**ecutwfc=200 结果**（此前测试，单位转换有误，仅参考趋势）：FD/AB = 1.018（1.82% 误差） + +**ecutwfc=400 结果**：待验证（需在对应分支上重新构建和测试）。 + +--- + +#### 方案 4：势中无 vlapl，vtxc 中加 e2·vlapl·∇²ρ + +**实现**： +- V_xc：不含 vlapl +- vtxc：含 e2·∫vlapl·∇²ρ·sgn +- stress：含 2·vlapl·Hess·e2 + +**问题**：V 无 vlapl 但 vtxc 有 → vtxc 与实际 V 不一致 → $-(E_{txc}-V_{txc})/\Omega$ 不正确 → 对角应力偏小。此前测试 FD/AB = 0.922（7.78% 误差）。 + +--- + +#### 方案 5：高斯衰减 + +**实现**： +- V_xc：e2·∇²(vlapl·sgn·damp(G))，damp = exp(-α·G²/G_max²) +- vtxc：从衰减势计算或从未衰减 vlapl·∇²ρ 计算 +- stress：含 2·vlapl·Hess·e2（无衰减） + +**此前测试结果**（ecutwfc=200 Ry，仅参考趋势）： + +| α | vtxc 来源 | FD/AB | +|---|----------|-------| +| 1.0 | 衰减势 | 0.789 | +| 2.0 | 衰减势 | 0.960 | +| 4.0 | 衰减势 | 0.918 | + +**问题**：衰减改变了势的物理形状，密度在错误势下自洽。不同 α 给出不同结果，无法确定哪个是"正确的"。本质上是调参拟合，没有物理依据。 + +--- + +#### 方案 6：绝热演化（ramp） + +**实现**： +- V_xc：e2·min(1, iter/ramp_steps)·∇²(vlapl·sgn) +- vtxc：e2·∫∇²(vlapl·sgn)·ρ·sgn（全强度） + +**结果**：**SCF 发散**。即使 ramp_steps=50 + mixing_beta=0.03，当 ramp_factor → 1 时仍开始振荡。 + +--- + +### 3.5 测试结果汇总 + +| 方案 | V 含 vlapl | vtxc 含 vlapl | stress 含 vlapl | SCF | FD/AB | 误差 | +|------|-----------|--------------|----------------|-----|-------|------| +| SCAN 基线 | N/A | N/A | N/A | 收敛 | 1.000 | 0.00% | +| 方案 1 (严格) | ✓ | ✓ | ✓ | **发散** | — | — | +| 方案 2 (sigma-as-lapl) | ✗ (假) | ✗ (假) | ✗ | 收敛 | 1.046 | 4.6% | +| 方案 3 (严格应力) | ✗ | ✗ | ✓ | 收敛 | ~1.02 | ~2% | +| 方案 4 (vtxc含vlapl) | ✗ | ✓ | ✓ | 收敛 | 0.922 | -7.8% | +| 方案 5 (衰减) | ✓(衰减) | ✓(衰减) | ✓ | 收敛 | ~0.9 | ~10% | +| 方案 6 (ramp) | ✓(渐增) | ✓ | ✓ | **发散** | — | — | + +注：方案 3、4、5 的 FD/AB 值来自 ecutwfc=200 Ry 的测试，由于此前存在单位转换错误，这些数值仅反映趋势，需用正确转换因子重新验证。 + +--- + +## 4. 模守恒赝势与 SCANL 的根本矛盾 + +### 4.1 ∇²(vlapl) 发散的 ecutwfc 依赖 + +vlapl 势在倒易空间为 V(G) = -e2·|G|²·vlapl(G)。|G|_max 由 ecutwfc 决定： + +$$|G|_{max} = \sqrt{2 \cdot E_{cut}}$$ + +|G|_max² 随 ecutwfc 线性增长，意味着 vlapl 势的高频放大随 ecutwfc 增强而加剧： + +- **低 ecutwfc**：|G|_max 小，放大范围有限，SCF 可能收敛，但基组不完备导致其他误差 +- **高 ecutwfc**：|G|_max 大，高频被剧烈放大，SCF 几乎必然发散 + +**严格实现方案 1 在任何实际使用的 ecutwfc 下都无法收敛**。 + +### 4.2 模守恒赝势需要高 ecutwfc + +模守恒赝势（NCPP）比 PAW 方法需要更高的 ecutwfc 才能收敛： + +| 赝势类型 | 典型 ecutwfc (Ry) | 典型 ecutwfc (eV) | +|---------|------------------|------------------| +| PAW/USPP | 30–50 | 400–700 | +| NCPP | 60–100+ | 800–1400+ | + +ABACUS 目前使用模守恒赝势。NCPP 要求的高 ecutwfc 与 vlapl 势的 |G|² 发散直接矛盾： + +- **低 ecutwfc**（< 50 Ry）：SCF 可能收敛，但 NCPP 的赝势不完备，结果不收敛 +- **高 ecutwfc**（> 60 Ry）：NCPP 收敛，但 vlapl 势发散 + +**这是一个不可调和的矛盾**。在 NCPP + 平面波框架下，不存在一个 ecutwfc 值能同时满足赝势收敛和 vlapl 势稳定。 + +### 4.3 VASP 的做法 + +VASP wiki [METAGGA 页面] 明确指出： + +> The results obtained with the meta-GGA functionals that depend on the Laplacian of the density ∇²n (e.g., SCAN-L) may not be reliable for large values of the energy cutoff ENCUT due to numerical instability. According to some tests, it is not recommended to use values of ENCUT above 800 eV. + +VASP 使用 PAW 方法，ecutwfc 要求较低（~400–700 eV），可以在 800 eV 以下的窗口中获得结果。但即便如此，VASP 也承认高截断能下结果不可靠。对于 NCPP，800 eV（约 60 Ry）可能刚达到赝势收敛的门槛，而更高截断又会触发不稳定——几乎不存在可用的 ecutwfc 窗口。 + +--- + +## 5. 结论与建议 + +### 5.1 不推荐使用 SCANL + +**在 ABACUS（NCPP + 平面波）框架下不推荐使用 SCANL 及其他 ∇²ρ 依赖泛函**。原因： + +1. **严格实现不可行**：vlapl 势 e2·∇²(vlapl) 在任何实际 ecutwfc 下导致 SCF 发散 +2. **NCPP 矛盾**：模守恒赝势需要 ≥60 Ry 的 ecutwfc，此范围内 vlapl 势必然不稳定 +3. **近似方案不可靠**：无论省略 vlapl 势（方案 3）还是用 sigma 伪装（方案 2），所得密度都不是 SCANL 自洽密度，结果与 FD 验证存在 2–5% 的偏差 +4. **VASP 也承认问题**：即使在 PAW 方法下，VASP 也建议不要使用超过 800 eV 的截断 + +### 5.2 推荐替代泛函 + +| 泛函 | 依赖变量 | 特点 | +|------|---------|------| +| **r²SCAN** | ρ, ∇ρ, τ | SCAN 的正则化版本，数值稳定，推荐首选 | +| **SCAN** | ρ, ∇ρ, τ | 与 SCANL 精度相近，无 ∇²ρ 依赖，数值稳定 | + +r²SCAN 是目前推荐的 meta-GGA 泛函，兼具 SCAN 的精度和更好的数值稳定性，且无 ∇²ρ 依赖问题。 + +### 5.3 代码实现建议 + +若仍需保留 SCANL 支持（例如用于与文献对比），建议采用**方案 3 + 运行时警告**： + +1. **V_xc 不含 vlapl**：e2·∇²(vlapl) 导致 |G|² 发散，无法加入 +2. **vtxc 不含 vlapl**：必须与 V 保持一致 +3. **etxc 含完整 exc**：LibXC 返回的能量密度已含 Laplacian 贡献 +4. **gradcorr stress 含 2·vlapl·Hess·e2**：从完整能量泛函严格推导 +5. **运行时打印 WARNING** + +在用户选择 ∇²ρ 依赖泛函时，应打印类似以下警告： + +``` +WARNING: The selected functional depends on the Laplacian of the density (∇²ρ). +The vlapl contribution to the Kohn-Sham potential is omitted because including +it (e2·∇²(vlapl)) causes SCF divergence due to |G|² amplification in reciprocal +space. The self-consistent density is therefore NOT the exact SCANL density. +The stress includes the vlapl term derived from the full energy functional. +Results may be unreliable, especially for norm-conserving pseudopotentials +which require high ecutwfc. Consider using r2SCAN or SCAN instead. +``` diff --git a/source/module_hamilt_general/module_xc/xc_functional_gradcorr.cpp b/source/module_hamilt_general/module_xc/xc_functional_gradcorr.cpp index edfedf326e1..b1a8f9501a9 100644 --- a/source/module_hamilt_general/module_xc/xc_functional_gradcorr.cpp +++ b/source/module_hamilt_general/module_xc/xc_functional_gradcorr.cpp @@ -367,15 +367,15 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v, double v3xc, vlaplc; double atau = chr->kin_r[0][ir]/2.0; XC_Functional_Libxc::tau_xc( func_id, arho, grho2a, lapl1[ir], atau, sxc, v1xc, v2xc, v3xc, vlaplc); - // add vlapl stress contribution: σ_αβ^vlapl = -2 ∫ vlapl ∂²ρ/∂r_α∂r_β dr - // Hessian components order: xx, yy, zz, xy, yz, zx + // vlapl stress: σ_αβ = (2/Ω) Σ vlapl·Hessian_αβ·e2 + // vlapl = ∂(ρε)/∂(∇²ρ) already includes ρ factor for(int l = 0; l < 3; l++) { for(int m = 0; m < l+1; m++) { int ind = l*3 + m; int ic = (l==0&&m==0) ? 0 : (l==1&&m==1) ? 1 : (l==2&&m==2) ? 2 : (l==0&&m==1) ? 3 : (l==1&&m==2) ? 4 : 5; - local_stress_gga[ind] -= 2.0 * vlaplc * hess1[ic * rhopw->nrxx + ir] * ModuleBase::e2; + local_stress_gga[ind] += 2.0 * vlaplc * hess1[ic * rhopw->nrxx + ir] * ModuleBase::e2; } } } @@ -442,7 +442,8 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v, rhotmp1[ir], rhotmp2[ir], gdr1[ir], gdr2[ir], lapl1[ir], lapl2[ir], atau1, atau2, sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud, v3xcup, v3xcdw, vlaplup, vlapldw); - // add vlapl stress contribution for both spin channels + // vlapl stress: σ_αβ = (2/Ω) Σ (vlapl_up·Hess_up + vlapl_dw·Hess_dw)·e2 + // vlapl = ∂(ρε)/∂(∇²ρ) already includes ρ factor if(is_stress) { for(int l = 0; l < 3; l++) @@ -451,7 +452,7 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v, { int ind = l*3 + m; int ic = (l==0&&m==0) ? 0 : (l==1&&m==1) ? 1 : (l==2&&m==2) ? 2 : (l==0&&m==1) ? 3 : (l==1&&m==2) ? 4 : 5; - local_stress_gga[ind] -= 2.0 * (vlaplup * hess1[ic * rhopw->nrxx + ir] + vlapldw * hess2[ic * rhopw->nrxx + ir]) * ModuleBase::e2; + local_stress_gga[ind] += 2.0 * (vlaplup * hess1[ic * rhopw->nrxx + ir] + vlapldw * hess2[ic * rhopw->nrxx + ir]) * ModuleBase::e2; } } } diff --git a/source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp b/source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp index 5eaaa347ef5..c5edc108ec9 100644 --- a/source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp +++ b/source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp @@ -392,45 +392,32 @@ std::tuple XC_Functional_Li } } - //process vlapl: ∇²(vlapl) contribution to potential - // vlapl = ∂ε/∂(∇²ρ), potential contribution = +∇²(vlapl) - std::vector vlapl_real(nrxx * nspin); - for( int is=0; is> vlapl_g(chr->rhopw->npw); - chr->rhopw->real2recip(&vlapl_real[is * nrxx], vlapl_g.data()); - - std::vector> vlapl_lapl_g(chr->rhopw->npw); -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 1024) -#endif - for(int ig=0; igrhopw->npw; ++ig) - vlapl_lapl_g[ig] = -vlapl_g[ig] * chr->rhopw->gg[ig] * tpiba2; - - std::vector> aux(chr->rhopw->nmaxgr); - chr->rhopw->recip2real(vlapl_lapl_g.data(), aux.data()); - - for( int ir=0; ir< nrxx; ++ir ) + static bool vlapl_warning_printed = false; + if (!vlapl_warning_printed) { -#ifdef __EXX - double vlapl_contrib = aux[ir].real(); - if (func.info->number == XC_MGGA_X_SCAN && XC_Functional::get_func_type() == 5) + 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-10) { - vlapl_contrib *= (1.0 - XC_Functional::get_hybrid_alpha()); + vlapl_warning_printed = true; + ModuleBase::WARNING("XC_Functional_Libxc::v_xc_meta", + "The selected functional depends on the Laplacian of the density. " + "The vlapl contribution to V_xc is omitted because including it " + "(e2*nabla^2(vlapl)) causes SCF divergence due to |G|^2 amplification. " + "The self-consistent density is NOT the exact functional density. " + "The stress includes the vlapl term from the full energy functional. " + "Results may be unreliable, especially for norm-conserving pseudopotentials. " + "Consider using r2SCAN or SCAN instead."); } - v(is,ir) += vlapl_contrib; -#else - v(is,ir) += aux[ir].real(); -#endif } } } From 5157d1c716951b2fd41c444ff5fea5cfa4f64bd4 Mon Sep 17 00:00:00 2001 From: dyzheng Date: Thu, 30 Apr 2026 18:20:24 +0800 Subject: [PATCH 3/5] Docs: update FD validation results with correct unit conversion Add ecutwfc 60/80/100 Ry PW and LCAO FD test results for scheme 3. All schemes that omit vlapl from V show 5-10% stress error due to etxc-vtxc mismatch. Update summary table accordingly. --- docs/scanl_laplacian_implementation_report.md | 24 +++++++++++++++---- 1 file changed, 19 insertions(+), 5 deletions(-) diff --git a/docs/scanl_laplacian_implementation_report.md b/docs/scanl_laplacian_implementation_report.md index dbb48fb9d4a..bd29f035745 100644 --- a/docs/scanl_laplacian_implementation_report.md +++ b/docs/scanl_laplacian_implementation_report.md @@ -141,7 +141,7 @@ FD/AB ≈ 1.04–1.05，存在系统性 4–5% 误差。SCANL 误差随 ecutwfc --- -#### 方案 3：V 和 vtxc 均不含 vlapl，stress 保留 vlapl 项 +#### 方案 3：V 和 vtxc 均不含 vlapl，stress 保留 vlapl 项（当前实现） **实现**： - `v_xc_libxc`：传真实 ∇²ρ 给 LibXC，获取正确的 exc、vrho、vsigma、vtau @@ -151,9 +151,23 @@ FD/AB ≈ 1.04–1.05，存在系统性 4–5% 误差。SCANL 误差随 ecutwfc **物理逻辑**：vlapl 势无法加入 V（会发散），因此 vtxc 也不含 vlapl 以保持一致。etxc 中包含完整的 exc（含 Laplacian 贡献）。gradcorr 中的 vlapl 应力项从完整能量泛函严格推导。 -**ecutwfc=200 结果**（此前测试，单位转换有误，仅参考趋势）：FD/AB = 1.018（1.82% 误差） +**FD 测试结果**： + +PW 基组： + +| ecutwfc (Ry) | FD (kbar) | AB (kbar) | FD/AB | 误差 | +|-------------|-----------|-----------|-------|------| +| 60 | 260.22 | 283.88 | 0.917 | -8.3% | +| 80 | 254.26 | 280.80 | 0.906 | -9.5% | +| 100 | 270.06 | 278.70 | 0.969 | -3.1% | + +LCAO 基组（ecutwfc=100 Ry）： + +| 基组 | FD (kbar) | AB (kbar) | FD/AB | 误差 | +|------|-----------|-----------|-------|------| +| LCAO | 248.30 | 276.09 | 0.899 | -10.1% | -**ecutwfc=400 结果**：待验证（需在对应分支上重新构建和测试）。 +**误差分析**：AB 应力系统性大于 FD 应力约 8–10%。根本原因是 V 不含 vlapl → 密度在近似势下自洽（接近 SCAN 密度）→ etxc 用 SCANL 的 exc 评估（含 Laplacian 贡献）→ etxc-vtxc 不匹配 → -(etxc-vtxc)/Ω 偏大 → 对角应力偏大。这是省略 vlapl 势的根本性限制，无法通过调整 vlapl stress 项修正。 --- @@ -204,12 +218,12 @@ FD/AB ≈ 1.04–1.05，存在系统性 4–5% 误差。SCANL 误差随 ecutwfc | SCAN 基线 | N/A | N/A | N/A | 收敛 | 1.000 | 0.00% | | 方案 1 (严格) | ✓ | ✓ | ✓ | **发散** | — | — | | 方案 2 (sigma-as-lapl) | ✗ (假) | ✗ (假) | ✗ | 收敛 | 1.046 | 4.6% | -| 方案 3 (严格应力) | ✗ | ✗ | ✓ | 收敛 | ~1.02 | ~2% | +| 方案 3 (严格应力) | ✗ | ✗ | ✓ | 收敛 | 0.91 | -9% | | 方案 4 (vtxc含vlapl) | ✗ | ✓ | ✓ | 收敛 | 0.922 | -7.8% | | 方案 5 (衰减) | ✓(衰减) | ✓(衰减) | ✓ | 收敛 | ~0.9 | ~10% | | 方案 6 (ramp) | ✓(渐增) | ✓ | ✓ | **发散** | — | — | -注：方案 3、4、5 的 FD/AB 值来自 ecutwfc=200 Ry 的测试，由于此前存在单位转换错误，这些数值仅反映趋势，需用正确转换因子重新验证。 +**结论**：所有能收敛的方案都存在 5–10% 的应力误差。这是省略 vlapl 势的根本性后果——密度不在 SCANL 自洽势下产生，导致 etxc-vtxc 不匹配，进而使应力偏大。方案 2（sigma-as-lapl）的误差最小（4.6%），但那是传假输入给 LibXC 导致的偶然偏差，并非物理上更正确。 --- From d5c2e814460621174a8b1ce9ff60730d447e30e6 Mon Sep 17 00:00:00 2001 From: dyzheng Date: Thu, 30 Apr 2026 21:46:44 +0800 Subject: [PATCH 4/5] Feature(module_xc): implement FD Laplacian kernel for SCANL vlapl potential MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Replace spectral Laplacian (|G|² amplification → SCF divergence) with finite-difference Laplacian kernel in G-space. The FD kernel matches the spectral kernel at low G but is bounded at high G, enabling SCF convergence while preserving density-potential consistency. Key changes: - xc_functional_libxc_vxc.cpp: compute ∇²(vlapl·sgn) using FD kernel in G-space with metric tensor GGT for non-orthogonal cells - xc_functional.cpp: register SCANL functional name (XC_MGGA_X_SCANL + XC_MGGA_C_SCANL) - Update vtxc to include vlapl contribution for stress consistency FD stress validation (Si2 FCC, Γ-point): SCAN: FD/AB = 1.0000 (all ecutwfc) SCANL: FD/AB = 0.995 (ecutwfc=80, -0.51% error) vs scheme 3 (no vlapl in V): FD/AB = 0.906 (-9.5% error) --- docs/scanl_laplacian_implementation_report.md | 152 ++++++++++++++---- .../module_xc/xc_functional.cpp | 9 +- .../module_xc/xc_functional_libxc_vxc.cpp | 108 ++++++++++--- 3 files changed, 214 insertions(+), 55 deletions(-) diff --git a/docs/scanl_laplacian_implementation_report.md b/docs/scanl_laplacian_implementation_report.md index bd29f035745..e869a043652 100644 --- a/docs/scanl_laplacian_implementation_report.md +++ b/docs/scanl_laplacian_implementation_report.md @@ -211,19 +211,112 @@ LCAO 基组（ecutwfc=100 Ry）： --- -### 3.5 测试结果汇总 +#### 方案 7：有限差分 Laplacian（FD kernel in G-space） + +**核心思想**：谱 Laplacian ∇²f(G) = -|G|²f(G) 中的 |G|² 是无界的，导致高频放大和 SCF 发散。有限差分（FD）Laplacian 在低 G 时与谱 Laplacian 一致，但在高 G 时自然衰减，最大放大因子有界。 + +**推导**： + +在分数坐标 u^α = n_α/N_α（α = 1,2,3）下，Laplacian 表示为： + +$$\nabla^2 = \sum_{\alpha,\beta} g^{\alpha\beta} \frac{\partial^2}{\partial u^\alpha \partial u^\beta}$$ + +其中 g^{αβ} = GGT[αβ]/lat0² 是逆变度规张量，GGT 是倒格子度规矩阵（ABACUS 中的 `rhopw->GGT`）。 + +对平面波 f(r) = exp(2πi Σ m_γ u^γ)，谱 Laplacian 给出： + +$$\frac{\partial^2 f}{\partial u^\alpha \partial u^\beta} = -(2\pi)^2 m_\alpha m_\beta \, f$$ + +FD 近似给出： + +- **对角项** (α = β)：中心差分 + +$$\left(\frac{\partial^2 f}{\partial (u^\alpha)^2}\right)_{\text{FD}} = -2N_\alpha^2 \left(1 - \cos\frac{2\pi m_\alpha}{N_\alpha}\right) \, f$$ + +- **交叉项** (α ≠ β)：四点差分 + +$$\left(\frac{\partial^2 f}{\partial u^\alpha \partial u^\beta}\right)_{\text{FD}} = -N_\alpha N_\beta \sin\frac{2\pi m_\alpha}{N_\alpha} \sin\frac{2\pi m_\beta}{N_\beta} \, f$$ + +对低频模式 (m_α ≪ N_α)，FD 退化为谱形式（1-cos(x) ≈ x²/2, sin(x) ≈ x）。对 Nyquist 频率 (m_α = N_α/2)，放大因子为有界值 N_α² · 4/(2π)² ≈ 0.405 N_α²，而非谱 Laplacian 的 (N_α/2)²。 + +组合所有项，FD Laplacian 的等效 G² 核为： + +$$\text{gg}_{\text{FD}} = \frac{1}{(2\pi)^2} \sum_{\alpha,\beta} \text{GGT}[\alpha][\beta] \cdot \text{FD}_{\alpha\beta}$$ + +其中： + +$$\text{FD}_{\alpha\alpha} = 2N_\alpha^2(1-\cos\frac{2\pi m_\alpha}{N_\alpha}), \quad \text{FD}_{\alpha\beta}^{(\alpha\neq\beta)} = N_\alpha N_\beta \sin\frac{2\pi m_\alpha}{N_\alpha} \sin\frac{2\pi m_\beta}{N_\beta}$$ + +最终，∇²_FD(vlapl·sgn) 在 G 空间中为： + +$$\widetilde{\nabla^2_{\text{FD}}(\text{vlapl}\cdot\text{sgn})}(\mathbf{G}) = -\text{gg}_{\text{FD}}(\mathbf{G}) \cdot \text{tpiba2} \cdot \widetilde{\text{vlapl}\cdot\text{sgn}}(\mathbf{G})$$ + +**实现**： +1. 计算 vlapl_sgn = vlapl·sgn（实空间） +2. FFT: vlapl_sgn → vlapl_g +3. 应用 FD kernel: vlapl_g_lapl[ig] = -vlapl_g[ig] · gg_FD[ig] · tpiba2 +4. IFFT: vlapl_g_lapl → vlapl_lapl（实空间） +5. 加入 V_xc 和 vtxc + +**交叉项的物理意义**：对于非正交晶胞（如 FCC 原胞），度规张量 GGT 的非对角元不为零。此时，沿一个晶格方向的位移会影响其他方向的梯度分量。FD kernel 的交叉项 sin·sin 正确反映了这一点。对于正交晶胞（GGT 对角），交叉项消失，FD kernel 退化为三个方向的独立贡献。 + +**数值验证**（FCC 原胞, a=10.2 Bohr, nx=ny=nz=36）： + +| Miller 指数 | gg_spectral | gg_FD | gg_FD/gg_spectral | +|-------------|-------------|-------|-------------------| +| (0,0,0) | 0 | 0 | 1.000 | +| (1,0,0) | 3.0 | 2.992 | 0.997 | +| (1,1,0) | 4.0 | 4.005 | 1.001 | +| (5,0,0) | 75.0 | 70.36 | 0.938 | +| (10,0,0) | 300.0 | 231.2 | 0.771 | +| (18,0,0) | 972.0 | 393.9 | 0.405 | + +FD kernel 在低 G 时与谱 kernel 完全一致（<1% 偏差），在高 G 时衰减到谱值的 ~40%，有界且不会导致发散。 + +**SCF 收敛性**：使用 Pulay mixing (beta=0.2)，SCANL 在 ecutwfc=80 Ry 下 SCF 正常收敛。 + +**FD 应力验证**（Si₂ FCC 原胞, a₀=10.2 Bohr, δ=0.001, Pulay mixing beta=0.1–0.2）： + +PW 基组： + +| 泛函 | ecutwfc (Ry) | σ_FD (kbar) | σ_AB (kbar) | FD/AB | 误差 | +|------|-------------|-------------|-------------|-------|------| +| SCAN | 60 | 434.29 | 434.30 | 1.0000 | 0.00% | +| SCAN | 80 | 434.47 | 434.48 | 1.0000 | 0.00% | +| SCAN | 100 | 434.54 | 434.55 | 1.0000 | 0.00% | +| SCANL | 60 | 431.89 | 435.60 | 0.9915 | -0.85% | +| SCANL | 80 | 432.37 | 434.58 | 0.9949 | -0.51% | + +**注意**：ecutwfc=100 时 SCANL 的 SCF 收敛较困难（能量振荡较大），FD 应力结果不稳定，未列入上表。建议 SCANL 使用 ecutwfc ≤ 80 Ry。 + +**与方案 3 的对比**： + +| 方案 | V 含 vlapl | vtxc 含 vlapl | SCF | FD/AB | 误差 | +|------|-----------|--------------|-----|-------|------| +| 方案 3 (无 vlapl 势) | ✗ | ✗ | 收敛 | 0.906 | -9.5% | +| 方案 7 (FD Laplacian) | ✓ (FD) | ✓ (FD) | 收敛 | 0.995 | -0.51% | + +**FD Laplacian 方案将应力误差从 9.5% 降至 0.51%，改善了约 20 倍。** + +残余 0.5% 误差的可能来源： +1. FD kernel 在高 G 区与谱 kernel 的偏差（已衰减但仍有差异） +2. 有限差分应力本身的高阶误差（δ=0.001 的中心差分为二阶精度） +3. vlapl stress 项中使用的 Hessian 仍为谱 Hessian + +--- | 方案 | V 含 vlapl | vtxc 含 vlapl | stress 含 vlapl | SCF | FD/AB | 误差 | |------|-----------|--------------|----------------|-----|-------|------| | SCAN 基线 | N/A | N/A | N/A | 收敛 | 1.000 | 0.00% | | 方案 1 (严格) | ✓ | ✓ | ✓ | **发散** | — | — | | 方案 2 (sigma-as-lapl) | ✗ (假) | ✗ (假) | ✗ | 收敛 | 1.046 | 4.6% | -| 方案 3 (严格应力) | ✗ | ✗ | ✓ | 收敛 | 0.91 | -9% | +| 方案 3 (严格应力) | ✗ | ✗ | ✓ | 收敛 | 0.906 | -9.5% | | 方案 4 (vtxc含vlapl) | ✗ | ✓ | ✓ | 收敛 | 0.922 | -7.8% | | 方案 5 (衰减) | ✓(衰减) | ✓(衰减) | ✓ | 收敛 | ~0.9 | ~10% | | 方案 6 (ramp) | ✓(渐增) | ✓ | ✓ | **发散** | — | — | +| **方案 7 (FD Laplacian)** | **✓ (FD)** | **✓ (FD)** | **✓** | **收敛** | **0.995** | **-0.5%** | -**结论**：所有能收敛的方案都存在 5–10% 的应力误差。这是省略 vlapl 势的根本性后果——密度不在 SCANL 自洽势下产生，导致 etxc-vtxc 不匹配，进而使应力偏大。方案 2（sigma-as-lapl）的误差最小（4.6%），但那是传假输入给 LibXC 导致的偶然偏差，并非物理上更正确。 +**结论**：方案 7（FD Laplacian）是唯一同时满足 SCF 收敛和应力精度的方案。它通过用有界的 FD kernel 替代无界的 |G|² 谱 Laplacian，既解决了 SCF 发散问题，又将 vlapl 势纳入 V_xc 使得密度-势自洽，从而大幅改善应力精度（从 9.5% 误差降至 0.5%）。 --- @@ -270,42 +363,37 @@ VASP 使用 PAW 方法，ecutwfc 要求较低（~400–700 eV），可以在 800 ## 5. 结论与建议 -### 5.1 不推荐使用 SCANL +### 5.1 推荐方案：FD Laplacian（方案 7） -**在 ABACUS（NCPP + 平面波）框架下不推荐使用 SCANL 及其他 ∇²ρ 依赖泛函**。原因： +**FD Laplacian 方案是当前最佳方案**，它： -1. **严格实现不可行**：vlapl 势 e2·∇²(vlapl) 在任何实际 ecutwfc 下导致 SCF 发散 -2. **NCPP 矛盾**：模守恒赝势需要 ≥60 Ry 的 ecutwfc，此范围内 vlapl 势必然不稳定 -3. **近似方案不可靠**：无论省略 vlapl 势（方案 3）还是用 sigma 伪装（方案 2），所得密度都不是 SCANL 自洽密度，结果与 FD 验证存在 2–5% 的偏差 -4. **VASP 也承认问题**：即使在 PAW 方法下，VASP 也建议不要使用超过 800 eV 的截断 +1. **SCF 收敛**：FD kernel 有界，不会导致 |G|² 发散 +2. **密度-势自洽**：V 包含 vlapl 贡献（通过 FD Laplacian），密度在近似自洽势下产生 +3. **应力精度高**：FD/AB = 0.995，误差仅 0.5%（远优于方案 3 的 9.5%） +4. **适用于非正交晶胞**：通过度规张量 GGT 正确处理交叉项 +5. **无需调参**：不依赖衰减因子或 ramp 步数等人为参数 -### 5.2 推荐替代泛函 +**残余 0.5% 误差**是 FD kernel 与谱 kernel 在高 G 区差异的代价，对于实际应用完全可以接受。 -| 泛函 | 依赖变量 | 特点 | -|------|---------|------| -| **r²SCAN** | ρ, ∇ρ, τ | SCAN 的正则化版本，数值稳定，推荐首选 | -| **SCAN** | ρ, ∇ρ, τ | 与 SCANL 精度相近，无 ∇²ρ 依赖，数值稳定 | +### 5.2 使用注意事项 + +1. **FD Laplacian 是对谱 Laplacian 的近似**：在 FFT 网格上，它是"正确"的 Laplacian（对网格函数的有限差分），但与连续 Laplacian 有差异 +2. **需要 mixing_tau = 1**：meta-GGA 泛函需要混合动能密度 +3. **建议使用 Pulay mixing**：mixing_beta=0.2，对 SCANL 效果较好 +4. **仍建议优先考虑 r²SCAN**：r²SCAN 无 ∇²ρ 依赖，数值更稳定 -r²SCAN 是目前推荐的 meta-GGA 泛函，兼具 SCAN 的精度和更好的数值稳定性，且无 ∇²ρ 依赖问题。 +### 5.3 代码实现 -### 5.3 代码实现建议 +SCANL 已注册为 `dft_functional scanl`，内部映射为 `XC_MGGA_X_SCANL + XC_MGGA_C_SCANL`。 -若仍需保留 SCANL 支持（例如用于与文献对比），建议采用**方案 3 + 运行时警告**： +vlapl 势项实现于 `xc_functional_libxc_vxc.cpp` 中的 `v_xc_meta` 函数，使用 G 空间 FD kernel 计算 ∇²(vlapl·sgn)。 -1. **V_xc 不含 vlapl**：e2·∇²(vlapl) 导致 |G|² 发散，无法加入 -2. **vtxc 不含 vlapl**：必须与 V 保持一致 -3. **etxc 含完整 exc**：LibXC 返回的能量密度已含 Laplacian 贡献 -4. **gradcorr stress 含 2·vlapl·Hess·e2**：从完整能量泛函严格推导 -5. **运行时打印 WARNING** +vlapl 应力项实现于 `xc_functional_gradcorr.cpp` 中的 `gradcorr` 函数。 -在用户选择 ∇²ρ 依赖泛函时，应打印类似以下警告： +### 5.4 替代泛函 -``` -WARNING: The selected functional depends on the Laplacian of the density (∇²ρ). -The vlapl contribution to the Kohn-Sham potential is omitted because including -it (e2·∇²(vlapl)) causes SCF divergence due to |G|² amplification in reciprocal -space. The self-consistent density is therefore NOT the exact SCANL density. -The stress includes the vlapl term derived from the full energy functional. -Results may be unreliable, especially for norm-conserving pseudopotentials -which require high ecutwfc. Consider using r2SCAN or SCAN instead. -``` +| 泛函 | 依赖变量 | 特点 | +|------|---------|------| +| **r²SCAN** | ρ, ∇ρ, τ | SCAN 的正则化版本，数值稳定，推荐首选 | +| **SCAN** | ρ, ∇ρ, τ | 无 ∇²ρ 依赖，数值稳定 | +| **SCANL** | ρ, ∇ρ, ∇²ρ | 需 FD Laplacian 近似，0.5% 应力误差 | diff --git a/source/module_hamilt_general/module_xc/xc_functional.cpp b/source/module_hamilt_general/module_xc/xc_functional.cpp index 2fc74ffbb78..e8758916d4f 100644 --- a/source/module_hamilt_general/module_xc/xc_functional.cpp +++ b/source/module_hamilt_general/module_xc/xc_functional.cpp @@ -199,6 +199,13 @@ void XC_Functional::set_xc_type(const std::string xc_func_in) func_type = 5; use_libxc = true; } + 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; + } else if( xc_func == "LC_PBE") { func_id.push_back(XC_HYB_GGA_XC_LC_PBEOP); @@ -350,7 +357,7 @@ void XC_Functional::set_xc_type(const std::string xc_func_in) #ifndef USE_LIBXC if(xc_func == "SCAN" || xc_func == "HSE" || xc_func == "SCAN0" - || xc_func == "MULLER" || xc_func == "POWER" || xc_func == "WP22" || xc_func == "CWP22" || + || xc_func == "SCANL" || xc_func == "MULLER" || xc_func == "POWER" || xc_func == "WP22" || xc_func == "CWP22" || xc_func == "LC_PBE" || xc_func == "LC_WPBE" || xc_func == "LRC_WPBE" || xc_func == "LRC_PBEH" || xc_func == "CAM_PBEH") { diff --git a/source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp b/source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp index c5edc108ec9..6c6558615bd 100644 --- a/source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp +++ b/source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp @@ -392,33 +392,97 @@ std::tuple XC_Functional_Li } } - //process vlapl: omitted from V_xc and vtxc for SCF stability - // V_xc^vlapl = e2·∇²(vlapl·sgn) is excluded because it causes - // |G|² amplification in reciprocal space leading to SCF divergence. - // vtxc_vlapl is also excluded for consistency with the actual potential. - // The vlapl energy is already included in etxc via exc. - // The vlapl stress (σ = 2·vlapl·Hess·e2) is computed in gradcorr. - // Note: LibXC vlapl = ∂(ρε)/∂(∇²ρ) already includes ρ factor. + //process vlapl: ∇²(vlapl) contribution to potential + // Use finite-difference Laplacian in G-space to avoid unbounded + // |G|² amplification from the spectral Laplacian. The FD kernel + // replaces m_α² with n_α²·2(1-cos(2πm_α/n_α)) and m_α·m_β + // with n_α·n_β·sin(2πm_α/n_α)·sin(2πm_β/n_β), which matches + // the spectral kernel at low G but is bounded at high G. + // Cross terms from the metric tensor (GGT) are included for + // non-orthogonal cells. + for( int is=0; is vlapl_sgn(nrxx); +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 1024) +#endif + for( int ir=0; ir> vlapl_g(chr->rhopw->npw); + chr->rhopw->real2recip(vlapl_sgn.data(), vlapl_g.data()); + + const int ngrid[3] = {chr->rhopw->nx, chr->rhopw->ny, chr->rhopw->nz}; + const ModuleBase::Matrix3& GGT_mat = chr->rhopw->GGT; + const double two_pi = ModuleBase::TWO_PI; + const double ggt[3][3] = { + {GGT_mat.e11, GGT_mat.e12, GGT_mat.e13}, + {GGT_mat.e21, GGT_mat.e22, GGT_mat.e23}, + {GGT_mat.e31, GGT_mat.e32, GGT_mat.e33} + }; + + const double four_pi_sq = ModuleBase::FOUR_PI * ModuleBase::PI; + std::vector> vlapl_g_lapl(chr->rhopw->npw); +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 256) +#endif + for(int ig=0; igrhopw->npw; ++ig) + { + double gg_fd_raw = 0.0; + const double m[3] = {chr->rhopw->gdirect[ig].x, + chr->rhopw->gdirect[ig].y, + chr->rhopw->gdirect[ig].z}; + for(int alpha=0; alpha<3; ++alpha) + { + const double phi_a = two_pi * m[alpha] / ngrid[alpha]; + for(int beta=0; beta<3; ++beta) + { + const double phi_b = two_pi * m[beta] / ngrid[beta]; + if(alpha == beta) + { + gg_fd_raw += ggt[alpha][alpha] * ngrid[alpha] * ngrid[alpha] * 2.0 * (1.0 - cos(phi_a)); + } + else + { + gg_fd_raw += ggt[alpha][beta] * ngrid[alpha] * ngrid[beta] * sin(phi_a) * sin(phi_b); + } + } + } + const double gg_fd = gg_fd_raw / four_pi_sq; + vlapl_g_lapl[ig] = -vlapl_g[ig] * gg_fd * tpiba2; + } + + std::vector> aux(chr->rhopw->nmaxgr); + chr->rhopw->recip2real(vlapl_g_lapl.data(), aux.data()); + +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 1024) +#endif + for( int ir=0; ir 1e-10) +#ifdef __EXX + double vlapl_contrib = aux[ir].real(); + if (func.info->number == XC_MGGA_X_SCAN && XC_Functional::get_func_type() == 5) { - vlapl_warning_printed = true; - ModuleBase::WARNING("XC_Functional_Libxc::v_xc_meta", - "The selected functional depends on the Laplacian of the density. " - "The vlapl contribution to V_xc is omitted because including it " - "(e2*nabla^2(vlapl)) causes SCF divergence due to |G|^2 amplification. " - "The self-consistent density is NOT the exact functional density. " - "The stress includes the vlapl term from the full energy functional. " - "Results may be unreliable, especially for norm-conserving pseudopotentials. " - "Consider using r2SCAN or SCAN instead."); + vlapl_contrib *= (1.0 - XC_Functional::get_hybrid_alpha()); } + v(is,ir) += ModuleBase::e2 * vlapl_contrib; +#else + v(is,ir) += ModuleBase::e2 * aux[ir].real(); +#endif + } + + double vtxc_vlapl = 0.0; +#ifdef _OPENMP +#pragma omp parallel for reduction(+:vtxc_vlapl) schedule(static, 1024) +#endif + for( int ir=0; irrho[is][ir] * sgn[ir*nspin+is]; } + vtxc += vtxc_vlapl; } } From e6915dbf4b6dcc7201653cfa5241abfd055a18b5 Mon Sep 17 00:00:00 2001 From: dyzheng Date: Thu, 30 Apr 2026 21:52:12 +0800 Subject: [PATCH 5/5] Docs: rewrite SCANL FD Laplacian report with detailed formulas and test data Add complete derivation of FD Laplacian kernel in G-space for general (non-orthogonal) cells, including metric tensor cross terms. Document numerical verification of gg_FD vs gg_spectral at multiple FFT grid sizes, full FD stress validation results with energy data, and error analysis. --- docs/scanl_laplacian_implementation_report.md | 452 ++++++++---------- 1 file changed, 193 insertions(+), 259 deletions(-) diff --git a/docs/scanl_laplacian_implementation_report.md b/docs/scanl_laplacian_implementation_report.md index e869a043652..35b5d493ffc 100644 --- a/docs/scanl_laplacian_implementation_report.md +++ b/docs/scanl_laplacian_implementation_report.md @@ -1,399 +1,333 @@ -# SCANL Laplacian 依赖项实现与验证报告 +# SCANL 有限差分 Laplacian 实现 -## 1. 背景 +## 1. 问题 -SCANL（SCAN-L）是 SCAN 的 deorbitalized 版本，将 SCAN 中对动能密度 τ 的依赖替换为对密度 Laplacian ∇²ρ 的依赖 [Mejia-Rodriguez & Trickey, Phys. Rev. B **98**, 115161 (2018)]。这导致其 Kohn-Sham 势和应力公式中出现了 ∇²(∂ε/∂(∇²ρ)) 项，在平面波基组下引发严重的数值稳定性问题。 +SCANL (deorbitalized SCAN) 将 SCAN 对动能密度 τ 的依赖替换为对密度 Laplacian ∇²ρ 的依赖。其 Kohn-Sham 势包含项 -本文档记录了 ABACUS 中 SCANL 泛函 Laplacian 相关项的实现方案、有限差分（FD）应力验证结果，以及最终建议。 +$$V_{xc}^{\nabla^2\rho}(\mathbf{r}) = e^2\,\nabla^2\!\left(\frac{\partial(\rho\varepsilon_{xc})}{\partial(\nabla^2\rho)}\right) \equiv e^2\,\nabla^2(\text{vlapl})$$ ---- - -## 2. 理论公式 - -### 2.1 meta-GGA 能量泛函 +在平面波基组下，倒易空间中 ∇² → −|G|²，因此 -meta-GGA 交换关联能量密度 ε_xc 依赖于四个变量： - -$$E_{xc} = \int \rho(\mathbf{r})\,\varepsilon_{xc}(\rho, |\nabla\rho|^2, \nabla^2\rho, \tau)\, d\mathbf{r}$$ - -对于 SCANL，∂ε/∂(∇²ρ) ≠ 0（与 SCAN 不同，SCAN 的此项为零）。 +$$\tilde{V}_{xc}^{\nabla^2\rho}(\mathbf{G}) = -e^2\,|\mathbf{G}|^2\,\widetilde{\text{vlapl}}(\mathbf{G})$$ -### 2.2 Kohn-Sham 势 +|G|² 无界增长导致 SCF 发散：vlapl(G) 的高频数值噪声被 |G|² 放大后反馈到密度，形成正反馈。此发散不依赖 mixing 策略，是 Laplace 逆变换的数学病态性。 -由变分原理，XC 势为： +--- -$$V_{xc} = e^2 \frac{\delta E_{xc}}{\delta \rho} = e^2 \left[ \frac{\partial(\rho\varepsilon)}{\partial\rho} - \nabla \cdot \frac{\partial(\rho\varepsilon)}{\partial(\nabla\rho)} + \nabla^2 \frac{\partial(\rho\varepsilon)}{\partial(\nabla^2\rho)} + \frac{\partial(\rho\varepsilon)}{\partial\tau} \hat{\tau} \right]$$ +## 2. 有限差分 Laplacian -其中 **vlapl 势项**为： +### 2.1 核心思想 -$$V_{xc}^{\nabla^2\rho} = e^2\,\nabla^2\!\left(\frac{\partial(\rho\varepsilon)}{\partial(\nabla^2\rho)}\right)$$ +FFT 网格上的密度是离散函数。对离散函数求 Laplacian，有限差分 (FD) 是自然且唯一自洽的方法。FD Laplacian 的 G 空间 kernel 在低 G 时与谱 kernel 一致，在高 G 时有界衰减，不会放大高频噪声。 -在倒易空间中，∇² → -|G|²，因此： +### 2.2 坐标系与度规 -$$\tilde{V}_{xc}^{\nabla^2\rho}(\mathbf{G}) = -e^2\,|\mathbf{G}|^2\,\widetilde{\text{vlapl}}(\mathbf{G})$$ +设晶格矢量 $\mathbf{a}_\alpha$（α = 1,2,3），以 lat0 为单位存储在 `rhopw->latvec` 中。分数坐标 -**关键问题**：|G|² 是高通滤波器，放大高频模式。vlapl(G) 中的数值噪声被 |G|² 放大后反馈到密度中，形成正反馈循环导致 SCF 发散。ecutwfc 越大，|G|_max 越大，放大效应越强，SCF 越容易发散。这并非 mixing 策略的问题，而是 Laplace 逆变换的数学病态性。 +$$u^\alpha = n_\alpha / N_\alpha, \quad n_\alpha = 0, 1, \ldots, N_\alpha - 1$$ -### 2.3 vtxc 与应力公式 +其中 $N_\alpha$ = `rhopw->nx, ny, nz`。实空间位置 -vtxc（势对密度的积分）在应力计算中至关重要： +$$\mathbf{r} = \text{lat0} \sum_\alpha u^\alpha \mathbf{a}_\alpha$$ -$$\text{vtxc} = \int V_{xc}(\mathbf{r})\,\rho(\mathbf{r})\, d\mathbf{r}$$ +Laplacian 在分数坐标下为 -应力中 XC 对角贡献为： +$$\nabla^2 = \sum_{\alpha,\beta} g^{\alpha\beta} \frac{\partial^2}{\partial u^\alpha \partial u^\beta}$$ -$$\sigma_{ii}^{xc} = -\frac{E_{txc} - V_{txc}}{\Omega}$$ +其中逆变度规张量 -其中 $E_{txc} = \int e^2\,\rho\,\varepsilon_{xc}\, d\mathbf{r}$。 +$$g^{\alpha\beta} = \frac{\text{GGT}[\alpha][\beta]}{\text{lat0}^2}$$ -**vtxc 必须与实际使用的 V_xc 一致**。若 V 中不含 vlapl 而 vtxc 中含 vlapl，则 $-(E_{txc}-V_{txc})/\Omega$ 不正确。 +GGT = G · Gᵀ 是倒格子度规矩阵（ABACUS 中的 `rhopw->GGT`），G = (latvec)⁻ᵀ。 -### 2.4 vlapl 应力项 +### 2.3 平面波的谱 Laplacian -由能量泛函对应变的变分推导，vlapl 对应力的贡献为： +对平面波 $f(\mathbf{r}) = \exp(2\pi i \sum_\gamma m_\gamma u^\gamma)$，其中 $m_\gamma$ 为 Miller 指数（`rhopw->gdirect[ig]`）： -$$\sigma_{\alpha\beta}^{\nabla^2\rho} = \frac{2}{\Omega} \sum_{\mathbf{r}} \frac{\partial(\rho\varepsilon)}{\partial(\nabla^2\rho)}\, H_{\alpha\beta}(\mathbf{r})\, e^2$$ +$$\frac{\partial^2 f}{\partial u^\alpha \partial u^\beta}\bigg|_{\text{spectral}} = -(2\pi)^2 m_\alpha m_\beta \, f$$ -其中 $H_{\alpha\beta} = \partial^2\rho / \partial r_\alpha \partial r_\beta$ 是密度的 Hessian 矩阵。 +组合得谱 Laplacian： -**注意**：LibXC 返回的 vlapl = ∂(ρε)/∂(∇²ρ) 已包含 ρ 因子，无需再乘 ρ。 +$$\nabla^2_{\text{spectral}} f = -\underbrace{\left[\sum_{\alpha,\beta} \text{GGT}[\alpha][\beta] \cdot m_\alpha m_\beta\right]}_{\text{gg}_{\text{spectral}}} \cdot \text{tpiba2} \cdot f$$ ---- +其中 tpiba2 = (2π/lat0)²。这就是 ABACUS 中 `gg[ig] * tpiba2` 的来源。 -## 3. 实现方案与测试结果 +### 2.4 平面波的 FD Laplacian -### 3.1 验证方法 +在分数坐标的均匀网格上对二阶导数做中心差分： -采用有限差分（FD）应力验证：对 Si2 FCC 原胞（Γ 点）在 a₀±δ 下计算总能量，由中心差分公式得到应力： +**对角项** (α = β)： -$$\sigma_{FD} = -\frac{E(a_0+\delta) - E(a_0-\delta)}{V(a_0+\delta) - V(a_0-\delta)}$$ +$$\frac{\partial^2 f}{\partial (u^\alpha)^2}\bigg|_{\text{FD}} = \frac{f(u^\alpha + 1/N_\alpha) + f(u^\alpha - 1/N_\alpha) - 2f(u^\alpha)}{(1/N_\alpha)^2}$$ -与解析应力 σ_AB 比较，FD/AB 越接近 1.0 越好。 +对平面波： -**单位转换**：ABACUS 总能量单位为 eV，应力原子单位为 Ry/Bohr³。FD 应力转换公式： +$$= -2N_\alpha^2 \left(1 - \cos\frac{2\pi m_\alpha}{N_\alpha}\right) f$$ -$$\sigma_{FD}(\text{kbar}) = \sigma_{FD}(\text{eV/Bohr}^3) \times \frac{2\,\text{Ry}}{27.211386\,\text{eV}} \times 147105\,\frac{\text{kbar}}{\text{Ry/Bohr}^3}$$ +**交叉项** (α ≠ β)： -### 3.2 测试参数 +$$\frac{\partial^2 f}{\partial u^\alpha \partial u^\beta}\bigg|_{\text{FD}} = \frac{1}{4 N_\alpha N_\beta}\Big[f\Big(\frac{+1}{N_\alpha}, \frac{+1}{N_\beta}\Big) - f\Big(\frac{+1}{N_\alpha}, \frac{-1}{N_\beta}\Big) - f\Big(\frac{-1}{N_\alpha}, \frac{+1}{N_\beta}\Big) + f\Big(\frac{-1}{N_\alpha}, \frac{-1}{N_\beta}\Big)\Big]$$ -| 参数 | 值 | -|------|-----| -| 体系 | Si₂ FCC 原胞 | -| a₀ | 10.2 Bohr | -| δ | 0.001 Bohr | -| ecutwfc | 60 / 80 / 100 / 400 Ry | -| smearing_sigma | 0.02 | -| ks_solver | dav_subspace | -| mixing_beta | 0.1 | -| mixing_tau | 1 | -| scf_thr | 1e-7 | -| pw_seed | 1 | +对平面波： -### 3.3 SCAN 基线验证 +$$= -N_\alpha N_\beta \sin\frac{2\pi m_\alpha}{N_\alpha} \sin\frac{2\pi m_\beta}{N_\beta} \, f$$ -SCAN 不依赖 ∇²ρ（LibXC 返回 vlapl=0），用于验证 FD 方法和代码的正确性： +### 2.5 FD kernel -| ecutwfc (Ry) | ecutwfc (eV) | FD (kbar) | AB (kbar) | FD/AB | 误差 | -|-------------|-------------|-----------|-----------|-------|------| -| 60 | 816 | 265.69 | 265.66 | 1.0001 | 0.01% | -| 80 | 1088 | 266.11 | 266.11 | 1.0000 | 0.00% | -| 100 | 1361 | 266.11 | 266.12 | 0.9999 | -0.01% | -| 400 | 5442 | 266.11 | 266.09 | 1.0001 | 0.01% | +定义 FD 因子矩阵： -**结论**：SCAN 在所有 ecutwfc 下 FD/AB ≈ 1.000，验证了 FD 方法和应力代码的正确性。 +$$\text{FD}_{\alpha\alpha} = 2N_\alpha^2\left(1 - \cos\frac{2\pi m_\alpha}{N_\alpha}\right)$$ -### 3.4 各方案详情与结果 +$$\text{FD}_{\alpha\beta} = N_\alpha N_\beta \sin\frac{2\pi m_\alpha}{N_\alpha} \sin\frac{2\pi m_\beta}{N_\beta} \quad (\alpha \neq \beta)$$ ---- +FD Laplacian 的等效 gg 为： -#### 方案 1：完整 e2·∇²(vlapl) 加入势和 vtxc（严格实现） +$$\text{gg}_{\text{FD}} = \frac{1}{(2\pi)^2} \sum_{\alpha,\beta} \text{GGT}[\alpha][\beta] \cdot \text{FD}_{\alpha\beta}$$ -**实现**： -- V_xc：含 e2·∇²(vlapl·sgn) -- vtxc：含 e2·∫∇²(vlapl·sgn)·ρ·sgn -- stress：含 2·vlapl·Hess·e2 +最终： -**结果**：**SCF 发散**。无论 mixing_beta（0.01–0.7）、Kerker mixing、绝热演化（50 步 ramp）均无法收敛。 +$$\nabla^2_{\text{FD}} f = -\text{gg}_{\text{FD}} \cdot \text{tpiba2} \cdot f$$ -**原因**：V(G) = -e2·|G|²·vlapl(G)，|G|² 高通滤波放大高频噪声，正反馈增益超过任何混合策略的阻尼能力。ecutwfc 越大，|G|_max 越大，发散越严重。 +**与谱 kernel 的关系**：当 $m_\alpha / N_\alpha \to 0$ 时， ---- +$$2(1-\cos\theta) \to \theta^2, \quad \sin\theta \to \theta$$ -#### 方案 2：sigma 伪装 lapl（原始代码方案） +因此 $\text{FD}_{\alpha\alpha} \to (2\pi m_\alpha)^2$，$\text{FD}_{\alpha\beta} \to (2\pi)^2 m_\alpha m_\beta$，$\text{gg}_{\text{FD}} \to \text{gg}_{\text{spectral}}$。 -**实现**： -- `v_xc_libxc`：传 sigma（= |∇ρ|²）给 LibXC 的 Laplacian 参数位置 -- `tau_xc`/`tau_xc_spin`：传 grho（= |∇ρ|²）给 LibXC 的 Laplacian 参数位置 -- gradcorr stress：无 vlapl 应力项 +**有界性**：gg_FD 的最大值有限（由 FFT 网格大小决定），而 gg_spectral 随 |m| 无界增长。 -**问题**：LibXC 在错误的 Laplacian 输入下计算 exc 和 vlapl，导致 etxc 和应力均不正确。 +### 2.6 正交晶胞简化 -**FD 测试结果**： +对正交晶胞，GGT 为对角矩阵（GGT[αβ] = 0 for α ≠ β），交叉项消失： -| ecutwfc (Ry) | FD (kbar) | AB (kbar) | FD/AB | 误差 | -|-------------|-----------|-----------|-------|------| -| 60 | 249.90 | 240.91 | 1.037 | 3.73% | -| 80 | 252.32 | 241.11 | 1.047 | 4.65% | -| 100 | 252.39 | 241.12 | 1.047 | 4.68% | -| 400 | 252.39 | 241.24 | 1.046 | 4.62% | +$$\text{gg}_{\text{FD}} = \frac{1}{(2\pi)^2} \sum_\alpha \text{GGT}[\alpha][\alpha] \cdot 2N_\alpha^2\left(1 - \cos\frac{2\pi m_\alpha}{N_\alpha}\right)$$ -FD/AB ≈ 1.04–1.05，存在系统性 4–5% 误差。SCANL 误差随 ecutwfc 趋于稳定（~4.6%），说明这是传假 Laplacian 导致的系统偏差，而非 ecutwfc 收敛问题。 +每个方向独立贡献。 --- -#### 方案 3：V 和 vtxc 均不含 vlapl，stress 保留 vlapl 项（当前实现） +## 3. 实现 -**实现**： -- `v_xc_libxc`：传真实 ∇²ρ 给 LibXC，获取正确的 exc、vrho、vsigma、vtau -- V_xc：不含 e2·∇²(vlapl) -- vtxc：不含 vlapl 贡献（与 V 保持一致） -- gradcorr stress：含 2·vlapl·Hess·e2 +### 3.1 vlapl 势 -**物理逻辑**：vlapl 势无法加入 V（会发散），因此 vtxc 也不含 vlapl 以保持一致。etxc 中包含完整的 exc（含 Laplacian 贡献）。gradcorr 中的 vlapl 应力项从完整能量泛函严格推导。 +在 `v_xc_meta` 函数中，对每个自旋通道 is： -**FD 测试结果**： +1. **构造 vlapl·sgn**（实空间）： -PW 基组： +$$\text{vlapl\_sgn}[ir] = \text{vlapl}[ir \cdot \text{nspin} + \text{is}] \cdot \text{sgn}[ir \cdot \text{nspin} + \text{is}]$$ -| ecutwfc (Ry) | FD (kbar) | AB (kbar) | FD/AB | 误差 | -|-------------|-----------|-----------|-------|------| -| 60 | 260.22 | 283.88 | 0.917 | -8.3% | -| 80 | 254.26 | 280.80 | 0.906 | -9.5% | -| 100 | 270.06 | 278.70 | 0.969 | -3.1% | +2. **FFT 到 G 空间**： -LCAO 基组（ecutwfc=100 Ry）： +$$\widetilde{\text{vlapl\_sgn}} = \text{FFT}(\text{vlapl\_sgn})$$ -| 基组 | FD (kbar) | AB (kbar) | FD/AB | 误差 | -|------|-----------|-----------|-------|------| -| LCAO | 248.30 | 276.09 | 0.899 | -10.1% | +3. **应用 FD kernel**：对每个 G 向量 ig，Miller 指数 m = gdirect[ig]， -**误差分析**：AB 应力系统性大于 FD 应力约 8–10%。根本原因是 V 不含 vlapl → 密度在近似势下自洽（接近 SCAN 密度）→ etxc 用 SCANL 的 exc 评估（含 Laplacian 贡献）→ etxc-vtxc 不匹配 → -(etxc-vtxc)/Ω 偏大 → 对角应力偏大。这是省略 vlapl 势的根本性限制，无法通过调整 vlapl stress 项修正。 +$$\widetilde{\nabla^2_{\text{FD}}(\text{vlapl\_sgn})}[ig] = -\text{gg}_{\text{FD}}[ig] \cdot \text{tpiba2} \cdot \widetilde{\text{vlapl\_sgn}}[ig]$$ ---- +其中 gg_FD[ig] 由 2.5 节公式计算。 -#### 方案 4：势中无 vlapl，vtxc 中加 e2·vlapl·∇²ρ +4. **IFFT 回实空间**： -**实现**： -- V_xc：不含 vlapl -- vtxc：含 e2·∫vlapl·∇²ρ·sgn -- stress：含 2·vlapl·Hess·e2 +$$\nabla^2_{\text{FD}}(\text{vlapl\_sgn}) = \text{IFFT}(\widetilde{\nabla^2_{\text{FD}}(\text{vlapl\_sgn})})$$ -**问题**：V 无 vlapl 但 vtxc 有 → vtxc 与实际 V 不一致 → $-(E_{txc}-V_{txc})/\Omega$ 不正确 → 对角应力偏小。此前测试 FD/AB = 0.922（7.78% 误差）。 +5. **加入 V_xc**： ---- +$$V_{xc}(\text{is}, ir) \mathrel{+}= e^2 \cdot \nabla^2_{\text{FD}}(\text{vlapl\_sgn})[ir]$$ -#### 方案 5：高斯衰减 +对于 HSE 等 hybrid 泛函 (func_type == 5)，乘以 (1 − hybrid_alpha)。 -**实现**： -- V_xc：e2·∇²(vlapl·sgn·damp(G))，damp = exp(-α·G²/G_max²) -- vtxc：从衰减势计算或从未衰减 vlapl·∇²ρ 计算 -- stress：含 2·vlapl·Hess·e2（无衰减） +6. **加入 vtxc**： -**此前测试结果**（ecutwfc=200 Ry，仅参考趋势）： +$$\text{vtxc} \mathrel{+}= e^2 \sum_{ir} \nabla^2_{\text{FD}}(\text{vlapl\_sgn})[ir] \cdot \rho[\text{is}][ir] \cdot \text{sgn}[ir \cdot \text{nspin} + \text{is}]$$ -| α | vtxc 来源 | FD/AB | -|---|----------|-------| -| 1.0 | 衰减势 | 0.789 | -| 2.0 | 衰减势 | 0.960 | -| 4.0 | 衰减势 | 0.918 | +### 3.2 vlapl 应力 -**问题**：衰减改变了势的物理形状，密度在错误势下自洽。不同 α 给出不同结果，无法确定哪个是"正确的"。本质上是调参拟合，没有物理依据。 +vlapl 对应力的贡献在 `gradcorr` 函数中计算，与谱 Laplacian 方案相同： ---- +$$\sigma_{\alpha\beta}^{\nabla^2\rho} = \frac{2}{\Omega} \sum_{\mathbf{r}} \text{vlapl}(\mathbf{r}) \cdot H_{\alpha\beta}(\mathbf{r}) \cdot e^2$$ -#### 方案 6：绝热演化（ramp） +其中 $H_{\alpha\beta} = \partial^2\rho / \partial r_\alpha \partial r_\beta$ 是密度的 Hessian 矩阵（谱方法计算）。LibXC 返回的 vlapl = ∂(ρε)/∂(∇²ρ) 已包含 ρ 因子。 -**实现**： -- V_xc：e2·min(1, iter/ramp_steps)·∇²(vlapl·sgn) -- vtxc：e2·∫∇²(vlapl·sgn)·ρ·sgn（全强度） +### 3.3 功能名注册 -**结果**：**SCF 发散**。即使 ramp_steps=50 + mixing_beta=0.03，当 ramp_factor → 1 时仍开始振荡。 +`dft_functional scanl` 映射为 LibXC 泛函 `XC_MGGA_X_SCANL` (ID=700) + `XC_MGGA_C_SCANL` (ID=702)。 --- -#### 方案 7：有限差分 Laplacian（FD kernel in G-space） - -**核心思想**：谱 Laplacian ∇²f(G) = -|G|²f(G) 中的 |G|² 是无界的，导致高频放大和 SCF 发散。有限差分（FD）Laplacian 在低 G 时与谱 Laplacian 一致，但在高 G 时自然衰减，最大放大因子有界。 +## 4. FD kernel 数值验证 -**推导**： +### 4.1 验证设置 -在分数坐标 u^α = n_α/N_α（α = 1,2,3）下，Laplacian 表示为： +FCC 原胞 (Si₂)，a = 10.2 Bohr。 -$$\nabla^2 = \sum_{\alpha,\beta} g^{\alpha\beta} \frac{\partial^2}{\partial u^\alpha \partial u^\beta}$$ - -其中 g^{αβ} = GGT[αβ]/lat0² 是逆变度规张量，GGT 是倒格子度规矩阵（ABACUS 中的 `rhopw->GGT`）。 - -对平面波 f(r) = exp(2πi Σ m_γ u^γ)，谱 Laplacian 给出： - -$$\frac{\partial^2 f}{\partial u^\alpha \partial u^\beta} = -(2\pi)^2 m_\alpha m_\beta \, f$$ - -FD 近似给出： - -- **对角项** (α = β)：中心差分 +晶格矢量（lat0 单位）： -$$\left(\frac{\partial^2 f}{\partial (u^\alpha)^2}\right)_{\text{FD}} = -2N_\alpha^2 \left(1 - \cos\frac{2\pi m_\alpha}{N_\alpha}\right) \, f$$ +$$\mathbf{a}_1 = (0, \tfrac{1}{2}, \tfrac{1}{2}), \quad \mathbf{a}_2 = (\tfrac{1}{2}, 0, \tfrac{1}{2}), \quad \mathbf{a}_3 = (\tfrac{1}{2}, \tfrac{1}{2}, 0)$$ -- **交叉项** (α ≠ β)：四点差分 +度规矩阵： -$$\left(\frac{\partial^2 f}{\partial u^\alpha \partial u^\beta}\right)_{\text{FD}} = -N_\alpha N_\beta \sin\frac{2\pi m_\alpha}{N_\alpha} \sin\frac{2\pi m_\beta}{N_\beta} \, f$$ +$$\text{GGT} = \begin{pmatrix} 3 & -1 & -1 \\ -1 & 3 & -1 \\ -1 & -1 & 3 \end{pmatrix}$$ -对低频模式 (m_α ≪ N_α)，FD 退化为谱形式（1-cos(x) ≈ x²/2, sin(x) ≈ x）。对 Nyquist 频率 (m_α = N_α/2)，放大因子为有界值 N_α² · 4/(2π)² ≈ 0.405 N_α²，而非谱 Laplacian 的 (N_α/2)²。 +### 4.2 gg_FD vs gg_spectral（网格 36×36×36, ecutwfc=60 Ry） -组合所有项，FD Laplacian 的等效 G² 核为： +| Miller (m₁,m₂,m₃) | gg_spectral | gg_FD | gg_FD / gg_spectral | +|---------------------|-------------|-------|---------------------| +| (0,0,0) | 0 | 0 | — | +| (1,0,0) | 3.000 | 2.992 | 0.9975 | +| (1,1,0) | 4.000 | 4.005 | 1.0013 | +| (2,0,0) | 12.00 | 11.88 | 0.9900 | +| (3,0,0) | 27.00 | 26.46 | 0.9800 | +| (5,0,0) | 75.00 | 70.36 | 0.9381 | +| (8,0,0) | 192.0 | 163.4 | 0.8511 | +| (10,0,0) | 300.0 | 231.2 | 0.7706 | +| (15,0,0) | 675.0 | 387.4 | 0.5740 | +| (18,0,0) | 972.0 | 393.9 | 0.4053 | -$$\text{gg}_{\text{FD}} = \frac{1}{(2\pi)^2} \sum_{\alpha,\beta} \text{GGT}[\alpha][\beta] \cdot \text{FD}_{\alpha\beta}$$ - -其中： - -$$\text{FD}_{\alpha\alpha} = 2N_\alpha^2(1-\cos\frac{2\pi m_\alpha}{N_\alpha}), \quad \text{FD}_{\alpha\beta}^{(\alpha\neq\beta)} = N_\alpha N_\beta \sin\frac{2\pi m_\alpha}{N_\alpha} \sin\frac{2\pi m_\beta}{N_\beta}$$ +**解读**： -最终，∇²_FD(vlapl·sgn) 在 G 空间中为： +- |m|/N ≤ 0.1 时：FD 与谱偏差 < 1%，物理上等价 +- |m|/N ≈ 0.2 时：偏差 ~6%，开始有差异 +- |m|/N ≈ 0.5（Nyquist）：FD = 0.405 × 谱，大幅衰减 -$$\widetilde{\nabla^2_{\text{FD}}(\text{vlapl}\cdot\text{sgn})}(\mathbf{G}) = -\text{gg}_{\text{FD}}(\mathbf{G}) \cdot \text{tpiba2} \cdot \widetilde{\text{vlapl}\cdot\text{sgn}}(\mathbf{G})$$ +### 4.3 gg_FD vs gg_spectral（网格 45×45×45, ecutwfc=80 Ry） -**实现**： -1. 计算 vlapl_sgn = vlapl·sgn（实空间） -2. FFT: vlapl_sgn → vlapl_g -3. 应用 FD kernel: vlapl_g_lapl[ig] = -vlapl_g[ig] · gg_FD[ig] · tpiba2 -4. IFFT: vlapl_g_lapl → vlapl_lapl（实空间） -5. 加入 V_xc 和 vtxc +| Miller (m₁,m₂,m₃) | gg_spectral | gg_FD | gg_FD / gg_spectral | +|---------------------|-------------|-------|---------------------| +| (1,0,0) | 3.000 | 2.994 | 0.9980 | +| (3,0,0) | 27.00 | 26.66 | 0.9875 | +| (5,0,0) | 75.00 | 72.55 | 0.9673 | +| (8,0,0) | 192.0 | 173.7 | 0.9048 | +| (9,0,0) | 243.0 | 211.7 | 0.8712 | -**交叉项的物理意义**：对于非正交晶胞（如 FCC 原胞），度规张量 GGT 的非对角元不为零。此时，沿一个晶格方向的位移会影响其他方向的梯度分量。FD kernel 的交叉项 sin·sin 正确反映了这一点。对于正交晶胞（GGT 对角），交叉项消失，FD kernel 退化为三个方向的独立贡献。 +ecutwfc=80 对应的最大 |m| ≈ 8.4，此时 FD/谱 ≈ 0.89。即 PW 截断处 FD kernel 衰减约 11%。 -**数值验证**（FCC 原胞, a=10.2 Bohr, nx=ny=nz=36）： +### 4.4 有界性验证 -| Miller 指数 | gg_spectral | gg_FD | gg_FD/gg_spectral | -|-------------|-------------|-------|-------------------| -| (0,0,0) | 0 | 0 | 1.000 | -| (1,0,0) | 3.0 | 2.992 | 0.997 | -| (1,1,0) | 4.0 | 4.005 | 1.001 | -| (5,0,0) | 75.0 | 70.36 | 0.938 | -| (10,0,0) | 300.0 | 231.2 | 0.771 | -| (18,0,0) | 972.0 | 393.9 | 0.405 | +Nyquist 频率 (m_α = N_α/2) 下 FD 因子的最大值： -FD kernel 在低 G 时与谱 kernel 完全一致（<1% 偏差），在高 G 时衰减到谱值的 ~40%，有界且不会导致发散。 +$$\text{FD}_{\alpha\alpha}^{\max} = 2N_\alpha^2 \cdot 2 = 4N_\alpha^2$$ -**SCF 收敛性**：使用 Pulay mixing (beta=0.2)，SCANL 在 ecutwfc=80 Ry 下 SCF 正常收敛。 +$$\text{gg}_{\text{FD}}^{\max} = \frac{1}{(2\pi)^2} \sum_\alpha \text{GGT}[\alpha][\alpha] \cdot 4N_\alpha^2 = \frac{4}{(2\pi)^2} \cdot 3 \cdot 3 N^2 = \frac{36 N^2}{4\pi^2}$$ -**FD 应力验证**（Si₂ FCC 原胞, a₀=10.2 Bohr, δ=0.001, Pulay mixing beta=0.1–0.2）： +对 N=36：gg_FD_max = 1181.8，对应 |G|_FD_max ≈ 21.2 Bohr⁻¹（有界）。 +对 N=45：gg_FD_max = 1846.5。 -PW 基组： +而谱 Laplacian 在同一点：gg_spectral = 9 · (N/2)² → 随 N² 无界增长。 -| 泛函 | ecutwfc (Ry) | σ_FD (kbar) | σ_AB (kbar) | FD/AB | 误差 | -|------|-------------|-------------|-------------|-------|------| -| SCAN | 60 | 434.29 | 434.30 | 1.0000 | 0.00% | -| SCAN | 80 | 434.47 | 434.48 | 1.0000 | 0.00% | -| SCAN | 100 | 434.54 | 434.55 | 1.0000 | 0.00% | -| SCANL | 60 | 431.89 | 435.60 | 0.9915 | -0.85% | -| SCANL | 80 | 432.37 | 434.58 | 0.9949 | -0.51% | - -**注意**：ecutwfc=100 时 SCANL 的 SCF 收敛较困难（能量振荡较大），FD 应力结果不稳定，未列入上表。建议 SCANL 使用 ecutwfc ≤ 80 Ry。 +--- -**与方案 3 的对比**： +## 5. 有限差分应力验证 -| 方案 | V 含 vlapl | vtxc 含 vlapl | SCF | FD/AB | 误差 | -|------|-----------|--------------|-----|-------|------| -| 方案 3 (无 vlapl 势) | ✗ | ✗ | 收敛 | 0.906 | -9.5% | -| 方案 7 (FD Laplacian) | ✓ (FD) | ✓ (FD) | 收敛 | 0.995 | -0.51% | +### 5.1 验证方法 -**FD Laplacian 方案将应力误差从 9.5% 降至 0.51%，改善了约 20 倍。** +对 Si₂ FCC 原胞（Γ 点）在 a₀ ± δ 下计算总能量 E，由差分得到应力： -残余 0.5% 误差的可能来源： -1. FD kernel 在高 G 区与谱 kernel 的偏差（已衰减但仍有差异） -2. 有限差分应力本身的高阶误差（δ=0.001 的中心差分为二阶精度） -3. vlapl stress 项中使用的 Hessian 仍为谱 Hessian +$$\sigma_{\text{FD}} = -\frac{E(a_0+\delta) - E(a_0-\delta)}{V(a_0+\delta) - V(a_0-\delta)}$$ ---- +FCC 原胞体积 $V = |\det(\text{latvec})| \cdot a^3 / 4 = a^3 / 4$。 -| 方案 | V 含 vlapl | vtxc 含 vlapl | stress 含 vlapl | SCF | FD/AB | 误差 | -|------|-----------|--------------|----------------|-----|-------|------| -| SCAN 基线 | N/A | N/A | N/A | 收敛 | 1.000 | 0.00% | -| 方案 1 (严格) | ✓ | ✓ | ✓ | **发散** | — | — | -| 方案 2 (sigma-as-lapl) | ✗ (假) | ✗ (假) | ✗ | 收敛 | 1.046 | 4.6% | -| 方案 3 (严格应力) | ✗ | ✗ | ✓ | 收敛 | 0.906 | -9.5% | -| 方案 4 (vtxc含vlapl) | ✗ | ✓ | ✓ | 收敛 | 0.922 | -7.8% | -| 方案 5 (衰减) | ✓(衰减) | ✓(衰减) | ✓ | 收敛 | ~0.9 | ~10% | -| 方案 6 (ramp) | ✓(渐增) | ✓ | ✓ | **发散** | — | — | -| **方案 7 (FD Laplacian)** | **✓ (FD)** | **✓ (FD)** | **✓** | **收敛** | **0.995** | **-0.5%** | +与 ABACUS 解析应力 σ_AB 比较。单位转换： -**结论**：方案 7（FD Laplacian）是唯一同时满足 SCF 收敛和应力精度的方案。它通过用有界的 FD kernel 替代无界的 |G|² 谱 Laplacian，既解决了 SCF 发散问题，又将 vlapl 势纳入 V_xc 使得密度-势自洽，从而大幅改善应力精度（从 9.5% 误差降至 0.5%）。 +$$\sigma(\text{kbar}) = \sigma(\text{eV/Bohr}^3) \times \frac{147105}{13.6057}$$ ---- +其中 147105 kbar/(Ry/Bohr³)，13.6057 eV/Ry。 -## 4. 模守恒赝势与 SCANL 的根本矛盾 +### 5.2 测试参数 -### 4.1 ∇²(vlapl) 发散的 ecutwfc 依赖 +| 参数 | 值 | +|------|-----| +| 体系 | Si₂ FCC 原胞 | +| a₀ | 10.2 Bohr | +| δ | 0.001 Bohr | +| 赝势 | Si_ONCV_PBE-1.0 (NCPP) | +| smearing | gauss, σ = 0.002 | +| mixing | Pulay, beta = 0.1–0.2, gg0 = 1.0 | +| mixing_tau | 1 | +| scf_thr | 1e-8 | +| pw_seed | 1 | -vlapl 势在倒易空间为 V(G) = -e2·|G|²·vlapl(G)。|G|_max 由 ecutwfc 决定： +### 5.3 SCAN 基线（无 ∇²ρ 依赖） -$$|G|_{max} = \sqrt{2 \cdot E_{cut}}$$ +SCAN 的 vlapl = 0，FD Laplacian 代码路径不激活，用于验证 FD 方法和应力代码。 -|G|_max² 随 ecutwfc 线性增长，意味着 vlapl 势的高频放大随 ecutwfc 增强而加剧： +| ecutwfc (Ry) | FFT grid | σ_FD (kbar) | σ_AB (kbar) | FD/AB | 误差 | +|-------------|----------|-------------|-------------|-------|------| +| 60 | 36³ | 434.29 | 434.30 | 1.0000 | 0.00% | +| 80 | 45³ | 434.47 | 434.48 | 1.0000 | 0.00% | +| 100 | 45³ | 434.54 | 434.55 | 1.0000 | 0.00% | -- **低 ecutwfc**：|G|_max 小，放大范围有限，SCF 可能收敛，但基组不完备导致其他误差 -- **高 ecutwfc**：|G|_max 大，高频被剧烈放大，SCF 几乎必然发散 +**结论**：SCAN 在所有 ecutwfc 下 FD/AB = 1.0000，验证了代码和方法的正确性。 -**严格实现方案 1 在任何实际使用的 ecutwfc 下都无法收敛**。 +### 5.4 SCANL（FD Laplacian 方案） -### 4.2 模守恒赝势需要高 ecutwfc +| ecutwfc (Ry) | FFT grid | σ_FD (kbar) | σ_AB (kbar) | FD/AB | 误差 | +|-------------|----------|-------------|-------------|-------|------| +| 60 | 36³ | 431.89 | 435.60 | 0.9915 | -0.85% | +| 80 | 45³ | 432.37 | 434.58 | 0.9949 | -0.51% | -模守恒赝势（NCPP）比 PAW 方法需要更高的 ecutwfc 才能收敛： +ecutwfc=100 时 SCANL 的 SCF 收敛困难（能量振荡），FD 结果不可靠，未列入。 -| 赝势类型 | 典型 ecutwfc (Ry) | 典型 ecutwfc (eV) | -|---------|------------------|------------------| -| PAW/USPP | 30–50 | 400–700 | -| NCPP | 60–100+ | 800–1400+ | +### 5.5 完整能量数据 -ABACUS 目前使用模守恒赝势。NCPP 要求的高 ecutwfc 与 vlapl 势的 |G|² 发散直接矛盾： +ecutwfc = 80 Ry： -- **低 ecutwfc**（< 50 Ry）：SCF 可能收敛，但 NCPP 的赝势不完备，结果不收敛 -- **高 ecutwfc**（> 60 Ry）：NCPP 收敛，但 vlapl 势发散 +| | E(a₀+δ) (eV) | E(a₀−δ) (eV) | ΔE (eV) | +|------|-------------|-------------|---------| +| SCAN | −197.22465 | −197.21838 | −0.006271 | +| SCANL | −198.70596 | −198.69972 | −0.006241 | -**这是一个不可调和的矛盾**。在 NCPP + 平面波框架下，不存在一个 ecutwfc 值能同时满足赝势收敛和 vlapl 势稳定。 +ecutwfc = 60 Ry： -### 4.3 VASP 的做法 +| | E(a₀+δ) (eV) | E(a₀−δ) (eV) | ΔE (eV) | +|------|-------------|-------------|---------| +| SCAN | −197.22344 | −197.21717 | −0.006269 | +| SCANL | −198.68669 | −198.68046 | −0.006234 | -VASP wiki [METAGGA 页面] 明确指出： +SCAN 的 ΔE 跨 ecutwfc 稳定（−0.00627 eV），SCANL 的 ΔE 在 ecut=60/80 时接近（−0.00623/−0.00624 eV），与 SCAN 的差异约 0.5%，与 FD 应力误差一致。 -> The results obtained with the meta-GGA functionals that depend on the Laplacian of the density ∇²n (e.g., SCAN-L) may not be reliable for large values of the energy cutoff ENCUT due to numerical instability. According to some tests, it is not recommended to use values of ENCUT above 800 eV. +### 5.6 与其他方案对比 -VASP 使用 PAW 方法，ecutwfc 要求较低（~400–700 eV），可以在 800 eV 以下的窗口中获得结果。但即便如此，VASP 也承认高截断能下结果不可靠。对于 NCPP，800 eV（约 60 Ry）可能刚达到赝势收敛的门槛，而更高截断又会触发不稳定——几乎不存在可用的 ecutwfc 窗口。 +| 方案 | V 含 vlapl | vtxc 含 vlapl | stress 含 vlapl | SCF | FD/AB (ecut=80) | 误差 | +|------|-----------|--------------|----------------|-----|-------|------| +| 严格谱 Laplacian | ✓ (谱) | ✓ (谱) | ✓ | **发散** | — | — | +| sigma-as-lapl | ✗ (假输入) | ✗ (假输入) | ✗ | 收敛 | 1.046 | +4.6% | +| V/vtxc 均省略 vlapl | ✗ | ✗ | ✓ | 收敛 | 0.906 | −9.5% | +| vtxc 含 vlapl | ✗ | ✓ | ✓ | 收敛 | 0.922 | −7.8% | +| 高斯衰减 | ✓ (衰减) | ✓ (衰减) | ✓ | 收敛 | ~0.9 | ~10% | +| 绝热 ramp | ✓ (渐增) | ✓ | ✓ | **发散** | — | — | +| **FD Laplacian** | **✓ (FD)** | **✓ (FD)** | **✓** | **收敛** | **0.995** | **−0.5%** | --- -## 5. 结论与建议 +## 6. 误差分析 + +### 6.1 残余 0.5% 应力误差来源 -### 5.1 推荐方案：FD Laplacian（方案 7） +1. **FD kernel 与谱 kernel 在高 G 的偏差**：在 PW 截断处（|m|/N ≈ 0.37 for ecut=80），FD kernel 衰减到谱值的 ~89%。vlapl 势在高 G 分量上偏弱，密度与严格 SCANL 密度有微小差异。 -**FD Laplacian 方案是当前最佳方案**，它： +2. **FD 应力本身的高阶误差**：δ = 0.001 Bohr 的中心差分为二阶精度，误差 O(δ²/a⁴) ≈ 10⁻⁸，对 0.5% 误差贡献可忽略。 -1. **SCF 收敛**：FD kernel 有界，不会导致 |G|² 发散 -2. **密度-势自洽**：V 包含 vlapl 贡献（通过 FD Laplacian），密度在近似自洽势下产生 -3. **应力精度高**：FD/AB = 0.995，误差仅 0.5%（远优于方案 3 的 9.5%） -4. **适用于非正交晶胞**：通过度规张量 GGT 正确处理交叉项 -5. **无需调参**：不依赖衰减因子或 ramp 步数等人为参数 +3. **vlapl stress 项中 Hessian 仍为谱 Hessian**：gradcorr 中的 vlapl 应力项使用谱方法计算 Hessian，而势用 FD Laplacian。此不一致性贡献未知但量级较小。 -**残余 0.5% 误差**是 FD kernel 与谱 kernel 在高 G 区差异的代价，对于实际应用完全可以接受。 +### 6.2 ecutwfc=100 的 SCF 困难 -### 5.2 使用注意事项 +ecutwfc=100 时 SCANL 的 ΔE = −0.00568 eV，显著偏离 SCAN 的 −0.00628 eV（偏差 ~10%），而 ecut=60/80 时偏差仅 ~0.5%。原因可能是： -1. **FD Laplacian 是对谱 Laplacian 的近似**：在 FFT 网格上，它是"正确"的 Laplacian（对网格函数的有限差分），但与连续 Laplacian 有差异 -2. **需要 mixing_tau = 1**：meta-GGA 泛函需要混合动能密度 -3. **建议使用 Pulay mixing**：mixing_beta=0.2，对 SCANL 效果较好 -4. **仍建议优先考虑 r²SCAN**：r²SCAN 无 ∇²ρ 依赖，数值更稳定 +- 45³ FFT 网格对 ecut=100 的 35749 个平面波不够大，某些高 G 分量的 FD kernel 在网格边界处行为异常 +- SCF 收敛到与 ecut=60/80 不同的密度 -### 5.3 代码实现 +**建议**：SCANL 使用 ecutwfc ≤ 80 Ry。若需更高截断，应确保 FFT 网格足够大。 -SCANL 已注册为 `dft_functional scanl`，内部映射为 `XC_MGGA_X_SCANL + XC_MGGA_C_SCANL`。 +### 6.3 与严格实现的关系 -vlapl 势项实现于 `xc_functional_libxc_vxc.cpp` 中的 `v_xc_meta` 函数，使用 G 空间 FD kernel 计算 ∇²(vlapl·sgn)。 +严格 SCANL 实现要求谱 Laplacian ∇² = −|G|²，但在任何实际 ecutwfc 下导致 SCF 发散。FD Laplacian 是在 FFT 网格上对离散函数求 Laplacian 的自洽方法——它是网格函数的"正确" Laplacian，但与连续 Laplacian 有差异。此差异随网格加密而减小（FD → spectral 当 N → ∞），但在实际网格上不可消除。 -vlapl 应力项实现于 `xc_functional_gradcorr.cpp` 中的 `gradcorr` 函数。 +--- -### 5.4 替代泛函 +## 7. 使用建议 -| 泛函 | 依赖变量 | 特点 | -|------|---------|------| -| **r²SCAN** | ρ, ∇ρ, τ | SCAN 的正则化版本，数值稳定，推荐首选 | -| **SCAN** | ρ, ∇ρ, τ | 无 ∇²ρ 依赖，数值稳定 | -| **SCANL** | ρ, ∇ρ, ∇²ρ | 需 FD Laplacian 近似，0.5% 应力误差 | +1. **dft_functional scanl**：已注册，映射到 XC_MGGA_X_SCANL + XC_MGGA_C_SCANL +2. **ecutwfc ≤ 80 Ry**：ecutwfc 过高时 SCF 收敛困难 +3. **mixing_tau = 1**：meta-GGA 必需 +4. **Pulay mixing, beta = 0.1–0.2**：对 SCANL 收敛效果较好 +5. **优先考虑 r²SCAN 或 SCAN**：无 ∇²ρ 依赖，数值更稳定