Skip to content

Commit b8e9264

Browse files
Flying-dragon-boxingzhangzh-pkuWHUweiqingzhoukirk0830haozhihan
authored
Feature: Support linear combination of coulomb_param for EXX PW (#6371)
* feat pexsi * fix : diag not completed * feat * feat: pexsi hsolver * CMake building implemented * Works * adapt to the new container * Turn off USE_PEXSI * Update LibRI to 553c91c * modify include files * namespace-ize * new inputs added * Configure Makefile Compiling, fix typos * Fix Makefile Intel toolchains compile errors * Fix even more PEXSI related Makefile compiling issues * Modify inputs and update to latest version (#2) * run INPUT.Default() in every process in InputParaTest (#3490) Co-authored-by: kirk0830 <[email protected]> * add blas support for FindLAPACK.cmake (#3497) * more unittest of QO: towards orbital selection (#3499) * Fix: fix bug in mulliken charge calculation (#3503) * fix phase * fix case test * Refactor: namespace Conv_Coulomb_Pot_K (#3446) * Refactor: namespace Conv_Coulomb_Pot_K * Refactor: namespace Conv_Coulomb_Pot_K --------- Co-authored-by: wqzhou <[email protected]> * enable the computation of all zeros in one function call (#3449) Co-authored-by: wqzhou <[email protected]> * replace ios.eof() by ios.good() to avoid meeting badbit and failbit in reading STRU (#3506) * Build: add ccache to accelerate the testing process (#3509) * Build: add ccache to accelerate the testing process * Update test.yml * Update test.yml * Update test.yml * Docs: to avoid the misunderstanding in docs (#3518) * to avoid the misunderstanding in docs * Update docs/quick_start/hands_on.md Co-authored-by: Chun Cai <[email protected]> --------- Co-authored-by: Chun Cai <[email protected]> * Docs: fix a missing depencency in conda build env (#3508) * Feature: Add ENABLE_RAPIDJSON option to control the output of abacus.json (#3519) Add ENABLE_RAPIDJSON option to control the output of abacus.json * Feature: add python wrapper for math sphbes (#3475) * recommit for review * add python wrapper * remove timer since performace tests add * Feature: support segment split in kline mode in KPT file and `out_band` band output precision control, `8` as default (#3493) * add precision control * correct serial version of nscf_band function * fix issue 3482 * update unit and integrated test * update document * correct unittest and make compatible with false and true * fix: bug in Autotest.sh when result.ref has no totaltimeref (#3523) * Fix : unit test of module_xc (#3524) * Fix: omit small magnetic moments to avoid numerical instability (#3530) * update deltalambda * avoid numerical error in orbMulP * add constrain on Mi * change case reference value * Fix: fix multiple compiler warnings (#3515) * Fix: add noreturn attribute to warning_quit * Add type conversion * fix string literal * fix small number trunctuation * Fix system call returned value not checked * fix missing braket * Refactor parameter_pool.cpp and parameter_pool.h * remove duplicated return statements * Change WARNING_QUIT occurances in tests * Add warning message to help debug UT * output the default precision flag (#3496) Co-authored-by: kirk0830 <[email protected]> * Build: Improving CMake performance for finding LibXC and ELPA (#3478) * Fix for finding LibXC and ELPA * For compatibility to previous routines * syntax fix for FindELPA.cmake * Update cmake/FindELPA.cmake Co-authored-by: Chun Cai <[email protected]> * Using CMake interface as default for finding LibXC * update docs * fix for FindLibxc: changing imcompatible if statement * fix for FindLibxc: changing imcompatible if statement * fix for FindLibxc: changing imcompatible if statement * update docs for installing pkg-config * Update FindLibxc.cmake * Update FindLibxc.cmake * remove previous LibXC routine in CMakeLists.txt Co-authored-by: Chun Cai <[email protected]> * Update easy_install.md with Makefile-built LibXC supported * Update easy_install.md to include different behavior in different version on finding ELPA --------- Co-authored-by: Chun Cai <[email protected]> * Docs: correct some docs about mp2 smearing method (#3533) * correct some docs about mp2 smearing method * add docs about mv method * Feature : printing band density (#3501) Co-authored-by: wenfei-li <[email protected]> Co-authored-by: wqzhou <[email protected]> * add some docs for PR#3501 (#3537) * Feature: enable restart charge density mixing during SCF (#3542) * add a new parameter mixing_restart * do not update rho if iter==mixing_restart * do not update rho if iter==mixing_restart-1 * reset mix and rho_mdata if iter==mixing_restart * fix SCF exit directly since drho=0 if iter=GlobalV::MIXING_RESTART * re-set_mixing in eachiterinit for PW and LCAO * enable SCF restarts in esolver_ks::RUN * add some UnitTests * add some Docs * new inputs added * Update input-main.md (#3551) Solve the format problem mentioned in issue 3543 * Build: fix compatibility issue against toolchain install (#3540) * Fix for finding LibXC and ELPA * For compatibility to previous routines * syntax fix for FindELPA.cmake * Update cmake/FindELPA.cmake Co-authored-by: Chun Cai <[email protected]> * Using CMake interface as default for finding LibXC * update docs * fix for FindLibxc: changing imcompatible if statement * fix for FindLibxc: changing imcompatible if statement * fix for FindLibxc: changing imcompatible if statement * update docs for installing pkg-config * Update FindLibxc.cmake * Update FindLibxc.cmake * remove previous LibXC routine in CMakeLists.txt Co-authored-by: Chun Cai <[email protected]> * Update easy_install.md with Makefile-built LibXC supported * Update easy_install.md to include different behavior in different version on finding ELPA * fix compatibility issue against toolchain * Change default ELPA install routine to old one --------- Co-authored-by: Chun Cai <[email protected]> * Test: Configure performance tests for math libraries (#3511) * add performace test of sphbes functions. * fix benchmark cmake errors * add dependencies for docker * update docs * add performance tests for sphbes * add google benchmark * rewrite benchmark tests in fixtures * disable internal testing in benchmark * merge benchmark into integration test --------- Co-authored-by: StarGrys <[email protected]> * Configure Makefile Compiling, fix typos * Fix Makefile Intel toolchains compile errors * Fix even more PEXSI related Makefile compiling issues * Update hsolver_pw.cpp (#3556) when use_uspp==false, overlap matrix should be E. * Fix: cuda build target (#3276) * Fix: cuda buid target * Update CMakeLists.txt --------- Co-authored-by: Denghui Lu <[email protected]> --------- Co-authored-by: wqzhou <[email protected]> Co-authored-by: kirk0830 <[email protected]> Co-authored-by: Haozhi Han <[email protected]> Co-authored-by: Zhao Tianqi <[email protected]> Co-authored-by: PeizeLin <[email protected]> Co-authored-by: jinzx10 <[email protected]> Co-authored-by: Chun Cai <[email protected]> Co-authored-by: Peng Xingliang <[email protected]> Co-authored-by: Jie Li <[email protected]> Co-authored-by: Wenfei Li <[email protected]> Co-authored-by: Denghui Lu <[email protected]> Co-authored-by: YI Zeping <[email protected]> Co-authored-by: wenfei-li <[email protected]> Co-authored-by: jingan-181 <[email protected]> Co-authored-by: StarGrys <[email protected]> Co-authored-by: Haozhi Han <[email protected]> * Revert "Modify inputs and update to latest version" * Update FindPEXSI.cmake to fix Comments * Fix CI errors * Fix CI Errors and Merge with Upstream * Resolve Pull Request Reviews * Fix parallel communication related issue * Fix vars in Makefile.vars, add input tests and comments for pexsi vars * Fix nspin > 1 cases * Improvement: take calculated mu as new initial guess, may slightly improve performance * Fix mistakes in the last commit * Fix: params and features - set default pexsi_temp - fix md in pexsi * fix empty lines * Fix: move params to pexsi_solver, rename USE_PEXSI to ENABLE_PEXSI * Docs: added docs for pexsi inputs * Fix unit test issues in input_conv * Change default pexsi_npole from 80 to 40 * Place pexsi_EDM in DensityMatrix, set size of pexsi_dm = 1 when GlobalV::NSPIN==4, and add comments for dmToRho * An unit test added for DiagoPexsi * modify for changed gint interface * correct nspin related behaviors * add efermi passthrough * Revert "add efermi passthrough" This reverts commit d7b402d. * commits to resolve conversations related to codes * DM and EDM pointers in pexsi now handled by diagopexsi, and copying h s matrices no longer needed * add pexsi examples * fix pexsi unit test (original version shouldn't run) * add building docs for pexsi * set cxx standard to c++14, which is required in make_unique * Fix: Fix typo related to pexsi * update to PPEXSIDFTDriver2 * default npoints to 1, so single core pexsi will work * Fix Compile errors * refactor to abandon `pdiagh` * Fix mu_buffer and nspin * Updates with latest * Refactor: in ESolver_KS_PW, calculate deband in iter_finish, not in hamilt2density * Fix: make files in consistent with upstream * Refactor * Refactor * Refactor * Refactor * Refactor * Refactor: fix unit test * Refactor: fix unit test * Refactor: fix unit test * Refactor: fix unit test * Refactor: Remove set kvec funcs in `K_Vectors` * Refactor: Remove final_scf * Refactor: Fix kvecc2d/d2c * Fix: Tests * Fix: Tests * Fix: Tests * Fix: Tests * Refactor: Final? * Fix * Fix * Fix * Fix * Fix: Compile Error on CUDA > 12.9 * Fix: Compile Error on CUDA > 12.9 * Feature: Support linear combination of coulomb_param for EXX PW * Fix: Fix compile issue --------- Co-authored-by: zhangzhihao <[email protected]> Co-authored-by: zhangzh-pku <[email protected]> Co-authored-by: wqzhou <[email protected]> Co-authored-by: kirk0830 <[email protected]> Co-authored-by: Haozhi Han <[email protected]> Co-authored-by: Zhao Tianqi <[email protected]> Co-authored-by: PeizeLin <[email protected]> Co-authored-by: jinzx10 <[email protected]> Co-authored-by: Chun Cai <[email protected]> Co-authored-by: Peng Xingliang <[email protected]> Co-authored-by: Jie Li <[email protected]> Co-authored-by: Wenfei Li <[email protected]> Co-authored-by: Denghui Lu <[email protected]> Co-authored-by: YI Zeping <[email protected]> Co-authored-by: wenfei-li <[email protected]> Co-authored-by: jingan-181 <[email protected]> Co-authored-by: StarGrys <[email protected]> Co-authored-by: Haozhi Han <[email protected]> Co-authored-by: Mohan Chen <[email protected]>
1 parent db41600 commit b8e9264

File tree

3 files changed

+145
-79
lines changed

3 files changed

+145
-79
lines changed

‎source/source_io/input_conv.cpp‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -410,6 +410,16 @@ void Input_Conv::Convert()
410410
{"singularity_correction", PARAM.inp.exx_singularity_correction} }};
411411
}
412412
}
413+
else if(PARAM.inp.basis_type == "pw")
414+
{
415+
GlobalC::exx_info.info_global.coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Erfc].resize(erfc_alpha.size());
416+
for(std::size_t i=0; i<erfc_alpha.size(); ++i)
417+
{
418+
GlobalC::exx_info.info_global.coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Erfc] = {{
419+
{"alpha", ModuleBase::GlobalFunc::TO_STRING(erfc_alpha[i])},
420+
{"omega", ModuleBase::GlobalFunc::TO_STRING(PARAM.inp.exx_erfc_omega[i])} }};
421+
}
422+
}
413423
}
414424
}
415425
#ifdef __EXX

