Skip to content

fix(xc): implement proper Laplacian calculation for SCAN-L functional #7457

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
mintleaf84 wants to merge 2 commits into
base: develop
Choose a base branch
from
Open

fix(xc): implement proper Laplacian calculation for SCAN-L functional #7457

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
0f94387
fix(xc): implement proper Laplacian calculation for SCAN-L functional
mintleaf84 Jun 9, 2026
4a546e6
fix: use SCAN-L instead of SCAN for laplacian sensitivity test
mintleaf84 Jun 9, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions source/source_hamilt/module_xc/libxc_abacus.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -198,6 +198,7 @@ namespace XC_Functional_Libxc
const std::vector<int> &func_id,
const double &rho,
const double &grho,
const double &lapl_rho,
const double &atau,
double &sxc,
double &v1xc,
Expand All @@ -211,6 +212,8 @@ namespace XC_Functional_Libxc
double rhodw,
ModuleBase::Vector3<double> gdr1,
ModuleBase::Vector3<double> gdr2,
double laplup,
double lapldw,
double tauup,
double taudw,
double &sxc,
Expand Down
6 changes: 4 additions & 2 deletions source/source_hamilt/module_xc/libxc_mgga_wrap.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@ void XC_Functional_Libxc::tau_xc(
const std::vector<int>& func_id,
const double& rho,
const double& grho,
const double& lapl_rho,
const double& atau,
double& sxc,
double& v1xc,
Expand All @@ -28,7 +29,6 @@ void XC_Functional_Libxc::tau_xc(
double v1 = 0.0;
double v2 = 0.0;
double v3 = 0.0;
double lapl_rho = grho;
double vlapl_rho = 0.0;
std::vector<xc_func_type> funcs = XC_Functional_Libxc::init_func(
/* func_id = */ func_id,
Expand Down Expand Up @@ -68,6 +68,8 @@ void XC_Functional_Libxc::tau_xc_spin(
double rhodw,
ModuleBase::Vector3<double> gdr1,
ModuleBase::Vector3<double> gdr2,
double laplup,
double lapldw,
double tauup,
double taudw,
double& sxc,
Expand Down Expand Up @@ -119,7 +121,7 @@ void XC_Functional_Libxc::tau_xc_spin(
double s = 0.0;
std::array<double, 2> v1xc = {0.0, 0.0};
std::array<double, 2> v3xc = {0.0, 0.0};
std::array<double, 2> lapl = {0.0, 0.0};
std::array<double, 2> lapl = {laplup, lapldw};
std::array<double, 2> vlapl = {0.0, 0.0};
std::array<double, 3> v2xc = {0.0, 0.0, 0.0};
// call Libxc function: xc_mgga_exc_vxc
Expand Down
24 changes: 23 additions & 1 deletion source/source_hamilt/module_xc/libxc_pot.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@
#include <xc.h>

#include <vector>
#include <complex>

std::tuple<double,double,ModuleBase::matrix> XC_Functional_Libxc::v_xc_libxc( // Peize Lin update for nspin==4 at 2023.01.14
const std::vector<int> &func_id,
Expand Down Expand Up @@ -237,6 +238,27 @@ std::tuple<double,double,ModuleBase::matrix,ModuleBase::matrix> XC_Functional_Li
= XC_Functional_Libxc::cal_gdr(nspin, nrxx, rho, tpiba, chr);
const std::vector<double> sigma = XC_Functional_Libxc::convert_sigma(gdr);

// compute laplacian for mGGA functionals
std::vector<double> lapl(nrxx * nspin, 0.0);
{
std::vector<std::complex<double>> rhog_tmp(chr->rhopw->npw);
std::vector<double> rhor(nrxx);
for( int is=0; is<nspin; ++is )
{
for( int ir=0; ir<nrxx; ++ir )
{
rhor[ir] = rho[ir*nspin+is];
}
chr->rhopw->real2recip(rhor.data(), rhog_tmp.data());
std::vector<double> lapl_spin(nrxx);
XC_Functional::laplacian_rho(rhog_tmp.data(), lapl_spin.data(), chr->rhopw, tpiba);
for( int ir=0; ir<nrxx; ++ir )
{
lapl[ir*nspin+is] = lapl_spin[ir];
}
}
}

//converting kin_r
std::vector<double> kin_r;
kin_r.resize(nrxx*nspin);
Expand Down Expand Up @@ -315,7 +337,7 @@ std::tuple<double,double,ModuleBase::matrix,ModuleBase::matrix> XC_Functional_Li
nrxx_thread,
rho.data() + ir_start * nspin,
sigma.data() + ir_start * ((1==nspin)?1:3),
sigma.data() + ir_start * ((1==nspin)?1:3),
lapl.data() + ir_start * nspin,
kin_r.data() + ir_start * nspin,
exc.data() + ir_start,
vrho.data() + ir_start * nspin,
Expand Down
13 changes: 13 additions & 0 deletions source/source_hamilt/module_xc/test/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -82,4 +82,17 @@ AddTest(
../../../source_base/module_fft/fft_bundle.cpp
../../../source_base/module_fft/fft_cpu.cpp
${FFT_SRC}
)

AddTest(
TARGET MODULE_HAMILT_XCTest_SCANL_LAPL
LIBS parameter MPI::MPI_CXX Libxc::xc
SOURCES test_xc6.cpp ../xc_functional.cpp ../xc_lda_wrap.cpp
../xc_gga_wrap.cpp
../libxc_setup.cpp
../libxc_lda_wrap.cpp
../libxc_gga_wrap.cpp
../libxc_mgga_wrap.cpp
../xc_gga_corr.cpp ../xc_lda_corr.cpp
../xc_gga_exch.cpp ../xc_lda_exch.cpp ../xc_hcth.cpp
)
3 changes: 2 additions & 1 deletion source/source_hamilt/module_xc/test/test_xc4.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,8 @@ class XCTest_SCAN : public XCTest
{
double e,v,v1,v2,v3;
double hybrid_alpha = 0.0;
XC_Functional_Libxc::tau_xc(XC_Functional::get_func_id(), rho[i],grho[i],tau[i],e,v1,v2,v3,hybrid_alpha);
double lapl_rho = grho[i];
XC_Functional_Libxc::tau_xc(XC_Functional::get_func_id(), rho[i],grho[i],lapl_rho,tau[i],e,v1,v2,v3,hybrid_alpha);
e_.push_back(e);
v1_.push_back(v1);
v2_.push_back(v2);
Expand Down
84 changes: 84 additions & 0 deletions source/source_hamilt/module_xc/test/test_xc6.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,84 @@
#include "../xc_functional.h"
#include "../libxc_abacus.h"
#include "gtest/gtest.h"
#include "xctest.h"
#include "../exx_info.h"
#include <cmath>
#include <iostream>
#include <iomanip>

namespace ModuleBase
{
void WARNING_QUIT(const std::string &file,const std::string &description) {exit(1);}
void TITLE(const std::string &class_function_name,bool disable){};
void TITLE(const std::string &class_name,const std::string &function_name,bool disable){};
}

namespace GlobalV
{
std::string BASIS_TYPE = "";
bool CAL_STRESS = false;
int CAL_FORCE = 0;
int NSPIN = 1;
}

namespace GlobalC
{
Exx_Info exx_info;
}

class XCTest_SCANL_Laplacian : public XCTest
{
protected:
double e_base, v1_base, v2_base, v3_base;
double e_modified, v1_modified, v2_modified, v3_modified;
double e_scaled, v1_scaled, v2_scaled, v3_scaled;

void SetUp()
{
XC_Functional::set_xc_type("MGGA_X_SCANL+MGGA_C_SCANL");

const double rho = 0.17E+01;
const double grho = 0.81E-11;
const double tau = 0.02403590412;
const double lapl_base = 0.15E+01;
double hybrid_alpha = 0.0;

XC_Functional_Libxc::tau_xc(
XC_Functional::get_func_id(),
rho, grho, lapl_base, tau,
e_base, v1_base, v2_base, v3_base, hybrid_alpha
);

XC_Functional_Libxc::tau_xc(
XC_Functional::get_func_id(),
rho, grho, lapl_base + 1.0, tau,
e_modified, v1_modified, v2_modified, v3_modified, hybrid_alpha
);

XC_Functional_Libxc::tau_xc(
XC_Functional::get_func_id(),
rho, grho, 2.0 * lapl_base, tau,
e_scaled, v1_scaled, v2_scaled, v3_scaled, hybrid_alpha
);
}
};

TEST_F(XCTest_SCANL_Laplacian, laplacian_affects_energy)
{
EXPECT_NE(e_base, e_modified);
EXPECT_NE(e_base, e_scaled);

std::cout << std::scientific << std::setprecision(15);
std::cout << "\n=== SCAN Laplacian Sensitivity Test ===" << std::endl;
std::cout << "Base Laplacian: " << 0.15E+01 << std::endl;
std::cout << " E_xc = " << e_base << std::endl;
std::cout << " dE/dlapl ≈ " << (e_modified - e_base) / 1.0 << std::endl;
std::cout << "Modified Laplacian:" << 0.15E+01 + 1.0 << std::endl;
std::cout << " E_xc = " << e_modified << std::endl;
std::cout << " Delta E = " << e_modified - e_base << std::endl;
std::cout << "Scaled Laplacian: " << 2.0 * 0.15E+01 << std::endl;
std::cout << " E_xc = " << e_scaled << std::endl;
std::cout << " Delta E = " << e_scaled - e_base << std::endl;
std::cout << "========================================" << std::endl;
}
6 changes: 6 additions & 0 deletions source/source_hamilt/module_xc/xc_functional.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -229,6 +229,12 @@ class XC_Functional
const ModulePW::PW_Basis* rho_basis,
const double tpiba);

static void laplacian_rho(
const std::complex<double>* rhog,
double* lapl,
const ModulePW::PW_Basis* rho_basis,
const double tpiba);

static void noncolin_rho(
double* rhoout1,
double* rhoout2,
Expand Down
78 changes: 74 additions & 4 deletions source/source_hamilt/module_xc/xc_grad.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -90,6 +90,8 @@ void XC_Functional::gradcorr(
double* neg = nullptr;
double** vsave = nullptr;
double** vgg = nullptr;
double* lapl1 = nullptr;
double* lapl2 = nullptr;

// for spin unpolarized case,
// calculate the gradient of (rho_core+rho) in reciprocal space.
Expand Down Expand Up @@ -118,6 +120,12 @@ void XC_Functional::gradcorr(

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

if(func_type == 3 || func_type == 5)
{
lapl1 = new double[rhopw->nrxx];
XC_Functional::laplacian_rho(rhogsum1, lapl1, rhopw, ucell->tpiba);
}

// for spin polarized case;
// calculate the gradient of (rho_core+rho) in reciprocal space.
if(PARAM.inp.nspin==2)
Expand Down Expand Up @@ -146,6 +154,12 @@ void XC_Functional::gradcorr(
}

XC_Functional::grad_rho( rhogsum2 , gdr2, rhopw, ucell->tpiba);

if(func_type == 3 || func_type == 5)
{
lapl2 = new double[rhopw->nrxx];
XC_Functional::laplacian_rho(rhogsum2, lapl2, rhopw, ucell->tpiba);
}
}

if(PARAM.inp.nspin == 4&&(PARAM.globalv.domag||PARAM.globalv.domag_z))
Expand Down Expand Up @@ -222,6 +236,20 @@ void XC_Functional::gradcorr(

XC_Functional::grad_rho( rhogsum1 , gdr1, rhopw, ucell->tpiba);
XC_Functional::grad_rho( rhogsum2 , gdr2, rhopw, ucell->tpiba);

if(func_type == 3 || func_type == 5)
{
if(lapl1 == nullptr)
{
lapl1 = new double[rhopw->nrxx];
}
XC_Functional::laplacian_rho(rhogsum1, lapl1, rhopw, ucell->tpiba);
if(lapl2 == nullptr)
{
lapl2 = new double[rhopw->nrxx];
}
XC_Functional::laplacian_rho(rhogsum2, lapl2, rhopw, ucell->tpiba);
}
}

const double epsr = 1.0e-6;
Expand Down Expand Up @@ -297,7 +325,8 @@ void XC_Functional::gradcorr(
#ifdef __EXX
hybrid_alpha = GlobalC::exx_info.info_global.hybrid_alpha;
#endif
XC_Functional_Libxc::tau_xc( func_id, arho, grho2a, atau, sxc, v1xc, v2xc, v3xc, hybrid_alpha);
double lapl_val = (lapl1 != nullptr) ? lapl1[ir] : 0.0;
XC_Functional_Libxc::tau_xc( func_id, arho, grho2a, lapl_val, atau, sxc, v1xc, v2xc, v3xc, hybrid_alpha);
}
else
{
Expand Down Expand Up @@ -366,10 +395,12 @@ void XC_Functional::gradcorr(
#ifdef __EXX
hybrid_alpha = GlobalC::exx_info.info_global.hybrid_alpha;
#endif
double laplup_val = (lapl1 != nullptr) ? lapl1[ir] : 0.0;
double lapldw_val = (lapl2 != nullptr) ? lapl2[ir] : 0.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, hybrid_alpha);
laplup_val, lapldw_val, atau1, atau2, sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud, v3xcup, v3xcdw, hybrid_alpha);
}
else
{
Expand Down Expand Up @@ -653,14 +684,18 @@ void XC_Functional::gradcorr(
delete[] rhotmp2;
delete[] rhogsum2;
delete[] gdr2;
delete[] lapl2;
if(!is_stress)
{
delete[] h2;
}
delete[] lapl1;
}
if(PARAM.inp.nspin == 4 && (PARAM.globalv.domag||PARAM.globalv.domag_z))
else if(PARAM.inp.nspin == 4 && (PARAM.globalv.domag||PARAM.globalv.domag_z))
{
delete[] neg;
delete[] lapl1;
delete[] lapl2;
if(!is_stress)
{
for(int i=0; i<nspin0; i++)
Expand All @@ -679,6 +714,10 @@ void XC_Functional::gradcorr(
delete[] rhogsum2;
delete[] gdr2;
}
else
{
delete[] lapl1;
}

return;
}
Expand Down Expand Up @@ -816,11 +855,42 @@ void XC_Functional::grad_dot(
dh[ir] = aux[ir].real() * tpiba;
}

delete[] aux;
delete[] aux;
delete[] gaux;
return;
}


void XC_Functional::laplacian_rho(
const std::complex<double>* rhog,
double* lapl,
const ModulePW::PW_Basis* rho_basis,
const double tpiba)
{
std::complex<double>* lapl_tmp = new std::complex<double>[rho_basis->nmaxgr];

for(int ir=0; ir<rho_basis->nrxx; ir++)
{
lapl[ir] = 0.0;
}

for(int i=0; i<3; i++)
{
for(int ig=0; ig<rho_basis->npw; ig++)
{
lapl_tmp[ig] = -rhog[ig] * rho_basis->gcar[ig][i] * rho_basis->gcar[ig][i];
}
rho_basis->recip2real(lapl_tmp, lapl_tmp);
for(int ir=0; ir<rho_basis->nrxx; ir++)
{
lapl[ir] += lapl_tmp[ir].real() * tpiba * tpiba;
}
}

delete[] lapl_tmp;
}


void XC_Functional::noncolin_rho(
double *rhoout1,
double *rhoout2,
Expand Down
Loading