Skip to content

Fix SCAN-L Laplacian: port to new module_xc structure + rebase #7652

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

Closed
mintleaf84 wants to merge 1 commit into from
Closed

Fix SCAN-L Laplacian: port to new module_xc structure + rebase #7652

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
31 changes: 30 additions & 1 deletion source/source_estate/module_pot/pot_xc.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
#include "pot_xc.h"

#include "source_base/timer.h"
#include "source_base/constants.h"
#include "source_hamilt/module_xc/xc_functional.h"
#include "source_io/module_parameter/parameter.h"

Expand Down Expand Up @@ -30,13 +31,41 @@ void PotXC::cal_v_eff(const Charge*const chg, const UnitCell*const ucell, Module
#else
const double hse_omega = 0.0;
#endif
const std::tuple<double, double, ModuleBase::matrix, ModuleBase::matrix> etxc_vtxc_v
const std::tuple<double, double, ModuleBase::matrix, ModuleBase::matrix, ModuleBase::matrix> etxc_vtxc_v
= XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), nrxx_current, ucell->omega, ucell->tpiba, chg,
PARAM.inp.nspin, hybrid_alpha, hse_omega);
*(this->etxc_) = std::get<0>(etxc_vtxc_v);
*(this->vtxc_) = std::get<1>(etxc_vtxc_v);
v_eff += std::get<2>(etxc_vtxc_v);
*(this->vofk) = std::get<3>(etxc_vtxc_v);

// Apply Laplacian potential correction using FD kernel
const ModuleBase::matrix& voflapl = std::get<4>(etxc_vtxc_v);
const int ng = chg->rhopw->npw;
const int nrxx = chg->rhopw->nrxx;
const double tpiba2 = ucell->tpiba * ucell->tpiba;
std::vector<double> gg_fd = XC_Functional::compute_fd_gg(chg->rhopw);

std::vector<std::complex<double>> lapl_tmp(chg->rhopw->nmaxgr);
for(int is = 0; is < voflapl.nr; is++)
{
for(int ir = 0; ir < nrxx; ir++)
lapl_tmp[ir] = std::complex<double>(voflapl(is, ir), 0.0);
for(int ig = ng; ig < chg->rhopw->nmaxgr; ig++)
lapl_tmp[ig] = std::complex<double>(0.0, 0.0);
chg->rhopw->real2recip(lapl_tmp.data(), lapl_tmp.data());
for(int ig = 0; ig < ng; ig++)
{
lapl_tmp[ig] *= -gg_fd[ig] * tpiba2;
}
chg->rhopw->recip2real(lapl_tmp.data(), lapl_tmp.data());
for(int ir = 0; ir < nrxx; ir++)
{
double vlapl_corr = ModuleBase::e2 * lapl_tmp[ir].real();
v_eff(is, ir) += vlapl_corr;
*(this->vtxc_) += vlapl_corr * chg->rho[is][ir] * ucell->omega / chg->rhopw->nxyz;
}
}
#else
ModuleBase::WARNING_QUIT("v_of_rho", "to use mGGA, compile with LIBXC");
#endif
Expand Down
8 changes: 7 additions & 1 deletion 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 @@ -70,7 +70,7 @@ namespace XC_Functional_Libxc
const double hse_omega);

// for mGGA functional
extern std::tuple<double, double, ModuleBase::matrix, ModuleBase::matrix> v_xc_meta(
extern std::tuple<double, double, ModuleBase::matrix, ModuleBase::matrix, ModuleBase::matrix> v_xc_meta(
const std::vector<int> &func_id,
const int &nrxx, // number of real-space grid
const double &omega, // volume of cell
Expand Down Expand Up @@ -214,11 +214,13 @@ 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,
double &v2xc,
double &v3xc,
double &vlaplxc,
const double &hybrid_alpha,
const double &hse_omega);

Expand All @@ -228,6 +230,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 All @@ -238,6 +242,8 @@ namespace XC_Functional_Libxc
double &v2xcud,
double &v3xcup,
double &v3xcdw,
double &vlaplxcup,
double &vlaplxcdw,
const double &hybrid_alpha,
const double &hse_omega);

