From 071ae7ad3725a54c30668b7c7ef960b3eb72d775 Mon Sep 17 00:00:00 2001 From: linpz Date: Wed, 13 Nov 2024 15:51:01 +0800 Subject: [PATCH] update code format of Exx_LIP --- source/Makefile.Objects | 1 - source/module_ri/CMakeLists.txt | 1 - source/module_ri/exx_lip.cpp | 716 -------------------------------- source/module_ri/exx_lip.h | 62 +-- source/module_ri/exx_lip.hpp | 611 +++++++++++++++++++++++++++ 5 files changed, 646 insertions(+), 745 deletions(-) delete mode 100644 source/module_ri/exx_lip.cpp create mode 100644 source/module_ri/exx_lip.hpp diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 3b95d05b35..f9cf8568c6 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -606,7 +606,6 @@ OBJS_MODULE_RI=conv_coulomb_pot_k.o\ ABFs_Construct-PCA.o \ exx_opt_orb.o \ exx_opt_orb-print.o \ - exx_lip.o\ Matrix_Orbs11.o\ Matrix_Orbs21.o\ Matrix_Orbs22.o\ diff --git a/source/module_ri/CMakeLists.txt b/source/module_ri/CMakeLists.txt index 10dc30d9a0..a201cab1a3 100644 --- a/source/module_ri/CMakeLists.txt +++ b/source/module_ri/CMakeLists.txt @@ -7,7 +7,6 @@ if (ENABLE_LIBRI) Matrix_Orbs21.cpp Matrix_Orbs22.cpp RI_2D_Comm.cpp - exx_lip.cpp Mix_DMk_2D.cpp Mix_Matrix.cpp ) diff --git a/source/module_ri/exx_lip.cpp b/source/module_ri/exx_lip.cpp deleted file mode 100644 index b62792310c..0000000000 --- a/source/module_ri/exx_lip.cpp +++ /dev/null @@ -1,716 +0,0 @@ -//========================================================== -// AUTHOR : Peize Lin - -// DATE : 2015-03-10 -//========================================================== - -#include "exx_lip.h" -#include "module_base/global_function.h" -#include "module_base/vector3.h" -#include "module_hamilt_pw/hamilt_pwdft/global.h" -#include "module_cell/klist.h" -#include "module_hamilt_pw/hamilt_pwdft/wavefunc.h" -#include "module_hamilt_lcao/hamilt_lcaodft/wavefunc_in_pw.h" -#include "module_base/lapack_connector.h" -#include -#include "module_base/parallel_global.h" -#include "module_parameter/parameter.h" -template -void Exx_Lip::cal_exx() -{ - ModuleBase::TITLE("Exx_Lip", "cal_exx"); - - wf_wg_cal(); - //cout_t("wf_wg_cal",my_time(t)); - psi_cal(); - //cout_t("psi_cal",my_time(t)); - for (int ik = 0; ik < k_pack->kv_ptr->get_nks(); ++ik) - { - phi_cal(k_pack, ik); - //t_phi_cal += my_time(t); - - judge_singularity(ik); - for (int iw_l = 0; iw_l < PARAM.globalv.nlocal; ++iw_l) { - for (int iw_r = 0; iw_r < PARAM.globalv.nlocal; ++iw_r) { - sum1[iw_l * PARAM.globalv.nlocal + iw_r] = T(0.0, 0.0); -} -} - if (Conv_Coulomb_Pot_K::Ccp_Type::Ccp == info.ccp_type || Conv_Coulomb_Pot_K::Ccp_Type::Hf == info.ccp_type) - { - sum2_factor = 0.0; - if (gzero_rank_in_pool == GlobalV::RANK_IN_POOL) { - for (int iw_l = 0; iw_l < PARAM.globalv.nlocal; ++iw_l) { - for (int iw_r = 0; iw_r < PARAM.globalv.nlocal; ++iw_r) { - sum3[iw_l][iw_r] = T(0.0, 0.0); -} -} -} - } - - for (int iq_tmp = iq_vecik; iq_tmp < iq_vecik + q_pack->kv_ptr->get_nks() / PARAM.inp.nspin; ++iq_tmp) // !!! k_point - // parallel incompleted.need to loop iq in other pool - { - int iq = (ik < (k_pack->kv_ptr->get_nks() / PARAM.inp.nspin)) ? (iq_tmp % (q_pack->kv_ptr->get_nks() / PARAM.inp.nspin)) : - (iq_tmp % (q_pack->kv_ptr->get_nks() / PARAM.inp.nspin) + (q_pack->kv_ptr->get_nks() / PARAM.inp.nspin)); - qkg2_exp(ik, iq); - //t_qkg2_exp += my_time(t); - for (int ib = 0; ib < PARAM.inp.nbands; ++ib) - { - b_cal(ik, iq, ib); - //t_b_cal += my_time(t); - if (Conv_Coulomb_Pot_K::Ccp_Type::Ccp == info.ccp_type || Conv_Coulomb_Pot_K::Ccp_Type::Hf == info.ccp_type) { - if (iq == iq_vecik) { - sum3_cal(iq, ib); -} -} - //t_sum3_cal += my_time(t); - b_sum(iq, ib); - //t_b_sum += my_time(t); - } - } - sum_all(ik); - //t_sum_all += my_time(t); - } - exx_energy_cal(); - - auto print_Hexxk = [&]() - { - static int istep = 1; - for (int ik = 0; ik != q_pack->kv_ptr->get_nks(); ++ik) - { - std::ofstream - ofs("Hexxk_" + ModuleBase::GlobalFunc::TO_STRING(istep++) + "_" + ModuleBase::GlobalFunc::TO_STRING(ik) + "_" + ModuleBase::GlobalFunc::TO_STRING(GlobalV::MY_RANK)); - for (int i = 0; i != PARAM.globalv.nlocal; ++i) - { - for (int j = 0; j != PARAM.globalv.nlocal; ++j) { - ofs << exx_matrix[ik][i][j] << "\t"; -} - ofs << std::endl; - } - }; - }; -} - -/* -void Exx_Lip::cal_exx() -{ - ModuleBase::TITLE("Exx_Lip","cal_exx"); - wf_wg_cal(); - psi_cal(); - for( int ik=0; ikkv_ptr->get_nks(); ++ik) - { - phi_cal(k_pack, ik); - - judge_singularity(ik); - for( int iw_l=0; iw_l(0.0,0.0); - if( Exx_Info::Hybrid_Type::HF==info.hybrid_type || Exx_Info::Hybrid_Type::PBE0==info.hybrid_type || Exx_Info::Hybrid_Type::SCAN0==info.hybrid_type ) - { - sum2_factor = 0.0; - if(gzero_rank_in_pool==GlobalV::RANK_IN_POOL) - for( int iw_l=0; iw_l(0.0,0.0); - } - - for( int iq_tmp=iq_vecik; iq_tmpkv_ptr->get_nks()/PARAM.inp.nspin; ++iq_tmp) // !!! k_point parallel incompleted. need to loop iq in other pool - { - int iq = (ik<(k_pack->kv_ptr->get_nks()/PARAM.inp.nspin)) ? (iq_tmp%(q_pack->kv_ptr->get_nks()/PARAM.inp.nspin)) : (iq_tmp%(q_pack->kv_ptr->get_nks()/PARAM.inp.nspin)+(q_pack->kv_ptr->get_nks()/PARAM.inp.nspin)); - qkg2_exp(ik, iq); - for( int ib=0; ib -Exx_Lip::Exx_Lip( - const Exx_Info::Exx_Info_Lip& info_in, - const ModuleSymmetry::Symmetry& symm, - K_Vectors* kv_ptr_in, - psi::WFInit* wf_ptr_in, - psi::Psi* kspw_psi_ptr_in, - // wavefunc* wf_ptr_in, - const ModulePW::PW_Basis_K* wfc_basis_in, - const ModulePW::PW_Basis* rho_basis_in, - const Structure_Factor& sf, - const UnitCell* ucell_ptr_in, - const elecstate::ElecState* pelec_in) : info(info_in) -{ - ModuleBase::TITLE("Exx_Lip", "init"); - try - { - k_pack = new k_package; - k_pack->kv_ptr = kv_ptr_in; - k_pack->wf_ptr = wf_ptr_in; - k_pack->pelec = pelec_in; - k_pack->kspw_psi_ptr = kspw_psi_ptr_in; - wfc_basis = wfc_basis_in; - rho_basis = rho_basis_in; - ucell_ptr = ucell_ptr_in; - - int gzero_judge(-1); - if (rho_basis->gg_uniq[0] < 1e-8) - { - gzero_judge = GlobalV::RANK_IN_POOL; - } -#ifdef __MPI - MPI_Allreduce(&gzero_judge, &gzero_rank_in_pool, 1, MPI_INT, MPI_MAX, POOL_WORLD); - #endif - k_pack->wf_wg.create(k_pack->kv_ptr->get_nks(),PARAM.inp.nbands); - - k_pack->hvec_array = new psi::Psi(k_pack->kv_ptr->get_nks(), PARAM.inp.nbands, PARAM.globalv.nlocal); - // k_pack->hvec_array = new ModuleBase::ComplexMatrix[k_pack->kv_ptr->get_nks()]; - // for( int ik=0; ikkv_ptr->get_nks(); ++ik) - // { - // k_pack->hvec_array[ik].create(PARAM.globalv.nlocal,PARAM.inp.nbands); - // } - - // if (PARAM.inp.init_chg=="atomic") - { - q_pack = k_pack; - } - // else if(PARAM.inp.init_chg=="file") - // { - // read_q_pack(symm, wfc_basis, sf); - // } - - phi.resize(PARAM.globalv.nlocal); - for (int iw = 0; iw < PARAM.globalv.nlocal; ++iw) { phi[iw].resize(rho_basis->nrxx); } - - psi.resize(q_pack->kv_ptr->get_nks()); - for (int iq = 0; iq < q_pack->kv_ptr->get_nks(); ++iq) - { - psi[iq].resize(PARAM.inp.nbands); - for (int ib = 0; ib < PARAM.inp.nbands; ++ib) { psi[iq][ib].resize(rho_basis->nrxx); } - } - - recip_qkg2.resize(rho_basis->npw); - - b.resize(PARAM.globalv.nlocal * rho_basis->npw); - - sum1.resize(PARAM.globalv.nlocal * PARAM.globalv.nlocal); - - if (Conv_Coulomb_Pot_K::Ccp_Type::Ccp == info.ccp_type || Conv_Coulomb_Pot_K::Ccp_Type::Hf == info.ccp_type) { - if (gzero_rank_in_pool == GlobalV::RANK_IN_POOL) - { - b0.resize(PARAM.globalv.nlocal); - sum3.resize(PARAM.globalv.nlocal); - for (int iw_l = 0; iw_l < PARAM.globalv.nlocal; ++iw_l) { sum3[iw_l].resize(PARAM.globalv.nlocal); } - } -} - - exx_matrix.resize(k_pack->kv_ptr->get_nks()); - for (int ik = 0; ik < k_pack->kv_ptr->get_nks(); ++ik) - { - exx_matrix[ik].resize(PARAM.globalv.nlocal); - for (int iw_l = 0; iw_l < PARAM.globalv.nlocal; ++iw_l) { exx_matrix[ik][iw_l].resize(PARAM.globalv.nlocal); } - } - } - catch (const std::bad_alloc& ex) - { - ModuleBase::WARNING_QUIT("exact_exchange", "Memory"); - } -} - -template -Exx_Lip::~Exx_Lip() -{ - if (k_pack) {delete k_pack->hvec_array; -} - delete k_pack; - - if (PARAM.inp.init_chg == "atomic") - { - q_pack = nullptr; - } - else if (PARAM.inp.init_chg == "file") - { - delete q_pack->kv_ptr; q_pack->kv_ptr = nullptr; - delete q_pack->wf_ptr; q_pack->wf_ptr = nullptr; - // delete[] q_pack->hvec_array; q_pack->hvec_array=NULL; - delete q_pack; q_pack = nullptr; - } -} - -template -void Exx_Lip::wf_wg_cal() -{ - ModuleBase::TITLE("Exx_Lip", "wf_wg_cal"); - if (PARAM.inp.nspin == 1) { - for (int ik = 0; ik < k_pack->kv_ptr->get_nks(); ++ik) { - for (int ib = 0; ib < PARAM.inp.nbands; ++ib) { - k_pack->wf_wg(ik, ib) = k_pack->pelec->wg(ik, ib) / 2; -} -} - } else if (PARAM.inp.nspin == 2) { - for (int ik = 0; ik < k_pack->kv_ptr->get_nks(); ++ik) { - for (int ib = 0; ib < PARAM.inp.nbands; ++ib) { - k_pack->wf_wg(ik, ib) = k_pack->pelec->wg(ik, ib); -} -} -} -} - -template -void Exx_Lip::phi_cal(k_package* kq_pack, int ikq) -{ - T* porter = new T[wfc_basis->nrxx]; - for (int iw = 0; iw < PARAM.globalv.nlocal; ++iw) - { - // wfc_basis->recip2real(&kq_pack->wf_ptr->wanf2[ikq](iw,0), porter, ikq); - wfc_basis->recip2real(&(kq_pack->wf_ptr->get_psig().lock()->operator()(ikq, iw, 0)), porter, ikq); - int ir(0); - for (int ix = 0; ix < rho_basis->nx; ++ix) - { - const Real phase_x = kq_pack->kv_ptr->kvec_d[ikq].x * ix / rho_basis->nx; - for (int iy = 0; iy < rho_basis->ny; ++iy) - { - const Real phase_xy = phase_x + kq_pack->kv_ptr->kvec_d[ikq].y * iy / rho_basis->ny; - for (int iz = rho_basis->startz_current; iz < rho_basis->startz_current + rho_basis->nplane; ++iz) - { - const Real phase_xyz = phase_xy + kq_pack->kv_ptr->kvec_d[ikq].z * iz / rho_basis->nz; - const T exp_tmp = exp(phase_xyz * two_pi_i); - phi[iw][ir] = porter[ir] * exp_tmp; - ++ir; - } - } - } - } - delete[] porter; -} -template -void Exx_Lip::psi_cal() -{ - ModuleBase::TITLE("Exx_Lip", "psi_cal"); - if (PARAM.inp.init_chg == "atomic") - { - T* porter = new T[wfc_basis->nrxx]; - for (int iq = 0; iq < q_pack->kv_ptr->get_nks(); ++iq) - { - for (int ib = 0; ib < PARAM.inp.nbands; ++ib) - { - wfc_basis->recip2real(&(q_pack->kspw_psi_ptr->operator()(iq, ib, 0)), porter, iq); - - int ir(0); - for (int ix = 0; ix < rho_basis->nx; ++ix) - { - const Real phase_x = q_pack->kv_ptr->kvec_d[iq].x * ix / rho_basis->nx; - for (int iy = 0; iy < rho_basis->ny; ++iy) - { - const Real phase_xy = phase_x + q_pack->kv_ptr->kvec_d[iq].y * iy / rho_basis->ny; - for (int iz = rho_basis->startz_current; iz < rho_basis->startz_current + rho_basis->nplane; ++iz) - { - const Real phase_xyz = phase_xy + q_pack->kv_ptr->kvec_d[iq].z * iz / rho_basis->nz; - const T exp_tmp = exp(phase_xyz * two_pi_i); - psi[iq][ib][ir] = porter[ir] * exp_tmp; - ++ir; - } - } - } - } - } - delete[] porter; - } - else if (PARAM.inp.init_chg == "file") - { - for (int iq = 0; iq < q_pack->kv_ptr->get_nks(); ++iq) - { - phi_cal(q_pack, iq); - for (int ib = 0; ib < PARAM.inp.nbands; ++ib) - { - ModuleBase::GlobalFunc::ZEROS(psi[iq][ib].data(), rho_basis->nrxx); - for (int iw = 0; iw < PARAM.globalv.nlocal; ++iw) - { - for (int ir = 0; ir < rho_basis->nrxx; ++ir) - { - psi[iq][ib][ir] += (*q_pack->hvec_array)(iq, ib, iw) * phi[iw][ir]; - } - } - } - } - } - /////////////////////////////////////////////////////////////////////////// - // !!! k_point parallel incompleted. need to loop iq in other pool - /////////////////////////////////////////////////////////////////////////// -} - -template -void Exx_Lip::judge_singularity(int ik) -{ - if (PARAM.inp.init_chg=="atomic") - { - iq_vecik = ik; - } - else if(PARAM.inp.init_chg=="file") - { - Real min_q_minus_k(std::numeric_limits::max()); - for( int iq=0; iqkv_ptr->get_nks(); ++iq) - { - const Real q_minus_k((q_pack->kv_ptr->kvec_c[iq] - k_pack->kv_ptr->kvec_c[ik]).norm2()); - if(q_minus_k < min_q_minus_k) - { - min_q_minus_k = q_minus_k; - iq_vecik = iq; - } - } - } -} - -template -void Exx_Lip::qkg2_exp(int ik, int iq) -{ - for( int ig=0; ignpw; ++ig) - { - const Real qkg2 = ((q_pack->kv_ptr->kvec_c[iq] - k_pack->kv_ptr->kvec_c[ik] + rho_basis->gcar[ig]) * (ModuleBase::TWO_PI / ucell_ptr->lat0)).norm2(); - if (Conv_Coulomb_Pot_K::Ccp_Type::Ccp == info.ccp_type || Conv_Coulomb_Pot_K::Ccp_Type::Hf == info.ccp_type) - { - if (std::abs(qkg2) < 1e-10) { - recip_qkg2[ig] = 0.0; // 0 to ignore bb/qkg2 when qkg2==0 - } else { - recip_qkg2[ig] = 1.0 / qkg2; -} - sum2_factor += recip_qkg2[ig] * exp(-info.lambda * qkg2); - recip_qkg2[ig] = sqrt(recip_qkg2[ig]); - } - else if (Conv_Coulomb_Pot_K::Ccp_Type::Hse == info.ccp_type) - { - if (std::abs(qkg2) < 1e-10) { - recip_qkg2[ig] = 1.0 / (2 * info.hse_omega); - } else { - recip_qkg2[ig] = sqrt((1 - exp(-qkg2 / (4 * info.hse_omega * info.hse_omega))) / qkg2); -} - } - } -} - -template -void Exx_Lip::b_cal(int ik, int iq, int ib) -{ - const ModuleBase::Vector3 q_minus_k = q_pack->kv_ptr->kvec_d[iq] - k_pack->kv_ptr->kvec_d[ik]; - std::vector mul_tmp(rho_basis->nrxx); - for( size_t ir=0,ix=0; ixnx; ++ix) - { - const Real phase_x = q_minus_k.x * ix / rho_basis->nx; - for( size_t iy=0; iyny; ++iy) - { - const Real phase_xy = phase_x + q_minus_k.y * iy / rho_basis->ny; - for( size_t iz=rho_basis->startz_current; izstartz_current+rho_basis->nplane; ++iz) - { - const Real phase_xyz = phase_xy + q_minus_k.z * iz / rho_basis->nz; - mul_tmp[ir] = exp(-phase_xyz * two_pi_i); - mul_tmp[ir] *= psi[iq][ib][ir]; - ++ir; - } - } - } - - T* const porter = new T[rho_basis->nrxx]; - - for(size_t iw=0; iw< PARAM.globalv.nlocal; ++iw) - { - auto& phi_w = phi[iw]; - for( size_t ir=0; irnrxx; ++ir) - { - porter[ir] = conj(phi_w[ir]) * mul_tmp[ir] ; -// porter[ir] = phi_w[ir] * psi_q_b[ir] *exp_tmp[ir] ; - } - T* const b_w = &b[iw * rho_basis->npw]; - rho_basis->real2recip( porter, b_w); - if (Conv_Coulomb_Pot_K::Ccp_Type::Ccp == info.ccp_type || Conv_Coulomb_Pot_K::Ccp_Type::Hf == info.ccp_type) { - if ((iq == iq_vecik) && (gzero_rank_in_pool == GlobalV::RANK_IN_POOL)) { /// need to check while use k_point parallel - b0[iw] = b_w[rho_basis->ig_gge0]; -} -} - - for (size_t ig = 0; ig < rho_basis->npw; ++ig) { - b_w[ig] *= recip_qkg2[ig]; -} - } - delete [] porter; -} - -template -void Exx_Lip::sum3_cal(int iq, int ib) -{ - if (gzero_rank_in_pool == GlobalV::RANK_IN_POOL) { - for (int iw_l = 0; iw_l < PARAM.globalv.nlocal; ++iw_l) { - for (int iw_r = 0; iw_r < PARAM.globalv.nlocal; ++iw_r) { - sum3[iw_l][iw_r] += b0[iw_l] * conj(b0[iw_r]) * (Real)q_pack->wf_wg(iq, ib); -} -} -} -} - -template -void Exx_Lip::b_sum(int iq, int ib) // Peize Lin change 2019-04-14 -{ - // sum1[iw_l,iw_r] += \sum_{ig} b[iw_l,ig] * conj(b[iw_r,ig]) * q_pack->wf_wg(iq,ib) - LapackConnector::herk( - 'U','N', - PARAM.globalv.nlocal, rho_basis->npw, - (Real)q_pack->wf_wg(iq, ib), b.data(), rho_basis->npw, - 1.0, sum1.data(), PARAM.globalv.nlocal); -// cblas_zherk( CblasRowMajor, CblasUpper, CblasNoTrans, -// PARAM.globalv.nlocal, rho_basis->npw, -// q_pack->wf_wg(iq,ib), static_cast(b), rho_basis->npw, -// 1.0, static_cast(sum1), PARAM.globalv.nlocal); -} - -template -void Exx_Lip::sum_all(int ik) -{ - Real sum2_factor_g(0.0); - Real fourpi_div_omega = 4 * (Real)(ModuleBase::PI / ucell_ptr->omega); - Real spin_fac = 2.0; -#ifdef __MPI - if (Conv_Coulomb_Pot_K::Ccp_Type::Ccp == info.ccp_type || Conv_Coulomb_Pot_K::Ccp_Type::Hf == info.ccp_type) { - MPI_Reduce(&sum2_factor, &sum2_factor_g, 1, MPI_DOUBLE, MPI_SUM, gzero_rank_in_pool, POOL_WORLD); -} -#endif - for (size_t iw_l = 1; iw_l < PARAM.globalv.nlocal; ++iw_l) { - for (size_t iw_r = 0; iw_r < iw_l; ++iw_r) { - sum1[iw_l * PARAM.globalv.nlocal + iw_r] = conj(sum1[iw_r * PARAM.globalv.nlocal + iw_l]); // Peize Lin add conj 2019-04-14 -} -} - - for (int iw_l = 0; iw_l < PARAM.globalv.nlocal; ++iw_l) - { - for( int iw_r=0; iw_rkv_ptr->get_nks() / PARAM.inp.nspin) * sum3[iw_l][iw_r]); - } -} - } - } -} - -template -void Exx_Lip::exx_energy_cal() -{ - ModuleBase::TITLE("Exx_Lip","exx_energy_cal"); - - Real exx_energy_tmp = 0.0; - - for( int ik=0; ikkv_ptr->get_nks(); ++ik) - { - for( int iw_l=0; iw_lhvec_array)(ik, ib, iw_l)) * (*k_pack->hvec_array)(ik, ib, iw_r)).real() * k_pack->wf_wg(ik, ib); - } - } - } - } -#ifdef __MPI - MPI_Allreduce( &exx_energy_tmp, &exx_energy, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); // !!! k_point parallel incompleted. different pools have different kv.set_nks(>) deadlock -#endif - exx_energy *= (PARAM.inp.nspin==1) ? 2 : 1; - exx_energy /= 2; // ETOT = E_band - 1/2 E_exx - - #if TEST_EXX==1 - { - std::ofstream ofs("exx_matrix.dat",std::ofstream::app); - static int istep=0; - ofs<<"istep:\t"<kv_ptr->get_nks(); ++ik) - { - ofs<<"ik:\t"<kv_ptr->get_nks(); ++ik) - { - ofs<<"ik:\t"<hvec_array)(ik, ib, iw_l)) * (*k_pack->hvec_array)(ik, ib, iw_r) * k_pack->wf_wg(ik, ib); - ofs< -void Exx_Lip::write_q_pack() const -{ - if (PARAM.inp.out_chg[0] == 0) { - return; -} - - if (!GlobalV::RANK_IN_POOL) - { - const std::string exx_q_pack = "exx_q_pack/"; - int return_value=0; - const std::string command_mkdir = "test -d " + PARAM.globalv.global_out_dir + exx_q_pack + " || mkdir " + PARAM.globalv.global_out_dir + exx_q_pack; - return_value = system(command_mkdir.c_str()); - assert(return_value == 0); - - const std::string command_kpoint = "test -f " + PARAM.globalv.global_out_dir + exx_q_pack + PARAM.inp.kpoint_file + " || cp " + PARAM.inp.kpoint_file + " " + PARAM.globalv.global_out_dir + exx_q_pack + PARAM.inp.kpoint_file; - return_value = system(command_kpoint.c_str()); - assert(return_value==0); - - std::stringstream ss_wf_wg; - ss_wf_wg << PARAM.globalv.global_out_dir << exx_q_pack << "wf_wg_" << GlobalV::MY_POOL; - std::ofstream ofs_wf_wg(ss_wf_wg.str().c_str()); - for( int iq = 0; iq < q_pack->kv_ptr->get_nks(); ++iq) - { - for( int ib=0; ibwf_wg(iq,ib)<<"\t"; - } - ofs_wf_wg<kv_ptr->get_nks(); ++iq) - { - for( int iw=0; iwhvec_array)(iq, ib, iw).real() << " " << (*q_pack->hvec_array)(iq, ib, iw).imag() << " "; - } - ofs_hvec<kv_ptr = new K_Vectors(); -// const std::string exx_kpoint_card = PARAM.globalv.global_out_dir + exx_q_pack + PARAM.inp.kpoint_file; -// q_pack->kv_ptr->set( symm, exx_kpoint_card, PARAM.inp.nspin, ucell_ptr->G, ucell_ptr->latvec, GlobalV::ofs_running ); - -// q_pack->wf_ptr = new wavefunc(); -// q_pack->wf_ptr->allocate(q_pack->kv_ptr->get_nkstot(), -// q_pack->kv_ptr->get_nks(), -// q_pack->kv_ptr->ngk.data(), -// wfc_basis->npwk_max); // mohan update 2021-02-25 -// // q_pack->wf_ptr->init(q_pack->kv_ptr->get_nks(),q_pack->kv_ptr,ucell_ptr,old_pwptr,&ppcell,&GlobalC::ORB,&hm,&Pkpoints); -// q_pack->wf_ptr->table_local.create(GlobalC::ucell.ntype, GlobalC::ucell.nmax_total, PARAM.globalv.nqx); -// // q_pack->wf_ptr->table_local.create(q_pack->wf_ptr->ucell_ptr->ntype, q_pack->wf_ptr->ucell_ptr->nmax_total, PARAM.globalv.nqx); -// #ifdef __LCAO -// Wavefunc_in_pw::make_table_q(GlobalC::ORB.orbital_file, q_pack->wf_ptr->table_local); -// // Wavefunc_in_pw::make_table_q(q_pack->wf_ptr->ORB_ptr->orbital_file, q_pack->wf_ptr->table_local, q_pack->wf_ptr); -// for(int iq=0; iqkv_ptr->get_nks(); ++iq) -// { -// Wavefunc_in_pw::produce_local_basis_in_pw(iq, -// wfc_basis, -// sf, -// q_pack->wf_ptr->wanf2[iq], -// q_pack->wf_ptr->table_local); -// // Wavefunc_in_pw::produce_local_basis_in_pw(iq, q_pack->wf_ptr->wanf2[iq], q_pack->wf_ptr->table_local, -// // q_pack->wf_ptr); -// } -// #endif -// q_pack->wf_wg.create(q_pack->kv_ptr->get_nks(),PARAM.inp.nbands); -// if(!GlobalV::RANK_IN_POOL) -// { -// std::stringstream ss_wf_wg; -// ss_wf_wg << PARAM.globalv.global_out_dir << exx_q_pack << "wf_wg_" << GlobalV::MY_POOL; -// std::ifstream ifs_wf_wg(ss_wf_wg.str().c_str()); -// for( int iq = 0; iq < q_pack->kv_ptr->get_nks(); ++iq) -// { -// for( int ib=0; ib>q_pack->wf_wg(iq,ib); -// } -// } -// ifs_wf_wg.close(); -// } -// #ifdef __MPI -// MPI_Bcast( q_pack->wf_wg.c, q_pack->kv_ptr->get_nks()*PARAM.inp.nbands, MPI_DOUBLE, 0, POOL_WORLD); -// #endif - -// q_pack->hvec_array = new ModuleBase::ComplexMatrix [q_pack->kv_ptr->get_nks()]; -// for( int iq=0; iqkv_ptr->get_nks(); ++iq) -// { -// q_pack->hvec_array[iq].create(PARAM.globalv.nlocal,PARAM.inp.nbands); -// } -// if(!GlobalV::RANK_IN_POOL) -// { -// std::stringstream ss_hvec; -// ss_hvec << PARAM.globalv.global_out_dir << exx_q_pack << "hvec_" << GlobalV::MY_POOL; -// std::ifstream ifs_hvec(ss_hvec.str().c_str()); -// for( int iq=0; iqkv_ptr->get_nks(); ++iq) -// { -// for( int iw=0; iw>a>>b; -// q_pack->hvec_array[iq](iw,ib) = {a,b}; -// } -// } -// } -// ifs_hvec.close(); -// } -// #ifdef __MPI -// for( int iq=0; iqkv_ptr->get_nks(); ++iq) -// { -// MPI_Bcast( q_pack->hvec_array[iq].c, PARAM.globalv.nlocal*PARAM.inp.nbands, MPI_DOUBLE_COMPLEX, 0, POOL_WORLD); -// } -// #endif - -// return; -// } - -template class Exx_Lip, base_device::DEVICE_CPU>; -template class Exx_Lip, base_device::DEVICE_CPU>; \ No newline at end of file diff --git a/source/module_ri/exx_lip.h b/source/module_ri/exx_lip.h index 0141d15c91..5d7916f888 100644 --- a/source/module_ri/exx_lip.h +++ b/source/module_ri/exx_lip.h @@ -5,19 +5,26 @@ #ifndef EXX_LIP_H #define EXX_LIP_H -#include "module_base/complexmatrix.h" -#include "module_base/vector3.h" #include "module_hamilt_general/module_xc/exx_info.h" -#include "module_basis/module_pw/pw_basis_k.h" -#include "module_elecstate/elecstate.h" -#include "module_cell/module_symmetry/symmetry.h" -#include "module_hamilt_pw/hamilt_pwdft/wfinit.h" -class K_Vectors; -class UnitCell; +#include "module_base/macros.h" +#include "module_base/matrix.h" + +#include +#include + + class K_Vectors; + class UnitCell; + class Structure_Factor; + namespace elecstate{ class ElecState; } + namespace ModulePW{ class PW_Basis_K; } + namespace ModulePW{ class PW_Basis; } + namespace ModuleSymmetry{ class Symmetry; } + namespace psi{ template class WFInit; } + template class Exx_Lip { - using Real = typename GetTypeReal::type; + using Treal = typename GetTypeReal::type; public: Exx_Lip(const Exx_Info::Exx_Info_Lip& info_in); ~Exx_Lip(); @@ -35,16 +42,16 @@ class Exx_Lip const Structure_Factor& sf, const UnitCell* ucell_ptr_in, const elecstate::ElecState* pelec_in); - + // void cal_exx(const int& nks); void cal_exx(); const std::vector>>& get_exx_matrix() const { - return exx_matrix; + return this->exx_matrix; } - Real get_exx_energy() const + Treal get_exx_energy() const { - return exx_energy; + return this->exx_energy; } void write_q_pack() const; @@ -55,7 +62,7 @@ class Exx_Lip } psi::Psi get_hvec() const { - return *k_pack->hvec_array; + return *this->k_pack->hvec_array; } private: @@ -80,32 +87,32 @@ class Exx_Lip std::vector> phi; std::vector>> psi; - std::vector recip_qkg2; - Real sum2_factor; + std::vector recip_qkg2; + Treal sum2_factor; std::vector b; std::vector b0; std::vector sum1; std::vector> sum3; std::vector>> exx_matrix; - Real exx_energy = 0.0; + Treal exx_energy = 0.0; - void wf_wg_cal(); - void phi_cal(k_package* kq_pack, int ikq); + void wf_wg_cal(); + void phi_cal(k_package* kq_pack, const int ikq); void psi_cal(); - void judge_singularity( int ik); - void qkg2_exp(int ik, int iq); - void b_cal(int ik, int iq, int ib); - void sum3_cal(int iq, int ib); - void b_sum(int iq, int ib); - void sum_all(int ik); - void exx_energy_cal(); + void judge_singularity(const int ik); + void qkg2_exp(const int ik, const int iq); + void b_cal(const int ik, int iq, const int ib); + void sum3_cal(const int iq, const int ib); + void b_sum(const int iq, const int ib); + void sum_all(const int ik); + void exx_energy_cal(); // void read_q_pack(const ModuleSymmetry::Symmetry& symm, // const ModulePW::PW_Basis_K* wfc_basis, // const Structure_Factor& sf); //2*pi*i - T two_pi_i = Real(ModuleBase::TWO_PI) * T(0.0, 1.0); + const T two_pi_i = Treal(ModuleBase::TWO_PI) * T(0.0, 1.0); public: const ModulePW::PW_Basis* rho_basis = nullptr; const ModulePW::PW_Basis_K* wfc_basis = nullptr; @@ -113,5 +120,6 @@ class Exx_Lip const UnitCell* ucell_ptr = nullptr; }; +#include "exx_lip.hpp" #endif \ No newline at end of file diff --git a/source/module_ri/exx_lip.hpp b/source/module_ri/exx_lip.hpp new file mode 100644 index 0000000000..e21bcea66e --- /dev/null +++ b/source/module_ri/exx_lip.hpp @@ -0,0 +1,611 @@ +//========================================================== +// AUTHOR : Peize Lin + +// DATE : 2015-03-10 +//========================================================== + +#ifndef EXX_LIP_HPP +#define EXX_LIP_HPP + +#include "exx_lip.h" +#include "module_base/vector3.h" +#include "module_base/global_function.h" +#include "module_base/vector3.h" +#include "module_hamilt_pw/hamilt_pwdft/global.h" +#include "module_cell/klist.h" +#include "module_hamilt_pw/hamilt_pwdft/wavefunc.h" +#include "module_hamilt_lcao/hamilt_lcaodft/wavefunc_in_pw.h" +#include "module_base/lapack_connector.h" +#include "module_base/parallel_global.h" +#include "module_parameter/parameter.h" +#include "module_elecstate/elecstate.h" +#include "module_basis/module_pw/pw_basis_k.h" +#include "module_cell/module_symmetry/symmetry.h" +#include "module_hamilt_pw/hamilt_pwdft/wfinit.h" +#include "module_hamilt_pw/hamilt_pwdft/structure_factor.h" +#include "module_base/tool_title.h" +#include "module_base/timer.h" + +#include + +template +void Exx_Lip::cal_exx() +{ + ModuleBase::TITLE("Exx_Lip", "cal_exx"); + ModuleBase::timer::tick("Exx_Lip", "cal_exx"); + + this->wf_wg_cal(); + this->psi_cal(); + for (int ik = 0; ik < this->k_pack->kv_ptr->get_nks(); ++ik) + { + this->phi_cal(this->k_pack, ik); + + this->judge_singularity(ik); + for (int iw_l = 0; iw_l < PARAM.globalv.nlocal; ++iw_l) { + for (int iw_r = 0; iw_r < PARAM.globalv.nlocal; ++iw_r) { + this->sum1[iw_l * PARAM.globalv.nlocal + iw_r] = T(0.0); + } } + if (Conv_Coulomb_Pot_K::Ccp_Type::Ccp == info.ccp_type || Conv_Coulomb_Pot_K::Ccp_Type::Hf == info.ccp_type) + { + this->sum2_factor = 0.0; + if (gzero_rank_in_pool == GlobalV::RANK_IN_POOL) { + for (int iw_l = 0; iw_l < PARAM.globalv.nlocal; ++iw_l) { + for (int iw_r = 0; iw_r < PARAM.globalv.nlocal; ++iw_r) { + this->sum3[iw_l][iw_r] = T(0.0); + } } } + } + + for (int iq_tmp = this->iq_vecik; iq_tmp < this->iq_vecik + this->q_pack->kv_ptr->get_nks() / PARAM.inp.nspin; ++iq_tmp) // !!! k_point + // parallel incompleted.need to loop iq in other pool + { + const int iq = + (ik < (this->k_pack->kv_ptr->get_nks() / PARAM.inp.nspin)) + ? (iq_tmp % (this->q_pack->kv_ptr->get_nks() / PARAM.inp.nspin)) + : (iq_tmp % (this->q_pack->kv_ptr->get_nks() / PARAM.inp.nspin) + (this->q_pack->kv_ptr->get_nks() / PARAM.inp.nspin)); + this->qkg2_exp(ik, iq); + for (int ib = 0; ib < PARAM.inp.nbands; ++ib) + { + this->b_cal(ik, iq, ib); + if (Conv_Coulomb_Pot_K::Ccp_Type::Ccp == info.ccp_type || Conv_Coulomb_Pot_K::Ccp_Type::Hf == info.ccp_type) { + if (iq == this->iq_vecik) { + this->sum3_cal(iq, ib); + } } + this->b_sum(iq, ib); + } + } + this->sum_all(ik); + } + this->exx_energy_cal(); + + ModuleBase::timer::tick("Exx_Lip", "cal_exx"); +} + +template +Exx_Lip::Exx_Lip( + const Exx_Info::Exx_Info_Lip& info_in, + const ModuleSymmetry::Symmetry& symm, + K_Vectors* kv_ptr_in, + psi::WFInit* wf_ptr_in, + psi::Psi* kspw_psi_ptr_in, + // wavefunc* wf_ptr_in, + const ModulePW::PW_Basis_K* wfc_basis_in, + const ModulePW::PW_Basis* rho_basis_in, + const Structure_Factor& sf, + const UnitCell* ucell_ptr_in, + const elecstate::ElecState* pelec_in) : info(info_in) +{ + ModuleBase::TITLE("Exx_Lip", "init"); + ModuleBase::timer::tick("Exx_Lip", "init"); + + this->k_pack = new k_package; + this->k_pack->kv_ptr = kv_ptr_in; + this->k_pack->wf_ptr = wf_ptr_in; + this->k_pack->pelec = pelec_in; + this->k_pack->kspw_psi_ptr = kspw_psi_ptr_in; + this->wfc_basis = wfc_basis_in; + this->rho_basis = rho_basis_in; + this->ucell_ptr = ucell_ptr_in; + + int gzero_judge = -1; + if (this->rho_basis->gg_uniq[0] < 1e-8) + { gzero_judge = GlobalV::RANK_IN_POOL; } + #ifdef __MPI + MPI_Allreduce(&gzero_judge, &gzero_rank_in_pool, 1, MPI_INT, MPI_MAX, POOL_WORLD); + #endif + this->k_pack->wf_wg.create(this->k_pack->kv_ptr->get_nks(),PARAM.inp.nbands); + + this->k_pack->hvec_array = new psi::Psi(this->k_pack->kv_ptr->get_nks(), PARAM.inp.nbands, PARAM.globalv.nlocal); + // this->k_pack->hvec_array = new ModuleBase::ComplexMatrix[this->k_pack->kv_ptr->get_nks()]; + // for( int ik=0; ikk_pack->kv_ptr->get_nks(); ++ik) + // { + // this->k_pack->hvec_array[ik].create(PARAM.globalv.nlocal,PARAM.inp.nbands); + // } + + // if (PARAM.inp.init_chg=="atomic") + { + this->q_pack = this->k_pack; + } + // else if(PARAM.inp.init_chg=="file") + // { + // read_q_pack(symm, this->wfc_basis, sf); + // } + + this->phi.resize(PARAM.globalv.nlocal); + for (int iw = 0; iw < PARAM.globalv.nlocal; ++iw) + { this->phi[iw].resize(this->rho_basis->nrxx); } + + this->psi.resize(this->q_pack->kv_ptr->get_nks()); + for (int iq = 0; iq < this->q_pack->kv_ptr->get_nks(); ++iq) + { + this->psi[iq].resize(PARAM.inp.nbands); + for (int ib = 0; ib < PARAM.inp.nbands; ++ib) + { this->psi[iq][ib].resize(this->rho_basis->nrxx); } + } + + this->recip_qkg2.resize(this->rho_basis->npw); + + this->b.resize(PARAM.globalv.nlocal * this->rho_basis->npw); + + this->sum1.resize(PARAM.globalv.nlocal * PARAM.globalv.nlocal); + + if (Conv_Coulomb_Pot_K::Ccp_Type::Ccp == info.ccp_type || Conv_Coulomb_Pot_K::Ccp_Type::Hf == info.ccp_type) + { + if (gzero_rank_in_pool == GlobalV::RANK_IN_POOL) + { + this->b0.resize(PARAM.globalv.nlocal); + this->sum3.resize(PARAM.globalv.nlocal); + for (int iw_l = 0; iw_l < PARAM.globalv.nlocal; ++iw_l) + { this->sum3[iw_l].resize(PARAM.globalv.nlocal); } + } + } + + this->exx_matrix.resize(this->k_pack->kv_ptr->get_nks()); + for (int ik = 0; ik < this->k_pack->kv_ptr->get_nks(); ++ik) + { + this->exx_matrix[ik].resize(PARAM.globalv.nlocal); + for (int iw_l = 0; iw_l < PARAM.globalv.nlocal; ++iw_l) + { this->exx_matrix[ik][iw_l].resize(PARAM.globalv.nlocal); } + } + + ModuleBase::timer::tick("Exx_Lip", "init"); +} + +template +Exx_Lip::~Exx_Lip() +{ + if (this->k_pack) + { delete this->k_pack->hvec_array; this->k_pack->hvec_array = nullptr; } + delete this->k_pack; this->k_pack = nullptr; + + if (PARAM.inp.init_chg == "atomic") + { + this->q_pack = nullptr; + } + else if (PARAM.inp.init_chg == "file") + { + delete this->q_pack->kv_ptr; this->q_pack->kv_ptr = nullptr; + delete this->q_pack->wf_ptr; this->q_pack->wf_ptr = nullptr; + // delete[] this->q_pack->hvec_array; this->q_pack->hvec_array=nullptr; + delete this->q_pack; this->q_pack = nullptr; + } +} + +template +void Exx_Lip::wf_wg_cal() +{ + ModuleBase::TITLE("Exx_Lip", "wf_wg_cal"); + ModuleBase::timer::tick("Exx_Lip", "wf_wg_cal"); + if (PARAM.inp.nspin == 1) { + for (int ik = 0; ik < this->k_pack->kv_ptr->get_nks(); ++ik) { + for (int ib = 0; ib < PARAM.inp.nbands; ++ib) { + this->k_pack->wf_wg(ik, ib) = this->k_pack->pelec->wg(ik, ib) / 2; + } } } + else if (PARAM.inp.nspin == 2) { + for (int ik = 0; ik < this->k_pack->kv_ptr->get_nks(); ++ik) { + for (int ib = 0; ib < PARAM.inp.nbands; ++ib) { + this->k_pack->wf_wg(ik, ib) = this->k_pack->pelec->wg(ik, ib); + } } } + ModuleBase::timer::tick("Exx_Lip", "wf_wg_cal"); +} + +template +void Exx_Lip::phi_cal(k_package* kq_pack, const int ikq) +{ + ModuleBase::timer::tick("Exx_Lip", "phi_cal"); + std::vector porter (this->wfc_basis->nrxx); + for (int iw = 0; iw < PARAM.globalv.nlocal; ++iw) + { + // this->wfc_basis->recip2real(&kq_pack->wf_ptr->wanf2[ikq](iw,0), porter.data(), ikq); + this->wfc_basis->recip2real(&(kq_pack->wf_ptr->get_psig().lock()->operator()(ikq, iw, 0)), porter.data(), ikq); + int ir = 0; + for (int ix = 0; ix < this->rho_basis->nx; ++ix) + { + const Treal phase_x = kq_pack->kv_ptr->kvec_d[ikq].x * ix / this->rho_basis->nx; + for (int iy = 0; iy < this->rho_basis->ny; ++iy) + { + const Treal phase_xy = phase_x + kq_pack->kv_ptr->kvec_d[ikq].y * iy / this->rho_basis->ny; + for (int iz = this->rho_basis->startz_current; iz < this->rho_basis->startz_current + this->rho_basis->nplane; ++iz) + { + const Treal phase_xyz = phase_xy + kq_pack->kv_ptr->kvec_d[ikq].z * iz / this->rho_basis->nz; + const T exp_tmp = std::exp(phase_xyz * this->two_pi_i); + this->phi[iw][ir] = porter[ir] * exp_tmp; + ++ir; + } + } + } + } + ModuleBase::timer::tick("Exx_Lip", "phi_cal"); +} + +template +void Exx_Lip::psi_cal() +{ + ModuleBase::TITLE("Exx_Lip", "psi_cal"); + ModuleBase::timer::tick("Exx_Lip", "psi_cal"); + if (PARAM.inp.init_chg == "atomic") + { + std::vector porter (this->wfc_basis->nrxx); + for (int iq = 0; iq < this->q_pack->kv_ptr->get_nks(); ++iq) + { + for (int ib = 0; ib < PARAM.inp.nbands; ++ib) + { + this->wfc_basis->recip2real(&(this->q_pack->kspw_psi_ptr->operator()(iq, ib, 0)), porter.data(), iq); + + int ir = 0; + for (int ix = 0; ix < this->rho_basis->nx; ++ix) + { + const Treal phase_x = this->q_pack->kv_ptr->kvec_d[iq].x * ix / this->rho_basis->nx; + for (int iy = 0; iy < this->rho_basis->ny; ++iy) + { + const Treal phase_xy = phase_x + this->q_pack->kv_ptr->kvec_d[iq].y * iy / this->rho_basis->ny; + for (int iz = this->rho_basis->startz_current; iz < this->rho_basis->startz_current + this->rho_basis->nplane; ++iz) + { + const Treal phase_xyz = phase_xy + this->q_pack->kv_ptr->kvec_d[iq].z * iz / this->rho_basis->nz; + const T exp_tmp = std::exp(phase_xyz * this->two_pi_i); + this->psi[iq][ib][ir] = porter[ir] * exp_tmp; + ++ir; + } + } + } + } + } + } + else if (PARAM.inp.init_chg == "file") + { + for (int iq = 0; iq < this->q_pack->kv_ptr->get_nks(); ++iq) + { + this->phi_cal(this->q_pack, iq); + for (int ib = 0; ib < PARAM.inp.nbands; ++ib) + { + ModuleBase::GlobalFunc::ZEROS(this->psi[iq][ib].data(), this->rho_basis->nrxx); + for (int iw = 0; iw < PARAM.globalv.nlocal; ++iw) + { + for (int ir = 0; ir < this->rho_basis->nrxx; ++ir) + { + this->psi[iq][ib][ir] += (*this->q_pack->hvec_array)(iq, ib, iw) * this->phi[iw][ir]; + } + } + } + } + } + /////////////////////////////////////////////////////////////////////////// + // !!! k_point parallel incompleted. need to loop iq in other pool + /////////////////////////////////////////////////////////////////////////// + ModuleBase::timer::tick("Exx_Lip", "psi_cal"); +} + +template +void Exx_Lip::judge_singularity(const int ik) +{ + ModuleBase::timer::tick("Exx_Lip", "judge_singularity"); + if (PARAM.inp.init_chg=="atomic") + { + this->iq_vecik = ik; + } + else if(PARAM.inp.init_chg=="file") + { + Treal min_q_minus_k(std::numeric_limits::max()); + for( int iq=0; iqq_pack->kv_ptr->get_nks(); ++iq) + { + const Treal q_minus_k((this->q_pack->kv_ptr->kvec_c[iq] - this->k_pack->kv_ptr->kvec_c[ik]).norm2()); + if(q_minus_k < min_q_minus_k) + { + min_q_minus_k = q_minus_k; + this->iq_vecik = iq; + } + } + } + ModuleBase::timer::tick("Exx_Lip", "judge_singularity"); +} + +template +void Exx_Lip::qkg2_exp(const int ik, const int iq) +{ + ModuleBase::timer::tick("Exx_Lip", "qkg2_exp"); + for( int ig=0; igrho_basis->npw; ++ig) + { + const Treal qkg2 = ((this->q_pack->kv_ptr->kvec_c[iq] - this->k_pack->kv_ptr->kvec_c[ik] + this->rho_basis->gcar[ig]) * (ModuleBase::TWO_PI / this->ucell_ptr->lat0)).norm2(); + if (Conv_Coulomb_Pot_K::Ccp_Type::Ccp == info.ccp_type || Conv_Coulomb_Pot_K::Ccp_Type::Hf == info.ccp_type) + { + if (std::abs(qkg2) < 1e-10) + { this->recip_qkg2[ig] = 0.0; } // 0 to ignore bb/qkg2 when qkg2==0 + else + { this->recip_qkg2[ig] = 1.0 / qkg2; } + this->sum2_factor += this->recip_qkg2[ig] * std::exp(-info.lambda * qkg2); + this->recip_qkg2[ig] = sqrt(this->recip_qkg2[ig]); + } + else if (Conv_Coulomb_Pot_K::Ccp_Type::Hse == info.ccp_type) + { + if (std::abs(qkg2) < 1e-10) + { this->recip_qkg2[ig] = 1.0 / (2 * info.hse_omega); } + else + { this->recip_qkg2[ig] = sqrt((1 - std::exp(-qkg2 / (4 * info.hse_omega * info.hse_omega))) / qkg2); } + } + } + ModuleBase::timer::tick("Exx_Lip", "qkg2_exp"); +} + +template +void Exx_Lip::b_cal(const int ik, const int iq, const int ib) +{ + ModuleBase::timer::tick("Exx_Lip", "b_cal"); + const ModuleBase::Vector3 q_minus_k = this->q_pack->kv_ptr->kvec_d[iq] - this->k_pack->kv_ptr->kvec_d[ik]; + std::vector mul_tmp(this->rho_basis->nrxx); + for( size_t ir=0,ix=0; ixrho_basis->nx; ++ix) + { + const Treal phase_x = q_minus_k.x * ix / this->rho_basis->nx; + for( size_t iy=0; iyrho_basis->ny; ++iy) + { + const Treal phase_xy = phase_x + q_minus_k.y * iy / this->rho_basis->ny; + for( size_t iz=this->rho_basis->startz_current; izrho_basis->startz_current+this->rho_basis->nplane; ++iz) + { + const Treal phase_xyz = phase_xy + q_minus_k.z * iz / this->rho_basis->nz; + mul_tmp[ir] = std::exp(-phase_xyz * this->two_pi_i); + mul_tmp[ir] *= this->psi[iq][ib][ir]; + ++ir; + } + } + } + + std::vector porter (this->rho_basis->nrxx); + for(size_t iw=0; iw< PARAM.globalv.nlocal; ++iw) + { + auto& phi_w = this->phi[iw]; + for( size_t ir=0; irrho_basis->nrxx; ++ir) + { + porter[ir] = conj(phi_w[ir]) * mul_tmp[ir] ; + // porter[ir] = phi_w[ir] * psi_q_b[ir] *exp_tmp[ir] ; + } + T* const b_w = &this->b[iw * this->rho_basis->npw]; + this->rho_basis->real2recip( porter.data(), b_w); + if (Conv_Coulomb_Pot_K::Ccp_Type::Ccp == info.ccp_type || Conv_Coulomb_Pot_K::Ccp_Type::Hf == info.ccp_type) { + if ((iq == this->iq_vecik) && (gzero_rank_in_pool == GlobalV::RANK_IN_POOL)) { /// need to check while use k_point parallel + this->b0[iw] = b_w[this->rho_basis->ig_gge0]; + } } + + for (size_t ig = 0; ig < this->rho_basis->npw; ++ig) + { b_w[ig] *= this->recip_qkg2[ig]; } + } + ModuleBase::timer::tick("Exx_Lip", "b_cal"); +} + +template +void Exx_Lip::sum3_cal(const int iq, const int ib) +{ + ModuleBase::timer::tick("Exx_Lip", "sum3_cal"); + if (gzero_rank_in_pool == GlobalV::RANK_IN_POOL) { + for (int iw_l = 0; iw_l < PARAM.globalv.nlocal; ++iw_l) { + for (int iw_r = 0; iw_r < PARAM.globalv.nlocal; ++iw_r) { + this->sum3[iw_l][iw_r] += this->b0[iw_l] * conj(this->b0[iw_r]) * (Treal)this->q_pack->wf_wg(iq, ib); + } } } + ModuleBase::timer::tick("Exx_Lip", "sum3_cal"); +} + +template +void Exx_Lip::b_sum(const int iq, const int ib) // Peize Lin change 2019-04-14 +{ + ModuleBase::timer::tick("Exx_Lip", "b_sum"); + // this->sum1[iw_l,iw_r] += \sum_{ig} this->b[iw_l,ig] * conj(this->b[iw_r,ig]) * this->q_pack->wf_wg(iq,ib) + LapackConnector::herk( + 'U','N', + PARAM.globalv.nlocal, this->rho_basis->npw, + (Treal)this->q_pack->wf_wg(iq, ib), this->b.data(), this->rho_basis->npw, + 1.0, this->sum1.data(), PARAM.globalv.nlocal); + // cblas_zherk( CblasRowMajor, CblasUpper, CblasNoTrans, + // PARAM.globalv.nlocal, this->rho_basis->npw, + // this->q_pack->wf_wg(iq,ib), static_cast(this->b), this->rho_basis->npw, + // 1.0, static_cast(this->sum1), PARAM.globalv.nlocal); + ModuleBase::timer::tick("Exx_Lip", "b_sum"); +} + +template +void Exx_Lip::sum_all(const int ik) +{ + ModuleBase::timer::tick("Exx_Lip", "sum_all"); + Treal sum2_factor_g = 0.0; + const Treal fourpi_div_omega = 4 * (Treal)(ModuleBase::PI / this->ucell_ptr->omega); + const Treal spin_fac = 2.0; + #ifdef __MPI + if (Conv_Coulomb_Pot_K::Ccp_Type::Ccp == info.ccp_type || Conv_Coulomb_Pot_K::Ccp_Type::Hf == info.ccp_type) + { MPI_Reduce(&this->sum2_factor, &sum2_factor_g, 1, MPI_DOUBLE, MPI_SUM, gzero_rank_in_pool, POOL_WORLD); } + #endif + for (size_t iw_l = 1; iw_l < PARAM.globalv.nlocal; ++iw_l) { + for (size_t iw_r = 0; iw_r < iw_l; ++iw_r) { + this->sum1[iw_l * PARAM.globalv.nlocal + iw_r] = conj(this->sum1[iw_r * PARAM.globalv.nlocal + iw_l]); // Peize Lin add conj 2019-04-14 + } } + + for (int iw_l = 0; iw_l < PARAM.globalv.nlocal; ++iw_l) + { + for( int iw_r=0; iw_rexx_matrix[ik][iw_l][iw_r] = -fourpi_div_omega * this->sum1[iw_l * PARAM.globalv.nlocal + iw_r] * spin_fac; + if (Conv_Coulomb_Pot_K::Ccp_Type::Ccp == info.ccp_type || Conv_Coulomb_Pot_K::Ccp_Type::Hf == info.ccp_type) + { + if (gzero_rank_in_pool == GlobalV::RANK_IN_POOL) + { + this->exx_matrix[ik][iw_l][iw_r] += spin_fac * (fourpi_div_omega * this->sum3[iw_l][iw_r] * sum2_factor_g); + this->exx_matrix[ik][iw_l][iw_r] += spin_fac * (-1 / (Treal)sqrt(info.lambda * ModuleBase::PI) * (Treal)(this->q_pack->kv_ptr->get_nks() / PARAM.inp.nspin) * this->sum3[iw_l][iw_r]); + } + } + } + } + ModuleBase::timer::tick("Exx_Lip", "sum_all"); +} + +template +void Exx_Lip::exx_energy_cal() +{ + ModuleBase::TITLE("Exx_Lip","exx_energy_cal"); + ModuleBase::timer::tick("Exx_Lip", "exx_energy_cal"); + + Treal exx_energy_tmp = 0.0; + + for( int ik=0; ikk_pack->kv_ptr->get_nks(); ++ik) { + for( int iw_l=0; iw_lexx_matrix[ik][iw_l][iw_r] * conj((*this->k_pack->hvec_array)(ik, ib, iw_l)) * (*this->k_pack->hvec_array)(ik, ib, iw_r)).real() * this->k_pack->wf_wg(ik, ib); + } } } } + #ifdef __MPI + MPI_Allreduce( &exx_energy_tmp, &this->exx_energy, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); // !!! k_point parallel incompleted. different pools have different kv.set_nks(>) deadlock + #endif + this->exx_energy *= (PARAM.inp.nspin==1) ? 2 : 1; + this->exx_energy /= 2; // ETOT = E_band - 1/2 E_exx + ModuleBase::timer::tick("Exx_Lip", "exx_energy_cal"); +} + +template +void Exx_Lip::write_q_pack() const +{ + ModuleBase::timer::tick("Exx_Lip", "write_q_pack"); + + if (PARAM.inp.out_chg[0] == 0) + { return; } + + if (!GlobalV::RANK_IN_POOL) + { + const std::string exx_q_pack = "exx_q_pack/"; + + const std::string command_mkdir = "test -d " + PARAM.globalv.global_out_dir + exx_q_pack + " || mkdir " + PARAM.globalv.global_out_dir + exx_q_pack; + assert( system(command_mkdir.c_str()) == 0); + + const std::string command_kpoint = "test -f " + PARAM.globalv.global_out_dir + exx_q_pack + PARAM.inp.kpoint_file + " || cp " + PARAM.inp.kpoint_file + " " + PARAM.globalv.global_out_dir + exx_q_pack + PARAM.inp.kpoint_file; + assert( system(command_kpoint.c_str()) == 0); + + std::stringstream ss_wf_wg; + ss_wf_wg << PARAM.globalv.global_out_dir << exx_q_pack << "wf_wg_" << GlobalV::MY_POOL; + std::ofstream ofs_wf_wg(ss_wf_wg.str().c_str()); + for( int iq = 0; iq < this->q_pack->kv_ptr->get_nks(); ++iq) + { + for( int ib=0; ibq_pack->wf_wg(iq,ib)<<"\t"; + } + ofs_wf_wg<q_pack->kv_ptr->get_nks(); ++iq) + { + for( int iw=0; iwq_pack->hvec_array)(iq, ib, iw).real() << " " << (*this->q_pack->hvec_array)(iq, ib, iw).imag() << " "; + } + ofs_hvec<wfc_basis, + const Structure_Factor& sf) +{ + const std::string exx_q_pack = "exx_q_pack/"; + this->q_pack = new k_package(); + this->q_pack->kv_ptr = new K_Vectors(); + const std::string exx_kpoint_card = PARAM.globalv.global_out_dir + exx_q_pack + PARAM.inp.kpoint_file; + this->q_pack->kv_ptr->set( symm, exx_kpoint_card, PARAM.inp.nspin, this->ucell_ptr->G, this->ucell_ptr->latvec, GlobalV::ofs_running ); + this->q_pack->wf_ptr = new wavefunc(); + this->q_pack->wf_ptr->allocate(this->q_pack->kv_ptr->get_nkstot(), + this->q_pack->kv_ptr->get_nks(), + this->q_pack->kv_ptr->ngk.data(), + this->wfc_basis->npwk_max); // mohan update 2021-02-25 + // this->q_pack->wf_ptr->init(this->q_pack->kv_ptr->get_nks(),this->q_pack->kv_ptr,this->ucell_ptr,old_pwptr,&ppcell,&GlobalC::ORB,&hm,&Pkpoints); + this->q_pack->wf_ptr->table_local.create(GlobalC::ucell.ntype, GlobalC::ucell.nmax_total, PARAM.globalv.nqx); + // this->q_pack->wf_ptr->table_local.create(this->q_pack->wf_ptr->this->ucell_ptr->ntype, this->q_pack->wf_ptr->this->ucell_ptr->nmax_total, PARAM.globalv.nqx); + #ifdef __LCAO + Wavefunc_in_pw::make_table_q(GlobalC::ORB.orbital_file, this->q_pack->wf_ptr->table_local); + // Wavefunc_in_pw::make_table_q(this->q_pack->wf_ptr->ORB_ptr->orbital_file, this->q_pack->wf_ptr->table_local, this->q_pack->wf_ptr); + for(int iq=0; iqq_pack->kv_ptr->get_nks(); ++iq) + { + Wavefunc_in_pw::produce_local_basis_in_pw(iq, + this->wfc_basis, + sf, + this->q_pack->wf_ptr->wanf2[iq], + this->q_pack->wf_ptr->table_local); + // Wavefunc_in_pw::produce_local_basis_in_pw(iq, this->q_pack->wf_ptr->wanf2[iq], this->q_pack->wf_ptr->table_local, + // this->q_pack->wf_ptr); + } + #endif + this->q_pack->wf_wg.create(this->q_pack->kv_ptr->get_nks(),PARAM.inp.nbands); + if(!GlobalV::RANK_IN_POOL) + { + std::stringstream ss_wf_wg; + ss_wf_wg << PARAM.globalv.global_out_dir << exx_q_pack << "wf_wg_" << GlobalV::MY_POOL; + std::ifstream ifs_wf_wg(ss_wf_wg.str().c_str()); + for( int iq = 0; iq < this->q_pack->kv_ptr->get_nks(); ++iq) + { + for( int ib=0; ib>this->q_pack->wf_wg(iq,ib); + } + } + ifs_wf_wg.close(); + } + #ifdef __MPI + MPI_Bcast( this->q_pack->wf_wg.c, this->q_pack->kv_ptr->get_nks()*PARAM.inp.nbands, MPI_DOUBLE, 0, POOL_WORLD); + #endif + this->q_pack->hvec_array = new ModuleBase::ComplexMatrix [this->q_pack->kv_ptr->get_nks()]; + for( int iq=0; iqq_pack->kv_ptr->get_nks(); ++iq) + { + this->q_pack->hvec_array[iq].create(PARAM.globalv.nlocal,PARAM.inp.nbands); + } + if(!GlobalV::RANK_IN_POOL) + { + std::stringstream ss_hvec; + ss_hvec << PARAM.globalv.global_out_dir << exx_q_pack << "hvec_" << GlobalV::MY_POOL; + std::ifstream ifs_hvec(ss_hvec.str().c_str()); + for( int iq=0; iqq_pack->kv_ptr->get_nks(); ++iq) + { + for( int iw=0; iwb; + ifs_hvec>>a>>this->b; + this->q_pack->hvec_array[iq](iw,ib) = {a,this->b}; + } + } + } + ifs_hvec.close(); + } + #ifdef __MPI + for( int iq=0; iqq_pack->kv_ptr->get_nks(); ++iq) + { + MPI_Bcast( this->q_pack->hvec_array[iq].c, PARAM.globalv.nlocal*PARAM.inp.nbands, MPI_DOUBLE_COMPLEX, 0, POOL_WORLD); + } + #endif + return; +} +*/ + +#endif \ No newline at end of file