diff --git a/source/source_cell/k_vector_utils.cpp b/source/source_cell/k_vector_utils.cpp index 0ea7cb77da..e8eb56397d 100644 --- a/source/source_cell/k_vector_utils.cpp +++ b/source/source_cell/k_vector_utils.cpp @@ -604,8 +604,8 @@ void kvec_ibz_kpoint(K_Vectors& kv, } return; }; - // for output in kpoints file - int ibz_index[kv.get_nkstot()]; + // update map k -> irreducible k + kv.ibz_index.assign( kv.get_nkstot_full(), -1); // -1 means not in ibz_kpoint list // search in all k-poins. for (int i = 0; i < kv.get_nkstot(); ++i) { @@ -667,7 +667,7 @@ void kvec_ibz_kpoint(K_Vectors& kv, // nkstot_ibz indicate the index of ibz kpoint. kvec_d_ibz[nkstot_ibz] = kv.kvec_d[i]; // output in kpoints file - ibz_index[i] = nkstot_ibz; + kv.ibz_index[i] = nkstot_ibz; // the weight should be averged k-point weight. wk_ibz[nkstot_ibz] = weight; @@ -688,7 +688,7 @@ void kvec_ibz_kpoint(K_Vectors& kv, double kmol_new = kv.kvec_d[i].norm2(); double kmol_old = kvec_d_ibz[exist_number].norm2(); - ibz_index[i] = exist_number; + kv.ibz_index[i] = exist_number; // std::cout << "\n kmol_new = " << kmol_new; // std::cout << "\n kmol_old = " << kmol_old; @@ -765,10 +765,10 @@ void kvec_ibz_kpoint(K_Vectors& kv, kv.kvec_d[i].x, kv.kvec_d[i].y, kv.kvec_d[i].z, - ibz_index[i] + 1, - kvec_d_ibz[ibz_index[i]].x, - kvec_d_ibz[ibz_index[i]].y, - kvec_d_ibz[ibz_index[i]].z); + kv.ibz_index[i] + 1, + kvec_d_ibz[kv.ibz_index[i]].x, + kvec_d_ibz[kv.ibz_index[i]].y, + kvec_d_ibz[kv.ibz_index[i]].z); } ss << table << std::endl; skpt = ss.str(); diff --git a/source/source_cell/klist.cpp b/source/source_cell/klist.cpp index 694411b3cf..139092b340 100644 --- a/source/source_cell/klist.cpp +++ b/source/source_cell/klist.cpp @@ -160,6 +160,13 @@ void K_Vectors::set(const UnitCell& ucell, // set the k vectors for the up and down spin this->set_kup_and_kdw(); + // initialize ibz_index + this->ibz_index.resize(this->nkstot_full); + for (int ik = 0; ik < this->nkstot_full; ik++) + { + this->ibz_index[ik] = ik; + } + // get ik2iktot this->cal_ik_global(); diff --git a/source/source_cell/klist.h b/source/source_cell/klist.h index 4fd2f41c1f..98f67918bb 100644 --- a/source/source_cell/klist.h +++ b/source/source_cell/klist.h @@ -126,6 +126,7 @@ class K_Vectors } std::vector ik2iktot; ///<[nks] map ik to the global index of k points + std::vector ibz_index; ///< map k points (before symmetry reduction) to irreducible k-points /** * @brief Updates the k-points to use the irreducible Brillouin zone (IBZ). diff --git a/source/source_cell/pseudo.h b/source/source_cell/pseudo.h index b654abcc8f..67322f78c7 100644 --- a/source/source_cell/pseudo.h +++ b/source/source_cell/pseudo.h @@ -1,6 +1,7 @@ #ifndef PSEUDO_H #define PSEUDO_H +#include #include "source_base/global_function.h" #include "source_io/output.h" diff --git a/source/source_lcao/module_ri/RPA_LRI.h b/source/source_lcao/module_ri/RPA_LRI.h index a6a7b275cb..75161a21a4 100644 --- a/source/source_lcao/module_ri/RPA_LRI.h +++ b/source/source_lcao/module_ri/RPA_LRI.h @@ -49,7 +49,7 @@ template class RPA_LRI const psi::Psi& psi, const elecstate::ElecState* pelec); void out_eigen_vector(const Parallel_Orbitals& parav, const psi::Psi& psi); - void out_struc(const ModuleBase::Matrix3& latvec, const ModuleBase::Matrix3& G); + void out_struc(const UnitCell &ucell); void out_bands(const elecstate::ElecState *pelec); void out_Cs(const UnitCell &ucell); diff --git a/source/source_lcao/module_ri/RPA_LRI.hpp b/source/source_lcao/module_ri/RPA_LRI.hpp index a505188b57..d17479b239 100644 --- a/source/source_lcao/module_ri/RPA_LRI.hpp +++ b/source/source_lcao/module_ri/RPA_LRI.hpp @@ -69,6 +69,7 @@ void RPA_LRI::cal_rpa_cv(const UnitCell& ucell) }); std::map>>& Cs = std::get<0>(Cs_dCs); this->Cs_period = RI::RI_Tools::cal_period(Cs, period); + this->Cs_period = exx_lri_rpa.exx_lri.post_2D.set_tensors_map2(this->Cs_period); } template @@ -112,6 +113,11 @@ void RPA_LRI::cal_postSCF_exx(const elecstate::DensityMatrix GlobalC::exx_info.info_global.ccp_type = Conv_Coulomb_Pot_K::Ccp_Type::Hf; GlobalC::exx_info.info_global.hybrid_alpha = 1; GlobalC::exx_info.info_ri.ccp_rmesh_times = PARAM.inp.rpa_ccp_rmesh_times; + GlobalC::exx_info.info_global.coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Fock].resize(1); + GlobalC::exx_info.info_global.coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Fock] = {{ + {"alpha", "1"}, + {"singularity_correction", "spencer"} }}; + exx_lri_rpa.init(mpi_comm_in, ucell, kv, orb); exx_lri_rpa.cal_exx_ions(ucell,PARAM.inp.out_ri_cv); @@ -120,8 +126,7 @@ void RPA_LRI::cal_postSCF_exx(const elecstate::DensityMatrix exx_lri_rpa.cal_exx_elec(Ds, ucell,*dm.get_paraV_pointer(), &symrot); } else { exx_lri_rpa.cal_exx_elec(Ds, ucell,*dm.get_paraV_pointer()); -} - // cout<<"postSCF_Eexx: "< @@ -133,7 +138,7 @@ void RPA_LRI::out_for_RPA(const UnitCell& ucell, ModuleBase::TITLE("DFT_RPA_interface", "out_for_RPA"); this->out_bands(pelec); this->out_eigen_vector(parav, psi); - this->out_struc(ucell.latvec, ucell.G); + this->out_struc(ucell); this->cal_rpa_cv(ucell); std::cout << "rpa_pca_threshold: " << this->info.pca_threshold << std::endl; @@ -217,7 +222,7 @@ void RPA_LRI::out_eigen_vector(const Parallel_Orbitals& parav, const p } template -void RPA_LRI::out_struc(const ModuleBase::Matrix3& latvec, const ModuleBase::Matrix3& G) +void RPA_LRI::out_struc(const UnitCell &ucell) { if (GlobalV::MY_RANK != 0) { @@ -225,14 +230,16 @@ void RPA_LRI::out_struc(const ModuleBase::Matrix3& latvec, const Modul } ModuleBase::TITLE("DFT_RPA_interface", "out_struc"); double TWOPI_Bohr2A = ModuleBase::TWO_PI * ModuleBase::BOHR_TO_A; - const int nks_tot = PARAM.inp.nspin == 2 ? (int)p_kv->get_nks() / 2 : p_kv->get_nks(); - ModuleBase::Matrix3 lat = latvec / ModuleBase::BOHR_TO_A; - ModuleBase::Matrix3 G_RPA = G * TWOPI_Bohr2A; + const int nks_tot = PARAM.inp.nspin == 2 ? (int)p_kv->get_nkstot() / 2 : p_kv->get_nkstot(); + const int nks_tot_full = p_kv->get_nkstot_full(); + const int natom = ucell.nat; + ModuleBase::Matrix3 lat = ucell.latvec / ModuleBase::BOHR_TO_A; + ModuleBase::Matrix3 G_RPA = ucell.G * TWOPI_Bohr2A; std::stringstream ss; ss << "stru_out"; std::ofstream ofs; ofs.open(ss.str().c_str(), std::ios::out); - ofs << lat.e11 << std::setw(15) << lat.e12 << std::setw(15) << lat.e13 << std::endl; + ofs << std::fixed << std::setprecision(9) << lat.e11 << std::setw(15) << lat.e12 << std::setw(15) << lat.e13 << std::endl; ofs << lat.e21 << std::setw(15) << lat.e22 << std::setw(15) << lat.e23 << std::endl; ofs << lat.e31 << std::setw(15) << lat.e32 << std::setw(15) << lat.e33 << std::endl; @@ -240,14 +247,25 @@ void RPA_LRI::out_struc(const ModuleBase::Matrix3& latvec, const Modul ofs << G_RPA.e21 << std::setw(15) << G_RPA.e22 << std::setw(15) << G_RPA.e23 << std::endl; ofs << G_RPA.e31 << std::setw(15) << G_RPA.e32 << std::setw(15) << G_RPA.e33 << std::endl; + ofs << natom << std::endl; + for(int iat=0; iat < natom; ++iat) { + int it = ucell.iat2it[iat]; + ModuleBase::Vector3 atom_pos = ucell.atoms[ it ].tau[ ucell.iat2ia[iat] ]/ ModuleBase::BOHR_TO_A; + ofs << atom_pos.x << std::setw(15) << atom_pos.y << std::setw(15) << atom_pos.z << std::setw(15) << (it + 1) << std::endl; + } + ofs << p_kv->nmp[0] << std::setw(6) << p_kv->nmp[1] << std::setw(6) << p_kv->nmp[2] << std::setw(6) << std::endl; for (int ik = 0; ik != nks_tot; ik++) { - ofs << std::setw(15) << std::fixed << std::setprecision(9) << p_kv->kvec_c[ik].x * TWOPI_Bohr2A << std::setw(15) - << std::fixed << std::setprecision(9) << p_kv->kvec_c[ik].y * TWOPI_Bohr2A << std::setw(15) << std::fixed - << std::setprecision(9) << p_kv->kvec_c[ik].z * TWOPI_Bohr2A << std::endl; + ofs << std::setw(15) << p_kv->kvec_c[ik].x * TWOPI_Bohr2A << std::setw(15) + << p_kv->kvec_c[ik].y * TWOPI_Bohr2A << std::setw(15) << p_kv->kvec_c[ik].z * TWOPI_Bohr2A << std::endl; + } + + for (int ik = 0; ik != nks_tot_full; ++ik){ + ofs << (p_kv->ibz_index[ik] + 1) << std::endl; } + ofs.close(); return; } @@ -381,11 +399,11 @@ void RPA_LRI::out_coulomb_k(const UnitCell &ucell) size_t nu_num = exx_lri_rpa.exx_objs[Conv_Coulomb_Pot_K::Coulomb_Method::Center2].cv.get_index_abfs_size(ucell.iat2it[iJ]); ofs << all_mu << " " << mu_shift[I] + 1 << " " << mu_shift[I] + mu_num << " " << mu_shift[iJ] + 1 << " " << mu_shift[iJ] + nu_num << std::endl; - ofs << ik + 1 << " " << p_kv->wk[ik] / 2.0 * PARAM.inp.nspin << std::endl; + ofs << ik + 1 << " " << std::fixed << std::setprecision(12) << p_kv->wk[ik] / 2.0 * PARAM.inp.nspin << std::endl; for (int i = 0; i != vq_J.data->size(); i++) { - ofs << std::setw(21) << std::fixed << std::setprecision(12) << (*vq_J.data)[i].real() - << std::setw(21) << std::fixed << std::setprecision(12) << (*vq_J.data)[i].imag() << std::endl; + ofs << std::setw(21) << (*vq_J.data)[i].real() + << std::setw(21) << (*vq_J.data)[i].imag() << std::endl; } } }