‎source/source_pw/module_pwdft/operator_pw/op_exx_pw.cpp‎

Lines changed: 126 additions & 72 deletions
Original file line numberDiff line numberDiff line change
@@ -88,7 +88,6 @@ OperatorEXXPW<T, Device>::OperatorEXXPW(const int* isk_in,
8888
tpiba = ucell->tpiba;
8989
Real tpiba2 = tpiba * tpiba;
9090
// calculate the exx_divergence
91-
exx_divergence();
9291

9392
}
9493

@@ -617,101 +616,156 @@ void OperatorEXXPW<T, Device>::get_potential() const
617616
Real nqs_half2 = 0.5 * kv->nmp[1];
618617
Real nqs_half3 = 0.5 * kv->nmp[2];
619618

619+
setmem_real_op()(pot, 0, rhopw->npw * wfcpw->nks * wfcpw->nks);
620620
int nks = wfcpw->nks, npw = rhopw->npw;
621621
double tpiba2 = tpiba * tpiba;
622-
// calculate the pot
623-
for (int ik = 0; ik < nks; ik++)
622+
// calculate Fock pot
623+
auto param_fock = GlobalC::exx_info.info_global.coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Fock];
624+
for (auto param : param_fock)
624625
{
625-
for (int iq = 0; iq < nks; iq++)
626+
double exx_div = exx_divergence(Conv_Coulomb_Pot_K::Coulomb_Type::Fock);
627+
double alpha = std::stod(param["alpha"]);
628+
for (int ik = 0; ik < nks; ik++)
626629
{
627-
const ModuleBase::Vector3<double> k_c = wfcpw->kvec_c[ik];
628-
const ModuleBase::Vector3<double> k_d = wfcpw->kvec_d[ik];
629-
const ModuleBase::Vector3<double> q_c = wfcpw->kvec_c[iq];
630-
const ModuleBase::Vector3<double> q_d = wfcpw->kvec_d[iq];
631-
632-
#ifdef _OPENMP
633-
#pragma omp parallel for schedule(static)
634-
#endif
635-
for (int ig = 0; ig < rhopw->npw; ig++)
630+
for (int iq = 0; iq < nks; iq++)
636631
{
637-
const ModuleBase::Vector3<double> g_d = rhopw->gdirect[ig];
638-
const ModuleBase::Vector3<double> kqg_d = k_d - q_d + g_d;
639-
// For gamma_extrapolation (https://doi.org/10.1103/PhysRevB.79.205114)
640-
// 7/8 of the points in the grid are "activated" and 1/8 are disabled.
641-
// grid_factor is designed for the 7/8 of the grid to function like all of the points
642-
Real grid_factor = 1;
643-
double extrapolate_grid = 8.0/7.0;
644-
if (gamma_extrapolation)
632+
const ModuleBase::Vector3<double> k_c = wfcpw->kvec_c[ik];
633+
const ModuleBase::Vector3<double> k_d = wfcpw->kvec_d[ik];
634+
const ModuleBase::Vector3<double> q_c = wfcpw->kvec_c[iq];
635+
const ModuleBase::Vector3<double> q_d = wfcpw->kvec_d[iq];
636+
637+
#ifdef _OPENMP
638+
#pragma omp parallel for schedule(static)
639+
#endif
640+
for (int ig = 0; ig < rhopw->npw; ig++)
645641
{
646-
// if isint(kqg_d[0] * nqs_half1) && isint(kqg_d[1] * nqs_half2) && isint(kqg_d[2] * nqs_half3)
647-
auto isint = [](double x)
642+
const ModuleBase::Vector3<double> g_d = rhopw->gdirect[ig];
643+
const ModuleBase::Vector3<double> kqg_d = k_d - q_d + g_d;
644+
// For gamma_extrapolation (https://doi.org/10.1103/PhysRevB.79.205114)
645+
// 7/8 of the points in the grid are "activated" and 1/8 are disabled.
646+
// grid_factor is designed for the 7/8 of the grid to function like all of the points
647+
Real grid_factor = 1;
648+
double extrapolate_grid = 8.0 / 7.0;
649+
if (gamma_extrapolation)
648650
{
649-
double epsilon = 1e-6; // this follows the isint judgement in q-e
650-
return std::abs(x - std::round(x)) < epsilon;
651-
};
652-
if (isint(kqg_d[0] * nqs_half1) &&
653-
isint(kqg_d[1] * nqs_half2) &&
654-
isint(kqg_d[2] * nqs_half3))
651+
// if isint(kqg_d[0] * nqs_half1) && isint(kqg_d[1] * nqs_half2) && isint(kqg_d[2] * nqs_half3)
652+
auto isint = [](double x) {
653+
double epsilon = 1e-6; // this follows the isint judgement in q-e
654+
return std::abs(x - std::round(x)) < epsilon;
655+
};
656+
if (isint(kqg_d[0] * nqs_half1) && isint(kqg_d[1] * nqs_half2) && isint(kqg_d[2] * nqs_half3))
657+
{
658+
grid_factor = 0;
659+
}
660+
else
661+
{
662+
grid_factor = extrapolate_grid;
663+
}
664+
}
665+
666+
const int nk_fac = PARAM.inp.nspin == 2 ? 2 : 1;
667+
const int nk = nks / nk_fac;
668+
const int ig_kq = ik * nks * npw + iq * npw + ig;
669+
670+
Real gg = (k_c - q_c + rhopw->gcar[ig]).norm2() * tpiba2;
671+
// if (kqgcar2 > 1e-12) // vasp uses 1/40 of the smallest (k spacing)**2
672+
if (gg >= 1e-8)
655673
{
656-
grid_factor = 0;
674+
Real fac = -ModuleBase::FOUR_PI * ModuleBase::e2 / gg;
675+
pot[ig_kq] += fac * grid_factor * alpha;
657676
}
677+
// }
658678
else
659679
{
660-
grid_factor = extrapolate_grid;
680+
pot[ig_kq] += exx_div * alpha;
661681
}
682+
// assert(is_finite(density_recip[ig]));
662683
}
684+
}
685+
}
686+
}
663687