Expand Down
22 changes: 18 additions & 4 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,19 +17,20 @@ 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,
double& v2xc,
double& v3xc,
double& vlaplxc,
const double& hybrid_alpha,
const double& hse_omega)
{
double s = 0.0;
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 All @@ -41,6 +42,7 @@ void XC_Functional_Libxc::tau_xc(
v1xc = 0.0;
v2xc = 0.0;
v3xc = 0.0;
vlaplxc = 0.0;

for (xc_func_type& func : funcs)
{
Expand All @@ -52,12 +54,14 @@ void XC_Functional_Libxc::tau_xc(
v1 *= (1.0 - hybrid_alpha);
v2 *= (1.0 - hybrid_alpha);
v3 *= (1.0 - hybrid_alpha);
vlapl_rho *= (1.0 - hybrid_alpha);
}
#endif
sxc += s * rho;
v2xc += v2 * 2.0;
v1xc += v1;
v3xc += v3;
vlaplxc += vlapl_rho;
}
XC_Functional_Libxc::finish_func(funcs);

Expand All @@ -71,6 +75,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 All @@ -81,6 +87,8 @@ void XC_Functional_Libxc::tau_xc_spin(
double& v2xcud,
double& v3xcup,
double& v3xcdw,
double& vlaplxcup,
double& vlaplxcdw,
const double& hybrid_alpha,
const double& hse_omega)
{
Expand All @@ -92,10 +100,13 @@ void XC_Functional_Libxc::tau_xc_spin(
v2xcud = 0.0;
v3xcup = 0.0;
v3xcdw = 0.0;
vlaplxcup = 0.0;
vlaplxcdw = 0.0;

const std::array<double, 2> rho = {rhoup, rhodw};
const std::array<double, 3> grho = {gdr1.norm2(), gdr1 * gdr2, gdr2.norm2()};
const std::array<double, 2> tau = {tauup, taudw};
const std::array<double, 2> lapl = {laplup, lapldw};

std::vector<xc_func_type> funcs = XC_Functional_Libxc::init_func(
/* func_id = */ func_id,
Expand Down Expand Up @@ -125,12 +136,11 @@ 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> vlapl = {0.0, 0.0};
std::array<double, 2> vlapl_out = {0.0, 0.0};
std::array<double, 3> v2xc = {0.0, 0.0, 0.0};
// call Libxc function: xc_mgga_exc_vxc
xc_mgga_exc_vxc(&func, 1, rho.data(), grho.data(), lapl.data(), tau.data(), &s,
v1xc.data(), v2xc.data(), vlapl.data(), v3xc.data());
v1xc.data(), v2xc.data(), vlapl_out.data(), v3xc.data());

#ifdef __EXX
if (func.info->number == XC_MGGA_X_SCAN && XC_Functional::get_func_type() == 5)
Expand All @@ -143,6 +153,8 @@ void XC_Functional_Libxc::tau_xc_spin(
v2xc[2] *= (1.0 - hybrid_alpha);
v3xc[0] *= (1.0 - hybrid_alpha);
v3xc[1] *= (1.0 - hybrid_alpha);
vlapl_out[0] *= (1.0 - hybrid_alpha);
vlapl_out[1] *= (1.0 - hybrid_alpha);
}
#endif

Expand All @@ -154,6 +166,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];
vlaplxcup += vlapl_out[0] * sgn[0];
vlaplxcdw += vlapl_out[1] * sgn[1];
}
}

Expand Down
35 changes: 30 additions & 5 deletions 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 @@ -15,6 +15,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 @@ -211,8 +212,8 @@ std::tuple<double,double,ModuleBase::matrix> XC_Functional_Libxc::v_xc_libxc( /

//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

// [etxc, vtxc, v, vofk] = XC_Functional::v_xc(...)
std::tuple<double,double,ModuleBase::matrix,ModuleBase::matrix> XC_Functional_Libxc::v_xc_meta(
// [etxc, vtxc, v, vofk, voflapl] = XC_Functional::v_xc(...)
std::tuple<double,double,ModuleBase::matrix,ModuleBase::matrix,ModuleBase::matrix> XC_Functional_Libxc::v_xc_meta(
const std::vector<int> &func_id,
const int &nrxx, // number of real-space grid
const double &omega, // volume of cell
Expand All @@ -232,6 +233,7 @@ std::tuple<double,double,ModuleBase::matrix,ModuleBase::matrix> XC_Functional_Li
double vtxc = 0.0;
ModuleBase::matrix v(nspin,nrxx);
ModuleBase::matrix vofk(nspin,nrxx);
ModuleBase::matrix voflapl(nspin,nrxx);

//----------------------------------------------------------
// xc_func_type is defined in Libxc package
Expand All @@ -250,6 +252,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 @@ -328,7 +351,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 Expand Up @@ -441,7 +464,7 @@ std::tuple<double,double,ModuleBase::matrix,ModuleBase::matrix> XC_Functional_Li
}
vtxc -= rvtxc;

//process vtau
//process vtau and vlapl
#ifdef _OPENMP
#pragma omp parallel for collapse(2) schedule(static, 1024)
#endif
Expand All @@ -453,9 +476,11 @@ std::tuple<double,double,ModuleBase::matrix,ModuleBase::matrix> XC_Functional_Li
if (func.info->number == XC_MGGA_X_SCAN && XC_Functional::get_func_type() == 5)
{
vtau[ir*nspin+is] *= (1.0 - XC_Functional::get_hybrid_alpha());
vlapl[ir*nspin+is] *= (1.0 - XC_Functional::get_hybrid_alpha());
}
#endif
vofk(is,ir) += vtau[ir*nspin+is] * sgn[ir*nspin+is];
voflapl(is,ir) += vlapl[ir*nspin+is] * sgn[ir*nspin+is];
}
}
}
Expand All @@ -474,7 +499,7 @@ std::tuple<double,double,ModuleBase::matrix,ModuleBase::matrix> XC_Functional_Li
XC_Functional_Libxc::finish_func(funcs);

