From 1dbe371e8f2ca601e686ab7f7f0a727f3476cc4b Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Mon, 8 Jul 2024 13:18:37 +0800 Subject: [PATCH 01/21] Add multi-k decomposed charge density calculation for LCAO (draft only) --- .../module_esolver/esolver_ks_lcao_elec.cpp | 3 +- source/module_io/istate_charge.cpp | 283 +++++++++++++++++- source/module_io/istate_charge.h | 53 +++- 3 files changed, 329 insertions(+), 10 deletions(-) diff --git a/source/module_esolver/esolver_ks_lcao_elec.cpp b/source/module_esolver/esolver_ks_lcao_elec.cpp index 92a39ccbea..8f1f20a49e 100644 --- a/source/module_esolver/esolver_ks_lcao_elec.cpp +++ b/source/module_esolver/esolver_ks_lcao_elec.cpp @@ -448,7 +448,8 @@ void ESolver_KS_LCAO::others(const int istep) { GlobalV::MY_RANK, GlobalV::ofs_warning, &GlobalC::ucell, - &GlobalC::GridD); + &GlobalC::GridD, + this->kv); } else if (cal_type == "get_wf") { IState_Envelope IEP(this->pelec); if (GlobalV::GAMMA_ONLY_LOCAL) { diff --git a/source/module_io/istate_charge.cpp b/source/module_io/istate_charge.cpp index f12cd9e903..29b0174970 100644 --- a/source/module_io/istate_charge.cpp +++ b/source/module_io/istate_charge.cpp @@ -16,6 +16,11 @@ IState_Charge::IState_Charge(psi::Psi* psi_gamma_in, const Parallel_Orbi { } +IState_Charge::IState_Charge(psi::Psi>* psi_k_in, const Parallel_Orbitals* ParaV_in) + : psi_k(psi_k_in), ParaV(ParaV_in) +{ +} + IState_Charge::~IState_Charge() { } @@ -43,7 +48,8 @@ void IState_Charge::begin(Gint_Gamma& gg, const int my_rank, std::ofstream& ofs_warning, const UnitCell* ucell_in, - Grid_Driver* GridD_in) + Grid_Driver* GridD_in, + const K_Vectors& kv) { ModuleBase::TITLE("IState_Charge", "begin"); @@ -206,7 +212,7 @@ void IState_Charge::begin(Gint_Gamma& gg, elecstate::DensityMatrix DM(this->ParaV, nspin); #ifdef __MPI - this->idmatrix(ib, nspin, nelec, nlocal, wg, DM); + this->idmatrix(ib, nspin, nelec, nlocal, wg, DM, kv); #else ModuleBase::WARNING_QUIT("IState_Charge::begin", "The `pchg` calculation is only available for MPI now!"); #endif @@ -281,13 +287,284 @@ void IState_Charge::begin(Gint_Gamma& gg, return; } +void IState_Charge::begin(Gint_k& gk, + double** rho, + const ModuleBase::matrix& wg, + const std::vector& ef_all_spin, + const int rhopw_nrxx, + const int rhopw_nplane, + const int rhopw_startz_current, + const int rhopw_nx, + const int rhopw_ny, + const int rhopw_nz, + const int bigpw_bz, + const int bigpw_nbz, + const bool gamma_only_local, + const int nbands_istate, + const std::vector& out_band_kb, + const int nbands, + const double nelec, + const int nspin, + const int nlocal, + const std::string& global_out_dir, + const int my_rank, + std::ofstream& ofs_warning, + const UnitCell* ucell_in, + Grid_Driver* GridD_in, + const K_Vectors& kv) +{ + ModuleBase::TITLE("IState_Charge", "begin"); + + std::cout << " Perform |psi(i)|^2 for selected bands (multi-k)." << std::endl; + + int mode = 0; + if (nbands_istate > 0 && static_cast(out_band_kb.size()) == 0) + { + mode = 1; + } + else if (static_cast(out_band_kb.size()) > 0) + { + mode = 2; + std::cout << " Notice: INPUT parameter `nbands_istate` overwritten by `bands_to_print`!" << std::endl; + } + else + { + mode = 3; + } + + int fermi_band = 0; + int bands_below = 0; + int bands_above = 0; + + this->bands_picked_.resize(nbands); + ModuleBase::GlobalFunc::ZEROS(bands_picked_.data(), nbands); + + std::cout << " number of electrons = " << nelec << std::endl; + fermi_band = static_cast((nelec + 1) / 2 + 1.0e-8); + std::cout << " number of occupied bands = " << fermi_band << std::endl; + + if (mode == 1) + { + bands_below = nbands_istate; + bands_above = nbands_istate; + + std::cout << " Plot band decomposed charge density below Fermi surface with " << bands_below << " bands." + << std::endl; + + std::cout << " Plot band decomposed charge density above Fermi surface with " << bands_above << " bands." + << std::endl; + + for (int ib = 0; ib < nbands; ++ib) + { + if (ib >= fermi_band - bands_below) + { + if (ib < fermi_band + bands_above) + { + bands_picked_[ib] = 1; + } + } + } + } + else if (mode == 2) + { + if (static_cast(out_band_kb.size()) > nbands) + { + ModuleBase::WARNING_QUIT( + "IState_Charge::begin", + "The number of bands specified by `bands_to_print` in the INPUT file exceeds `nbands`!"); + } + for (int value: out_band_kb) + { + if (value != 0 && value != 1) + { + ModuleBase::WARNING_QUIT( + "IState_Charge::begin", + "The elements of `bands_to_print` must be either 0 or 1. Invalid values found!"); + } + } + int length = std::min(static_cast(out_band_kb.size()), nbands); + for (int i = 0; i < length; ++i) + { + bands_picked_[i] = out_band_kb[i]; + } + + std::cout << " Plot band decomposed charge density below the Fermi surface: band "; + for (int i = 0; i + 1 <= fermi_band; ++i) + { + if (bands_picked_[i] == 1) + { + std::cout << i + 1 << " "; + } + } + std::cout << std::endl; + std::cout << " Plot band decomposed charge density above the Fermi surface: band "; + for (int i = fermi_band; i < nbands; ++i) + { + if (bands_picked_[i] == 1) + { + std::cout << i + 1 << " "; + } + } + std::cout << std::endl; + } + else if (mode == 3) + { + bool stop = false; + std::stringstream ss; + ss << global_out_dir << "istate.info"; + std::cout << " Open the file : " << ss.str() << std::endl; + if (my_rank == 0) + { + std::ifstream ifs(ss.str().c_str()); + if (!ifs) + { + stop = true; + } + else + { + for (int ib = 0; ib < nbands; ++ib) + { + ModuleBase::GlobalFunc::READ_VALUE(ifs, bands_picked_[ib]); + } + } + } + +#ifdef __MPI + Parallel_Common::bcast_bool(stop); + Parallel_Common::bcast_int(bands_picked_.data(), nbands); +#endif + if (stop) + { + ofs_warning << " Can't find the file : " << ss.str() << std::endl; + ModuleBase::WARNING_QUIT("IState_Charge::begin", "can't find the istate file."); + } + } + + for (int ib = 0; ib < nbands; ++ib) + { + if (bands_picked_[ib]) + { + std::cout << " Perform band decomposed charge density for band " << ib + 1 << std::endl; + + elecstate::DensityMatrix, double> DM(this->ParaV, nspin); + +#ifdef __MPI + this->idmatrix(ib, nspin, nelec, nlocal, wg, DM, kv); +#else + ModuleBase::WARNING_QUIT("IState_Charge::begin", "The `pchg` calculation is only available for MPI now!"); +#endif + + for (int is = 0; is < nspin; ++is) + { + ModuleBase::GlobalFunc::ZEROS(rho[is], rhopw_nrxx); + } + + DM.init_DMR(GridD_in, ucell_in); + DM.cal_DMR(); + + gk.initialize_pvpR(*ucell_in, GridD_in); + + gk.transfer_DM2DtoGrid(DM.get_DMR_vector()); + + Gint_inout inout((double***)nullptr, rho, Gint_Tools::job_type::rho); + gk.cal_gint(&inout); + + double** rho_save = new double*[nspin]; + for (int is = 0; is < nspin; is++) + { + rho_save[is] = new double[rhopw_nrxx]; + ModuleBase::GlobalFunc::DCOPY(rho[is], rho_save[is], rhopw_nrxx); + } + + std::stringstream ssc; + ssc << global_out_dir << "BAND" << ib + 1; + for (int is = 0; is < nspin; ++is) + { + ssc << "_SPIN" << is << "_CHG.cube"; + + double ef_spin = ef_all_spin[is]; + ModuleIO::write_rho( +#ifdef __MPI + bigpw_bz, + bigpw_nbz, + rhopw_nplane, + rhopw_startz_current, +#endif + rho_save[is], + is, + nspin, + 0, + ssc.str(), + rhopw_nx, + rhopw_ny, + rhopw_nz, + ef_spin, + ucell_in); + } + + for (int is = 0; is < nspin; is++) + { + delete[] rho_save[is]; + } + delete[] rho_save; + } + } + + return; +} + +#ifdef __MPI +void IState_Charge::idmatrix(const int& ib, + const int nspin, + const double& nelec, + const int nlocal, + const ModuleBase::matrix& wg, + elecstate::DensityMatrix, double>& DM, + const K_Vectors& kv) +{ + ModuleBase::TITLE("IState_Charge", "idmatrix"); + assert(wg.nr == nspin); + + int fermi_band = static_cast((nelec + 1) / 2 + 1.0e-8); + + for (int is = 0; is < nspin; ++is) + { + std::vector wg_local(this->ParaV->ncol, 0.0); + const int ib_local = this->ParaV->global2local_col(ib); + + if (ib_local >= 0) + { + wg_local[ib_local] = (ib < fermi_band) ? wg(is, ib) : wg(is, fermi_band - 1); + } + + this->psi_k->fix_k(is); + psi::Psi> wg_wfc(*this->psi_k, 1); + + for (int ir = 0; ir < wg_wfc.get_nbands(); ++ir) + { + BlasConnector::scal(wg_wfc.get_nbasis(), wg_local[ir], wg_wfc.get_pointer() + ir * wg_wfc.get_nbasis(), 1); + } + + for (int ik = 0; ik < kv.get_nks(); ++ik) + { + elecstate::psiMulPsiMpi(wg_wfc, + *(this->psi_k), + DM.get_DMK_pointer(ik), + this->ParaV->desc_wfc, + this->ParaV->desc); + } + } +} +#endif + #ifdef __MPI void IState_Charge::idmatrix(const int& ib, const int nspin, const double& nelec, const int nlocal, const ModuleBase::matrix& wg, - elecstate::DensityMatrix& DM) + elecstate::DensityMatrix& DM, + const K_Vectors& kv) { ModuleBase::TITLE("IState_Charge", "idmatrix"); assert(wg.nr == nspin); diff --git a/source/module_io/istate_charge.h b/source/module_io/istate_charge.h index 579491d510..57445adc49 100644 --- a/source/module_io/istate_charge.h +++ b/source/module_io/istate_charge.h @@ -5,6 +5,7 @@ #include "module_hamilt_lcao/hamilt_lcaodft/local_orbital_charge.h" #include "module_hamilt_lcao/module_gint/gint.h" #include "module_hamilt_lcao/module_gint/gint_gamma.h" +#include "module_hamilt_lcao/module_gint/gint_k.h" #include "module_psi/psi.h" #include @@ -24,13 +25,14 @@ class IState_Charge { public: IState_Charge(psi::Psi* psi_gamma_in, const Parallel_Orbitals* ParaV_in); - IState_Charge(psi::Psi>* psi_k_in, const Parallel_Orbitals* ParaV_in) - { - throw std::logic_error("IState_Charge for multi-k is not implemented."); - }; + IState_Charge(psi::Psi>* psi_k_in, const Parallel_Orbitals* ParaV_in); + // { + // throw std::logic_error("IState_Charge for multi-k is not implemented."); + // }; ~IState_Charge(); + // for gamma_only void begin(Gint_Gamma& gg, double** rho, const ModuleBase::matrix& wg, @@ -54,7 +56,35 @@ class IState_Charge const int my_rank, std::ofstream& ofs_warning, const UnitCell* ucell_in, - Grid_Driver* GridD_in); + Grid_Driver* GridD_in, + const K_Vectors& kv); + + // For multi-k + void begin(Gint_k& gk, + double** rho, + const ModuleBase::matrix& wg, + const std::vector& ef_all_spin, + const int rhopw_nrxx, + const int rhopw_nplane, + const int rhopw_startz_current, + const int rhopw_nx, + const int rhopw_ny, + const int rhopw_nz, + const int bigpw_bz, + const int bigpw_nbz, + const bool gamma_only_local, + const int nbands_istate, + const std::vector& out_band_kb, + const int nbands, + const double nelec, + const int nspin, + const int nlocal, + const std::string& global_out_dir, + const int my_rank, + std::ofstream& ofs_warning, + const UnitCell* ucell_in, + Grid_Driver* GridD_in, + const K_Vectors& kv); private: std::vector bands_picked_; @@ -79,9 +109,20 @@ class IState_Charge const double& nelec, const int nlocal, const ModuleBase::matrix& wg, - elecstate::DensityMatrix& DM); + elecstate::DensityMatrix& DM, + const K_Vectors& kv); + + void idmatrix(const int& ib, + const int nspin, + const double& nelec, + const int nlocal, + const ModuleBase::matrix& wg, + elecstate::DensityMatrix, double>& DM, + const K_Vectors& kv); + #endif psi::Psi* psi_gamma; + psi::Psi>* psi_k; const Parallel_Orbitals* ParaV; }; #endif From cdaa1cdcaf884df667f03b17b25129d3f6feb268 Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Tue, 9 Jul 2024 13:17:26 +0800 Subject: [PATCH 02/21] Add multi-k decomposed charge density calculation for LCAO (draft 2) --- .../module_esolver/esolver_ks_lcao_elec.cpp | 78 +++++++++++++------ source/module_io/istate_charge.cpp | 30 +++++-- 2 files changed, 77 insertions(+), 31 deletions(-) diff --git a/source/module_esolver/esolver_ks_lcao_elec.cpp b/source/module_esolver/esolver_ks_lcao_elec.cpp index 534ff1ecac..f96e73ca5e 100644 --- a/source/module_esolver/esolver_ks_lcao_elec.cpp +++ b/source/module_esolver/esolver_ks_lcao_elec.cpp @@ -432,31 +432,59 @@ void ESolver_KS_LCAO::others(const int istep) { } else if (cal_type == "get_pchg") { std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "getting partial charge"); IState_Charge ISC(this->psi, &(this->orb_con.ParaV)); - ISC.begin(this->GG, - this->pelec->charge->rho, - this->pelec->wg, - this->pelec->eferm.get_all_ef(), - this->pw_rho->nrxx, - this->pw_rho->nplane, - this->pw_rho->startz_current, - this->pw_rho->nx, - this->pw_rho->ny, - this->pw_rho->nz, - this->pw_big->bz, - this->pw_big->nbz, - GlobalV::GAMMA_ONLY_LOCAL, - GlobalV::NBANDS_ISTATE, - INPUT.get_out_band_kb(), - GlobalV::NBANDS, - GlobalV::nelec, - GlobalV::NSPIN, - GlobalV::NLOCAL, - GlobalV::global_out_dir, - GlobalV::MY_RANK, - GlobalV::ofs_warning, - &GlobalC::ucell, - &GlobalC::GridD, - this->kv); + if (GlobalV::GAMMA_ONLY_LOCAL) { + ISC.begin(this->GG, + this->pelec->charge->rho, + this->pelec->wg, + this->pelec->eferm.get_all_ef(), + this->pw_rho->nrxx, + this->pw_rho->nplane, + this->pw_rho->startz_current, + this->pw_rho->nx, + this->pw_rho->ny, + this->pw_rho->nz, + this->pw_big->bz, + this->pw_big->nbz, + GlobalV::GAMMA_ONLY_LOCAL, + GlobalV::NBANDS_ISTATE, + INPUT.get_out_band_kb(), + GlobalV::NBANDS, + GlobalV::nelec, + GlobalV::NSPIN, + GlobalV::NLOCAL, + GlobalV::global_out_dir, + GlobalV::MY_RANK, + GlobalV::ofs_warning, + &GlobalC::ucell, + &GlobalC::GridD, + this->kv); + } else { + ISC.begin(this->GK, + this->pelec->charge->rho, + this->pelec->wg, + this->pelec->eferm.get_all_ef(), + this->pw_rho->nrxx, + this->pw_rho->nplane, + this->pw_rho->startz_current, + this->pw_rho->nx, + this->pw_rho->ny, + this->pw_rho->nz, + this->pw_big->bz, + this->pw_big->nbz, + GlobalV::GAMMA_ONLY_LOCAL, + GlobalV::NBANDS_ISTATE, + INPUT.get_out_band_kb(), + GlobalV::NBANDS, + GlobalV::nelec, + GlobalV::NSPIN, + GlobalV::NLOCAL, + GlobalV::global_out_dir, + GlobalV::MY_RANK, + GlobalV::ofs_warning, + &GlobalC::ucell, + &GlobalC::GridD, + this->kv); + } std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "getting partial charge"); } else if (cal_type == "get_wf") { std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "getting wave function"); diff --git a/source/module_io/istate_charge.cpp b/source/module_io/istate_charge.cpp index 29b0174970..6ce7afea58 100644 --- a/source/module_io/istate_charge.cpp +++ b/source/module_io/istate_charge.cpp @@ -444,9 +444,11 @@ void IState_Charge::begin(Gint_k& gk, { if (bands_picked_[ib]) { - std::cout << " Perform band decomposed charge density for band " << ib + 1 << std::endl; + // std::cout << " Perform band decomposed charge density for kpoint " << ik << ", band" << ib + 1 + // << std::endl; - elecstate::DensityMatrix, double> DM(this->ParaV, nspin); + // elecstate::DensityMatrix, double> DM(this->ParaV, nspin); + elecstate::DensityMatrix, double> DM(&kv, this->ParaV, nspin); #ifdef __MPI this->idmatrix(ib, nspin, nelec, nlocal, wg, DM, kv); @@ -523,10 +525,20 @@ void IState_Charge::idmatrix(const int& ib, const K_Vectors& kv) { ModuleBase::TITLE("IState_Charge", "idmatrix"); - assert(wg.nr == nspin); + // assert(wg.nr == nspin); + std::cout << "wg.nr = " << wg.nr << std::endl; int fermi_band = static_cast((nelec + 1) / 2 + 1.0e-8); + for (int ik = 0; ik < kv.get_nks(); ++ik) + { + std::cout << " Perform band decomposed charge density for kpoint " << ik << ", band" << ib + 1 << std::endl; + + const int ispin = kv.isk[ik]; + + + } + for (int is = 0; is < nspin; ++is) { std::vector wg_local(this->ParaV->ncol, 0.0); @@ -544,14 +556,16 @@ void IState_Charge::idmatrix(const int& ib, { BlasConnector::scal(wg_wfc.get_nbasis(), wg_local[ir], wg_wfc.get_pointer() + ir * wg_wfc.get_nbasis(), 1); } - + std::cout << "Here is Okay." << std::endl; for (int ik = 0; ik < kv.get_nks(); ++ik) { + std::cout << "Here is Okay." << std::endl; elecstate::psiMulPsiMpi(wg_wfc, *(this->psi_k), DM.get_DMK_pointer(ik), this->ParaV->desc_wfc, this->ParaV->desc); + std::cout << "Here is Okay." << std::endl; } } } @@ -569,6 +583,10 @@ void IState_Charge::idmatrix(const int& ib, ModuleBase::TITLE("IState_Charge", "idmatrix"); assert(wg.nr == nspin); + std::cout << "wg.nr = " << wg.nr << std::endl; + std::cout << "wg.nc = " << wg.nc << std::endl; + std::cout << "nspin = " << nspin << std::endl; + int fermi_band = static_cast((nelec + 1) / 2 + 1.0e-8); for (int is = 0; is < nspin; ++is) @@ -590,11 +608,11 @@ void IState_Charge::idmatrix(const int& ib, { BlasConnector::scal(wg_wfc.get_nbasis(), wg_local[ir], wg_wfc.get_pointer() + ir * wg_wfc.get_nbasis(), 1); } - const int ik = 0; // Gamma point only + // const int ik = 0; // Gamma point only elecstate::psiMulPsiMpi(wg_wfc, *(this->psi_gamma), - DM.get_DMK_pointer(ik), + DM.get_DMK_pointer(is), this->ParaV->desc_wfc, this->ParaV->desc); } From 1df5013f412b7544e833bbeda09f89857333cd68 Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Tue, 9 Jul 2024 13:56:27 +0800 Subject: [PATCH 03/21] Fix wrong cube file name for get_pchg calculation when nspin=2 --- source/module_io/istate_charge.cpp | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/source/module_io/istate_charge.cpp b/source/module_io/istate_charge.cpp index 6ce7afea58..c40a25a013 100644 --- a/source/module_io/istate_charge.cpp +++ b/source/module_io/istate_charge.cpp @@ -246,12 +246,12 @@ void IState_Charge::begin(Gint_Gamma& gg, rhopw_nrxx); // Copy data after allocation } - std::stringstream ssc; - ssc << global_out_dir << "BAND" << ib + 1; // 0 means definitely output charge density. for (int is = 0; is < nspin; ++is) { - ssc << "_SPIN" << is << "_CHG.cube"; + // ssc should be inside the inner loop to reset the string stream each time + std::stringstream ssc; + ssc << global_out_dir << "BAND" << ib + 1 << "_SPIN" << is << "_CHG.cube"; // Use a const vector to store efermi for all spins, replace the original implementation: // const double ef_tmp = pelec->eferm.get_efval(is); @@ -478,11 +478,11 @@ void IState_Charge::begin(Gint_k& gk, ModuleBase::GlobalFunc::DCOPY(rho[is], rho_save[is], rhopw_nrxx); } - std::stringstream ssc; - ssc << global_out_dir << "BAND" << ib + 1; for (int is = 0; is < nspin; ++is) { - ssc << "_SPIN" << is << "_CHG.cube"; + // ssc should be inside the inner loop to reset the string stream each time + std::stringstream ssc; + ssc << global_out_dir << "BAND" << ib + 1 << "_SPIN" << is << "_CHG.cube"; double ef_spin = ef_all_spin[is]; ModuleIO::write_rho( @@ -535,8 +535,6 @@ void IState_Charge::idmatrix(const int& ib, std::cout << " Perform band decomposed charge density for kpoint " << ik << ", band" << ib + 1 << std::endl; const int ispin = kv.isk[ik]; - - } for (int is = 0; is < nspin; ++is) From 1b71e2172a2b1405f4f2c7b3afbaacaf090d9694 Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Tue, 9 Jul 2024 23:35:21 +0800 Subject: [PATCH 04/21] Add multi-k decomposed charge density calculation for LCAO (draft 3) --- .../module_dm/density_matrix.cpp | 271 ++++++++++++------ .../module_dm/density_matrix.h | 5 + source/module_io/istate_charge.cpp | 127 ++++---- 3 files changed, 265 insertions(+), 138 deletions(-) diff --git a/source/module_elecstate/module_dm/density_matrix.cpp b/source/module_elecstate/module_dm/density_matrix.cpp index c80151caeb..e981a0e45d 100644 --- a/source/module_elecstate/module_dm/density_matrix.cpp +++ b/source/module_elecstate/module_dm/density_matrix.cpp @@ -1,9 +1,9 @@ #include "density_matrix.h" #include "module_base/libm/libm.h" -#include "module_base/tool_title.h" #include "module_base/memory.h" #include "module_base/timer.h" +#include "module_base/tool_title.h" #include "module_cell/klist.h" namespace elecstate @@ -126,7 +126,7 @@ void DensityMatrix::init_DMR(Grid_Driver* GridD_in, const UnitCell* ucel } } // allocate the memory of BaseMatrix in SR, and set the new values to zero - if(std::is_same::value) + if (std::is_same::value) { tmp_DMR->fix_gamma(); } @@ -171,12 +171,17 @@ void DensityMatrix::init_DMR(Record_adj& ra, const UnitCell* ucell) { continue; } - hamilt::AtomPair tmp_ap(iat1, iat2, ra.info[iat1][ad][0], ra.info[iat1][ad][1], ra.info[iat1][ad][2], this->_paraV); + hamilt::AtomPair tmp_ap(iat1, + iat2, + ra.info[iat1][ad][0], + ra.info[iat1][ad][1], + ra.info[iat1][ad][2], + this->_paraV); tmp_DMR->insert_pair(tmp_ap); } } // allocate the memory of BaseMatrix in SR, and set the new values to zero - if(std::is_same::value) + if (std::is_same::value) { tmp_DMR->fix_gamma(); } @@ -227,15 +232,15 @@ void DensityMatrix::init_DMR(const hamilt::HContainer& DMR_in) this->_DMR.clear(); // set up a HContainer using another one int size_ap = DMR_in.size_atom_pairs(); - if(size_ap > 0) - { + if (size_ap > 0) + { const Parallel_Orbitals* paraV_ = DMR_in.get_atom_pair(0).get_paraV(); hamilt::HContainer* tmp_DMR = new hamilt::HContainer(paraV_); - for(int iap=0;iap R_index = DMR_in.get_atom_pair(iap).get_R_index(ir); hamilt::AtomPair tmp_ap(iat1, iat2, R_index, paraV_); @@ -244,7 +249,7 @@ void DensityMatrix::init_DMR(const hamilt::HContainer& DMR_in) } tmp_DMR->allocate(nullptr, true); this->_DMR.push_back(tmp_DMR); - if(this->_nspin == 2) + if (this->_nspin == 2) { hamilt::HContainer* tmp_DMR1 = new hamilt::HContainer(*tmp_DMR); this->_DMR.push_back(tmp_DMR1); @@ -301,7 +306,8 @@ void DensityMatrix::set_DMK_zero() { for (int ik = 0; ik < _nspin * _nks; ik++) { - ModuleBase::GlobalFunc::ZEROS(this->_DMK[ik].data(), this->_paraV->get_row_size() * this->_paraV->get_col_size()); + ModuleBase::GlobalFunc::ZEROS(this->_DMK[ik].data(), + this->_paraV->get_row_size() * this->_paraV->get_col_size()); } } @@ -359,29 +365,29 @@ void DensityMatrix::save_DMR() { ModuleBase::TITLE("DensityMatrix", "save_DMR"); ModuleBase::timer::tick("DensityMatrix", "save_DMR"); - + const int nnr = this->_DMR[0]->get_nnr(); - // allocate if _DMR_save is empty - if(_DMR_save.size() == 0) + // allocate if _DMR_save is empty + if (_DMR_save.size() == 0) { _DMR_save.resize(this->_DMR.size()); } // resize if _DMR_save[is].size is not equal to _DMR.size - for(int is = 0; is < _DMR_save.size(); is++) + for (int is = 0; is < _DMR_save.size(); is++) { - if(_DMR_save[is].size() != nnr) + if (_DMR_save[is].size() != nnr) { _DMR_save[is].resize(nnr); } } // save _DMR to _DMR_save - for(int is=0;is_DMR.size();is++) + for (int is = 0; is < this->_DMR.size(); is++) { TR* DMR_pointer = this->_DMR[is]->get_wrapper(); TR* DMR_save_pointer = _DMR_save[is].data(); // set to zero ModuleBase::GlobalFunc::ZEROS(DMR_save_pointer, nnr); - for(int i=0;i::save_DMR() // calculate DMR from DMK using add_element template -void DensityMatrix::cal_DMR_test() +void DensityMatrix::cal_DMR_test() { for (int is = 1; is <= this->_nspin; ++is) { - int ik_begin = this->_nks*(is-1); // jump this->_nks for spin_down if nspin==2 - hamilt::HContainer* tmp_DMR = this->_DMR[is-1]; + int ik_begin = this->_nks * (is - 1); // jump this->_nks for spin_down if nspin==2 + hamilt::HContainer* tmp_DMR = this->_DMR[is - 1]; // set zero since this function is called in every scf step tmp_DMR->set_zero(); #ifdef _OPENMP - #pragma omp parallel for +#pragma omp parallel for #endif for (int iap = 0; iap < tmp_DMR->size_atom_pairs(); ++iap) { @@ -420,11 +426,11 @@ void DensityMatrix::cal_DMR_test() const ModuleBase::Vector3 r_index = tmp_ap.get_R_index(ir); hamilt::BaseMatrix* tmp_matrix = tmp_ap.find_matrix(r_index); #ifdef __DEBUG - if (tmp_matrix == nullptr) - { - std::cout << "tmp_matrix is nullptr" << std::endl; - continue; - } + if (tmp_matrix == nullptr) + { + std::cout << "tmp_matrix is nullptr" << std::endl; + continue; + } #endif std::complex tmp_res; // loop over k-points @@ -443,7 +449,8 @@ void DensityMatrix::cal_DMR_test() for (int j = 0; j < this->_paraV->get_col_size(iat2); ++j) { // since DMK is column-major, we need to transpose it col=>row - tmp_res = kphase * this->_DMK[ik_begin+ik][(col_ap+j)*this->_paraV->nrow+row_ap+i]; + tmp_res + = kphase * this->_DMK[ik_begin + ik][(col_ap + j) * this->_paraV->nrow + row_ap + i]; tmp_matrix->add_element(i, j, tmp_res.real()); } } @@ -463,7 +470,7 @@ void DensityMatrix, double>::cal_DMR() #ifdef __DEBUG assert(!this->_DMR.empty() && "DMR has not been initialized!"); #endif - + ModuleBase::timer::tick("DensityMatrix", "cal_DMR"); int ld_hk = this->_paraV->nrow; int ld_hk2 = 2 * ld_hk; @@ -474,7 +481,7 @@ void DensityMatrix, double>::cal_DMR() // set zero since this function is called in every scf step tmp_DMR->set_zero(); #ifdef _OPENMP - #pragma omp parallel for +#pragma omp parallel for #endif for (int i = 0; i < tmp_DMR->size_atom_pairs(); ++i) { @@ -500,52 +507,54 @@ void DensityMatrix, double>::cal_DMR() } #endif // loop over k-points - if(GlobalV::NSPIN !=4 ) - { - for (int ik = 0; ik < this->_nks; ++ik) - { - // cal k_phase - // if TK==std::complex, kphase is e^{ikR} - const ModuleBase::Vector3 dR(r_index[0], r_index[1], r_index[2]); - const double arg = (this->_kv->kvec_d[ik] * dR) * ModuleBase::TWO_PI; - double sinp, cosp; - ModuleBase::libm::sincos(arg, &sinp, &cosp); - std::complex kphase = std::complex(cosp, sinp); - // set DMR element - double* tmp_DMR_pointer = tmp_matrix->get_pointer(); - std::complex* tmp_DMK_pointer = this->_DMK[ik + ik_begin].data(); - double* DMK_real_pointer = nullptr; - double* DMK_imag_pointer = nullptr; - // jump DMK to fill DMR - // DMR is row-major, DMK is column-major - tmp_DMK_pointer += col_ap * this->_paraV->nrow + row_ap; - for (int mu = 0; mu < this->_paraV->get_row_size(iat1); ++mu) - { - DMK_real_pointer = (double*)tmp_DMK_pointer; - DMK_imag_pointer = DMK_real_pointer + 1; - BlasConnector::axpy(this->_paraV->get_col_size(iat2), - kphase.real(), - DMK_real_pointer, - ld_hk2, - tmp_DMR_pointer, - 1); - // "-" since i^2 = -1 - BlasConnector::axpy(this->_paraV->get_col_size(iat2), - -kphase.imag(), - DMK_imag_pointer, - ld_hk2, - tmp_DMR_pointer, - 1); - tmp_DMK_pointer += 1; - tmp_DMR_pointer += this->_paraV->get_col_size(iat2); - } - } - } + if (GlobalV::NSPIN != 4) + { + for (int ik = 0; ik < this->_nks; ++ik) + { + // cal k_phase + // if TK==std::complex, kphase is e^{ikR} + const ModuleBase::Vector3 dR(r_index[0], r_index[1], r_index[2]); + const double arg = (this->_kv->kvec_d[ik] * dR) * ModuleBase::TWO_PI; + double sinp, cosp; + ModuleBase::libm::sincos(arg, &sinp, &cosp); + std::complex kphase = std::complex(cosp, sinp); + // set DMR element + double* tmp_DMR_pointer = tmp_matrix->get_pointer(); + std::complex* tmp_DMK_pointer = this->_DMK[ik + ik_begin].data(); + double* DMK_real_pointer = nullptr; + double* DMK_imag_pointer = nullptr; + // jump DMK to fill DMR + // DMR is row-major, DMK is column-major + tmp_DMK_pointer += col_ap * this->_paraV->nrow + row_ap; + for (int mu = 0; mu < this->_paraV->get_row_size(iat1); ++mu) + { + DMK_real_pointer = (double*)tmp_DMK_pointer; + DMK_imag_pointer = DMK_real_pointer + 1; + BlasConnector::axpy(this->_paraV->get_col_size(iat2), + kphase.real(), + DMK_real_pointer, + ld_hk2, + tmp_DMR_pointer, + 1); + // "-" since i^2 = -1 + BlasConnector::axpy(this->_paraV->get_col_size(iat2), + -kphase.imag(), + DMK_imag_pointer, + ld_hk2, + tmp_DMR_pointer, + 1); + tmp_DMK_pointer += 1; + tmp_DMR_pointer += this->_paraV->get_col_size(iat2); + } + } + } // treat DMR as pauli matrix when NSPIN=4 - if(GlobalV::NSPIN==4) + if (GlobalV::NSPIN == 4) { - std::vector> tmp_DMR(this->_paraV->get_col_size() * this->_paraV->get_row_size(), std::complex(0.0, 0.0)); + std::vector> tmp_DMR(this->_paraV->get_col_size() + * this->_paraV->get_row_size(), + std::complex(0.0, 0.0)); for (int ik = 0; ik < this->_nks; ++ik) { // cal k_phase @@ -583,15 +592,15 @@ void DensityMatrix, double>::cal_DMR() for (int is2 = 0; is2 < npol; is2++) { step_trace[is * npol + is2] = this->_paraV->get_col_size(iat2) * is + is2; - //step_trace[is + is2 * npol] = this->_paraV->get_col_size(iat2) * is + is2; + // step_trace[is + is2 * npol] = this->_paraV->get_col_size(iat2) * is + is2; } } std::complex tmp[4]; double* target_DMR = tmp_matrix->get_pointer(); std::complex* tmp_DMR_pointer = tmp_DMR.data(); - for(int irow=0;irow_paraV->get_row_size(iat1);irow += 2) + for (int irow = 0; irow < this->_paraV->get_row_size(iat1); irow += 2) { - for(int icol=0;icol_paraV->get_col_size(iat2);icol += 2) + for (int icol = 0; icol < this->_paraV->get_col_size(iat2); icol += 2) { // catch the 4 spin component value of one orbital pair tmp[0] = tmp_DMR_pointer[icol + step_trace[0]]; @@ -602,7 +611,8 @@ void DensityMatrix, double>::cal_DMR() // save them back to the tmp_matrix target_DMR[icol + step_trace[0]] = tmp[0].real() + tmp[3].real(); target_DMR[icol + step_trace[1]] = tmp[1].real() + tmp[2].real(); - target_DMR[icol + step_trace[2]] = - tmp[1].imag() + tmp[2].imag();// (i * (rho_updown - rho_downup)).real() + target_DMR[icol + step_trace[2]] + = -tmp[1].imag() + tmp[2].imag(); // (i * (rho_updown - rho_downup)).real() target_DMR[icol + step_trace[3]] = tmp[0].real() - tmp[3].real(); } tmp_DMR_pointer += this->_paraV->get_col_size(iat2) * 2; @@ -615,6 +625,96 @@ void DensityMatrix, double>::cal_DMR() ModuleBase::timer::tick("DensityMatrix", "cal_DMR"); } +// calculate DMR from DMK using blas for multi-k calculation, without summing over k-points +template <> +void DensityMatrix, double>::cal_DMR(const int ik) +{ + ModuleBase::TITLE("DensityMatrix", "cal_DMR"); + + // To check whether DMR has been initialized +#ifdef __DEBUG + assert(!this->_DMR.empty() && "DMR has not been initialized!"); +#endif + + ModuleBase::timer::tick("DensityMatrix", "cal_DMR"); + int ld_hk = this->_paraV->nrow; + int ld_hk2 = 2 * ld_hk; + for (int is = 1; is <= this->_nspin; ++is) + { + int ik_begin = this->_nks * (is - 1); // jump this->_nks for spin_down if nspin==2 + hamilt::HContainer* tmp_DMR = this->_DMR[is - 1]; + // set zero since this function is called in every scf step + tmp_DMR->set_zero(); +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (int i = 0; i < tmp_DMR->size_atom_pairs(); ++i) + { + hamilt::AtomPair& tmp_ap = tmp_DMR->get_atom_pair(i); + int iat1 = tmp_ap.get_atom_i(); + int iat2 = tmp_ap.get_atom_j(); + // get global indexes of whole matrix for each atom in this process + int row_ap = this->_paraV->atom_begin_row[iat1]; + int col_ap = this->_paraV->atom_begin_col[iat2]; + if (row_ap == -1 || col_ap == -1) + { + throw std::string("Atom-pair not belong this process"); + } + for (int ir = 0; ir < tmp_ap.get_R_size(); ++ir) + { + const ModuleBase::Vector3 r_index = tmp_ap.get_R_index(ir); + hamilt::BaseMatrix* tmp_matrix = tmp_ap.find_matrix(r_index); +#ifdef __DEBUG + if (tmp_matrix == nullptr) + { + std::cout << "tmp_matrix is nullptr" << std::endl; + continue; + } +#endif + // Remove loop over k-points and directly use the provided ik + if (GlobalV::NSPIN != 4) + { + // cal k_phase + // if TK==std::complex, kphase is e^{ikR} + const ModuleBase::Vector3 dR(r_index[0], r_index[1], r_index[2]); + const double arg = (this->_kv->kvec_d[ik] * dR) * ModuleBase::TWO_PI; + double sinp, cosp; + ModuleBase::libm::sincos(arg, &sinp, &cosp); + std::complex kphase = std::complex(cosp, sinp); + // set DMR element + double* tmp_DMR_pointer = tmp_matrix->get_pointer(); + std::complex* tmp_DMK_pointer = this->_DMK[ik + ik_begin].data(); + double* DMK_real_pointer = nullptr; + double* DMK_imag_pointer = nullptr; + // jump DMK to fill DMR + // DMR is row-major, DMK is column-major + tmp_DMK_pointer += col_ap * this->_paraV->nrow + row_ap; + for (int mu = 0; mu < this->_paraV->get_row_size(iat1); ++mu) + { + DMK_real_pointer = (double*)tmp_DMK_pointer; + DMK_imag_pointer = DMK_real_pointer + 1; + BlasConnector::axpy(this->_paraV->get_col_size(iat2), + kphase.real(), + DMK_real_pointer, + ld_hk2, + tmp_DMR_pointer, + 1); + // "-" since i^2 = -1 + BlasConnector::axpy(this->_paraV->get_col_size(iat2), + -kphase.imag(), + DMK_imag_pointer, + ld_hk2, + tmp_DMR_pointer, + 1); + tmp_DMK_pointer += 1; + tmp_DMR_pointer += this->_paraV->get_col_size(iat2); + } + } + } + } + } + ModuleBase::timer::tick("DensityMatrix", "cal_DMR"); +} // calculate DMR from DMK using blas for gamma-only calculation template <> @@ -635,13 +735,13 @@ void DensityMatrix::cal_DMR() hamilt::HContainer* tmp_DMR = this->_DMR[is - 1]; // set zero since this function is called in every scf step tmp_DMR->set_zero(); - + #ifdef __DEBUG - //assert(tmp_DMR->is_gamma_only() == true); + // assert(tmp_DMR->is_gamma_only() == true); assert(this->_nks == 1); #endif #ifdef _OPENMP - #pragma omp parallel for +#pragma omp parallel for #endif for (int i = 0; i < tmp_DMR->size_atom_pairs(); ++i) { @@ -678,7 +778,12 @@ void DensityMatrix::cal_DMR() tmp_DMK_pointer += col_ap * this->_paraV->nrow + row_ap; for (int mu = 0; mu < this->_paraV->get_row_size(iat1); ++mu) { - BlasConnector::axpy(this->_paraV->get_col_size(iat2), kphase, tmp_DMK_pointer, ld_hk, tmp_DMR_pointer, 1); + BlasConnector::axpy(this->_paraV->get_col_size(iat2), + kphase, + tmp_DMK_pointer, + ld_hk, + tmp_DMR_pointer, + 1); tmp_DMK_pointer += 1; tmp_DMR_pointer += this->_paraV->get_col_size(iat2); } @@ -795,7 +900,9 @@ void DensityMatrix::write_DMK(const std::string directory, const for (int j = 0; j < this->_paraV->ncol; ++j) { if (j % 8 == 0) + { ofs << "\n"; + } ofs << " " << this->_DMK[ik + this->_nks * (ispin - 1)][i * this->_paraV->ncol + j]; } } @@ -830,7 +937,9 @@ void DensityMatrix, double>::write_DMK(const std::string di for (int j = 0; j < this->_paraV->ncol; ++j) { if (j % 8 == 0) + { ofs << "\n"; + } ofs << " " << this->_DMK[ik + this->_nks * (ispin - 1)][i * this->_paraV->ncol + j].real(); } } @@ -841,6 +950,6 @@ void DensityMatrix, double>::write_DMK(const std::string di // T of HContainer can be double or complex template class DensityMatrix; // Gamma-Only case template class DensityMatrix, double>; // Multi-k case -//template class DensityMatrix, std::complex>; // For EXX in future +// template class DensityMatrix, std::complex>; // For EXX in future } // namespace elecstate diff --git a/source/module_elecstate/module_dm/density_matrix.h b/source/module_elecstate/module_dm/density_matrix.h index 90f0db4ae7..0bc28c6cad 100644 --- a/source/module_elecstate/module_dm/density_matrix.h +++ b/source/module_elecstate/module_dm/density_matrix.h @@ -176,6 +176,11 @@ class DensityMatrix */ void cal_DMR(); + /** + * @brief calculate density matrix DMR from dm(k) using blas::axpy for multi-k calculation, without summing over k-points + */ + void cal_DMR(const int ik); + /** * @brief calculate density matrix DMR from dm(k) using base_matrix->add_element() */ diff --git a/source/module_io/istate_charge.cpp b/source/module_io/istate_charge.cpp index c40a25a013..23d35b96b1 100644 --- a/source/module_io/istate_charge.cpp +++ b/source/module_io/istate_charge.cpp @@ -456,59 +456,78 @@ void IState_Charge::begin(Gint_k& gk, ModuleBase::WARNING_QUIT("IState_Charge::begin", "The `pchg` calculation is only available for MPI now!"); #endif - for (int is = 0; is < nspin; ++is) + for (int ik = 0; ik < kv.get_nks() / nspin; ++ik) { - ModuleBase::GlobalFunc::ZEROS(rho[is], rhopw_nrxx); - } + const int ispin = kv.isk[ik]; + std::cout << "ispin = " << ispin << std::endl; - DM.init_DMR(GridD_in, ucell_in); - DM.cal_DMR(); + for (int is = 0; is < nspin; ++is) + { + ModuleBase::GlobalFunc::ZEROS(rho[is], rhopw_nrxx); + } - gk.initialize_pvpR(*ucell_in, GridD_in); + std::cout << "Here is Okay." << std::endl; - gk.transfer_DM2DtoGrid(DM.get_DMR_vector()); + DM.init_DMR(GridD_in, ucell_in); - Gint_inout inout((double***)nullptr, rho, Gint_Tools::job_type::rho); - gk.cal_gint(&inout); + std::cout << "Here is Okay." << std::endl; - double** rho_save = new double*[nspin]; - for (int is = 0; is < nspin; is++) - { - rho_save[is] = new double[rhopw_nrxx]; - ModuleBase::GlobalFunc::DCOPY(rho[is], rho_save[is], rhopw_nrxx); - } + DM.cal_DMR(ik); + // DM.cal_DMR(); - for (int is = 0; is < nspin; ++is) - { - // ssc should be inside the inner loop to reset the string stream each time - std::stringstream ssc; - ssc << global_out_dir << "BAND" << ib + 1 << "_SPIN" << is << "_CHG.cube"; + std::cout << "Here is Okay." << std::endl; - double ef_spin = ef_all_spin[is]; - ModuleIO::write_rho( + gk.initialize_pvpR(*ucell_in, GridD_in); + + gk.transfer_DM2DtoGrid(DM.get_DMR_vector()); + + Gint_inout inout(rho, Gint_Tools::job_type::rho); + gk.cal_gint(&inout); + + double** rho_save = new double*[nspin]; // Initialize an array of pointers + for (int is = 0; is < nspin; is++) + { + rho_save[is] = new double[rhopw_nrxx]; // Allocate memory for each internal array + ModuleBase::GlobalFunc::DCOPY(rho[is], rho_save[is], + rhopw_nrxx); // Copy data after allocation + } + + // 0 means definitely output charge density. + for (int is = 0; is < nspin; ++is) + { + // ssc should be inside the inner loop to reset the string stream each time + std::stringstream ssc; + ssc << global_out_dir << "BAND" << ib + 1 << "_K" << ik + 1 << "_SPIN" << is + 1 << "_CHG.cube"; + + // Use a const vector to store efermi for all spins, replace the original implementation: + // const double ef_tmp = pelec->eferm.get_efval(is); + double ef_spin = ef_all_spin[is]; + ModuleIO::write_rho( #ifdef __MPI - bigpw_bz, - bigpw_nbz, - rhopw_nplane, - rhopw_startz_current, + bigpw_bz, + bigpw_nbz, + rhopw_nplane, + rhopw_startz_current, #endif - rho_save[is], - is, - nspin, - 0, - ssc.str(), - rhopw_nx, - rhopw_ny, - rhopw_nz, - ef_spin, - ucell_in); - } + rho_save[is], + is, + nspin, + 0, + ssc.str(), + rhopw_nx, + rhopw_ny, + rhopw_nz, + ef_spin, + ucell_in); + } - for (int is = 0; is < nspin; is++) - { - delete[] rho_save[is]; + // Release memory of rho_save + for (int is = 0; is < nspin; is++) + { + delete[] rho_save[is]; // Release memory of each internal array + } + delete[] rho_save; // Release memory of the array of pointers } - delete[] rho_save; } } @@ -535,36 +554,30 @@ void IState_Charge::idmatrix(const int& ib, std::cout << " Perform band decomposed charge density for kpoint " << ik << ", band" << ib + 1 << std::endl; const int ispin = kv.isk[ik]; - } - for (int is = 0; is < nspin; ++is) - { + std::cout << "ispin = " << ispin << std::endl; + std::vector wg_local(this->ParaV->ncol, 0.0); const int ib_local = this->ParaV->global2local_col(ib); if (ib_local >= 0) { - wg_local[ib_local] = (ib < fermi_band) ? wg(is, ib) : wg(is, fermi_band - 1); + wg_local[ib_local] = (ib < fermi_band) ? wg(ik, ib) : wg(ik, fermi_band - 1); } - this->psi_k->fix_k(is); + this->psi_k->fix_k(ik); psi::Psi> wg_wfc(*this->psi_k, 1); for (int ir = 0; ir < wg_wfc.get_nbands(); ++ir) { BlasConnector::scal(wg_wfc.get_nbasis(), wg_local[ir], wg_wfc.get_pointer() + ir * wg_wfc.get_nbasis(), 1); } - std::cout << "Here is Okay." << std::endl; - for (int ik = 0; ik < kv.get_nks(); ++ik) - { - std::cout << "Here is Okay." << std::endl; - elecstate::psiMulPsiMpi(wg_wfc, - *(this->psi_k), - DM.get_DMK_pointer(ik), - this->ParaV->desc_wfc, - this->ParaV->desc); - std::cout << "Here is Okay." << std::endl; - } + + elecstate::psiMulPsiMpi(wg_wfc, + *(this->psi_k), + DM.get_DMK_pointer(ik), + this->ParaV->desc_wfc, + this->ParaV->desc); } } #endif From 3b19e248e6f17e89fcc8ea5bc07c1d47b0139387 Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Sun, 14 Jul 2024 20:50:40 +0800 Subject: [PATCH 05/21] Fix the issue of parameter bands_to_print not working after the refactoring of INPUT --- .../module_esolver/esolver_ks_lcao_elec.cpp | 4 ++-- source/module_io/read_input_item_general.cpp | 23 ++++++++++++------- source/module_parameter/input_parameter.h | 3 +++ 3 files changed, 20 insertions(+), 10 deletions(-) diff --git a/source/module_esolver/esolver_ks_lcao_elec.cpp b/source/module_esolver/esolver_ks_lcao_elec.cpp index a8c8892031..422a2c1293 100644 --- a/source/module_esolver/esolver_ks_lcao_elec.cpp +++ b/source/module_esolver/esolver_ks_lcao_elec.cpp @@ -456,7 +456,7 @@ void ESolver_KS_LCAO::others(const int istep) this->pw_big->nbz, GlobalV::GAMMA_ONLY_LOCAL, GlobalV::NBANDS_ISTATE, - INPUT.get_out_band_kb(), + PARAM.inp.sup.out_band_kb, GlobalV::NBANDS, GlobalV::nelec, GlobalV::NSPIN, @@ -482,7 +482,7 @@ void ESolver_KS_LCAO::others(const int istep) this->pw_big->nbz, GlobalV::GAMMA_ONLY_LOCAL, GlobalV::NBANDS_ISTATE, - INPUT.get_out_band_kb(), + PARAM.inp.sup.out_band_kb, GlobalV::NBANDS, GlobalV::nelec, GlobalV::NSPIN, diff --git a/source/module_io/read_input_item_general.cpp b/source/module_io/read_input_item_general.cpp index bcdc818bfd..53cc8e63f9 100644 --- a/source/module_io/read_input_item_general.cpp +++ b/source/module_io/read_input_item_general.cpp @@ -1,12 +1,12 @@ -#include - -#include - +#include "input_conv.h" #include "module_base/global_function.h" #include "module_base/tool_quit.h" #include "read_input.h" #include "read_input_tool.h" +#include +#include + namespace ModuleIO { // There are some examples: @@ -174,8 +174,8 @@ void ReadInput::item_general() if (!find_str(esolver_types, para.input.esolver_type)) { ModuleBase::WARNING_QUIT("ReadInput", - "esolver_type should be ksdft, sdft, " - "ofdft, tddft, lr, ks-lr, lj or dp."); + "esolver_type should be ksdft, sdft, " + "ofdft, tddft, lr, ks-lr, lj or dp."); } if (para.input.esolver_type == "dp") { @@ -277,17 +277,24 @@ void ReadInput::item_general() } { Input_Item item("nbands_istate"); - item.annotation = "number of bands around Fermi level for get_pchg calulation"; + item.annotation = "number of bands around Fermi level for get_wf and get_pchg calulation"; read_sync_int(nbands_istate); this->add_item(item); } { Input_Item item("bands_to_print"); - item.annotation = "specify the bands to be calculated in the get_pchg calculation"; + item.annotation = "specify the bands to be calculated in get_wf and get_pchg calculation"; item.read_value = [](const Input_Item& item, Parameter& para) { para.input.bands_to_print = longstring(item.str_values, item.get_size()); }; sync_string(bands_to_print); + item.reset_value = [](const Input_Item& item, Parameter& para) { + Input_Conv::parse_expression(para.input.bands_to_print, para.input.sup.out_band_kb); + para.input.sup.out_band_kb_size = para.input.sup.out_band_kb.size(); + }; + // We need to broadcast the size of out_band_kb first, then broadcast the content of out_band_kb + add_int_bcast(sup.out_band_kb_size); + add_intvec_bcast(sup.out_band_kb, para.input.sup.out_band_kb_size, 0); this->add_item(item); } { diff --git a/source/module_parameter/input_parameter.h b/source/module_parameter/input_parameter.h index 36d9cc2ead..f0115a492e 100644 --- a/source/module_parameter/input_parameter.h +++ b/source/module_parameter/input_parameter.h @@ -32,6 +32,9 @@ struct Input_supplement double uramping = -10.0 / 13.6; /// U-Ramping method (Ry) std::vector hubbard_u = {}; ///< Hubbard Coulomb interaction parameter U (Ry) std::string global_calculation = "scf"; ///< global calculation type decided by "calculation" + + std::vector out_band_kb = {}; ///< Return parsed bands_to_print as a vector of integers + int out_band_kb_size = 0; ///< Size of out_band_kb }; // It stores all input parameters both defined in INPUT file and not defined in From 6bb5b939b3258d09c142afec264f0d003374df06 Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Sun, 14 Jul 2024 21:24:46 +0800 Subject: [PATCH 06/21] Refactor istate_charge begin function into several individual functions (stage 1) --- source/module_io/istate_charge.cpp | 314 +++++++++-------------------- source/module_io/istate_charge.h | 12 ++ 2 files changed, 107 insertions(+), 219 deletions(-) diff --git a/source/module_io/istate_charge.cpp b/source/module_io/istate_charge.cpp index 23d35b96b1..3ea6d6269b 100644 --- a/source/module_io/istate_charge.cpp +++ b/source/module_io/istate_charge.cpp @@ -55,11 +55,6 @@ void IState_Charge::begin(Gint_Gamma& gg, std::cout << " Perform |psi(i)|^2 for selected bands." << std::endl; - if (!gamma_only_local) - { - ModuleBase::WARNING_QUIT("IState_Charge::begin", "Only available for GAMMA_ONLY_LOCAL now."); - } - int mode = 0; if (nbands_istate > 0 && static_cast(out_band_kb.size()) == 0) { @@ -76,129 +71,14 @@ void IState_Charge::begin(Gint_Gamma& gg, mode = 3; } - int fermi_band = 0; - int bands_below = 0; - int bands_above = 0; - - // (2) cicle: - // (2.1) calculate the selected density matrix from wave functions. - // (2.2) carry out the grid integration to get the charge density. - this->bands_picked_.resize(nbands); - ModuleBase::GlobalFunc::ZEROS(bands_picked_.data(), nbands); - - // (1) - // (1.2) read in WFC_NAO_GAMMA1.dat - std::cout << " number of electrons = " << nelec << std::endl; - - // mohan update 2011-03-21 // if ucell is odd, it's correct, // if ucell is even, it's also correct. // +1.0e-8 in case like (2.999999999+1)/2 - fermi_band = static_cast((nelec + 1) / 2 + 1.0e-8); + int fermi_band = static_cast((nelec + 1) / 2 + 1.0e-8); + std::cout << " number of electrons = " << nelec << std::endl; std::cout << " number of occupied bands = " << fermi_band << std::endl; - if (mode == 1) - { - bands_below = nbands_istate; - bands_above = nbands_istate; - - std::cout << " Plot band decomposed charge density below Fermi surface with " << bands_below << " bands." - << std::endl; - - std::cout << " Plot band decomposed charge density above Fermi surface with " << bands_above << " bands." - << std::endl; - - for (int ib = 0; ib < nbands; ++ib) - { - if (ib >= fermi_band - bands_below) - { - if (ib < fermi_band + bands_above) - { - bands_picked_[ib] = 1; - } - } - } - } - else if (mode == 2) - { - // Check if length of out_band_kb is valid - if (static_cast(out_band_kb.size()) > nbands) - { - ModuleBase::WARNING_QUIT( - "IState_Charge::begin", - "The number of bands specified by `bands_to_print` in the INPUT file exceeds `nbands`!"); - } - // Check if all elements in bands_picked_ are 0 or 1 - for (int value: out_band_kb) - { - if (value != 0 && value != 1) - { - ModuleBase::WARNING_QUIT( - "IState_Charge::begin", - "The elements of `bands_to_print` must be either 0 or 1. Invalid values found!"); - } - } - // Fill bands_picked_ with values from out_band_kb - // Remaining bands are already set to 0 - int length = std::min(static_cast(out_band_kb.size()), nbands); - for (int i = 0; i < length; ++i) - { - // out_band_kb rely on function parse_expression from input_conv.cpp - bands_picked_[i] = out_band_kb[i]; - } - - std::cout << " Plot band decomposed charge density below the Fermi surface: band "; - for (int i = 0; i + 1 <= fermi_band; ++i) - { - if (bands_picked_[i] == 1) - { - std::cout << i + 1 << " "; - } - } - std::cout << std::endl; - std::cout << " Plot band decomposed charge density above the Fermi surface: band "; - for (int i = fermi_band; i < nbands; ++i) - { - if (bands_picked_[i] == 1) - { - std::cout << i + 1 << " "; - } - } - std::cout << std::endl; - } - else if (mode == 3) - { - bool stop = false; - std::stringstream ss; - ss << global_out_dir << "istate.info"; - std::cout << " Open the file : " << ss.str() << std::endl; - if (my_rank == 0) - { - std::ifstream ifs(ss.str().c_str()); - if (!ifs) - { - stop = true; - } - else - { - // int band_index; - for (int ib = 0; ib < nbands; ++ib) - { - ModuleBase::GlobalFunc::READ_VALUE(ifs, bands_picked_[ib]); - } - } - } - -#ifdef __MPI - Parallel_Common::bcast_bool(stop); - Parallel_Common::bcast_int(bands_picked_.data(), nbands); -#endif - if (stop) - { - ofs_warning << " Can't find the file : " << ss.str() << std::endl; - ModuleBase::WARNING_QUIT("IState_Charge::begin", "can't find the istate file."); - } - } + select_bands(nbands_istate, out_band_kb, nbands, nelec, mode, fermi_band); for (int ib = 0; ib < nbands; ++ib) { @@ -206,8 +86,6 @@ void IState_Charge::begin(Gint_Gamma& gg, { std::cout << " Perform band decomposed charge density for band " << ib + 1 << std::endl; - // (1) calculate the density matrix for a partuclar band, whenever it is occupied or not. - // Using new density matrix inplementation elecstate::DensityMatrix DM(this->ParaV, nspin); @@ -217,14 +95,11 @@ void IState_Charge::begin(Gint_Gamma& gg, ModuleBase::WARNING_QUIT("IState_Charge::begin", "The `pchg` calculation is only available for MPI now!"); #endif - // (2) zero out of charge density array. for (int is = 0; is < nspin; ++is) { ModuleBase::GlobalFunc::ZEROS(rho[is], rhopw_nrxx); } - // (3) calculate charge density for a particular band. - DM.init_DMR(GridD_in, ucell_in); DM.cal_DMR(); @@ -332,82 +207,13 @@ void IState_Charge::begin(Gint_k& gk, mode = 3; } - int fermi_band = 0; - int bands_below = 0; - int bands_above = 0; - - this->bands_picked_.resize(nbands); - ModuleBase::GlobalFunc::ZEROS(bands_picked_.data(), nbands); - + int fermi_band = static_cast((nelec + 1) / 2 + 1.0e-8); std::cout << " number of electrons = " << nelec << std::endl; - fermi_band = static_cast((nelec + 1) / 2 + 1.0e-8); std::cout << " number of occupied bands = " << fermi_band << std::endl; - if (mode == 1) - { - bands_below = nbands_istate; - bands_above = nbands_istate; - - std::cout << " Plot band decomposed charge density below Fermi surface with " << bands_below << " bands." - << std::endl; - - std::cout << " Plot band decomposed charge density above Fermi surface with " << bands_above << " bands." - << std::endl; + select_bands(nbands_istate, out_band_kb, nbands, nelec, mode, fermi_band); - for (int ib = 0; ib < nbands; ++ib) - { - if (ib >= fermi_band - bands_below) - { - if (ib < fermi_band + bands_above) - { - bands_picked_[ib] = 1; - } - } - } - } - else if (mode == 2) - { - if (static_cast(out_band_kb.size()) > nbands) - { - ModuleBase::WARNING_QUIT( - "IState_Charge::begin", - "The number of bands specified by `bands_to_print` in the INPUT file exceeds `nbands`!"); - } - for (int value: out_band_kb) - { - if (value != 0 && value != 1) - { - ModuleBase::WARNING_QUIT( - "IState_Charge::begin", - "The elements of `bands_to_print` must be either 0 or 1. Invalid values found!"); - } - } - int length = std::min(static_cast(out_band_kb.size()), nbands); - for (int i = 0; i < length; ++i) - { - bands_picked_[i] = out_band_kb[i]; - } - - std::cout << " Plot band decomposed charge density below the Fermi surface: band "; - for (int i = 0; i + 1 <= fermi_band; ++i) - { - if (bands_picked_[i] == 1) - { - std::cout << i + 1 << " "; - } - } - std::cout << std::endl; - std::cout << " Plot band decomposed charge density above the Fermi surface: band "; - for (int i = fermi_band; i < nbands; ++i) - { - if (bands_picked_[i] == 1) - { - std::cout << i + 1 << " "; - } - } - std::cout << std::endl; - } - else if (mode == 3) + if (mode == 3) { bool stop = false; std::stringstream ss; @@ -458,25 +264,16 @@ void IState_Charge::begin(Gint_k& gk, for (int ik = 0; ik < kv.get_nks() / nspin; ++ik) { - const int ispin = kv.isk[ik]; - std::cout << "ispin = " << ispin << std::endl; - for (int is = 0; is < nspin; ++is) { ModuleBase::GlobalFunc::ZEROS(rho[is], rhopw_nrxx); } - std::cout << "Here is Okay." << std::endl; - DM.init_DMR(GridD_in, ucell_in); - std::cout << "Here is Okay." << std::endl; - DM.cal_DMR(ik); // DM.cal_DMR(); - std::cout << "Here is Okay." << std::endl; - gk.initialize_pvpR(*ucell_in, GridD_in); gk.transfer_DM2DtoGrid(DM.get_DMR_vector()); @@ -534,6 +331,94 @@ void IState_Charge::begin(Gint_k& gk, return; } +void IState_Charge::select_bands(const int nbands_istate, + const std::vector& out_band_kb, + const int nbands, + const double nelec, + const int mode, + const int fermi_band) +{ + int bands_below = 0; + int bands_above = 0; + + this->bands_picked_.resize(nbands); + ModuleBase::GlobalFunc::ZEROS(bands_picked_.data(), nbands); + + if (mode == 1) + { + bands_below = nbands_istate; + bands_above = nbands_istate; + + std::cout << " Plot band decomposed charge density below Fermi surface with " << bands_below << " bands." + << std::endl; + + std::cout << " Plot band decomposed charge density above Fermi surface with " << bands_above << " bands." + << std::endl; + + for (int ib = 0; ib < nbands; ++ib) + { + if (ib >= fermi_band - bands_below) + { + if (ib < fermi_band + bands_above) + { + bands_picked_[ib] = 1; + } + } + } + } + else if (mode == 2) + { + // Check if length of out_band_kb is valid + if (static_cast(out_band_kb.size()) > nbands) + { + ModuleBase::WARNING_QUIT( + "IState_Charge::select_bands", + "The number of bands specified by `bands_to_print` in the INPUT file exceeds `nbands`!"); + } + // Check if all elements in out_band_kb are 0 or 1 + for (int value: out_band_kb) + { + if (value != 0 && value != 1) + { + ModuleBase::WARNING_QUIT( + "IState_Charge::select_bands", + "The elements of `bands_to_print` must be either 0 or 1. Invalid values found!"); + } + } + // Fill bands_picked_ with values from out_band_kb + // Remaining bands are already set to 0 + int length = std::min(static_cast(out_band_kb.size()), nbands); + for (int i = 0; i < length; ++i) + { + // out_band_kb rely on function parse_expression from input_conv.cpp + bands_picked_[i] = out_band_kb[i]; + } + + std::cout << " Plot band decomposed charge density below the Fermi surface: band "; + for (int i = 0; i + 1 <= fermi_band; ++i) + { + if (bands_picked_[i] == 1) + { + std::cout << i + 1 << " "; + } + } + std::cout << std::endl; + std::cout << " Plot band decomposed charge density above the Fermi surface: band "; + for (int i = fermi_band; i < nbands; ++i) + { + if (bands_picked_[i] == 1) + { + std::cout << i + 1 << " "; + } + } + std::cout << std::endl; + } + else + { + ModuleBase::WARNING_QUIT("IState_Charge::select_bands", "Invalid mode! Please check the code."); + } +} + #ifdef __MPI void IState_Charge::idmatrix(const int& ib, const int nspin, @@ -544,8 +429,7 @@ void IState_Charge::idmatrix(const int& ib, const K_Vectors& kv) { ModuleBase::TITLE("IState_Charge", "idmatrix"); - // assert(wg.nr == nspin); - std::cout << "wg.nr = " << wg.nr << std::endl; + assert(wg.nr == kv.get_nks()); int fermi_band = static_cast((nelec + 1) / 2 + 1.0e-8); @@ -553,10 +437,6 @@ void IState_Charge::idmatrix(const int& ib, { std::cout << " Perform band decomposed charge density for kpoint " << ik << ", band" << ib + 1 << std::endl; - const int ispin = kv.isk[ik]; - - std::cout << "ispin = " << ispin << std::endl; - std::vector wg_local(this->ParaV->ncol, 0.0); const int ib_local = this->ParaV->global2local_col(ib); @@ -594,10 +474,6 @@ void IState_Charge::idmatrix(const int& ib, ModuleBase::TITLE("IState_Charge", "idmatrix"); assert(wg.nr == nspin); - std::cout << "wg.nr = " << wg.nr << std::endl; - std::cout << "wg.nc = " << wg.nc << std::endl; - std::cout << "nspin = " << nspin << std::endl; - int fermi_band = static_cast((nelec + 1) / 2 + 1.0e-8); for (int is = 0; is < nspin; ++is) diff --git a/source/module_io/istate_charge.h b/source/module_io/istate_charge.h index dffb07ba22..406a7c62c4 100644 --- a/source/module_io/istate_charge.h +++ b/source/module_io/istate_charge.h @@ -88,6 +88,18 @@ class IState_Charge private: std::vector bands_picked_; + void select_bands(const int nbands_istate, + const std::vector& out_band_kb, + const int nbands, + const double nelec, + const int mode, + const int fermi_band); + + void read_istate_file(const std::string& global_out_dir, + const int nbands, + const int my_rank, + std::ofstream& ofs_warning); + #ifdef __MPI /** * @brief Calculates the density matrix for a given band and spin. From 4711afa4b5ee3edfe6814fa975f27b4491aff7d3 Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Sun, 14 Jul 2024 22:35:35 +0800 Subject: [PATCH 07/21] Fix the issue of parameter bands_to_print not working after the refactoring of INPUT (stage 2) --- source/module_esolver/esolver_ks_lcao_elec.cpp | 4 ++-- source/module_io/bcast_globalv.cpp | 4 ++++ source/module_io/read_input_item_general.cpp | 7 ++----- source/module_parameter/system_parameter.h | 3 +++ 4 files changed, 11 insertions(+), 7 deletions(-) diff --git a/source/module_esolver/esolver_ks_lcao_elec.cpp b/source/module_esolver/esolver_ks_lcao_elec.cpp index 131b7f42d0..cb6e4e80b5 100644 --- a/source/module_esolver/esolver_ks_lcao_elec.cpp +++ b/source/module_esolver/esolver_ks_lcao_elec.cpp @@ -456,7 +456,7 @@ void ESolver_KS_LCAO::others(const int istep) this->pw_big->nbz, GlobalV::GAMMA_ONLY_LOCAL, GlobalV::NBANDS_ISTATE, - PARAM.inp.sup.out_band_kb, + PARAM.globalv.out_band_kb, GlobalV::NBANDS, GlobalV::nelec, GlobalV::NSPIN, @@ -482,7 +482,7 @@ void ESolver_KS_LCAO::others(const int istep) this->pw_big->nbz, GlobalV::GAMMA_ONLY_LOCAL, GlobalV::NBANDS_ISTATE, - PARAM.inp.sup.out_band_kb, + PARAM.globalv.out_band_kb, GlobalV::NBANDS, GlobalV::nelec, GlobalV::NSPIN, diff --git a/source/module_io/bcast_globalv.cpp b/source/module_io/bcast_globalv.cpp index 88c90e18a9..a17fa78a69 100644 --- a/source/module_io/bcast_globalv.cpp +++ b/source/module_io/bcast_globalv.cpp @@ -18,5 +18,9 @@ void ReadInput::set_globalv_bcast() add_bool_bcast(sys.double_grid); add_double_bcast(sys.uramping); add_string_bcast(sys.global_calculation); + + // We need to broadcast the size of out_band_kb first, then broadcast the content of out_band_kb + add_int_bcast(sys.out_band_kb_size); + add_intvec_bcast(sys.out_band_kb, para.sys.out_band_kb_size, 0); } } // namespace ModuleIO \ No newline at end of file diff --git a/source/module_io/read_input_item_general.cpp b/source/module_io/read_input_item_general.cpp index 19037fb6f6..766e8e19b2 100644 --- a/source/module_io/read_input_item_general.cpp +++ b/source/module_io/read_input_item_general.cpp @@ -288,12 +288,9 @@ void ReadInput::item_general() }; sync_string(input.bands_to_print); item.reset_value = [](const Input_Item& item, Parameter& para) { - Input_Conv::parse_expression(para.input.bands_to_print, para.input.sup.out_band_kb); - para.input.sup.out_band_kb_size = para.input.sup.out_band_kb.size(); + Input_Conv::parse_expression(para.input.bands_to_print, para.sys.out_band_kb); + para.sys.out_band_kb_size = para.sys.out_band_kb.size(); }; - // We need to broadcast the size of out_band_kb first, then broadcast the content of out_band_kb - add_int_bcast(sup.out_band_kb_size); - add_intvec_bcast(sup.out_band_kb, para.input.sup.out_band_kb_size, 0); this->add_item(item); } { diff --git a/source/module_parameter/system_parameter.h b/source/module_parameter/system_parameter.h index 2c40b87e05..04b0a36ec2 100644 --- a/source/module_parameter/system_parameter.h +++ b/source/module_parameter/system_parameter.h @@ -37,5 +37,8 @@ struct System_para double uramping = -10.0 / 13.6; /// U-Ramping method (Ry) std::vector hubbard_u = {}; ///< Hubbard Coulomb interaction parameter U (Ry) std::string global_calculation = "scf"; ///< global calculation type decided by "calculation" + + std::vector out_band_kb = {}; ///< return parsed bands_to_print as a vector of integers + int out_band_kb_size = 0; ///< size of out_band_kb }; #endif \ No newline at end of file From 950bd3ffd40e776d635589fdd5cfcd7f6bd49c11 Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Mon, 15 Jul 2024 13:34:33 +0800 Subject: [PATCH 08/21] Some slight modifications, and added some doxygen-style annotation for istate_charge.cpp --- source/module_io/istate_charge.cpp | 163 ++++++++++++++--------------- source/module_io/istate_charge.h | 39 +++---- 2 files changed, 99 insertions(+), 103 deletions(-) diff --git a/source/module_io/istate_charge.cpp b/source/module_io/istate_charge.cpp index 3ea6d6269b..0dd000417e 100644 --- a/source/module_io/istate_charge.cpp +++ b/source/module_io/istate_charge.cpp @@ -53,7 +53,7 @@ void IState_Charge::begin(Gint_Gamma& gg, { ModuleBase::TITLE("IState_Charge", "begin"); - std::cout << " Perform |psi(i)|^2 for selected bands." << std::endl; + std::cout << " Perform |psi(i)|^2 for selected bands (band-decomposed charge densities, gamma only)." << std::endl; int mode = 0; if (nbands_istate > 0 && static_cast(out_band_kb.size()) == 0) @@ -78,15 +78,14 @@ void IState_Charge::begin(Gint_Gamma& gg, std::cout << " number of electrons = " << nelec << std::endl; std::cout << " number of occupied bands = " << fermi_band << std::endl; + // Set this->bands_picked_ according to the mode select_bands(nbands_istate, out_band_kb, nbands, nelec, mode, fermi_band); for (int ib = 0; ib < nbands; ++ib) { if (bands_picked_[ib]) { - std::cout << " Perform band decomposed charge density for band " << ib + 1 << std::endl; - - // Using new density matrix inplementation + // Using new density matrix inplementation (gamma only) elecstate::DensityMatrix DM(this->ParaV, nspin); #ifdef __MPI @@ -100,14 +99,12 @@ void IState_Charge::begin(Gint_Gamma& gg, ModuleBase::GlobalFunc::ZEROS(rho[is], rhopw_nrxx); } + std::cout << " Performing grid integral over real space grid for band " << ib + 1 << "..." << std::endl; + DM.init_DMR(GridD_in, ucell_in); DM.cal_DMR(); - - // gg.DMRGint.resize(nspin); gg.initialize_pvpR(*ucell_in, GridD_in); - gg.transfer_DM2DtoGrid(DM.get_DMR_vector()); - Gint_inout inout((double***)nullptr, rho, Gint_Tools::job_type::rho); gg.cal_gint(&inout); @@ -121,12 +118,14 @@ void IState_Charge::begin(Gint_Gamma& gg, rhopw_nrxx); // Copy data after allocation } + std::cout << " Writting cube files..."; + // 0 means definitely output charge density. for (int is = 0; is < nspin; ++is) { // ssc should be inside the inner loop to reset the string stream each time std::stringstream ssc; - ssc << global_out_dir << "BAND" << ib + 1 << "_SPIN" << is << "_CHG.cube"; + ssc << global_out_dir << "BAND" << ib + 1 << "_GAMMA" << "_SPIN" << is + 1 << "_CHG.cube"; // Use a const vector to store efermi for all spins, replace the original implementation: // const double ef_tmp = pelec->eferm.get_efval(is); @@ -150,6 +149,8 @@ void IState_Charge::begin(Gint_Gamma& gg, ucell_in); } + std::cout << " Complete!" << std::endl; + // Release memory of rho_save for (int is = 0; is < nspin; is++) { @@ -190,7 +191,7 @@ void IState_Charge::begin(Gint_k& gk, { ModuleBase::TITLE("IState_Charge", "begin"); - std::cout << " Perform |psi(i)|^2 for selected bands (multi-k)." << std::endl; + std::cout << " Perform |psi(i)|^2 for selected bands (band-decomposed charge densities, multi-k)." << std::endl; int mode = 0; if (nbands_istate > 0 && static_cast(out_band_kb.size()) == 0) @@ -199,6 +200,7 @@ void IState_Charge::begin(Gint_k& gk, } else if (static_cast(out_band_kb.size()) > 0) { + // If out_band_kb (bands_to_print) is not empty, set mode to 2 mode = 2; std::cout << " Notice: INPUT parameter `nbands_istate` overwritten by `bands_to_print`!" << std::endl; } @@ -211,49 +213,14 @@ void IState_Charge::begin(Gint_k& gk, std::cout << " number of electrons = " << nelec << std::endl; std::cout << " number of occupied bands = " << fermi_band << std::endl; + // Set this->bands_picked_ according to the mode select_bands(nbands_istate, out_band_kb, nbands, nelec, mode, fermi_band); - if (mode == 3) - { - bool stop = false; - std::stringstream ss; - ss << global_out_dir << "istate.info"; - std::cout << " Open the file : " << ss.str() << std::endl; - if (my_rank == 0) - { - std::ifstream ifs(ss.str().c_str()); - if (!ifs) - { - stop = true; - } - else - { - for (int ib = 0; ib < nbands; ++ib) - { - ModuleBase::GlobalFunc::READ_VALUE(ifs, bands_picked_[ib]); - } - } - } - -#ifdef __MPI - Parallel_Common::bcast_bool(stop); - Parallel_Common::bcast_int(bands_picked_.data(), nbands); -#endif - if (stop) - { - ofs_warning << " Can't find the file : " << ss.str() << std::endl; - ModuleBase::WARNING_QUIT("IState_Charge::begin", "can't find the istate file."); - } - } - for (int ib = 0; ib < nbands; ++ib) { if (bands_picked_[ib]) { - // std::cout << " Perform band decomposed charge density for kpoint " << ik << ", band" << ib + 1 - // << std::endl; - - // elecstate::DensityMatrix, double> DM(this->ParaV, nspin); + // Using new density matrix inplementation (multi-k) elecstate::DensityMatrix, double> DM(&kv, this->ParaV, nspin); #ifdef __MPI @@ -262,6 +229,7 @@ void IState_Charge::begin(Gint_k& gk, ModuleBase::WARNING_QUIT("IState_Charge::begin", "The `pchg` calculation is only available for MPI now!"); #endif + // For multi-k, loop over all real k-points for (int ik = 0; ik < kv.get_nks() / nspin; ++ik) { for (int is = 0; is < nspin; ++is) @@ -269,15 +237,13 @@ void IState_Charge::begin(Gint_k& gk, ModuleBase::GlobalFunc::ZEROS(rho[is], rhopw_nrxx); } - DM.init_DMR(GridD_in, ucell_in); + std::cout << " Performing grid integral over real space grid for band " << ib + 1 << ", k-point " + << ik + 1 << "..." << std::endl; + DM.init_DMR(GridD_in, ucell_in); DM.cal_DMR(ik); - // DM.cal_DMR(); - gk.initialize_pvpR(*ucell_in, GridD_in); - gk.transfer_DM2DtoGrid(DM.get_DMR_vector()); - Gint_inout inout(rho, Gint_Tools::job_type::rho); gk.cal_gint(&inout); @@ -289,6 +255,8 @@ void IState_Charge::begin(Gint_k& gk, rhopw_nrxx); // Copy data after allocation } + std::cout << " Writting cube files..."; + // 0 means definitely output charge density. for (int is = 0; is < nspin; ++is) { @@ -296,8 +264,6 @@ void IState_Charge::begin(Gint_k& gk, std::stringstream ssc; ssc << global_out_dir << "BAND" << ib + 1 << "_K" << ik + 1 << "_SPIN" << is + 1 << "_CHG.cube"; - // Use a const vector to store efermi for all spins, replace the original implementation: - // const double ef_tmp = pelec->eferm.get_efval(is); double ef_spin = ef_all_spin[is]; ModuleIO::write_rho( #ifdef __MPI @@ -318,6 +284,8 @@ void IState_Charge::begin(Gint_k& gk, ucell_in); } + std::cout << " Complete!" << std::endl; + // Release memory of rho_save for (int is = 0; is < nspin; is++) { @@ -349,10 +317,10 @@ void IState_Charge::select_bands(const int nbands_istate, bands_below = nbands_istate; bands_above = nbands_istate; - std::cout << " Plot band decomposed charge density below Fermi surface with " << bands_below << " bands." + std::cout << " Plot band-decomposed charge density below Fermi surface with " << bands_below << " bands." << std::endl; - std::cout << " Plot band decomposed charge density above Fermi surface with " << bands_above << " bands." + std::cout << " Plot band-decomposed charge density above Fermi surface with " << bands_above << " bands." << std::endl; for (int ib = 0; ib < nbands; ++ib) @@ -394,24 +362,51 @@ void IState_Charge::select_bands(const int nbands_istate, bands_picked_[i] = out_band_kb[i]; } - std::cout << " Plot band decomposed charge density below the Fermi surface: band "; + // Check if there are selected bands below the Fermi surface + bool has_below = false; for (int i = 0; i + 1 <= fermi_band; ++i) { if (bands_picked_[i] == 1) { - std::cout << i + 1 << " "; + has_below = true; + break; } } - std::cout << std::endl; - std::cout << " Plot band decomposed charge density above the Fermi surface: band "; + if (has_below) + { + std::cout << " Plot band-decomposed charge density below the Fermi surface: band "; + for (int i = 0; i + 1 <= fermi_band; ++i) + { + if (bands_picked_[i] == 1) + { + std::cout << i + 1 << " "; + } + } + std::cout << std::endl; + } + + // Check if there are selected bands above the Fermi surface + bool has_above = false; for (int i = fermi_band; i < nbands; ++i) { if (bands_picked_[i] == 1) { - std::cout << i + 1 << " "; + has_above = true; + break; + } + } + if (has_above) + { + std::cout << " Plot band-decomposed charge density above the Fermi surface: band "; + for (int i = fermi_band; i < nbands; ++i) + { + if (bands_picked_[i] == 1) + { + std::cout << i + 1 << " "; + } } + std::cout << std::endl; } - std::cout << std::endl; } else { @@ -425,28 +420,30 @@ void IState_Charge::idmatrix(const int& ib, const double& nelec, const int nlocal, const ModuleBase::matrix& wg, - elecstate::DensityMatrix, double>& DM, + elecstate::DensityMatrix& DM, const K_Vectors& kv) { ModuleBase::TITLE("IState_Charge", "idmatrix"); - assert(wg.nr == kv.get_nks()); + assert(wg.nr == nspin); int fermi_band = static_cast((nelec + 1) / 2 + 1.0e-8); - for (int ik = 0; ik < kv.get_nks(); ++ik) + for (int is = 0; is < nspin; ++is) { - std::cout << " Perform band decomposed charge density for kpoint " << ik << ", band" << ib + 1 << std::endl; + std::cout << " Calculating density matrix for band " << ib + 1 << ", spin " << is + 1 << std::endl; std::vector wg_local(this->ParaV->ncol, 0.0); const int ib_local = this->ParaV->global2local_col(ib); if (ib_local >= 0) { - wg_local[ib_local] = (ib < fermi_band) ? wg(ik, ib) : wg(ik, fermi_band - 1); + // For unoccupied bands, use occupation of HOMO + wg_local[ib_local] = (ib < fermi_band) ? wg(is, ib) : wg(is, fermi_band - 1); } - this->psi_k->fix_k(ik); - psi::Psi> wg_wfc(*this->psi_k, 1); + // wg_wfc(ib,iw) = wg[ib] * wfc(ib,iw); + this->psi_gamma->fix_k(is); + psi::Psi wg_wfc(*this->psi_gamma, 1); for (int ir = 0; ir < wg_wfc.get_nbands(); ++ir) { @@ -454,52 +451,50 @@ void IState_Charge::idmatrix(const int& ib, } elecstate::psiMulPsiMpi(wg_wfc, - *(this->psi_k), - DM.get_DMK_pointer(ik), + *(this->psi_gamma), + DM.get_DMK_pointer(is), this->ParaV->desc_wfc, this->ParaV->desc); } } -#endif -#ifdef __MPI void IState_Charge::idmatrix(const int& ib, const int nspin, const double& nelec, const int nlocal, const ModuleBase::matrix& wg, - elecstate::DensityMatrix& DM, + elecstate::DensityMatrix, double>& DM, const K_Vectors& kv) { ModuleBase::TITLE("IState_Charge", "idmatrix"); - assert(wg.nr == nspin); + assert(wg.nr == kv.get_nks()); int fermi_band = static_cast((nelec + 1) / 2 + 1.0e-8); - for (int is = 0; is < nspin; ++is) + for (int ik = 0; ik < kv.get_nks(); ++ik) { + std::cout << " Calculating density matrix for band " << ib + 1 << ", k-point " + << ik % (kv.get_nks() / nspin) + 1 << ", spin " << kv.isk[ik] + 1 << std::endl; + std::vector wg_local(this->ParaV->ncol, 0.0); const int ib_local = this->ParaV->global2local_col(ib); if (ib_local >= 0) { - // For unoccupied bands, use occupation of HOMO - wg_local[ib_local] = (ib < fermi_band) ? wg(is, ib) : wg(is, fermi_band - 1); + wg_local[ib_local] = (ib < fermi_band) ? wg(ik, ib) : wg(ik, fermi_band - 1); } - // wg_wfc(ib,iw) = wg[ib] * wfc(ib,iw); - this->psi_gamma->fix_k(is); - psi::Psi wg_wfc(*this->psi_gamma, 1); + this->psi_k->fix_k(ik); + psi::Psi> wg_wfc(*this->psi_k, 1); for (int ir = 0; ir < wg_wfc.get_nbands(); ++ir) { BlasConnector::scal(wg_wfc.get_nbasis(), wg_local[ir], wg_wfc.get_pointer() + ir * wg_wfc.get_nbasis(), 1); } - // const int ik = 0; // Gamma point only elecstate::psiMulPsiMpi(wg_wfc, - *(this->psi_gamma), - DM.get_DMK_pointer(is), + *(this->psi_k), + DM.get_DMK_pointer(ik), this->ParaV->desc_wfc, this->ParaV->desc); } diff --git a/source/module_io/istate_charge.h b/source/module_io/istate_charge.h index 406a7c62c4..17baa1c7f0 100644 --- a/source/module_io/istate_charge.h +++ b/source/module_io/istate_charge.h @@ -13,25 +13,21 @@ #include /** - * @brief Manages the computation of the charge density for different bands. + * @brief Manages the computation of the charge densities for different bands (band-decomposed charge densities). * * This class is responsible for initializing and managing the * charge state computation process, offering functionality to - * calculate and plot the decomposed charge density below and above - * the Fermi surface based on specified bands. + * calculate and plot the decomposed charge density for specified bands. */ class IState_Charge { public: IState_Charge(psi::Psi* psi_gamma_in, const Parallel_Orbitals* ParaV_in); IState_Charge(psi::Psi>* psi_k_in, const Parallel_Orbitals* ParaV_in); - // { - // throw std::logic_error("IState_Charge for multi-k is not implemented."); - // }; ~IState_Charge(); - // for gamma_only + // for gamma only void begin(Gint_Gamma& gg, double** rho, const ModuleBase::matrix& wg, @@ -86,8 +82,16 @@ class IState_Charge const K_Vectors& kv); private: - std::vector bands_picked_; - + /** + * @brief Set this->bands_picked_ according to the mode, and process an error if the mode is not recognized. + * + * @param nbands_istate INPUT parameter nbands_istate. + * @param out_band_kb Calculated from INPUT parameter bands_to_print, vector. + * @param nbands INPUT parameter nbands. + * @param nelec Total number of electrons. + * @param mode Selected mode. + * @param fermi_band Calculated Fermi band. + */ void select_bands(const int nbands_istate, const std::vector& out_band_kb, const int nbands, @@ -95,25 +99,20 @@ class IState_Charge const int mode, const int fermi_band); - void read_istate_file(const std::string& global_out_dir, - const int nbands, - const int my_rank, - std::ofstream& ofs_warning); - #ifdef __MPI /** - * @brief Calculates the density matrix for a given band and spin. + * @brief Calculates the density matrix for a given band. * - * This method calculates the density matrix for a given band and spin using the wave function coefficients. - * It adjusts the coefficients based on the Fermi energy and performs a matrix multiplication to produce the density - * matrix. + * This method calculates the density matrix for a given band using the wave function coefficients. + * It performs a matrix multiplication to produce the density matrix. * * @param ib Band index. * @param nspin Number of spin channels. * @param nelec Total number of electrons. * @param nlocal Number of local orbitals. - * @param wg Weight matrix for bands and spins. + * @param wg Weight matrix for bands and spins (k-points). * @param DM Density matrix to be calculated. + * @param kv K-vectors. */ void idmatrix(const int& ib, const int nspin, @@ -123,6 +122,7 @@ class IState_Charge elecstate::DensityMatrix& DM, const K_Vectors& kv); + // For multi-k void idmatrix(const int& ib, const int nspin, const double& nelec, @@ -132,6 +132,7 @@ class IState_Charge const K_Vectors& kv); #endif + std::vector bands_picked_; psi::Psi* psi_gamma; psi::Psi>* psi_k; const Parallel_Orbitals* ParaV; From 744c1668c32a6102f8175504936e37a78b01ed68 Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Mon, 15 Jul 2024 14:26:50 +0800 Subject: [PATCH 09/21] Add INPUT parameter if_sep_k --- docs/advanced/input_files/input-main.md | 16 +++- .../module_esolver/esolver_ks_lcao_elec.cpp | 3 +- source/module_io/istate_charge.cpp | 89 ++++++++++++++++--- source/module_io/istate_charge.h | 23 ++--- source/module_io/read_input_item_general.cpp | 6 ++ source/module_parameter/input_parameter.h | 1 + 6 files changed, 110 insertions(+), 28 deletions(-) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 53052e87af..68f5890d98 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -158,6 +158,7 @@ - [rpa](#rpa) - [nbands\_istate](#nbands_istate) - [bands\_to\_print](#bands_to_print) + - [if\_separate\_k](#if_separate_k) - [Density of states](#density-of-states) - [dos\_edelta\_ev](#dos_edelta_ev) - [dos\_sigma](#dos_sigma) @@ -449,8 +450,8 @@ These variables are used to control general system parameters. - **relax**: perform structure relaxation calculations, the `relax_nmax` parameter depicts the maximal number of ionic iterations - **cell-relax**: perform cell relaxation calculations - **md**: perform molecular dynamics simulations - - **get_pchg**: obtain partial charge density (for LCAO basis only). See `nbands_istate` and `bands_to_print` for more information - - **get_wf**: obtain wave functions (for LCAO basis only). See `nbands_istate` for more information + - **get_pchg**: obtain partial (band-decomposed) charge densities (for LCAO basis only). See `nbands_istate` and `bands_to_print` for more information + - **get_wf**: obtain wave functions (for LCAO basis only). See `nbands_istate` and `bands_to_print` for more information - **get_S** : obtain the overlap matrix formed by localized orbitals (for LCAO basis with multiple k points). the file name is `SR.csr` with file format being the same as that generated by [out_mat_hs2](#out_mat_hs2) - **gen_bessel** : generates projectors, i.e., a series of Bessel functions, for the DeePKS method (for LCAO basis only); see also keywords `bessel_descriptor_lmax`, `bessel_descriptor_rcut` and `bessel_descriptor_tolerence`. A file named `jle.orb` will be generated which contains the projectors. An example is provided in examples/H2O-deepks-pw - **test_memory** : obtain a rough estimation of memory consuption for the calculation @@ -1755,10 +1756,17 @@ The band (KS orbital) energy for each (k-point, spin, band) will be printed in t ### bands_to_print - **Type**: String -- **Availability**: For both PW and LCAO. When `basis_type = lcao`, only used when `calculation = get_pchg`. -- **Description**: Specifies the bands to calculate the charge density for, using a space-separated string of 0s and 1s, providing a more flexible selection compared to `nbands_istate`. Each digit in the string corresponds to a band, starting from the first band. A `1` indicates that the charge density should be calculated for that band, while a `0` means the band will be ignored. The parameter allows a compact and flexible notation (similar to [`ocp_set`](#ocp_set)), for example the syntax `1 4*0 5*1 0` is used to denote the selection of bands: `1` means calculate for the first band, `4*0` skips the next four bands, `5*1` means calculate for the following five bands, and the final `0` skips the next band. It's essential that the total count of bands does not exceed the total number of bands (`nbands`); otherwise, it results in an error, and the process exits. The input string must contain only numbers and the asterisk (`*`) for repetition, ensuring correct format and intention of band selection. +- **Availability**: For both PW and LCAO. When `basis_type = lcao`, used when `calculation = get_wf` or `calculation = get_pchg`. +- **Description**: Specifies the bands to calculate the wave functions/charge densities for, using a space-separated string of 0s and 1s, providing a more flexible selection compared to `nbands_istate`. Each digit in the string corresponds to a band, starting from the first band. A `1` indicates that the charge density should be calculated for that band, while a `0` means the band will be ignored. The parameter allows a compact and flexible notation (similar to [`ocp_set`](#ocp_set)), for example the syntax `1 4*0 5*1 0` is used to denote the selection of bands: `1` means calculate for the first band, `4*0` skips the next four bands, `5*1` means calculate for the following five bands, and the final `0` skips the next band. It's essential that the total count of bands does not exceed the total number of bands (`nbands`); otherwise, it results in an error, and the process exits. The input string must contain only numbers and the asterisk (`*`) for repetition, ensuring correct format and intention of band selection. - **Default**: none +### if_separate_k + +- **Type**: Boolean +- **Availability**: Only for LCAO, used only when `calculation = get_pchg` and `gamma_only` is turned off. +- **Description**: Specifies whether to write the partial charge densities for all k-points to individual files or merge them. **Warning**: Enabling symmetry may produce incorrect results due to incorrect k-point weights. Therefore, when calculating partial charge densities, it is strongly recommended to set `symmetry = -1`. +- **Default**: false + [back to top](#full-list-of-input-keywords) ## Density of states diff --git a/source/module_esolver/esolver_ks_lcao_elec.cpp b/source/module_esolver/esolver_ks_lcao_elec.cpp index cb6e4e80b5..9eefd6d680 100644 --- a/source/module_esolver/esolver_ks_lcao_elec.cpp +++ b/source/module_esolver/esolver_ks_lcao_elec.cpp @@ -492,7 +492,8 @@ void ESolver_KS_LCAO::others(const int istep) GlobalV::ofs_warning, &GlobalC::ucell, &GlobalC::GridD, - this->kv); + this->kv, + PARAM.inp.if_separate_k); } std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "getting partial charge"); } else if (cal_type == "get_wf") { diff --git a/source/module_io/istate_charge.cpp b/source/module_io/istate_charge.cpp index 0dd000417e..6f61c908a9 100644 --- a/source/module_io/istate_charge.cpp +++ b/source/module_io/istate_charge.cpp @@ -187,7 +187,8 @@ void IState_Charge::begin(Gint_k& gk, std::ofstream& ofs_warning, const UnitCell* ucell_in, Grid_Driver* GridD_in, - const K_Vectors& kv) + const K_Vectors& kv, + const bool if_separate_k) { ModuleBase::TITLE("IState_Charge", "begin"); @@ -228,20 +229,84 @@ void IState_Charge::begin(Gint_k& gk, #else ModuleBase::WARNING_QUIT("IState_Charge::begin", "The `pchg` calculation is only available for MPI now!"); #endif - - // For multi-k, loop over all real k-points - for (int ik = 0; ik < kv.get_nks() / nspin; ++ik) + if (if_separate_k) + { + // For multi-k, loop over all real k-points + for (int ik = 0; ik < kv.get_nks() / nspin; ++ik) + { + for (int is = 0; is < nspin; ++is) + { + ModuleBase::GlobalFunc::ZEROS(rho[is], rhopw_nrxx); + } + + std::cout << " Performing grid integral over real space grid for band " << ib + 1 << ", k-point " + << ik + 1 << "..." << std::endl; + + DM.init_DMR(GridD_in, ucell_in); + DM.cal_DMR(ik); + gk.initialize_pvpR(*ucell_in, GridD_in); + gk.transfer_DM2DtoGrid(DM.get_DMR_vector()); + Gint_inout inout(rho, Gint_Tools::job_type::rho); + gk.cal_gint(&inout); + + double** rho_save = new double*[nspin]; // Initialize an array of pointers + for (int is = 0; is < nspin; is++) + { + rho_save[is] = new double[rhopw_nrxx]; // Allocate memory for each internal array + ModuleBase::GlobalFunc::DCOPY(rho[is], rho_save[is], + rhopw_nrxx); // Copy data after allocation + } + + std::cout << " Writting cube files..."; + + // 0 means definitely output charge density. + for (int is = 0; is < nspin; ++is) + { + // ssc should be inside the inner loop to reset the string stream each time + std::stringstream ssc; + ssc << global_out_dir << "BAND" << ib + 1 << "_K" << ik + 1 << "_SPIN" << is + 1 << "_CHG.cube"; + + double ef_spin = ef_all_spin[is]; + ModuleIO::write_rho( +#ifdef __MPI + bigpw_bz, + bigpw_nbz, + rhopw_nplane, + rhopw_startz_current, +#endif + rho_save[is], + is, + nspin, + 0, + ssc.str(), + rhopw_nx, + rhopw_ny, + rhopw_nz, + ef_spin, + ucell_in); + } + + std::cout << " Complete!" << std::endl; + + // Release memory of rho_save + for (int is = 0; is < nspin; is++) + { + delete[] rho_save[is]; // Release memory of each internal array + } + delete[] rho_save; // Release memory of the array of pointers + } + } + else { for (int is = 0; is < nspin; ++is) { ModuleBase::GlobalFunc::ZEROS(rho[is], rhopw_nrxx); } - std::cout << " Performing grid integral over real space grid for band " << ib + 1 << ", k-point " - << ik + 1 << "..." << std::endl; + std::cout << " Performing grid integral over real space grid for band " << ib + 1 << "..." << std::endl; DM.init_DMR(GridD_in, ucell_in); - DM.cal_DMR(ik); + DM.cal_DMR(); gk.initialize_pvpR(*ucell_in, GridD_in); gk.transfer_DM2DtoGrid(DM.get_DMR_vector()); Gint_inout inout(rho, Gint_Tools::job_type::rho); @@ -262,7 +327,7 @@ void IState_Charge::begin(Gint_k& gk, { // ssc should be inside the inner loop to reset the string stream each time std::stringstream ssc; - ssc << global_out_dir << "BAND" << ib + 1 << "_K" << ik + 1 << "_SPIN" << is + 1 << "_CHG.cube"; + ssc << global_out_dir << "BAND" << ib + 1 << "_SPIN" << is + 1 << "_CHG.cube"; double ef_spin = ef_all_spin[is]; ModuleIO::write_rho( @@ -317,10 +382,10 @@ void IState_Charge::select_bands(const int nbands_istate, bands_below = nbands_istate; bands_above = nbands_istate; - std::cout << " Plot band-decomposed charge density below Fermi surface with " << bands_below << " bands." + std::cout << " Plot band-decomposed charge densities below Fermi surface with " << bands_below << " bands." << std::endl; - std::cout << " Plot band-decomposed charge density above Fermi surface with " << bands_above << " bands." + std::cout << " Plot band-decomposed charge densities above Fermi surface with " << bands_above << " bands." << std::endl; for (int ib = 0; ib < nbands; ++ib) @@ -374,7 +439,7 @@ void IState_Charge::select_bands(const int nbands_istate, } if (has_below) { - std::cout << " Plot band-decomposed charge density below the Fermi surface: band "; + std::cout << " Plot band-decomposed charge densities below the Fermi surface: band "; for (int i = 0; i + 1 <= fermi_band; ++i) { if (bands_picked_[i] == 1) @@ -397,7 +462,7 @@ void IState_Charge::select_bands(const int nbands_istate, } if (has_above) { - std::cout << " Plot band-decomposed charge density above the Fermi surface: band "; + std::cout << " Plot band-decomposed charge densities above the Fermi surface: band "; for (int i = fermi_band; i < nbands; ++i) { if (bands_picked_[i] == 1) diff --git a/source/module_io/istate_charge.h b/source/module_io/istate_charge.h index 17baa1c7f0..78c3529a4c 100644 --- a/source/module_io/istate_charge.h +++ b/source/module_io/istate_charge.h @@ -79,19 +79,20 @@ class IState_Charge std::ofstream& ofs_warning, const UnitCell* ucell_in, Grid_Driver* GridD_in, - const K_Vectors& kv); + const K_Vectors& kv, + const bool if_separate_k); private: - /** - * @brief Set this->bands_picked_ according to the mode, and process an error if the mode is not recognized. - * - * @param nbands_istate INPUT parameter nbands_istate. - * @param out_band_kb Calculated from INPUT parameter bands_to_print, vector. - * @param nbands INPUT parameter nbands. - * @param nelec Total number of electrons. - * @param mode Selected mode. - * @param fermi_band Calculated Fermi band. - */ + /** + * @brief Set this->bands_picked_ according to the mode, and process an error if the mode is not recognized. + * + * @param nbands_istate INPUT parameter nbands_istate. + * @param out_band_kb Calculated from INPUT parameter bands_to_print, vector. + * @param nbands INPUT parameter nbands. + * @param nelec Total number of electrons. + * @param mode Selected mode. + * @param fermi_band Calculated Fermi band. + */ void select_bands(const int nbands_istate, const std::vector& out_band_kb, const int nbands, diff --git a/source/module_io/read_input_item_general.cpp b/source/module_io/read_input_item_general.cpp index 766e8e19b2..3ecc44b152 100644 --- a/source/module_io/read_input_item_general.cpp +++ b/source/module_io/read_input_item_general.cpp @@ -293,6 +293,12 @@ void ReadInput::item_general() }; this->add_item(item); } + { + Input_Item item("if_separate_k"); + item.annotation = "specify whether to write the partial charge densities for all k-points to individual files or merge them"; + read_sync_bool(input.if_separate_k); + this->add_item(item); + } { Input_Item item("symmetry"); item.annotation = "the control of symmetry"; diff --git a/source/module_parameter/input_parameter.h b/source/module_parameter/input_parameter.h index 764f3ca2f7..52ebe51d4d 100644 --- a/source/module_parameter/input_parameter.h +++ b/source/module_parameter/input_parameter.h @@ -39,6 +39,7 @@ struct Input_para int nbands_istate = 5; ///< number of bands around fermi level for get_pchg calculation. std::string bands_to_print = ""; ///< specify the bands to be calculated in the get_pchg ///< calculation, formalism similar to ocp_set. + bool if_separate_k = false; ///< whether to write partial charge for all k-points to individual files or merge them /* symmetry level: -1, no symmetry at all; 0, only basic time reversal would be considered; From 27a716110ff8c434892f086336ccc7c5582db845 Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Mon, 15 Jul 2024 15:32:06 +0800 Subject: [PATCH 10/21] And INPUT test for if_separate_k --- source/module_io/test/read_input_ptest.cpp | 1 + source/module_io/test/support/INPUT | 1 + 2 files changed, 2 insertions(+) diff --git a/source/module_io/test/read_input_ptest.cpp b/source/module_io/test/read_input_ptest.cpp index d2ed8b1ddc..665da1badc 100644 --- a/source/module_io/test/read_input_ptest.cpp +++ b/source/module_io/test/read_input_ptest.cpp @@ -50,6 +50,7 @@ TEST_F(InputParaTest, ParaRead) EXPECT_EQ(param.inp.nbands, 8); EXPECT_EQ(param.inp.nbands_sto, 256); EXPECT_EQ(param.inp.nbands_istate, 5); + EXPECT_FALSE(param.inp.if_separate_k) EXPECT_EQ(param.inp.pw_seed, 1); EXPECT_EQ(param.inp.emin_sto, 0.0); EXPECT_EQ(param.inp.emax_sto, 0.0); diff --git a/source/module_io/test/support/INPUT b/source/module_io/test/support/INPUT index 7b47adecca..482c87fb15 100644 --- a/source/module_io/test/support/INPUT +++ b/source/module_io/test/support/INPUT @@ -19,6 +19,7 @@ kspacing 0 #unit in 1/bohr, should be > 0, default is 0 wh min_dist_coef 0.2 #factor related to the allowed minimum distance between two atoms nbands 8 #number of bands nbands_istate 5 #number of bands around Fermi level for get_pchg calulation +if_separate_k false #specify whether to write the partial charge densities for all k-points to individual files or merge them symmetry 1 #the control of symmetry init_vel False #read velocity from STRU or not symmetry_prec 1e-06 #accuracy for symmetry From 1e7ca06e145c0b8e8f2f1f23b5bbf4d83ff5d9d0 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci-lite[bot]" <117423508+pre-commit-ci-lite[bot]@users.noreply.github.com> Date: Mon, 15 Jul 2024 07:41:13 +0000 Subject: [PATCH 11/21] [pre-commit.ci lite] apply automatic fixes --- source/module_io/test/read_input_ptest.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/module_io/test/read_input_ptest.cpp b/source/module_io/test/read_input_ptest.cpp index 665da1badc..65bfdee02a 100644 --- a/source/module_io/test/read_input_ptest.cpp +++ b/source/module_io/test/read_input_ptest.cpp @@ -50,7 +50,7 @@ TEST_F(InputParaTest, ParaRead) EXPECT_EQ(param.inp.nbands, 8); EXPECT_EQ(param.inp.nbands_sto, 256); EXPECT_EQ(param.inp.nbands_istate, 5); - EXPECT_FALSE(param.inp.if_separate_k) + EXPECT_FALSE(param.inp.if_separate_k); EXPECT_EQ(param.inp.pw_seed, 1); EXPECT_EQ(param.inp.emin_sto, 0.0); EXPECT_EQ(param.inp.emax_sto, 0.0); From 92a58c2b82821156facd26fabbb4fbef644c494f Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Mon, 15 Jul 2024 15:46:47 +0800 Subject: [PATCH 12/21] Explicitly instantiate the template function Input_Conv::parse_expression --- source/module_io/input_conv.cpp | 57 ++++++++++++++++++++++----------- 1 file changed, 39 insertions(+), 18 deletions(-) diff --git a/source/module_io/input_conv.cpp b/source/module_io/input_conv.cpp index 11a28193bf..f1eda9fb4c 100644 --- a/source/module_io/input_conv.cpp +++ b/source/module_io/input_conv.cpp @@ -1,7 +1,5 @@ #include "module_io/input_conv.h" -#include - #include "module_base/global_function.h" #include "module_base/global_variable.h" #include "module_cell/module_symmetry/symmetry.h" @@ -15,6 +13,8 @@ #include "module_relax/relax_old/ions_move_basic.h" #include "module_relax/relax_old/lattice_change_basic.h" +#include + #ifdef __EXX #include "module_ri/exx_abfs-jle.h" #endif @@ -109,9 +109,10 @@ void Input_Conv::parse_expression(const std::string& fn, std::vector& vec) // const std::vector::iterator dest = vec.begin() + count; // copy(ocp_temp.begin(), ocp_temp.end(), dest); // count += num; - for (size_t k = 0; k != num; k++) { + for (size_t k = 0; k != num; k++) + { vec.emplace_back(occ); -} + } } else { @@ -128,6 +129,9 @@ void Input_Conv::parse_expression(const std::string& fn, std::vector& vec) regfree(®); } +// Explicitly instantiate the template function +template void Input_Conv::parse_expression(const std::string& fn, std::vector& vec); + #ifdef __LCAO std::vector Input_Conv::convert_units(std::string params, double c) { @@ -313,7 +317,11 @@ void Input_Conv::Convert() ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "pseudo_dir", GlobalV::global_pseudo_dir); ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "orbital_dir", GlobalV::global_orbital_dir); // GlobalV::global_pseudo_type = PARAM.inp.pseudo_type; - GlobalC::ucell.setup(PARAM.inp.latname, PARAM.inp.ntype, PARAM.inp.lmaxmax, PARAM.inp.init_vel, PARAM.inp.fixed_axes); + GlobalC::ucell.setup(PARAM.inp.latname, + PARAM.inp.ntype, + PARAM.inp.lmaxmax, + PARAM.inp.init_vel, + PARAM.inp.fixed_axes); if (PARAM.inp.calculation == "relax" || PARAM.inp.calculation == "cell-relax") { @@ -564,7 +572,10 @@ void Input_Conv::Convert() if (PARAM.inp.restart_save) { std::string dft_functional_lower = PARAM.inp.dft_functional; - std::transform(PARAM.inp.dft_functional.begin(), PARAM.inp.dft_functional.end(), dft_functional_lower.begin(), tolower); + std::transform(PARAM.inp.dft_functional.begin(), + PARAM.inp.dft_functional.end(), + dft_functional_lower.begin(), + tolower); GlobalC::restart.folder = GlobalV::global_readin_dir + "restart/"; ModuleBase::GlobalFunc::MAKE_DIR(GlobalC::restart.folder); if (dft_functional_lower == "hf" || dft_functional_lower == "pbe0" || dft_functional_lower == "hse" @@ -581,7 +592,10 @@ void Input_Conv::Convert() if (PARAM.inp.restart_load) { std::string dft_functional_lower = PARAM.inp.dft_functional; - std::transform(PARAM.inp.dft_functional.begin(), PARAM.inp.dft_functional.end(), dft_functional_lower.begin(), tolower); + std::transform(PARAM.inp.dft_functional.begin(), + PARAM.inp.dft_functional.end(), + dft_functional_lower.begin(), + tolower); GlobalC::restart.folder = GlobalV::global_readin_dir + "restart/"; if (dft_functional_lower == "hf" || dft_functional_lower == "pbe0" || dft_functional_lower == "hse" || dft_functional_lower == "opt_orb" || dft_functional_lower == "scan0") @@ -602,7 +616,10 @@ void Input_Conv::Convert() #ifdef __LCAO std::string dft_functional_lower = PARAM.inp.dft_functional; - std::transform(PARAM.inp.dft_functional.begin(), PARAM.inp.dft_functional.end(), dft_functional_lower.begin(), tolower); + std::transform(PARAM.inp.dft_functional.begin(), + PARAM.inp.dft_functional.end(), + dft_functional_lower.begin(), + tolower); if (dft_functional_lower == "hf" || dft_functional_lower == "pbe0" || dft_functional_lower == "scan0") { GlobalC::exx_info.info_global.cal_exx = true; @@ -655,8 +672,8 @@ void Input_Conv::Convert() if (PARAM.inp.calculation != "nscf" && PARAM.inp.symmetry == "1") ModuleSymmetry::Symmetry::symm_flag = 0; } -#endif // __LCAO -#endif // __EXX +#endif // __LCAO +#endif // __EXX GlobalC::ppcell.cell_factor = PARAM.inp.cell_factor; // LiuXh add 20180619 //---------------------------------------------------------- @@ -778,19 +795,23 @@ void Input_Conv::Convert() { GlobalV::deepks_out_labels = true; GlobalV::deepks_scf = true; - if (GlobalV::NPROC > 1) { + if (GlobalV::NPROC > 1) + { ModuleBase::WARNING_QUIT("Input_conv", "generate deepks unittest with only 1 processor"); -} - if (GlobalV::CAL_FORCE != 1) { + } + if (GlobalV::CAL_FORCE != 1) + { ModuleBase::WARNING_QUIT("Input_conv", "force is required in generating deepks unittest"); -} - if (GlobalV::CAL_STRESS != 1) { + } + if (GlobalV::CAL_STRESS != 1) + { ModuleBase::WARNING_QUIT("Input_conv", "stress is required in generating deepks unittest"); -} + } } - if (GlobalV::deepks_scf || GlobalV::deepks_out_labels) { + if (GlobalV::deepks_scf || GlobalV::deepks_out_labels) + { GlobalV::deepks_setorb = true; -} + } #else if (PARAM.inp.deepks_scf || PARAM.inp.deepks_out_labels || PARAM.inp.deepks_bandgap || PARAM.inp.deepks_v_delta) { From ad2e258f13c90115a1ff6f8fe127034efdfe7002 Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Mon, 15 Jul 2024 15:49:26 +0800 Subject: [PATCH 13/21] Fix a ; problem --- source/module_io/test/read_input_ptest.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/module_io/test/read_input_ptest.cpp b/source/module_io/test/read_input_ptest.cpp index 665da1badc..65bfdee02a 100644 --- a/source/module_io/test/read_input_ptest.cpp +++ b/source/module_io/test/read_input_ptest.cpp @@ -50,7 +50,7 @@ TEST_F(InputParaTest, ParaRead) EXPECT_EQ(param.inp.nbands, 8); EXPECT_EQ(param.inp.nbands_sto, 256); EXPECT_EQ(param.inp.nbands_istate, 5); - EXPECT_FALSE(param.inp.if_separate_k) + EXPECT_FALSE(param.inp.if_separate_k); EXPECT_EQ(param.inp.pw_seed, 1); EXPECT_EQ(param.inp.emin_sto, 0.0); EXPECT_EQ(param.inp.emax_sto, 0.0); From d7438164ee8a0d214ac2cf475a36181d1f899d29 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci-lite[bot]" <117423508+pre-commit-ci-lite[bot]@users.noreply.github.com> Date: Mon, 15 Jul 2024 07:52:13 +0000 Subject: [PATCH 14/21] [pre-commit.ci lite] apply automatic fixes --- source/module_io/input_conv.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/source/module_io/input_conv.cpp b/source/module_io/input_conv.cpp index f1eda9fb4c..0500694711 100644 --- a/source/module_io/input_conv.cpp +++ b/source/module_io/input_conv.cpp @@ -669,8 +669,9 @@ void Input_Conv::Convert() Exx_Abfs::Jle::tolerence = PARAM.inp.exx_opt_orb_tolerence; // EXX does not support symmetry=1 - if (PARAM.inp.calculation != "nscf" && PARAM.inp.symmetry == "1") + if (PARAM.inp.calculation != "nscf" && PARAM.inp.symmetry == "1") { ModuleSymmetry::Symmetry::symm_flag = 0; +} } #endif // __LCAO #endif // __EXX From 8403566e29180d4cef656b24b3c07d103bc97a76 Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Mon, 15 Jul 2024 16:11:04 +0800 Subject: [PATCH 15/21] Explicitly instantiate the template function Input_Conv::parse_expression 2 --- source/module_io/input_conv.h | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/source/module_io/input_conv.h b/source/module_io/input_conv.h index 1aadcb55c3..77ed591294 100644 --- a/source/module_io/input_conv.h +++ b/source/module_io/input_conv.h @@ -5,13 +5,12 @@ #ifndef INPUT_CONVERT_H #define INPUT_CONVERT_H -#include -#include -#include - #include #include #include +#include +#include +#include #include #include @@ -42,6 +41,9 @@ void Convert(void); template void parse_expression(const std::string& fn, std::vector& arr); +// Explicitly instantiate the template function +extern template void parse_expression(const std::string& fn, std::vector& arr); + #ifdef __LCAO /** * @brief convert units of different parameters From 391f44328870dc8fee8d20e7ab55a3b5800e89a4 Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Mon, 15 Jul 2024 16:37:17 +0800 Subject: [PATCH 16/21] Move definition of Input_Conv::parse_expression into the header file input_conv.h --- source/module_io/input_conv.cpp | 94 ++------------------------------- source/module_io/input_conv.h | 93 ++++++++++++++++++++++++++++++-- 2 files changed, 93 insertions(+), 94 deletions(-) diff --git a/source/module_io/input_conv.cpp b/source/module_io/input_conv.cpp index 0500694711..db9880fd60 100644 --- a/source/module_io/input_conv.cpp +++ b/source/module_io/input_conv.cpp @@ -43,94 +43,9 @@ #include "module_hsolver/hsolver_pw.h" #include "module_md/md_func.h" -template -void Input_Conv::parse_expression(const std::string& fn, std::vector& vec) -{ - ModuleBase::TITLE("Input_Conv", "parse_expression"); - int count = 0; - std::string pattern("([0-9]+\\*[0-9.]+|[0-9,.]+)"); - std::vector str; - std::stringstream ss(fn); - std::string section; - while (ss >> section) - { - int index = 0; - if (str.empty()) - { - while (index < section.size() && std::isspace(section[index])) - { - index++; - } - } - section.erase(0, index); - str.push_back(section); - } - // std::string::size_type pos1, pos2; - // std::string c = " "; - // pos2 = fn.find(c); - // pos1 = 0; - // while (std::string::npos != pos2) - // { - // str.push_back(fn.substr(pos1, pos2 - pos1)); - // pos1 = pos2 + c.size(); - // pos2 = fn.find(c, pos1); - // } - // if (pos1 != fn.length()) - // { - // str.push_back(fn.substr(pos1)); - // } - regex_t reg; - regcomp(®, pattern.c_str(), REG_EXTENDED); - regmatch_t pmatch[1]; - const size_t nmatch = 1; - for (size_t i = 0; i < str.size(); ++i) - { - if (str[i] == "") - { - continue; - } - int status = regexec(®, str[i].c_str(), nmatch, pmatch, 0); - std::string sub_str = ""; - for (size_t j = pmatch[0].rm_so; j != pmatch[0].rm_eo; ++j) - { - sub_str += str[i][j]; - } - std::string sub_pattern("\\*"); - regex_t sub_reg; - regcomp(&sub_reg, sub_pattern.c_str(), REG_EXTENDED); - regmatch_t sub_pmatch[1]; - const size_t sub_nmatch = 1; - if (regexec(&sub_reg, sub_str.c_str(), sub_nmatch, sub_pmatch, 0) == 0) - { - int pos = sub_str.find("*"); - int num = stoi(sub_str.substr(0, pos)); - T occ = stof(sub_str.substr(pos + 1, sub_str.size())); - // std::vector ocp_temp(num, occ); - // const std::vector::iterator dest = vec.begin() + count; - // copy(ocp_temp.begin(), ocp_temp.end(), dest); - // count += num; - for (size_t k = 0; k != num; k++) - { - vec.emplace_back(occ); - } - } - else - { - // vec[count] = stof(sub_str); - // count += 1; - std::stringstream convert; - convert << sub_str; - T occ; - convert >> occ; - vec.emplace_back(occ); - } - regfree(&sub_reg); - } - regfree(®); -} - // Explicitly instantiate the template function -template void Input_Conv::parse_expression(const std::string& fn, std::vector& vec); +template void Input_Conv::parse_expression(const std::string& fn, std::vector& arr); +template void Input_Conv::parse_expression(const std::string& fn, std::vector& arr); #ifdef __LCAO std::vector Input_Conv::convert_units(std::string params, double c) @@ -669,9 +584,10 @@ void Input_Conv::Convert() Exx_Abfs::Jle::tolerence = PARAM.inp.exx_opt_orb_tolerence; // EXX does not support symmetry=1 - if (PARAM.inp.calculation != "nscf" && PARAM.inp.symmetry == "1") { + if (PARAM.inp.calculation != "nscf" && PARAM.inp.symmetry == "1") + { ModuleSymmetry::Symmetry::symm_flag = 0; -} + } } #endif // __LCAO #endif // __EXX diff --git a/source/module_io/input_conv.h b/source/module_io/input_conv.h index 77ed591294..45226ca64f 100644 --- a/source/module_io/input_conv.h +++ b/source/module_io/input_conv.h @@ -5,6 +5,9 @@ #ifndef INPUT_CONVERT_H #define INPUT_CONVERT_H +#include "module_base/global_function.h" +#include "module_base/global_variable.h" + #include #include #include @@ -34,15 +37,95 @@ void Convert(void); * * @tparam T * @param fn (string): expressions such as "3*1 0 2*0.5 3*0" - * @param arr (vector): stores parsing results, + * @param vec (vector): stores parsing results, * for example, "3*1 0 2*0.5 1*1.5" can be parsed as * [1, 1, 1, 0, 0.5, 0.5, 1.5] */ template -void parse_expression(const std::string& fn, std::vector& arr); - -// Explicitly instantiate the template function -extern template void parse_expression(const std::string& fn, std::vector& arr); +void parse_expression(const std::string& fn, std::vector& vec) +{ + ModuleBase::TITLE("Input_Conv", "parse_expression"); + int count = 0; + std::string pattern("([0-9]+\\*[0-9.]+|[0-9,.]+)"); + std::vector str; + std::stringstream ss(fn); + std::string section; + while (ss >> section) + { + int index = 0; + if (str.empty()) + { + while (index < section.size() && std::isspace(section[index])) + { + index++; + } + } + section.erase(0, index); + str.push_back(section); + } + // std::string::size_type pos1, pos2; + // std::string c = " "; + // pos2 = fn.find(c); + // pos1 = 0; + // while (std::string::npos != pos2) + // { + // str.push_back(fn.substr(pos1, pos2 - pos1)); + // pos1 = pos2 + c.size(); + // pos2 = fn.find(c, pos1); + // } + // if (pos1 != fn.length()) + // { + // str.push_back(fn.substr(pos1)); + // } + regex_t reg; + regcomp(®, pattern.c_str(), REG_EXTENDED); + regmatch_t pmatch[1]; + const size_t nmatch = 1; + for (size_t i = 0; i < str.size(); ++i) + { + if (str[i] == "") + { + continue; + } + int status = regexec(®, str[i].c_str(), nmatch, pmatch, 0); + std::string sub_str = ""; + for (size_t j = pmatch[0].rm_so; j != pmatch[0].rm_eo; ++j) + { + sub_str += str[i][j]; + } + std::string sub_pattern("\\*"); + regex_t sub_reg; + regcomp(&sub_reg, sub_pattern.c_str(), REG_EXTENDED); + regmatch_t sub_pmatch[1]; + const size_t sub_nmatch = 1; + if (regexec(&sub_reg, sub_str.c_str(), sub_nmatch, sub_pmatch, 0) == 0) + { + int pos = sub_str.find("*"); + int num = stoi(sub_str.substr(0, pos)); + T occ = stof(sub_str.substr(pos + 1, sub_str.size())); + // std::vector ocp_temp(num, occ); + // const std::vector::iterator dest = vec.begin() + count; + // copy(ocp_temp.begin(), ocp_temp.end(), dest); + // count += num; + for (size_t k = 0; k != num; k++) + { + vec.emplace_back(occ); + } + } + else + { + // vec[count] = stof(sub_str); + // count += 1; + std::stringstream convert; + convert << sub_str; + T occ; + convert >> occ; + vec.emplace_back(occ); + } + regfree(&sub_reg); + } + regfree(®); +} #ifdef __LCAO /** From 9fdc80ae4c30275f78745bc0fb8f8572a31ba71e Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Mon, 15 Jul 2024 18:59:34 +0800 Subject: [PATCH 17/21] Delete out_band_kb_size from system parameter due to latest INPUT refactoring --- source/module_io/bcast_globalv.cpp | 5 +---- source/module_io/input_conv.cpp | 4 ---- source/module_io/read_input_item_general.cpp | 5 +---- source/module_parameter/system_parameter.h | 1 - 4 files changed, 2 insertions(+), 13 deletions(-) diff --git a/source/module_io/bcast_globalv.cpp b/source/module_io/bcast_globalv.cpp index a17fa78a69..b7f522e238 100644 --- a/source/module_io/bcast_globalv.cpp +++ b/source/module_io/bcast_globalv.cpp @@ -18,9 +18,6 @@ void ReadInput::set_globalv_bcast() add_bool_bcast(sys.double_grid); add_double_bcast(sys.uramping); add_string_bcast(sys.global_calculation); - - // We need to broadcast the size of out_band_kb first, then broadcast the content of out_band_kb - add_int_bcast(sys.out_band_kb_size); - add_intvec_bcast(sys.out_band_kb, para.sys.out_band_kb_size, 0); + add_intvec_bcast(sys.out_band_kb, para.sys.out_band_kb.size(), 0); } } // namespace ModuleIO \ No newline at end of file diff --git a/source/module_io/input_conv.cpp b/source/module_io/input_conv.cpp index db9880fd60..f23e789efc 100644 --- a/source/module_io/input_conv.cpp +++ b/source/module_io/input_conv.cpp @@ -43,10 +43,6 @@ #include "module_hsolver/hsolver_pw.h" #include "module_md/md_func.h" -// Explicitly instantiate the template function -template void Input_Conv::parse_expression(const std::string& fn, std::vector& arr); -template void Input_Conv::parse_expression(const std::string& fn, std::vector& arr); - #ifdef __LCAO std::vector Input_Conv::convert_units(std::string params, double c) { diff --git a/source/module_io/read_input_item_general.cpp b/source/module_io/read_input_item_general.cpp index 3ecc44b152..ced6c33c22 100644 --- a/source/module_io/read_input_item_general.cpp +++ b/source/module_io/read_input_item_general.cpp @@ -285,12 +285,9 @@ void ReadInput::item_general() item.annotation = "specify the bands to be calculated in get_wf and get_pchg calculation"; item.read_value = [](const Input_Item& item, Parameter& para) { para.input.bands_to_print = longstring(item.str_values, item.get_size()); - }; - sync_string(input.bands_to_print); - item.reset_value = [](const Input_Item& item, Parameter& para) { Input_Conv::parse_expression(para.input.bands_to_print, para.sys.out_band_kb); - para.sys.out_band_kb_size = para.sys.out_band_kb.size(); }; + sync_string(input.bands_to_print); this->add_item(item); } { diff --git a/source/module_parameter/system_parameter.h b/source/module_parameter/system_parameter.h index 04b0a36ec2..9624d29bb5 100644 --- a/source/module_parameter/system_parameter.h +++ b/source/module_parameter/system_parameter.h @@ -39,6 +39,5 @@ struct System_para std::string global_calculation = "scf"; ///< global calculation type decided by "calculation" std::vector out_band_kb = {}; ///< return parsed bands_to_print as a vector of integers - int out_band_kb_size = 0; ///< size of out_band_kb }; #endif \ No newline at end of file From 8efe88ef65bbef09ad9a7fe7d2e8c644ce65490e Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Tue, 16 Jul 2024 21:04:12 +0800 Subject: [PATCH 18/21] Fix a testing problem related to get_pchg due to output file name change --- tests/integrate/312_NO_GO_wfc_istate/result.ref | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/integrate/312_NO_GO_wfc_istate/result.ref b/tests/integrate/312_NO_GO_wfc_istate/result.ref index e022a789de..66f747f6b9 100644 --- a/tests/integrate/312_NO_GO_wfc_istate/result.ref +++ b/tests/integrate/312_NO_GO_wfc_istate/result.ref @@ -1,5 +1,5 @@ -BAND1_SPIN0_CHG.cube 2.06078 -BAND2_SPIN0_CHG.cube 2.07625 -BAND3_SPIN0_CHG.cube 2.07027 -BAND4_SPIN0_CHG.cube 2.06983 +BAND1_GAMMA_SPIN1_CHG.cube 2.06078 +BAND2_GAMMA_SPIN1_CHG.cube 2.07625 +BAND3_GAMMA_SPIN1_CHG.cube 2.07027 +BAND4_GAMMA_SPIN1_CHG.cube 2.06983 totaltimeref 0.16856 From 81a341ad643c39e0053286a0d8700f4637dd75ff Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Tue, 16 Jul 2024 21:32:43 +0800 Subject: [PATCH 19/21] Replace the pointer array double** rho_save with std::vector> rho_save in istate_charge.cpp --- source/module_base/global_function.h | 7 ++++ source/module_io/istate_charge.cpp | 57 +++++++++------------------- 2 files changed, 25 insertions(+), 39 deletions(-) diff --git a/source/module_base/global_function.h b/source/module_base/global_function.h index 46b8871a17..54d7dd61ec 100644 --- a/source/module_base/global_function.h +++ b/source/module_base/global_function.h @@ -161,6 +161,13 @@ static inline void DCOPY(const T& a, T& b, const int& dim) } } +template +inline void DCOPY(const T* a, T* b, const int& dim) { + for (int i = 0; i < dim; ++i) { + b[i] = a[i]; + } +} + template inline void COPYARRAY(const T* a, T* b, const long dim); diff --git a/source/module_io/istate_charge.cpp b/source/module_io/istate_charge.cpp index 6f61c908a9..a33f0b3052 100644 --- a/source/module_io/istate_charge.cpp +++ b/source/module_io/istate_charge.cpp @@ -110,12 +110,12 @@ void IState_Charge::begin(Gint_Gamma& gg, // A solution to replace the original implementation of the following code: // pelec->charge->save_rho_before_sum_band(); - double** rho_save = new double*[nspin]; // Initialize an array of pointers - for (int is = 0; is < nspin; is++) + // Using std::vector to replace the original double** rho_save + std::vector> rho_save(nspin, std::vector(rhopw_nrxx)); + + for (int is = 0; is < nspin; ++is) { - rho_save[is] = new double[rhopw_nrxx]; // Allocate memory for each internal array - ModuleBase::GlobalFunc::DCOPY(rho[is], rho_save[is], - rhopw_nrxx); // Copy data after allocation + ModuleBase::GlobalFunc::DCOPY(rho[is], rho_save[is].data(), rhopw_nrxx); // Copy data } std::cout << " Writting cube files..."; @@ -137,7 +137,7 @@ void IState_Charge::begin(Gint_Gamma& gg, rhopw_nplane, rhopw_startz_current, #endif - rho_save[is], + rho_save[is].data(), is, nspin, 0, @@ -150,13 +150,6 @@ void IState_Charge::begin(Gint_Gamma& gg, } std::cout << " Complete!" << std::endl; - - // Release memory of rho_save - for (int is = 0; is < nspin; is++) - { - delete[] rho_save[is]; // Release memory of each internal array - } - delete[] rho_save; // Release memory of the array of pointers } } @@ -249,12 +242,12 @@ void IState_Charge::begin(Gint_k& gk, Gint_inout inout(rho, Gint_Tools::job_type::rho); gk.cal_gint(&inout); - double** rho_save = new double*[nspin]; // Initialize an array of pointers - for (int is = 0; is < nspin; is++) + // Using std::vector to replace the original double** rho_save + std::vector> rho_save(nspin, std::vector(rhopw_nrxx)); + + for (int is = 0; is < nspin; ++is) { - rho_save[is] = new double[rhopw_nrxx]; // Allocate memory for each internal array - ModuleBase::GlobalFunc::DCOPY(rho[is], rho_save[is], - rhopw_nrxx); // Copy data after allocation + ModuleBase::GlobalFunc::DCOPY(rho[is], rho_save[is].data(), rhopw_nrxx); // Copy data } std::cout << " Writting cube files..."; @@ -274,7 +267,7 @@ void IState_Charge::begin(Gint_k& gk, rhopw_nplane, rhopw_startz_current, #endif - rho_save[is], + rho_save[is].data(), is, nspin, 0, @@ -287,13 +280,6 @@ void IState_Charge::begin(Gint_k& gk, } std::cout << " Complete!" << std::endl; - - // Release memory of rho_save - for (int is = 0; is < nspin; is++) - { - delete[] rho_save[is]; // Release memory of each internal array - } - delete[] rho_save; // Release memory of the array of pointers } } else @@ -312,12 +298,12 @@ void IState_Charge::begin(Gint_k& gk, Gint_inout inout(rho, Gint_Tools::job_type::rho); gk.cal_gint(&inout); - double** rho_save = new double*[nspin]; // Initialize an array of pointers - for (int is = 0; is < nspin; is++) + // Using std::vector to replace the original double** rho_save + std::vector> rho_save(nspin, std::vector(rhopw_nrxx)); + + for (int is = 0; is < nspin; ++is) { - rho_save[is] = new double[rhopw_nrxx]; // Allocate memory for each internal array - ModuleBase::GlobalFunc::DCOPY(rho[is], rho_save[is], - rhopw_nrxx); // Copy data after allocation + ModuleBase::GlobalFunc::DCOPY(rho[is], rho_save[is].data(), rhopw_nrxx); // Copy data } std::cout << " Writting cube files..."; @@ -337,7 +323,7 @@ void IState_Charge::begin(Gint_k& gk, rhopw_nplane, rhopw_startz_current, #endif - rho_save[is], + rho_save[is].data(), is, nspin, 0, @@ -350,13 +336,6 @@ void IState_Charge::begin(Gint_k& gk, } std::cout << " Complete!" << std::endl; - - // Release memory of rho_save - for (int is = 0; is < nspin; is++) - { - delete[] rho_save[is]; // Release memory of each internal array - } - delete[] rho_save; // Release memory of the array of pointers } } } From 8a6f7325ab01d1a13a416afc8eaa85cb5827d5b4 Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Tue, 16 Jul 2024 21:35:56 +0800 Subject: [PATCH 20/21] Add the const keyword before some variables --- source/module_io/istate_charge.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/source/module_io/istate_charge.cpp b/source/module_io/istate_charge.cpp index a33f0b3052..ab7a65b175 100644 --- a/source/module_io/istate_charge.cpp +++ b/source/module_io/istate_charge.cpp @@ -74,7 +74,7 @@ void IState_Charge::begin(Gint_Gamma& gg, // if ucell is odd, it's correct, // if ucell is even, it's also correct. // +1.0e-8 in case like (2.999999999+1)/2 - int fermi_band = static_cast((nelec + 1) / 2 + 1.0e-8); + const int fermi_band = static_cast((nelec + 1) / 2 + 1.0e-8); std::cout << " number of electrons = " << nelec << std::endl; std::cout << " number of occupied bands = " << fermi_band << std::endl; @@ -203,7 +203,7 @@ void IState_Charge::begin(Gint_k& gk, mode = 3; } - int fermi_band = static_cast((nelec + 1) / 2 + 1.0e-8); + const int fermi_band = static_cast((nelec + 1) / 2 + 1.0e-8); std::cout << " number of electrons = " << nelec << std::endl; std::cout << " number of occupied bands = " << fermi_band << std::endl; @@ -399,7 +399,7 @@ void IState_Charge::select_bands(const int nbands_istate, } // Fill bands_picked_ with values from out_band_kb // Remaining bands are already set to 0 - int length = std::min(static_cast(out_band_kb.size()), nbands); + const int length = std::min(static_cast(out_band_kb.size()), nbands); for (int i = 0; i < length; ++i) { // out_band_kb rely on function parse_expression from input_conv.cpp @@ -470,7 +470,7 @@ void IState_Charge::idmatrix(const int& ib, ModuleBase::TITLE("IState_Charge", "idmatrix"); assert(wg.nr == nspin); - int fermi_band = static_cast((nelec + 1) / 2 + 1.0e-8); + const int fermi_band = static_cast((nelec + 1) / 2 + 1.0e-8); for (int is = 0; is < nspin; ++is) { @@ -513,7 +513,7 @@ void IState_Charge::idmatrix(const int& ib, ModuleBase::TITLE("IState_Charge", "idmatrix"); assert(wg.nr == kv.get_nks()); - int fermi_band = static_cast((nelec + 1) / 2 + 1.0e-8); + const int fermi_band = static_cast((nelec + 1) / 2 + 1.0e-8); for (int ik = 0; ik < kv.get_nks(); ++ik) { From 740eec6c265e3698134b0fe0292bc1ba96386cfe Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Wed, 17 Jul 2024 10:16:28 +0800 Subject: [PATCH 21/21] Did nothing, just to restart GitHub auto check --- source/module_io/istate_charge.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/source/module_io/istate_charge.cpp b/source/module_io/istate_charge.cpp index ab7a65b175..87116b8fa7 100644 --- a/source/module_io/istate_charge.cpp +++ b/source/module_io/istate_charge.cpp @@ -120,7 +120,6 @@ void IState_Charge::begin(Gint_Gamma& gg, std::cout << " Writting cube files..."; - // 0 means definitely output charge density. for (int is = 0; is < nspin; ++is) { // ssc should be inside the inner loop to reset the string stream each time @@ -252,7 +251,6 @@ void IState_Charge::begin(Gint_k& gk, std::cout << " Writting cube files..."; - // 0 means definitely output charge density. for (int is = 0; is < nspin; ++is) { // ssc should be inside the inner loop to reset the string stream each time @@ -308,7 +306,6 @@ void IState_Charge::begin(Gint_k& gk, std::cout << " Writting cube files..."; - // 0 means definitely output charge density. for (int is = 0; is < nspin; ++is) { // ssc should be inside the inner loop to reset the string stream each time