664-
const int nk_fac = PARAM.inp.nspin == 2 ? 2 : 1;
665-
const int nk = nks / nk_fac;
666-
const int ig_kq = ik * nks * npw + iq * npw + ig;
688+
// calculate erfc pot
689+
auto param_erfc = GlobalC::exx_info.info_global.coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Erfc];
690+
for (auto param : param_erfc)
691+
{
692+
double erfc_omega = std::stod(param["omega"]);
693+
double erfc_omega2 = erfc_omega * erfc_omega;
694+
double alpha = std::stod(param["alpha"]);
695+
double exx_div = exx_divergence(Conv_Coulomb_Pot_K::Coulomb_Type::Erfc, erfc_omega);
696+
for (int ik = 0; ik < nks; ik++)
697+
{
698+
for (int iq = 0; iq < nks; iq++)
699+
{
700+
const ModuleBase::Vector3<double> k_c = wfcpw->kvec_c[ik];
701+
const ModuleBase::Vector3<double> k_d = wfcpw->kvec_d[ik];
702+
const ModuleBase::Vector3<double> q_c = wfcpw->kvec_c[iq];
703+
const ModuleBase::Vector3<double> q_d = wfcpw->kvec_d[iq];
667704

668-
Real gg = (k_c - q_c + rhopw->gcar[ig]).norm2() * tpiba2;
669-
Real hse_omega2 = GlobalC::exx_info.info_global.hse_omega * GlobalC::exx_info.info_global.hse_omega;
670-
// if (kqgcar2 > 1e-12) // vasp uses 1/40 of the smallest (k spacing)**2
671-
if (gg >= 1e-8)
705+
#ifdef _OPENMP
706+
#pragma omp parallel for schedule(static)
707+
#endif
708+
for (int ig = 0; ig < rhopw->npw; ig++)
672709
{
673-
Real fac = -ModuleBase::FOUR_PI * ModuleBase::e2 / gg;
674-
// if (PARAM.inp.dft_functional == "hse")
675-
if (GlobalC::exx_info.info_global.ccp_type == Conv_Coulomb_Pot_K::Ccp_Type::Erfc)
710+
const ModuleBase::Vector3<double> g_d = rhopw->gdirect[ig];
711+
const ModuleBase::Vector3<double> kqg_d = k_d - q_d + g_d;
712+
// For gamma_extrapolation (https://doi.org/10.1103/PhysRevB.79.205114)
713+
// 7/8 of the points in the grid are "activated" and 1/8 are disabled.
714+
// grid_factor is designed for the 7/8 of the grid to function like all of the points
715+
Real grid_factor = 1;
716+
double extrapolate_grid = 8.0 / 7.0;
717+
if (gamma_extrapolation)
676718
{
677-
pot[ig_kq] = fac * (1.0 - std::exp(-gg / 4.0 / hse_omega2)) * grid_factor;
719+
// if isint(kqg_d[0] * nqs_half1) && isint(kqg_d[1] * nqs_half2) && isint(kqg_d[2] * nqs_half3)
720+
auto isint = [](double x) {
721+
double epsilon = 1e-6; // this follows the isint judgement in q-e
722+
return std::abs(x - std::round(x)) < epsilon;
723+
};
724+
if (isint(kqg_d[0] * nqs_half1) && isint(kqg_d[1] * nqs_half2) && isint(kqg_d[2] * nqs_half3))
725+
{
726+
grid_factor = 0;
727+
}
728+
else
729+
{
730+
grid_factor = extrapolate_grid;
731+
}
678732
}
679-
else if (GlobalC::exx_info.info_global.ccp_type == Conv_Coulomb_Pot_K::Ccp_Type::Erf)
680-
{
681-
pot[ig_kq] = fac * (std::exp(-gg / 4.0 / hse_omega2)) * grid_factor;
682-
}
683-
else
684-
{
685-
pot[ig_kq] = fac * grid_factor;
686-
}
687-
}
688-
// }
689-
else
690-
{
691-
// if (PARAM.inp.dft_functional == "hse")
692-
if (GlobalC::exx_info.info_global.ccp_type == Conv_Coulomb_Pot_K::Ccp_Type::Erfc &&
693-
!gamma_extrapolation)
733+
734+
const int nk_fac = PARAM.inp.nspin == 2 ? 2 : 1;
735+
const int nk = nks / nk_fac;
736+
const int ig_kq = ik * nks * npw + iq * npw + ig;
737+
738+
Real gg = (k_c - q_c + rhopw->gcar[ig]).norm2() * tpiba2;
739+
// if (kqgcar2 > 1e-12) // vasp uses 1/40 of the smallest (k spacing)**2
740+
if (gg >= 1e-8)
694741
{
695-
pot[ig_kq] = exx_div - ModuleBase::PI * ModuleBase::e2 / hse_omega2;
742+
Real fac = -ModuleBase::FOUR_PI * ModuleBase::e2 / gg;
743+
pot[ig_kq] += fac * (1.0 - std::exp(-gg / 4.0 / erfc_omega2)) * grid_factor * alpha;
696744
}
745+
// }
697746
else
698747
{
699-
pot[ig_kq] = exx_div;
748+
// if (PARAM.inp.dft_functional == "hse")
749+
if (!gamma_extrapolation)
750+
{
751+
pot[ig_kq] += (exx_div - ModuleBase::PI * ModuleBase::e2 / erfc_omega2) * alpha;
752+
}
753+
else
754+
{
755+
pot[ig_kq] += exx_div * alpha;
756+
}
700757
}
758+
// assert(is_finite(density_recip[ig]));
701759
}
702-
// assert(is_finite(density_recip[ig]));
703760
}
704761
}
705762
}
706763
}
707764

