Skip to content

Commit c6f0d57

Browse files
dyzhengdyzheng
committed
fix(xc): improve SCAN-L Laplacian implementation with bug fixes and tests
Bug fixes: - libxc_mgga_wrap: add hybrid alpha scaling for vlapl in tau_xc/tau_xc_spin - xc_grad: replace hardcoded func_type check with XC_FLAGS_NEEDS_LAPLACIAN flag to support future Laplacian-dependent meta-GGA functionals generically Efficiency: - libxc_tools: add cal_gdr_and_lapl() to share real2recip FFT between gradient and Laplacian computation (one FFT saved per spin channel) - libxc_pot: use cal_gdr_and_lapl in v_xc_meta instead of separate cal_gdr + cal_lapl calls New features: - libxc_setup: add generic XC_FLAGS_NEEDS_LAPLACIAN warning for all Laplacian-dependent functionals, not just SCANL Tests: - test_xc6.cpp: SCANL Laplacian sensitivity test with libxc reference values (MGGA_X_SCANL, MGGA_C_SCANL, tau_xc wrapper) - test_xc7.cpp: analytical Laplacian tests for laplacian_rho() (zero input, imaginary rhog, linearity, single plane wave) - Updated test CMakeLists.txt with new test targets
1 parent 89dcd28 commit c6f0d57

9 files changed

Lines changed: 386 additions & 16 deletions

File tree

‎source/source_hamilt/module_xc/libxc_abacus.h‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -105,6 +105,15 @@ namespace XC_Functional_Libxc
105105
const double tpiba,
106106
const Charge* const chr);
107107

108+
extern void cal_gdr_and_lapl(
109+
const int nspin,
110+
const std::size_t nrxx,
111+
const std::vector<double> &rho,
112+
const double tpiba,
113+
const Charge* const chr,
114+
std::vector<std::vector<ModuleBase::Vector3<double>>> &gdr,
115+
std::vector<double> &lapl);
116+
108117
// converting grho (abacus=>libxc)
109118
extern std::vector<double> convert_sigma(
110119
const std::vector<std::vector<ModuleBase::Vector3<double>>> &gdr);

‎source/source_hamilt/module_xc/libxc_mgga_wrap.cpp‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -54,6 +54,7 @@ void XC_Functional_Libxc::tau_xc(
5454
v1 *= (1.0 - hybrid_alpha);
5555
v2 *= (1.0 - hybrid_alpha);
5656
v3 *= (1.0 - hybrid_alpha);
57+
vlapl_rho *= (1.0 - hybrid_alpha);
5758
}
5859
#endif
5960
sxc += s * rho;
@@ -151,6 +152,8 @@ void XC_Functional_Libxc::tau_xc_spin(
151152
v2xc[2] *= (1.0 - hybrid_alpha);
152153
v3xc[0] *= (1.0 - hybrid_alpha);
153154
v3xc[1] *= (1.0 - hybrid_alpha);
155+
vlapl[0] *= (1.0 - hybrid_alpha);
156+
vlapl[1] *= (1.0 - hybrid_alpha);
154157
}
155158
#endif
156159

‎source/source_hamilt/module_xc/libxc_pot.cpp‎

Lines changed: 3 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -384,13 +384,11 @@ std::tuple<double,double,ModuleBase::matrix,ModuleBase::matrix> XC_Functional_Li
384384
/* hse_omega = */ hse_omega);
385385

386386
const std::vector<double> rho = XC_Functional_Libxc::convert_rho(nspin, nrxx, chr);
387-
const std::vector<std::vector<ModuleBase::Vector3<double>>> gdr
388-
= XC_Functional_Libxc::cal_gdr(nspin, nrxx, rho, tpiba, chr);
387+
std::vector<std::vector<ModuleBase::Vector3<double>>> gdr;
388+
std::vector<double> lapl;
389+
XC_Functional_Libxc::cal_gdr_and_lapl(nspin, nrxx, rho, tpiba, chr, gdr, lapl);
389390
const std::vector<double> sigma = XC_Functional_Libxc::convert_sigma(gdr);
390391

391-
// compute laplacian for mGGA functionals
392-
const std::vector<double> lapl = XC_Functional_Libxc::cal_lapl(nspin, nrxx, rho, tpiba, chr);
393-
394392
//converting kin_r
395393
std::vector<double> kin_r;
396394
kin_r.resize(nrxx*nspin);

‎source/source_hamilt/module_xc/libxc_setup.cpp‎

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -180,6 +180,22 @@ XC_Functional_Libxc::set_xc_type_libxc(const std::string& xc_func_in)
180180
ModuleBase::WARNING_QUIT("XC_Functional::set_xc_type_libxc", message);
181181
}
182182

