diff --git a/source/source_io/input_conv.cpp b/source/source_io/input_conv.cpp index 2107fa3d45..8867346787 100644 --- a/source/source_io/input_conv.cpp +++ b/source/source_io/input_conv.cpp @@ -409,6 +409,16 @@ void Input_Conv::Convert() {"singularity_correction", PARAM.inp.exx_singularity_correction} }}; } } + else if(PARAM.inp.basis_type == "pw") + { + GlobalC::exx_info.info_global.coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Erfc].resize(erfc_alpha.size()); + for(std::size_t i=0; i::OperatorEXXPW(const int* isk_in, tpiba = ucell->tpiba; Real tpiba2 = tpiba * tpiba; // calculate the exx_divergence - exx_divergence(); } @@ -617,101 +616,156 @@ void OperatorEXXPW::get_potential() const Real nqs_half2 = 0.5 * kv->nmp[1]; Real nqs_half3 = 0.5 * kv->nmp[2]; + setmem_real_op()(pot, 0, rhopw->npw * wfcpw->nks * wfcpw->nks); int nks = wfcpw->nks, npw = rhopw->npw; double tpiba2 = tpiba * tpiba; - // calculate the pot - for (int ik = 0; ik < nks; ik++) + // calculate Fock pot + auto param_fock = GlobalC::exx_info.info_global.coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Fock]; + for (auto param : param_fock) { - for (int iq = 0; iq < nks; iq++) + double exx_div = exx_divergence(Conv_Coulomb_Pot_K::Coulomb_Type::Fock); + double alpha = std::stod(param["alpha"]); + for (int ik = 0; ik < nks; ik++) { - const ModuleBase::Vector3 k_c = wfcpw->kvec_c[ik]; - const ModuleBase::Vector3 k_d = wfcpw->kvec_d[ik]; - const ModuleBase::Vector3 q_c = wfcpw->kvec_c[iq]; - const ModuleBase::Vector3 q_d = wfcpw->kvec_d[iq]; - - #ifdef _OPENMP - #pragma omp parallel for schedule(static) - #endif - for (int ig = 0; ig < rhopw->npw; ig++) + for (int iq = 0; iq < nks; iq++) { - const ModuleBase::Vector3 g_d = rhopw->gdirect[ig]; - const ModuleBase::Vector3 kqg_d = k_d - q_d + g_d; - // For gamma_extrapolation (https://doi.org/10.1103/PhysRevB.79.205114) - // 7/8 of the points in the grid are "activated" and 1/8 are disabled. - // grid_factor is designed for the 7/8 of the grid to function like all of the points - Real grid_factor = 1; - double extrapolate_grid = 8.0/7.0; - if (gamma_extrapolation) + const ModuleBase::Vector3 k_c = wfcpw->kvec_c[ik]; + const ModuleBase::Vector3 k_d = wfcpw->kvec_d[ik]; + const ModuleBase::Vector3 q_c = wfcpw->kvec_c[iq]; + const ModuleBase::Vector3 q_d = wfcpw->kvec_d[iq]; + +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif + for (int ig = 0; ig < rhopw->npw; ig++) { - // if isint(kqg_d[0] * nqs_half1) && isint(kqg_d[1] * nqs_half2) && isint(kqg_d[2] * nqs_half3) - auto isint = [](double x) + const ModuleBase::Vector3 g_d = rhopw->gdirect[ig]; + const ModuleBase::Vector3 kqg_d = k_d - q_d + g_d; + // For gamma_extrapolation (https://doi.org/10.1103/PhysRevB.79.205114) + // 7/8 of the points in the grid are "activated" and 1/8 are disabled. + // grid_factor is designed for the 7/8 of the grid to function like all of the points + Real grid_factor = 1; + double extrapolate_grid = 8.0 / 7.0; + if (gamma_extrapolation) { - double epsilon = 1e-6; // this follows the isint judgement in q-e - return std::abs(x - std::round(x)) < epsilon; - }; - if (isint(kqg_d[0] * nqs_half1) && - isint(kqg_d[1] * nqs_half2) && - isint(kqg_d[2] * nqs_half3)) + // if isint(kqg_d[0] * nqs_half1) && isint(kqg_d[1] * nqs_half2) && isint(kqg_d[2] * nqs_half3) + auto isint = [](double x) { + double epsilon = 1e-6; // this follows the isint judgement in q-e + return std::abs(x - std::round(x)) < epsilon; + }; + if (isint(kqg_d[0] * nqs_half1) && isint(kqg_d[1] * nqs_half2) && isint(kqg_d[2] * nqs_half3)) + { + grid_factor = 0; + } + else + { + grid_factor = extrapolate_grid; + } + } + + const int nk_fac = PARAM.inp.nspin == 2 ? 2 : 1; + const int nk = nks / nk_fac; + const int ig_kq = ik * nks * npw + iq * npw + ig; + + Real gg = (k_c - q_c + rhopw->gcar[ig]).norm2() * tpiba2; + // if (kqgcar2 > 1e-12) // vasp uses 1/40 of the smallest (k spacing)**2 + if (gg >= 1e-8) { - grid_factor = 0; + Real fac = -ModuleBase::FOUR_PI * ModuleBase::e2 / gg; + pot[ig_kq] += fac * grid_factor * alpha; } + // } else { - grid_factor = extrapolate_grid; + pot[ig_kq] += exx_div * alpha; } + // assert(is_finite(density_recip[ig])); } + } + } + } - const int nk_fac = PARAM.inp.nspin == 2 ? 2 : 1; - const int nk = nks / nk_fac; - const int ig_kq = ik * nks * npw + iq * npw + ig; + // calculate erfc pot + auto param_erfc = GlobalC::exx_info.info_global.coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Erfc]; + for (auto param : param_erfc) + { + double erfc_omega = std::stod(param["omega"]); + double erfc_omega2 = erfc_omega * erfc_omega; + double alpha = std::stod(param["alpha"]); + double exx_div = exx_divergence(Conv_Coulomb_Pot_K::Coulomb_Type::Erfc, erfc_omega); + for (int ik = 0; ik < nks; ik++) + { + for (int iq = 0; iq < nks; iq++) + { + const ModuleBase::Vector3 k_c = wfcpw->kvec_c[ik]; + const ModuleBase::Vector3 k_d = wfcpw->kvec_d[ik]; + const ModuleBase::Vector3 q_c = wfcpw->kvec_c[iq]; + const ModuleBase::Vector3 q_d = wfcpw->kvec_d[iq]; - Real gg = (k_c - q_c + rhopw->gcar[ig]).norm2() * tpiba2; - Real hse_omega2 = GlobalC::exx_info.info_global.hse_omega * GlobalC::exx_info.info_global.hse_omega; - // if (kqgcar2 > 1e-12) // vasp uses 1/40 of the smallest (k spacing)**2 - if (gg >= 1e-8) +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif + for (int ig = 0; ig < rhopw->npw; ig++) { - Real fac = -ModuleBase::FOUR_PI * ModuleBase::e2 / gg; - // if (PARAM.inp.dft_functional == "hse") - if (GlobalC::exx_info.info_global.ccp_type == Conv_Coulomb_Pot_K::Ccp_Type::Erfc) + const ModuleBase::Vector3 g_d = rhopw->gdirect[ig]; + const ModuleBase::Vector3 kqg_d = k_d - q_d + g_d; + // For gamma_extrapolation (https://doi.org/10.1103/PhysRevB.79.205114) + // 7/8 of the points in the grid are "activated" and 1/8 are disabled. + // grid_factor is designed for the 7/8 of the grid to function like all of the points + Real grid_factor = 1; + double extrapolate_grid = 8.0 / 7.0; + if (gamma_extrapolation) { - pot[ig_kq] = fac * (1.0 - std::exp(-gg / 4.0 / hse_omega2)) * grid_factor; + // if isint(kqg_d[0] * nqs_half1) && isint(kqg_d[1] * nqs_half2) && isint(kqg_d[2] * nqs_half3) + auto isint = [](double x) { + double epsilon = 1e-6; // this follows the isint judgement in q-e + return std::abs(x - std::round(x)) < epsilon; + }; + if (isint(kqg_d[0] * nqs_half1) && isint(kqg_d[1] * nqs_half2) && isint(kqg_d[2] * nqs_half3)) + { + grid_factor = 0; + } + else + { + grid_factor = extrapolate_grid; + } } - else if (GlobalC::exx_info.info_global.ccp_type == Conv_Coulomb_Pot_K::Ccp_Type::Erf) - { - pot[ig_kq] = fac * (std::exp(-gg / 4.0 / hse_omega2)) * grid_factor; - } - else - { - pot[ig_kq] = fac * grid_factor; - } - } - // } - else - { - // if (PARAM.inp.dft_functional == "hse") - if (GlobalC::exx_info.info_global.ccp_type == Conv_Coulomb_Pot_K::Ccp_Type::Erfc && - !gamma_extrapolation) + + const int nk_fac = PARAM.inp.nspin == 2 ? 2 : 1; + const int nk = nks / nk_fac; + const int ig_kq = ik * nks * npw + iq * npw + ig; + + Real gg = (k_c - q_c + rhopw->gcar[ig]).norm2() * tpiba2; + // if (kqgcar2 > 1e-12) // vasp uses 1/40 of the smallest (k spacing)**2 + if (gg >= 1e-8) { - pot[ig_kq] = exx_div - ModuleBase::PI * ModuleBase::e2 / hse_omega2; + Real fac = -ModuleBase::FOUR_PI * ModuleBase::e2 / gg; + pot[ig_kq] += fac * (1.0 - std::exp(-gg / 4.0 / erfc_omega2)) * grid_factor * alpha; } + // } else { - pot[ig_kq] = exx_div; + // if (PARAM.inp.dft_functional == "hse") + if (!gamma_extrapolation) + { + pot[ig_kq] += (exx_div - ModuleBase::PI * ModuleBase::e2 / erfc_omega2) * alpha; + } + else + { + pot[ig_kq] += exx_div * alpha; + } } + // assert(is_finite(density_recip[ig])); } - // assert(is_finite(density_recip[ig])); } } } } template -void OperatorEXXPW::exx_divergence() +double OperatorEXXPW::exx_divergence(Conv_Coulomb_Pot_K::Coulomb_Type coulomb_type, double erfc_omega) const { - if (GlobalC::exx_info.info_lip.lambda == 0.0) - { - return; - } + double exx_div = 0; Real nqs_half1 = 0.5 * kv->nmp[0]; Real nqs_half2 = 0.5 * kv->nmp[1]; @@ -764,9 +818,9 @@ void OperatorEXXPW::exx_divergence() if (qq <= 1e-8) continue; // else if (PARAM.inp.dft_functional == "hse") - else if (GlobalC::exx_info.info_global.ccp_type == Conv_Coulomb_Pot_K::Ccp_Type::Erfc) + else if (coulomb_type == Conv_Coulomb_Pot_K::Coulomb_Type::Erfc) { - double omega = GlobalC::exx_info.info_global.hse_omega; + double omega = erfc_omega; double omega2 = omega * omega; div += std::exp(-alpha * qq) / qq * (1.0 - std::exp(-qq*tpiba2 / 4.0 / omega2)) * grid_factor; } @@ -783,9 +837,9 @@ void OperatorEXXPW::exx_divergence() // if (PARAM.inp.dft_functional == "hse") if (!gamma_extrapolation) { - if (GlobalC::exx_info.info_global.ccp_type == Conv_Coulomb_Pot_K::Ccp_Type::Erfc) + if (coulomb_type == Conv_Coulomb_Pot_K::Coulomb_Type::Erfc) { - double omega = GlobalC::exx_info.info_global.hse_omega; + double omega = erfc_omega; div += tpiba2 / 4.0 / omega / omega; // compensate for the finite value when qq = 0 } else @@ -805,9 +859,9 @@ void OperatorEXXPW::exx_divergence() double dq = 5.0 / std::sqrt(alpha) / nqq; double aa = 0.0; // if (PARAM.inp.dft_functional == "hse") - if (GlobalC::exx_info.info_global.ccp_type == Conv_Coulomb_Pot_K::Ccp_Type::Erfc) + if (coulomb_type == Conv_Coulomb_Pot_K::Coulomb_Type::Erfc) { - double omega = GlobalC::exx_info.info_global.hse_omega; + double omega = erfc_omega; double omega2 = omega * omega; #ifdef _OPENMP #pragma omp parallel for reduction(+:aa) @@ -828,7 +882,7 @@ void OperatorEXXPW::exx_divergence() // exx_div = 0; // std::cout << "EXX divergence: " << exx_div << std::endl; - return; + return exx_div; } template diff --git a/source/source_pw/module_pwdft/operator_pw/op_exx_pw.h b/source/source_pw/module_pwdft/operator_pw/op_exx_pw.h index 76f5cdf3dc..b5da2f8355 100644 --- a/source/source_pw/module_pwdft/operator_pw/op_exx_pw.h +++ b/source/source_pw/module_pwdft/operator_pw/op_exx_pw.h @@ -1,15 +1,16 @@ #ifndef OPEXXPW_H #define OPEXXPW_H +#include "operator_pw.h" +#include "source_base/blas_connector.h" +#include "source_base/kernels/math_kernel_op.h" +#include "source_base/macros.h" #include "source_base/matrix.h" #include "source_basis/module_pw/pw_basis.h" +#include "source_basis/module_pw/pw_basis_k.h" #include "source_cell/klist.h" +#include "source_lcao/module_ri/conv_coulomb_pot_k.h" #include "source_psi/psi.h" -#include "operator_pw.h" -#include "source_basis/module_pw/pw_basis_k.h" -#include "source_base/macros.h" -#include "source_base/kernels/math_kernel_op.h" -#include "source_base/blas_connector.h" #include #include @@ -59,7 +60,7 @@ class OperatorEXXPW : public OperatorPW const ModulePW::PW_Basis_K *wfcpw = nullptr; const ModulePW::PW_Basis *rhopw = nullptr; const UnitCell *ucell = nullptr; - Real exx_div = 0; +// Real exx_div = 0; Real tpiba = 0; std::vector get_q_points(const int ik) const; @@ -67,7 +68,7 @@ class OperatorEXXPW : public OperatorPW void multiply_potential(T *density_recip, int ik, int iq) const; - void exx_divergence(); + double exx_divergence(Conv_Coulomb_Pot_K::Coulomb_Type coulomb_type, double erfc_omega = 0) const; void get_potential() const; @@ -135,6 +136,7 @@ class OperatorEXXPW : public OperatorPW base_device::AbacusDevice_t device = {}; using setmem_complex_op = base_device::memory::set_memory_op; + using setmem_real_op = base_device::memory::set_memory_op; using resmem_complex_op = base_device::memory::resize_memory_op; using delmem_complex_op = base_device::memory::delete_memory_op; using syncmem_complex_op = base_device::memory::synchronize_memory_op;