708765
template <typename T, typename Device>
709-
void OperatorEXXPW<T, Device>::exx_divergence()
766+
double OperatorEXXPW<T, Device>::exx_divergence(Conv_Coulomb_Pot_K::Coulomb_Type coulomb_type, double erfc_omega) const
710767
{
711-
if (GlobalC::exx_info.info_lip.lambda == 0.0)
712-
{
713-
return;
714-
}
768+
double exx_div = 0;
715769

716770
Real nqs_half1 = 0.5 * kv->nmp[0];
717771
Real nqs_half2 = 0.5 * kv->nmp[1];
@@ -764,9 +818,9 @@ void OperatorEXXPW<T, Device>::exx_divergence()
764818

765819
if (qq <= 1e-8) continue;
766820
// else if (PARAM.inp.dft_functional == "hse")
767-
else if (GlobalC::exx_info.info_global.ccp_type == Conv_Coulomb_Pot_K::Ccp_Type::Erfc)
821+
else if (coulomb_type == Conv_Coulomb_Pot_K::Coulomb_Type::Erfc)
768822
{
769-
double omega = GlobalC::exx_info.info_global.hse_omega;
823+
double omega = erfc_omega;
770824
double omega2 = omega * omega;
771825
div += std::exp(-alpha * qq) / qq * (1.0 - std::exp(-qq*tpiba2 / 4.0 / omega2)) * grid_factor;
772826
}
@@ -783,9 +837,9 @@ void OperatorEXXPW<T, Device>::exx_divergence()
783837
// if (PARAM.inp.dft_functional == "hse")
784838
if (!gamma_extrapolation)
785839
{
786-
if (GlobalC::exx_info.info_global.ccp_type == Conv_Coulomb_Pot_K::Ccp_Type::Erfc)
840+
if (coulomb_type == Conv_Coulomb_Pot_K::Coulomb_Type::Erfc)
787841
{
788-
double omega = GlobalC::exx_info.info_global.hse_omega;
842+
double omega = erfc_omega;
789843
div += tpiba2 / 4.0 / omega / omega; // compensate for the finite value when qq = 0
790844
}
791845
else
@@ -805,9 +859,9 @@ void OperatorEXXPW<T, Device>::exx_divergence()
805859
double dq = 5.0 / std::sqrt(alpha) / nqq;
806860
double aa = 0.0;
807861
// if (PARAM.inp.dft_functional == "hse")
808-
if (GlobalC::exx_info.info_global.ccp_type == Conv_Coulomb_Pot_K::Ccp_Type::Erfc)
862+
if (coulomb_type == Conv_Coulomb_Pot_K::Coulomb_Type::Erfc)
809863
{
810-
double omega = GlobalC::exx_info.info_global.hse_omega;
864+
double omega = erfc_omega;
811865
double omega2 = omega * omega;
812866
#ifdef _OPENMP
813867
#pragma omp parallel for reduction(+:aa)
@@ -828,7 +882,7 @@ void OperatorEXXPW<T, Device>::exx_divergence()
828882
// exx_div = 0;
829883
// std::cout << "EXX divergence: " << exx_div << std::endl;
830884

831-
return;
885+
return exx_div;
832886
}
833887

834888
template <typename T, typename Device>

‎source/source_pw/module_pwdft/operator_pw/op_exx_pw.h‎

Lines changed: 9 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -1,15 +1,16 @@
11
#ifndef OPEXXPW_H
22
#define OPEXXPW_H
33

4+
#include "operator_pw.h"
5+
#include "source_base/blas_connector.h"
6+
#include "source_base/kernels/math_kernel_op.h"
7+
#include "source_base/macros.h"
48
#include "source_base/matrix.h"
59
#include "source_basis/module_pw/pw_basis.h"
10+
#include "source_basis/module_pw/pw_basis_k.h"
611
#include "source_cell/klist.h"
12+
#include "source_lcao/module_ri/conv_coulomb_pot_k.h"
713
#include "source_psi/psi.h"
8-
#include "operator_pw.h"
9-
#include "source_basis/module_pw/pw_basis_k.h"
10-
#include "source_base/macros.h"
11-
#include "source_base/kernels/math_kernel_op.h"
12-
#include "source_base/blas_connector.h"
1314

1415
#include <memory>
1516
#include <utility>
@@ -59,15 +60,15 @@ class OperatorEXXPW : public OperatorPW<T, Device>
5960
const ModulePW::PW_Basis_K *wfcpw = nullptr;
6061
const ModulePW::PW_Basis *rhopw = nullptr;
6162
const UnitCell *ucell = nullptr;
62-
Real exx_div = 0;
63+
// Real exx_div = 0;
6364
Real tpiba = 0;
6465

6566
std::vector<int> get_q_points(const int ik) const;
6667
const T *get_pw(const int m, const int iq) const;
6768

6869
void multiply_potential(T *density_recip, int ik, int iq) const;
6970

70-
void exx_divergence();
71+
double exx_divergence(Conv_Coulomb_Pot_K::Coulomb_Type coulomb_type, double erfc_omega = 0) const;
7172

7273
void get_potential() const;
7374

@@ -135,6 +136,7 @@ class OperatorEXXPW : public OperatorPW<T, Device>
135136
base_device::AbacusDevice_t device = {};
136137

137138
using setmem_complex_op = base_device::memory::set_memory_op<T, Device>;
139+
using setmem_real_op = base_device::memory::set_memory_op<Real, Device>;
138140
using resmem_complex_op = base_device::memory::resize_memory_op<T, Device>;
139141
using delmem_complex_op = base_device::memory::delete_memory_op<T, Device>;
140142
using syncmem_complex_op = base_device::memory::synchronize_memory_op<T, Device, Device>;

0 commit comments

Comments
 (0)