183+
{
184+
std::vector<xc_func_type> tmp_funcs = XC_Functional_Libxc::init_func(func_id, XC_UNPOLARIZED, 0.0, 0.0);
185+
for (auto& f : tmp_funcs)
186+
{
187+
if (f.info->flags & XC_FLAGS_NEEDS_LAPLACIAN)
188+
{
189+
std::cout << " WARNING: XC functional \"" << f.info->name
190+
<< "\" requires Laplacian of density (nabla^2 rho)."
191+
<< " This may require a higher energy cutoff for numerical stability."
192+
<< std::endl;
193+
break;
194+
}
195+
}
196+
XC_Functional_Libxc::finish_func(tmp_funcs);
197+
}
198+
183199
// return
184200
return std::make_pair(func_type, func_id);
185201
}

‎source/source_hamilt/module_xc/libxc_tools.cpp‎

Lines changed: 39 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -88,6 +88,45 @@ XC_Functional_Libxc::cal_gdr(
8888
return gdr;
8989
}
9090

91+
void XC_Functional_Libxc::cal_gdr_and_lapl(
92+
const int nspin,
93+
const std::size_t nrxx,
94+
const std::vector<double> &rho,
95+
const double tpiba,
96+
const Charge* const chr,
97+
std::vector<std::vector<ModuleBase::Vector3<double>>> &gdr,
98+
std::vector<double> &lapl)
99+
{
100+
gdr.resize(nspin);
101+
lapl.assign(nrxx * nspin, 0.0);
102+
for( int is=0; is!=nspin; ++is )
103+
{
104+
std::vector<double> rhor(nrxx);
105+
#ifdef _OPENMP
106+
#pragma omp parallel for schedule(static, 1024)
107+
#endif
108+
for(std::size_t ir=0; ir<nrxx; ++ir)
109+
{
110+
rhor[ir] = rho[ir*nspin+is];
111+
}
112+
std::vector<std::complex<double>> rhog(chr->rhopw->npw);
113+
chr->rhopw->real2recip(rhor.data(), rhog.data());
114+
115+
gdr[is].resize(nrxx);
116+
XC_Functional::grad_rho(rhog.data(), gdr[is].data(), chr->rhopw, tpiba);
117+
118+
std::vector<double> lapl_spin(nrxx);
119+
XC_Functional::laplacian_rho(rhog.data(), lapl_spin.data(), chr->rhopw, tpiba);
120+
#ifdef _OPENMP
121+
#pragma omp parallel for schedule(static, 1024)
122+
#endif
123+
for(std::size_t ir=0; ir<nrxx; ++ir)
124+
{
125+
lapl[ir*nspin+is] = lapl_spin[ir];
126+
}
127+
}
128+
}
129+
91130
// converting grho (abacus=>libxc)
92131
std::vector<double> XC_Functional_Libxc::convert_sigma(
93132
const std::vector<std::vector<ModuleBase::Vector3<double>>> &gdr)

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

Lines changed: 41 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -98,3 +98,44 @@ AddTest(
9898
../../../source_base/module_fft/fft_cpu.cpp
9999
${FFT_SRC}
100100
)
101+
102+
AddTest(
103+
TARGET MODULE_HAMILT_XCTest_SCANL_LAPL
104+
LIBS parameter MPI::MPI_CXX Libxc::xc
105+
SOURCES test_xc6.cpp ../xc_functional.cpp ../xc_lda_wrap.cpp
106+
../xc_gga_wrap.cpp
107+
../libxc_setup.cpp
108+
../libxc_lda_wrap.cpp
109+
../libxc_gga_wrap.cpp
110+
../libxc_mgga_wrap.cpp
111+
../xc_gga_corr.cpp ../xc_lda_corr.cpp
112+
../xc_gga_exch.cpp ../xc_lda_exch.cpp ../xc_hcth.cpp
113+
)
114+
115+
if (USE_CUDA)
116+
list(APPEND FFT_SRC ../../../source_base/module_fft/fft_cuda.cpp)
117+
endif()
118+
if (USE_ROCM)
119+
list(APPEND FFT_SRC ../../../source_base/module_fft/fft_rocm.cpp)
120+
endif()
121+
AddTest(
122+
TARGET MODULE_HAMILT_XCTest_LAPL
123+
LIBS parameter MPI::MPI_CXX Libxc::xc psi device container
124+
SOURCES test_xc7.cpp ../xc_grad.cpp ../xc_functional.cpp
125+
../xc_lda_wrap.cpp ../xc_gga_wrap.cpp
126+
../libxc_setup.cpp
127+
../libxc_lda_wrap.cpp
128+
../libxc_gga_wrap.cpp
129+
../libxc_mgga_wrap.cpp
130+
../libxc_tools.cpp
131+
../xc_gga_corr.cpp ../xc_lda_corr.cpp ../xc_gga_exch.cpp
132+
../xc_lda_exch.cpp ../xc_hcth.cpp
133+
../../../source_base/matrix.cpp
134+
../../../source_base/memory_recorder.cpp
135+
../../../source_base/libm/branred.cpp
136+
../../../source_base/libm/sincos.cpp
137+
../../../source_base/module_external/blas_connector_base.cpp ../../../source_base/module_external/blas_connector_vector.cpp ../../../source_base/module_external/blas_connector_matrix.cpp
138+
../../../source_base/module_fft/fft_bundle.cpp
139+
../../../source_base/module_fft/fft_cpu.cpp
140+
${FFT_SRC}
141+
)

