Skip to content

Commit 812061e

Browse files
linpeizedyzheng
linpeize
authored and
dyzheng
committed
Fix(module_xc): compute real density laplacian and vlapl stress for meta-GGA functionals
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
1 parent 4ef86d8 commit 812061e

11 files changed

Lines changed: 394 additions & 42 deletions

‎source/module_esolver/esolver_ks.cpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -185,6 +185,7 @@ void ESolver_KS<T, Device>::before_all_runners(UnitCell& ucell, const Input_para
185185
{
186186
XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func);
187187
}
188+
GlobalV::ofs_running<<XC_Functional::output_info()<<std::endl;
188189
ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SETUP UNITCELL");
189190

190191
//! 4) setup the charge mixing parameters

‎source/module_hamilt_general/module_xc/test/CMakeLists.txt‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,7 @@ AddTest(
1111
TARGET XCTest_HSE
1212
LIBS parameter MPI::MPI_CXX Libxc::xc # required by global.h; for details, `remove_definitions(-D__MPI)`.
1313
SOURCES test_xc1.cpp ../xc_functional.cpp ../xc_functional_libxc.cpp
14+
../xc_functional_libxc_tools.cpp
1415
)
1516

1617

@@ -36,6 +37,7 @@ AddTest(
3637
../xc_functional_libxc_wrapper_xc.cpp
3738
../xc_functional_libxc_wrapper_gcxc.cpp
3839
../xc_functional_libxc_wrapper_tauxc.cpp
40+
../xc_functional_libxc_tools.cpp
3941
../xc_funct_corr_gga.cpp ../xc_funct_corr_lda.cpp ../xc_funct_exch_gga.cpp
4042
../xc_funct_exch_lda.cpp ../xc_funct_hcth.cpp
4143
../../../module_base/matrix.cpp
@@ -57,6 +59,7 @@ AddTest(
5759
../xc_functional_libxc_wrapper_xc.cpp
5860
../xc_functional_libxc_wrapper_gcxc.cpp
5961
../xc_functional_libxc_wrapper_tauxc.cpp
62+
../xc_functional_libxc_tools.cpp
6063
../xc_funct_corr_gga.cpp ../xc_funct_corr_lda.cpp
6164
../xc_funct_exch_gga.cpp ../xc_funct_exch_lda.cpp ../xc_funct_hcth.cpp
6265
)

‎source/module_hamilt_general/module_xc/test/test_xc4.cpp‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -43,8 +43,9 @@ class XCTest_SCAN : public XCTest
4343

4444
for(int i=0;i<5;i++)
4545
{
46-
double e,v,v1,v2,v3;
47-
XC_Functional_Libxc::tau_xc(XC_Functional::get_func_id(), rho[i],grho[i],tau[i],e,v1,v2,v3);
46+
double e,v,v1,v2,v3,v4;
47+
// SCAN does not depend on laplacian, pass 0.0
48+
XC_Functional_Libxc::tau_xc(XC_Functional::get_func_id(), rho[i],grho[i],0.0,tau[i],e,v1,v2,v3,v4);
4849
e_.push_back(e);
4950
v1_.push_back(v1);
5051
v2_.push_back(v2);

‎source/module_hamilt_general/module_xc/xc_functional.cpp‎

Lines changed: 60 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -34,24 +34,37 @@ int XC_Functional::get_func_type()
3434
{
3535
return func_type;
3636
}
37+
3738
void XC_Functional::set_xc_first_loop(const UnitCell& ucell)
3839
{
40+
ModuleBase::TITLE("XC_Functional", "set_xc_first_loop");
3941
/** In the special "two-level" calculation case,
4042
the first scf iteration only calculate the functional without exact
4143
exchange. but in "nscf" calculation, there is no need of "two-level"
4244
method. */
43-
if (ucell.atoms[0].ncpp.xc_func == "HF" || ucell.atoms[0].ncpp.xc_func == "HSE"
44-
|| ucell.atoms[0].ncpp.xc_func == "PBE0"|| ucell.atoms[0].ncpp.xc_func == "LC_PBE"
45-
|| ucell.atoms[0].ncpp.xc_func == "LC_WPBE" || ucell.atoms[0].ncpp.xc_func == "LRC_WPBEH"
46-
|| ucell.atoms[0].ncpp.xc_func == "CAM_PBEH") {
45+
if (ucell.atoms[0].ncpp.xc_func == "HF" || ucell.atoms[0].ncpp.xc_func == "PBE0" || ucell.atoms[0].ncpp.xc_func == "HSE")
46+
{
4747
XC_Functional::set_xc_type("pbe");
4848
}
49-
else if (ucell.atoms[0].ncpp.xc_func == "SCAN0") {
50-
XC_Functional::set_xc_type("scan");
49+
else if ( ucell.atoms[0].ncpp.xc_func == "LC_PBE" || ucell.atoms[0].ncpp.xc_func == "LC_WPBE"
50+
|| ucell.atoms[0].ncpp.xc_func == "LRC_WPBEH" || ucell.atoms[0].ncpp.xc_func == "CAM_PBEH" )
51+
{
52+
XC_Functional::set_xc_type("pbe");
5153
}
52-
else if (ucell.atoms[0].ncpp.xc_func == "B3LYP") {
54+
// added by jghan, 2024-07-07
55+
else if ( ucell.atoms[0].ncpp.xc_func == "MULLER" || ucell.atoms[0].ncpp.xc_func == "POWER"
56+
|| ucell.atoms[0].ncpp.xc_func == "WP22" || ucell.atoms[0].ncpp.xc_func == "CWP22" )
57+
{
58+
XC_Functional::set_xc_type("pbe");
59+
}
60+
else if (ucell.atoms[0].ncpp.xc_func == "B3LYP")
61+
{
5362
XC_Functional::set_xc_type("blyp");
5463
}
64+
else if (ucell.atoms[0].ncpp.xc_func == "SCAN0")
65+
{
66+
XC_Functional::set_xc_type("scan");
67+
}
5568
}
5669

5770
// The setting values of functional id according to the index in LIBXC
@@ -346,4 +359,44 @@ void XC_Functional::set_xc_type(const std::string xc_func_in)
346359
use_libxc = false;
347360
#endif
348361

362+
}
363+
364+
365+
std::string XC_Functional::output_info()
366+
{
367+
#ifdef USE_LIBXC
368+
if(use_libxc)
369+
{
370+
std::stringstream ss;
371+
ss<<" Libxc v"<<xc_version_string()<<std::endl;
372+
ss<<"\t"<<xc_reference()<<std::endl;
373+
374+
std::vector<xc_func_type> funcs = XC_Functional_Libxc::init_func(func_id, XC_UNPOLARIZED);
375+
for(const auto &func : funcs)
376+
{
377+
const xc_func_info_type *info = xc_func_get_info(&func);
378+
ss<<" XC: "<<xc_func_info_get_name(info)<<std::endl;
379+
for(int i=0; i<XC_MAX_REFERENCES; ++i)
380+
{
381+
const func_reference_type *ref = xc_func_info_get_references(func.info, i);
382+
if(ref)
383+
ss<<"\t"<<xc_func_reference_get_ref(ref)<<std::endl;
384+
}
385+
}
386+
XC_Functional_Libxc::finish_func(funcs);
387+
return ss.str();
388+
}
389+
else
390+
{
391+
std::string s = " XC:\t";
392+
for(const auto &id: func_id)
393+
s += std::string(xc_functional_get_name(id))+"\t";
394+
return s;
395+
}
396+
#else
397+
std::string s = " XC:\t";
398+
for(const auto &id: func_id)
399+
s += std::to_string(id)+"\t";
400+
return s;
401+
#endif
349402
}

‎source/module_hamilt_general/module_xc/xc_functional.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -88,6 +88,7 @@ class XC_Functional
8888

8989
public:
9090
static std::vector<int> get_func_id() { return func_id; }
91+
static std::string output_info();
9192

9293
//-------------------
9394
// xc_functional_wrapper_xc.cpp

‎source/module_hamilt_general/module_xc/xc_functional_gradcorr.cpp‎

Lines changed: 151 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -67,6 +67,10 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v,
6767
ModuleBase::Vector3<double>* gdr2 = nullptr;
6868
ModuleBase::Vector3<double>* h1 = nullptr;
6969
ModuleBase::Vector3<double>* h2 = nullptr;
70+
double* lapl1 = nullptr;
71+
double* lapl2 = nullptr;
72+
double* hess1 = nullptr; // nspin=1 Hessian (6 components)
73+
double* hess2 = nullptr; // nspin=1 Hessian for spin-down channel
7074
double* neg = nullptr;
7175
double** vsave = nullptr;
7276
double** vgg = nullptr;
@@ -95,6 +99,46 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v,
9599

96100
XC_Functional::grad_rho( rhogsum1 , gdr1, rhopw, ucell->tpiba);
97101

102+
XC_Functional::grad_rho( rhogsum1 , gdr1, rhopw, ucell->tpiba);
103+
104+
// compute laplacian of total density (rho[0] + rho_core) for nspin=1,2
105+
// and as the "up" channel equivalent for nspin=4 noncollinear case
106+
{
107+
std::vector<double> rho_total(rhopw->nrxx);
108+
#ifdef _OPENMP
109+
#pragma omp parallel for schedule(static, 1024)
110+
#endif
111+
for(int ir=0; ir<rhopw->nrxx; ++ir)
112+
rho_total[ir] = rhotmp1[ir];
113+
std::vector<double> lapl_all = XC_Functional_Libxc::cal_lapl(
114+
1, rhopw->nrxx, rho_total, ucell->tpiba2, chr);
115+
lapl1 = new double[rhopw->nrxx];
116+
#ifdef _OPENMP
117+
#pragma omp parallel for schedule(static, 1024)
118+
#endif
119+
for(int ir=0; ir<rhopw->nrxx; ++ir)
120+
lapl1[ir] = lapl_all[ir];
121+
}
122+
123+
// compute Hessian of total density for meta-GGA stress
124+
if(is_stress && (func_type == 3 || func_type == 5))
125+
{
126+
std::vector<double> rho_total(rhopw->nrxx);
127+
#ifdef _OPENMP
128+
#pragma omp parallel for schedule(static, 1024)
129+
#endif
130+
for(int ir=0; ir<rhopw->nrxx; ++ir)
131+
rho_total[ir] = rhotmp1[ir];
132+
std::vector<double> hess_all = XC_Functional_Libxc::cal_rho_hessian(
133+
1, rhopw->nrxx, rho_total, chr);
134+
hess1 = new double[rhopw->nrxx * 6];
135+
#ifdef _OPENMP
136+
#pragma omp parallel for schedule(static, 1024)
137+
#endif
138+
for(int i=0; i<rhopw->nrxx * 6; ++i)
139+
hess1[i] = hess_all[i];
140+
}
141+
98142
// for spin polarized case;
99143
// calculate the gradient of (rho_core+rho) in reciprocal space.
100144
if(PARAM.inp.nspin==2)
@@ -120,6 +164,43 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v,
120164
if(!is_stress) { h2 = new ModuleBase::Vector3<double>[rhopw->nrxx]; }
121165

122166
XC_Functional::grad_rho( rhogsum2 , gdr2, rhopw, ucell->tpiba);
167+
168+
// compute laplacian of spin-down density (rho[1] + rho_core)
169+
{
170+
std::vector<double> rho_dw(rhopw->nrxx);
171+
#ifdef _OPENMP
172+
#pragma omp parallel for schedule(static, 1024)
173+
#endif
174+
for(int ir=0; ir<rhopw->nrxx; ++ir)
175+
rho_dw[ir] = rhotmp2[ir];
176+
std::vector<double> lapl_all = XC_Functional_Libxc::cal_lapl(
177+
1, rhopw->nrxx, rho_dw, ucell->tpiba2, chr);
178+
lapl2 = new double[rhopw->nrxx];
179+
#ifdef _OPENMP
180+
#pragma omp parallel for schedule(static, 1024)
181+
#endif
182+
for(int ir=0; ir<rhopw->nrxx; ++ir)
183+
lapl2[ir] = lapl_all[ir];
184+
}
185+
186+
// compute Hessian of spin-down density for meta-GGA stress
187+
if(is_stress && (func_type == 3 || func_type == 5))
188+
{
189+
std::vector<double> rho_dw(rhopw->nrxx);
190+
#ifdef _OPENMP
191+
#pragma omp parallel for schedule(static, 1024)
192+
#endif
193+
for(int ir=0; ir<rhopw->nrxx; ++ir)
194+
rho_dw[ir] = rhotmp2[ir];
195+
std::vector<double> hess_all = XC_Functional_Libxc::cal_rho_hessian(
196+
1, rhopw->nrxx, rho_dw, chr);
197+
hess2 = new double[rhopw->nrxx * 6];
198+
#ifdef _OPENMP
199+
#pragma omp parallel for schedule(static, 1024)
200+
#endif
201+
for(int i=0; i<rhopw->nrxx * 6; ++i)
202+
hess2[i] = hess_all[i];
203+
}
123204
}
124205

125206
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,
189270
XC_Functional::grad_rho( rhogsum1 , gdr1, rhopw, ucell->tpiba);
190271
XC_Functional::grad_rho( rhogsum2 , gdr2, rhopw, ucell->tpiba);
191272

273+
// re-compute laplacians after noncollinear rho update
274+
{
275+
std::vector<double> rho_up(rhopw->nrxx);
276+
#ifdef _OPENMP
277+
#pragma omp parallel for schedule(static, 1024)
278+
#endif
279+
for(int ir=0; ir<rhopw->nrxx; ++ir)
280+
rho_up[ir] = rhotmp1[ir];
281+
std::vector<double> lapl_up = XC_Functional_Libxc::cal_lapl(
282+
1, rhopw->nrxx, rho_up, ucell->tpiba2, chr);
283+
lapl1 = new double[rhopw->nrxx];
284+
#ifdef _OPENMP
285+
#pragma omp parallel for schedule(static, 1024)
286+
#endif
287+
for(int ir=0; ir<rhopw->nrxx; ++ir)
288+
lapl1[ir] = lapl_up[ir];
289+
}
290+
{
291+
std::vector<double> rho_dw(rhopw->nrxx);
292+
#ifdef _OPENMP
293+
#pragma omp parallel for schedule(static, 1024)
294+
#endif
295+
for(int ir=0; ir<rhopw->nrxx; ++ir)
296+
rho_dw[ir] = rhotmp2[ir];
297+
std::vector<double> lapl_dw = XC_Functional_Libxc::cal_lapl(
298+
1, rhopw->nrxx, rho_dw, ucell->tpiba2, chr);
299+
lapl2 = new double[rhopw->nrxx];
300+
#ifdef _OPENMP
301+
#pragma omp parallel for schedule(static, 1024)
302+
#endif
303+
for(int ir=0; ir<rhopw->nrxx; ++ir)
304+
lapl2[ir] = lapl_dw[ir];
305+
}
306+
192307
}
193308

194309
const double epsr = 1.0e-6;
@@ -249,9 +364,20 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v,
249364
#ifdef USE_LIBXC
250365
if(func_type == 3 || func_type == 5) //the gradcorr part to stress of mGGA
251366
{
252-
double v3xc;
367+
double v3xc, vlaplc;
253368
double atau = chr->kin_r[0][ir]/2.0;
254-
XC_Functional_Libxc::tau_xc( func_id, arho, grho2a, atau, sxc, v1xc, v2xc, v3xc);
369+
XC_Functional_Libxc::tau_xc( func_id, arho, grho2a, lapl1[ir], atau, sxc, v1xc, v2xc, v3xc, vlaplc);
370+
// add vlapl stress contribution: σ_αβ^vlapl = -2 ∫ vlapl ∂²ρ/∂r_α∂r_β dr
371+
// Hessian components order: xx, yy, zz, xy, yz, zx
372+
for(int l = 0; l < 3; l++)
373+
{
374+
for(int m = 0; m < l+1; m++)
375+
{
376+
int ind = l*3 + m;
377+
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;
378+
local_stress_gga[ind] -= 2.0 * vlaplc * hess1[ic * rhopw->nrxx + ir] * ModuleBase::e2;
379+
}
380+
}
255381
}
256382
else
257383
{
@@ -308,13 +434,27 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v,
308434
double sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud;
309435
if(func_type == 3 || func_type == 5) //the gradcorr part to stress of mGGA
310436
{
311-
double v3xcup, v3xcdw;
437+
double v3xcup, v3xcdw, vlaplup, vlapldw;
312438
double atau1 = chr->kin_r[0][ir]/2.0;
313439
double atau2 = chr->kin_r[1][ir]/2.0;
314440
XC_Functional_Libxc::tau_xc_spin(
315441
func_id,
316-
rhotmp1[ir], rhotmp2[ir], gdr1[ir], gdr2[ir],
317-
atau1, atau2, sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud, v3xcup, v3xcdw);
442+
rhotmp1[ir], rhotmp2[ir], gdr1[ir], gdr2[ir],
443+
lapl1[ir], lapl2[ir],
444+
atau1, atau2, sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud, v3xcup, v3xcdw, vlaplup, vlapldw);
445+
// add vlapl stress contribution for both spin channels
446+
if(is_stress)
447+
{
448+
for(int l = 0; l < 3; l++)
449+
{
450+
for(int m = 0; m < l+1; m++)
451+
{
452+
int ind = l*3 + m;
453+
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;
454+
local_stress_gga[ind] -= 2.0 * (vlaplup * hess1[ic * rhopw->nrxx + ir] + vlapldw * hess2[ic * rhopw->nrxx + ir]) * ModuleBase::e2;
455+
}
456+
}
457+
}
318458
}
319459
else
320460
{
@@ -582,6 +722,8 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v,
582722
delete[] rhotmp1;
583723
delete[] rhogsum1;
584724
delete[] gdr1;
725+
delete[] lapl1;
726+
delete[] hess1;
585727
if(!is_stress) { delete[] h1;
586728
}
587729

@@ -590,6 +732,8 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v,
590732
delete[] rhotmp2;
591733
delete[] rhogsum2;
592734
delete[] gdr2;
735+
delete[] lapl2;
736+
delete[] hess2;
593737
if(!is_stress) { delete[] h2;
594738
}
595739
}
@@ -609,6 +753,8 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v,
609753
delete[] rhotmp2;
610754
delete[] rhogsum2;
611755
delete[] gdr2;
756+
delete[] lapl2;
757+
delete[] hess2;
612758
}
613759

614760
return;

0 commit comments

Comments
 (0)