ModuleBase::timer::end("XC_Functional_Libxc","v_xc_meta");
return std::make_tuple( etxc, vtxc, std::move(v), std::move(vofk) );
return std::make_tuple( etxc, vtxc, std::move(v), std::move(vofk), std::move(voflapl) );
}

#endif
17 changes: 17 additions & 0 deletions source/source_hamilt/module_xc/libxc_setup.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -180,6 +180,23 @@ XC_Functional_Libxc::set_xc_type_libxc(const std::string& xc_func_in)
ModuleBase::WARNING_QUIT("XC_Functional::set_xc_type_libxc", message);
}

// warn if any functional needs Laplacian of density
{
std::vector<xc_func_type> tmp_funcs = XC_Functional_Libxc::init_func(func_id, XC_UNPOLARIZED, 0.0, 0.0);
for (auto& f : tmp_funcs)
{
if (f.info->flags & XC_FLAGS_NEEDS_LAPLACIAN)
{
std::cout << " WARNING: XC functional \"" << f.info->name
<< "\" requires Laplacian of density (nabla^2 rho)."
<< " This may require a higher energy cutoff for numerical stability."
<< std::endl;
break;
}
}
XC_Functional_Libxc::finish_func(tmp_funcs);
}

// return
return std::make_pair(func_type, func_id);
}
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 @@ -83,3 +83,16 @@ AddTest(
../../../source_base/module_fft/fft_cpu.cpp
${FFT_SRC}
)

AddTest(
TARGET MODULE_HAMILT_XCTest_SCANL_LAPL
LIBS parameter MPI::MPI_CXX Libxc::xc
SOURCES test_xc6.cpp ../xc_functional.cpp ../xc_lda_wrap.cpp
../xc_gga_wrap.cpp
../libxc_setup.cpp
../libxc_lda_wrap.cpp
../libxc_gga_wrap.cpp
../libxc_mgga_wrap.cpp
../xc_gga_corr.cpp ../xc_lda_corr.cpp
../xc_gga_exch.cpp ../xc_lda_exch.cpp ../xc_hcth.cpp
)
5 changes: 3 additions & 2 deletions 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 @@ -45,10 +45,11 @@ class XCTest_SCAN : public XCTest

for(int i=0;i<5;i++)
{
double e,v,v1,v2,v3;
double e,v,v1,v2,v3,vlapl;
double hybrid_alpha = 0.0;
double hse_omega = 0.0;
XC_Functional_Libxc::tau_xc(XC_Functional::get_func_id(), rho[i],grho[i],tau[i],e,v1,v2,v3,hybrid_alpha, hse_omega);
double lapl_rho = 0.0; // SCAN does not use Laplacian; value has no effect on result
XC_Functional_Libxc::tau_xc(XC_Functional::get_func_id(), rho[i],grho[i],lapl_rho,tau[i],e,v1,v2,v3,vlapl,hybrid_alpha, hse_omega);
e_.push_back(e);
v1_.push_back(v1);
v2_.push_back(v2);
Expand Down
4 changes: 3 additions & 1 deletion source/source_hamilt/module_xc/test/test_xc5.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -296,19 +296,21 @@ class XCTest_VXC_meta : public XCTest

const double hybrid_alpha = XC_Functional::get_hybrid_alpha();
const double hse_omega = XC_Functional::get_hse_omega();
std::tuple<double, double, ModuleBase::matrix, ModuleBase::matrix> etxc_vtxc_v
std::tuple<double, double, ModuleBase::matrix, ModuleBase::matrix, ModuleBase::matrix> etxc_vtxc_v
= XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), rhopw.nrxx,ucell.omega,ucell.tpiba,&chr,nspin1, hybrid_alpha, hse_omega);
et1 = std::get<0>(etxc_vtxc_v);
vt1 = std::get<1>(etxc_vtxc_v);
v1 = std::get<2>(etxc_vtxc_v);
vtau1 = std::get<3>(etxc_vtxc_v);
ModuleBase::matrix vlapl1 = std::get<4>(etxc_vtxc_v);

etxc_vtxc_v
= XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), rhopw.nrxx,ucell.omega,ucell.tpiba,&chr,nspin2, hybrid_alpha, hse_omega);
et2 = std::get<0>(etxc_vtxc_v);
vt2 = std::get<1>(etxc_vtxc_v);
v2 = std::get<2>(etxc_vtxc_v);
vtau2 = std::get<3>(etxc_vtxc_v);
ModuleBase::matrix vlapl2 = std::get<4>(etxc_vtxc_v);
}
};

Expand Down
Loading
Loading