‎source/source_hamilt/module_xc/test/test_xc6.cpp‎

Lines changed: 155 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,155 @@
1+
#include "../xc_functional.h"
2+
#include "../libxc_abacus.h"
3+
#include "gtest/gtest.h"
4+
#include "xctest.h"
5+
#include "../exx_info.h"
6+
#include <cmath>
7+
#include <iostream>
8+
#include <iomanip>
9+
10+
namespace ModuleBase
11+
{
12+
void WARNING_QUIT(const std::string &file,const std::string &description) {exit(1);}
13+
void TITLE(const std::string &class_function_name,bool disable){};
14+
void TITLE(const std::string &class_name,const std::string &function_name,bool disable){};
15+
}
16+
17+
namespace GlobalV
18+
{
19+
std::string BASIS_TYPE = "";
20+
bool CAL_STRESS = false;
21+
int CAL_FORCE = 0;
22+
int NSPIN = 1;
23+
}
24+
25+
namespace GlobalC
26+
{
27+
Exx_Info exx_info;
28+
}
29+
30+
class XCTest_SCANL_Laplacian : public XCTest
31+
{
32+
protected:
33+
double e_base, v1_base, v2_base, v3_base, vlapl_base;
34+
double e_modified, v1_modified, v2_modified, v3_modified, vlapl_modified;
35+
double e_scaled, v1_scaled, v2_scaled, v3_scaled, vlapl_scaled;
36+
37+
void SetUp()
38+
{
39+
XC_Functional::set_xc_type("MGGA_X_SCANL+MGGA_C_SCANL");
40+
41+
const double rho = 0.17E+01;
42+
const double grho = 0.81E-11;
43+
const double tau = 0.02403590412;
44+
const double lapl_base = 0.15E+01;
45+
double hybrid_alpha = 0.0;
46+
double hse_omega = 0.0;
47+
48+
XC_Functional_Libxc::tau_xc(
49+
XC_Functional::get_func_id(),
50+
rho, grho, lapl_base, tau,
51+
e_base, v1_base, v2_base, v3_base, vlapl_base, hybrid_alpha, hse_omega
52+
);
53+
54+
XC_Functional_Libxc::tau_xc(
55+
XC_Functional::get_func_id(),
56+
rho, grho, lapl_base + 1.0, tau,
57+
e_modified, v1_modified, v2_modified, v3_modified, vlapl_modified, hybrid_alpha, hse_omega
58+
);
59+
60+
XC_Functional_Libxc::tau_xc(
61+
XC_Functional::get_func_id(),
62+
rho, grho, 2.0 * lapl_base, tau,
63+
e_scaled, v1_scaled, v2_scaled, v3_scaled, vlapl_scaled, hybrid_alpha, hse_omega
64+
);
65+
}
66+
};
67+
68+
TEST_F(XCTest_SCANL_Laplacian, laplacian_affects_energy)
69+
{
70+
EXPECT_NE(e_base, e_modified);
71+
EXPECT_NE(e_base, e_scaled);
72+
73+
std::cout << std::scientific << std::setprecision(15);
74+
std::cout << "\n=== SCAN-L Laplacian Sensitivity Test ===" << std::endl;
75+
std::cout << "Base Laplacian: " << 0.15E+01 << std::endl;
76+
std::cout << " E_xc = " << e_base << std::endl;
77+
std::cout << " dE/dlapl ~= " << (e_modified - e_base) / 1.0 << std::endl;
78+
std::cout << "Modified Laplacian:" << 0.15E+01 + 1.0 << std::endl;
79+
std::cout << " E_xc = " << e_modified << std::endl;
80+
std::cout << " Delta E = " << e_modified - e_base << std::endl;
81+
std::cout << "Scaled Laplacian: " << 2.0 * 0.15E+01 << std::endl;
82+
std::cout << " E_xc = " << e_scaled << std::endl;
83+
std::cout << " Delta E = " << e_scaled - e_base << std::endl;
84+
std::cout << "=========================================" << std::endl;
85+
}
86+
87+
TEST(XC_ScanL_Reference, MGGA_X_SCANL_direct_libxc)
88+
{
89+
const double rho = 35.536521214608185;
90+
const double sigma = 1.149382202334535e+05;
91+
const double lapl = 8.411855859277239e+02;
92+
const double tau = 1.389887953757970e+02;
93+
94+
std::vector<int> func_id = {XC_MGGA_X_SCANL};
95+
std::vector<xc_func_type> funcs = XC_Functional_Libxc::init_func(func_id, XC_UNPOLARIZED, 0.0, 0.0);
96+
97+
double exc = 0.0, vrho = 0.0, vsigma = 0.0, vlapl = 0.0, vtau = 0.0;
98+
xc_mgga_exc_vxc(&funcs[0], 1, &rho, &sigma, &lapl, &tau,
99+
&exc, &vrho, &vsigma, &vlapl, &vtau);
100+
101+
EXPECT_NEAR(exc, -2.569101943337681e+00, 5.0e-04);
102+
EXPECT_NEAR(vrho, -2.912514021641506e+00, 1.0e-03);
103+
EXPECT_NEAR(vsigma, -8.289754258579972e-05, 1.0e-12);
104+
EXPECT_NEAR(vlapl, 3.971745124876924e-03, 1.0e-12);
105+
EXPECT_EQ(vtau, 0.0);
106+
107+
XC_Functional_Libxc::finish_func(funcs);
108+
}
109+
110+
TEST(XC_ScanL_Reference, MGGA_C_SCANL_direct_libxc)
111+
{
112+
const double rho = 35.536521214608185;
113+
const double sigma = 1.149382202334535e+05;
114+
const double lapl = 8.411855859277239e+02;
115+
const double tau = 1.389887953757970e+02;
116+
117+
std::vector<int> func_id = {XC_MGGA_C_SCANL};
118+
std::vector<xc_func_type> funcs = XC_Functional_Libxc::init_func(func_id, XC_UNPOLARIZED, 0.0, 0.0);
119+
120+
double exc = 0.0, vrho = 0.0, vsigma = 0.0, vlapl = 0.0, vtau = 0.0;
121+
xc_mgga_exc_vxc(&funcs[0], 1, &rho, &sigma, &lapl, &tau,
122+
&exc, &vrho, &vsigma, &vlapl, &vtau);
123+
124+
EXPECT_NEAR(exc, -5.250261832501450e-02, 1.0e-06);
125+
EXPECT_NEAR(vrho, -1.323740339176898e-01, 1.0e-05);
126+
EXPECT_NEAR(vsigma, 1.138021340988220e-05, 1.0e-10);
127+
EXPECT_NEAR(vlapl, -3.867564402989276e-04, 1.0e-10);
128+
EXPECT_EQ(vtau, 0.0);
129+
130+
XC_Functional_Libxc::finish_func(funcs);
131+
}
132+
133+
TEST(XC_ScanL_Reference, tau_xc_wrapper_libxc)
134+
{
135+
XC_Functional::set_xc_type("MGGA_X_SCANL+MGGA_C_SCANL");
136+
137+
const double rho = 35.536521214608185;
138+
const double grho = 1.149382202334535e+05;
139+
const double lapl = 8.411855859277239e+02;
140+
const double tau = 1.389887953757970e+02;
141+
142+
double sxc, v1xc, v2xc, v3xc, vlaplxc;
143+
double hybrid_alpha = 0.0;
144+
double hse_omega = 0.0;
145+
XC_Functional_Libxc::tau_xc(XC_Functional::get_func_id(), rho, grho, lapl, tau,
146+
sxc, v1xc, v2xc, v3xc, vlaplxc, hybrid_alpha, hse_omega);
147+
148+
double exc_x_ref = -2.569101943337681e+00;
149+
double exc_c_ref = -5.250261832501450e-02;
150+
double sxc_ref = (exc_x_ref + exc_c_ref) * rho;
151+
152+
EXPECT_NEAR(sxc, sxc_ref, 0.05);
153+
EXPECT_NE(v1xc, 0.0);
154+
EXPECT_NE(vlaplxc, 0.0);
155+
}

0 commit comments

Comments
 (0)