diff --git a/source/module_hamilt_pw/hamilt_pwdft/operator_pw/op_exx_pw.cpp b/source/module_hamilt_pw/hamilt_pwdft/operator_pw/op_exx_pw.cpp index fc7733370c..022d439fc5 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/operator_pw/op_exx_pw.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/operator_pw/op_exx_pw.cpp @@ -80,8 +80,9 @@ OperatorEXXPW::OperatorEXXPW(const int* isk_in, // allocate h_psi recip space memory resmem_complex_op()(h_psi_recip, wfcpw->npwk_max); // resmem_complex_op()(this->ctx, psi_all_real, wfcpw->nrxx * GlobalV::NBANDS); + int nks = wfcpw->nks; -// std::cout << "nks: " << nks << std::endl; + int nk_fac = PARAM.inp.nspin == 2 ? 2 : 1; resmem_real_op()(pot, rhopw->npw * nks * nks); tpiba = ucell->tpiba; @@ -168,6 +169,12 @@ void OperatorEXXPW::act_op(const int nbands, const int ngk_ik, const bool is_first_node) const { +// std::cout << "nbands: " << nbands +// << " nbasis: " << nbasis +// << " npol: " << npol +// << " ngk_ik: " << ngk_ik +// << " is_first_node: " << is_first_node +// << std::endl; if (!potential_got) { get_potential(); @@ -199,6 +206,7 @@ void OperatorEXXPW::act_op(const int nbands, Real nqs = q_points.size(); for (int iq: q_points) { +// std::cout << "ik" << this->ik << " iq" << iq << std::endl; for (int m_iband = 0; m_iband < psi.get_nbands(); m_iband++) { // double wg_mqb_real = GlobalC::exx_helper.wg(iq, m_iband); @@ -518,7 +526,23 @@ std::vector OperatorEXXPW::get_q_points(const int ik) const { for (int iq = 0; iq < wfcpw->nks; iq++) { - q_points_ik.push_back(iq); + if (PARAM.inp.nspin ==1 ) + { + q_points_ik.push_back(iq); + } + else if (PARAM.inp.nspin == 2) + { + int nk_fac = 2; + int nk = wfcpw->nks / nk_fac; + if (iq / nk == ik / nk) + { + q_points_ik.push_back(iq); + } + } + else + { + ModuleBase::WARNING_QUIT("OperatorEXXPW", "nspin == 4 not supported"); + } } } // else @@ -539,6 +563,8 @@ void OperatorEXXPW::multiply_potential(T *density_recip, int ik, int ModuleBase::timer::tick("OperatorEXXPW", "multiply_potential"); int npw = rhopw->npw; int nks = wfcpw->nks; + int nk_fac = PARAM.inp.nspin == 2 ? 2 : 1; + int nk = nks / nk_fac; #ifdef _OPENMP #pragma omp parallel for schedule(static) @@ -635,6 +661,8 @@ void OperatorEXXPW::get_potential() const } } + 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; @@ -689,6 +717,8 @@ void OperatorEXXPW::exx_divergence() Real nqs_half2 = 0.5 * kv->nmp[1]; Real nqs_half3 = 0.5 * kv->nmp[2]; + int nk_fac = PARAM.inp.nspin == 2 ? 2 : 1; + // here we follow the exx_divergence subroutine in q-e (PW/src/exx_base.f90) double alpha = 10.0 / wfcpw->gk_ecut; double tpiba2 = tpiba * tpiba; @@ -766,6 +796,7 @@ void OperatorEXXPW::exx_divergence() } div *= ModuleBase::e2 * ModuleBase::FOUR_PI / tpiba2 / wfcpw->nks; +// std::cout << "div: " << div << std::endl; // numerically value the mean value of F(q) in the reciprocal space // This means we need to calculate the average of F(q) in the first brillouin zone @@ -793,8 +824,9 @@ void OperatorEXXPW::exx_divergence() // printf("ucell: %p\n", ucell); double omega = ucell->omega; div -= ModuleBase::e2 * omega * aa; - exx_div = div * wfcpw->nks; - // std::cout << "EXX divergence: " << exx_div << std::endl; + exx_div = div * wfcpw->nks / nk_fac; +// exx_div = 0; +// std::cout << "EXX divergence: " << exx_div << std::endl; return; } @@ -862,8 +894,7 @@ double OperatorEXXPW::cal_exx_energy_op(psi::Psi *ppsi_) c T* density_recip = new T[rhopw->npw]; if (wg == nullptr) return 0.0; - // evaluate the Eexx - // T Eexx_ik = 0.0; + const int nk_fac = PARAM.inp.nspin == 2 ? 2 : 1; double Eexx_ik_real = 0.0; for (int ik = 0; ik < wfcpw->nks; ik++) { @@ -884,8 +915,6 @@ double OperatorEXXPW::cal_exx_energy_op(psi::Psi *ppsi_) c continue; } - // std::cout << "ik = " << ik << " nb = " << n_iband << " wg_ikb = " << wg_ikb_real << std::endl; - // const T *psi_nk = get_pw(n_iband, ik); psi.fix_kb(ik, n_iband); const T* psi_nk = psi.get_pointer(); @@ -901,7 +930,6 @@ double OperatorEXXPW::cal_exx_energy_op(psi::Psi *ppsi_) c } double nqs = q_points.size(); - // std::cout << "ik = " << ik << " ib = " << n_iband << " wg_kb = " << wg_ikb_real << " wk_ik = " << kv->wk[ik] << std::endl; for (int iq: q_points) { for (int m_iband = 0; m_iband < psi.get_nbands(); m_iband++) @@ -914,8 +942,6 @@ double OperatorEXXPW::cal_exx_energy_op(psi::Psi *ppsi_) c continue; } - // std::cout << "iq = " << iq << " mb = " << m_iband << " wg_iqb = " << wg_iqb_real << std::endl; - psi_.fix_kb(iq, m_iband); const T* psi_mq = psi_.get_pointer(); // const T* psi_mq = get_pw(m_iband, iq); @@ -945,6 +971,7 @@ double OperatorEXXPW::cal_exx_energy_op(psi::Psi *ppsi_) c { int nks = wfcpw->nks; int npw = rhopw->npw; + int nk = nks / nk_fac; Real Fac = pot[ik * nks * npw + iq * npw + ig]; Eexx_ik_real += Fac * (density_recip[ig] * std::conj(density_recip[ig])).real() * wg_iqb_real / nqs * wg_ikb_real / kv->wk[ik]; diff --git a/source/module_io/input_conv.cpp b/source/module_io/input_conv.cpp index d0b5bf5512..00dc63e15a 100644 --- a/source/module_io/input_conv.cpp +++ b/source/module_io/input_conv.cpp @@ -419,9 +419,9 @@ void Input_Conv::Convert() ModuleSymmetry::Symmetry::symm_flag = -1; } - if (PARAM.inp.nspin != 1) + if (PARAM.inp.nspin != 1 && PARAM.inp.nspin != 2) { - ModuleBase::WARNING_QUIT("Input_Conv", "EXX PW works only with nspin=1"); + ModuleBase::WARNING_QUIT("Input_Conv", "EXX PW works only with nspin=1 and 2"); } if (PARAM.inp.device != "cpu")