From bc0499df80ce08ea18a3c7d3556489b3b255604f Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Fri, 19 Jul 2024 19:02:56 +0800 Subject: [PATCH 1/8] Change some annotations for bands_to_print due to INPUT refactoring --- source/module_esolver/esolver_ks_pw.cpp | 6 +++--- source/module_io/istate_charge.cpp | 2 +- source/module_io/istate_envelope.cpp | 4 ++-- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index 969c8deb3e..139b9ec348 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -752,7 +752,7 @@ void ESolver_KS_PW::after_scf(const int istep) { ModuleBase::GlobalFunc::ZEROS(bands_picked.data(), this->kspw_psi->get_nbands()); - // Check if length of out_band_kb is valid + // Check if length of bands_to_print is valid if (static_cast(bands_to_print.size()) > this->kspw_psi->get_nbands()) { ModuleBase::WARNING_QUIT( @@ -771,12 +771,12 @@ void ESolver_KS_PW::after_scf(const int istep) { } } - // Fill bands_picked with values from out_band_kb, converting to int + // Fill bands_picked with values from bands_to_print // Remaining bands are already set to 0 int length = std::min(static_cast(bands_to_print.size()), this->kspw_psi->get_nbands()); for (int i = 0; i < length; ++i) { - // out_band_kb rely on function parse_expression from input_conv.cpp + // bands_to_print rely on function parse_expression // Initially designed for ocp_set, which can be double bands_picked[i] = static_cast(bands_to_print[i]); } diff --git a/source/module_io/istate_charge.cpp b/source/module_io/istate_charge.cpp index e699c61f04..49e5a20828 100644 --- a/source/module_io/istate_charge.cpp +++ b/source/module_io/istate_charge.cpp @@ -401,7 +401,7 @@ void IState_Charge::select_bands(const int nbands_istate, 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 + // out_band_kb rely on function parse_expression bands_picked_[i] = out_band_kb[i]; } diff --git a/source/module_io/istate_envelope.cpp b/source/module_io/istate_envelope.cpp index b537c20275..5d7d224d24 100644 --- a/source/module_io/istate_envelope.cpp +++ b/source/module_io/istate_envelope.cpp @@ -112,7 +112,7 @@ void IState_Envelope::begin(const psi::Psi* psid, 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 + // out_band_kb rely on function parse_expression bands_picked_[i] = out_band_kb[i]; } @@ -347,7 +347,7 @@ void IState_Envelope::begin(const psi::Psi>* psi, 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 + // out_band_kb rely on function parse_expression bands_picked_[i] = out_band_kb[i]; } From 545de0786b15033510940a36883c79466da664cd Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Fri, 26 Jul 2024 14:28:17 +0800 Subject: [PATCH 2/8] Fix the bug of pchg incorrect calculation due to symmetry=1 --- .../module_charge/symmetry_rho.cpp | 199 +++++++++++------- .../module_charge/symmetry_rho.h | 54 +++-- .../module_esolver/esolver_ks_lcao_elec.cpp | 5 +- source/module_io/istate_charge.cpp | 26 ++- source/module_io/istate_charge.h | 5 +- 5 files changed, 192 insertions(+), 97 deletions(-) diff --git a/source/module_elecstate/module_charge/symmetry_rho.cpp b/source/module_elecstate/module_charge/symmetry_rho.cpp index 224c344751..8e9bd38be4 100644 --- a/source/module_elecstate/module_charge/symmetry_rho.cpp +++ b/source/module_elecstate/module_charge/symmetry_rho.cpp @@ -1,14 +1,13 @@ #include "symmetry_rho.h" + #include "module_hamilt_general/module_xc/xc_functional.h" Symmetry_rho::Symmetry_rho() { - } Symmetry_rho::~Symmetry_rho() { - } void Symmetry_rho::begin(const int& spin_now, @@ -17,96 +16,142 @@ void Symmetry_rho::begin(const int& spin_now, Parallel_Grid& Pgrid, ModuleSymmetry::Symmetry& symm) const { - assert(spin_now < 4);//added by zhengdy-soc + assert(spin_now < 4); // added by zhengdy-soc - if(ModuleSymmetry::Symmetry::symm_flag != 1) return; - // both parallel and serial + if (ModuleSymmetry::Symmetry::symm_flag != 1) + return; + // both parallel and serial // if(symm.nrot==symm.nrotk) //pure point-group, do rho_symm in real space // { // psymm(CHR.rho[spin_now], rho_basis, Pgrid, symm); - // if(XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) psymm(CHR.kin_r[spin_now], rho_basis,Pgrid,symm); + // if(XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) psymm(CHR.kin_r[spin_now], + // rho_basis,Pgrid,symm); // } // else //space group, do rho_symm in reciprocal space - { - rho_basis->real2recip(CHR.rho[spin_now], CHR.rhog[spin_now]); - psymmg(CHR.rhog[spin_now], rho_basis, Pgrid, symm); //need to modify - rho_basis->recip2real(CHR.rhog[spin_now], CHR.rho[spin_now]); + { + rho_basis->real2recip(CHR.rho[spin_now], CHR.rhog[spin_now]); + psymmg(CHR.rhog[spin_now], rho_basis, Pgrid, symm); // need to modify + rho_basis->recip2real(CHR.rhog[spin_now], CHR.rho[spin_now]); - if(XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) - { - std::complex* kin_g = new std::complex[CHR.ngmc]; - rho_basis->real2recip(CHR.kin_r[spin_now], kin_g); - psymmg(kin_g, rho_basis,Pgrid,symm); - rho_basis->recip2real(kin_g, CHR.kin_r[spin_now]); - delete[] kin_g; - } - } - return; + if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) + { + std::complex* kin_g = new std::complex[CHR.ngmc]; + rho_basis->real2recip(CHR.kin_r[spin_now], kin_g); + psymmg(kin_g, rho_basis, Pgrid, symm); + rho_basis->recip2real(kin_g, CHR.kin_r[spin_now]); + delete[] kin_g; + } + } + return; } -void Symmetry_rho::psymm(double* rho_part, const ModulePW::PW_Basis *rho_basis, Parallel_Grid &Pgrid, ModuleSymmetry::Symmetry &symm) const +void Symmetry_rho::begin(const int& spin_now, + double** rho, + std::complex** rhog, + int ngmc, + double** kin_r, + const ModulePW::PW_Basis* rho_basis, + Parallel_Grid& Pgrid, + ModuleSymmetry::Symmetry& symm) const +{ + assert(spin_now < 4); // added by zhengdy-soc + + if (ModuleSymmetry::Symmetry::symm_flag != 1) + { + return; + } + // both parallel and serial + // if(symm.nrot==symm.nrotk) //pure point-group, do rho_symm in real space + // { + // psymm(CHR.rho[spin_now], rho_basis, Pgrid, symm); + // if(XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) psymm(CHR.kin_r[spin_now], + // rho_basis,Pgrid,symm); + // } + // else //space group, do rho_symm in reciprocal space + { + rho_basis->real2recip(rho[spin_now], rhog[spin_now]); + psymmg(rhog[spin_now], rho_basis, Pgrid, symm); + rho_basis->recip2real(rhog[spin_now], rho[spin_now]); + + if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) + { + std::complex* kin_g = new std::complex[ngmc]; + rho_basis->real2recip(kin_r[spin_now], kin_g); + psymmg(kin_g, rho_basis, Pgrid, symm); + rho_basis->recip2real(kin_g, kin_r[spin_now]); + delete[] kin_g; + } + } + return; +} + +void Symmetry_rho::psymm(double* rho_part, + const ModulePW::PW_Basis* rho_basis, + Parallel_Grid& Pgrid, + ModuleSymmetry::Symmetry& symm) const { #ifdef __MPI - // (1) reduce all rho from the first pool. - double* rhotot; - if(GlobalV::MY_RANK == 0) - { - rhotot = new double[rho_basis->nxyz]; - ModuleBase::GlobalFunc::ZEROS(rhotot, rho_basis->nxyz); - } - Pgrid.reduce_to_fullrho(rhotot, rho_part); + // (1) reduce all rho from the first pool. + double* rhotot; + if (GlobalV::MY_RANK == 0) + { + rhotot = new double[rho_basis->nxyz]; + ModuleBase::GlobalFunc::ZEROS(rhotot, rho_basis->nxyz); + } + Pgrid.reduce_to_fullrho(rhotot, rho_part); - // (2) - if(GlobalV::MY_RANK==0) - { - symm.rho_symmetry(rhotot, rho_basis->nx, rho_basis->ny, rho_basis->nz); + // (2) + if (GlobalV::MY_RANK == 0) + { + symm.rho_symmetry(rhotot, rho_basis->nx, rho_basis->ny, rho_basis->nz); #else - symm.rho_symmetry(rho_part, rho_basis->nx, rho_basis->ny, rho_basis->nz); + symm.rho_symmetry(rho_part, rho_basis->nx, rho_basis->ny, rho_basis->nz); #endif - /* - int count = 0; - GlobalV::ofs_running << scientific; - for(int iz=0; iznz; iz++) - { - GlobalV::ofs_running << "\n iz=" << iz; - for(int iy=0; iyny; iy++) - { - for(int ix=0; ixnx; ix++) - { - if(count%5==0) GlobalV::ofs_running << "\n"; - ++count; - GlobalV::ofs_running << " " << rhotot[ix*rho_basis->ny*rho_basis->nz+iy*rho_basis->nz+iz]; - } - } - } - */ - #ifdef __MPI - } + /* + int count = 0; + GlobalV::ofs_running << scientific; + for(int iz=0; iznz; iz++) + { + GlobalV::ofs_running << "\n iz=" << iz; + for(int iy=0; iyny; iy++) + { + for(int ix=0; ixnx; ix++) + { + if(count%5==0) GlobalV::ofs_running << "\n"; + ++count; + GlobalV::ofs_running << " " << rhotot[ix*rho_basis->ny*rho_basis->nz+iy*rho_basis->nz+iz]; + } + } + } + */ +#ifdef __MPI + } - // (3) - const int ncxy = rho_basis->nx * rho_basis->ny; - double* zpiece = new double[ncxy]; - for(int iz=0; iznz; iz++) - { - //GlobalV::ofs_running << "\n iz=" << iz; - ModuleBase::GlobalFunc::ZEROS(zpiece, ncxy); - if(GlobalV::MY_RANK==0) - { - for(int ix=0; ixnx; ix++) - { - for(int iy=0; iyny; iy++) - { - const int ir = ix * rho_basis->ny + iy; - zpiece[ir] = rhotot[ix * rho_basis->ny * rho_basis->nz + iy * rho_basis->nz + iz]; - //rho[ir*nczp+znow] = zpiece[ir]; - } - } - } - Pgrid.zpiece_to_all(zpiece,iz, rho_part); - } + // (3) + const int ncxy = rho_basis->nx * rho_basis->ny; + double* zpiece = new double[ncxy]; + for (int iz = 0; iz < rho_basis->nz; iz++) + { + // GlobalV::ofs_running << "\n iz=" << iz; + ModuleBase::GlobalFunc::ZEROS(zpiece, ncxy); + if (GlobalV::MY_RANK == 0) + { + for (int ix = 0; ix < rho_basis->nx; ix++) + { + for (int iy = 0; iy < rho_basis->ny; iy++) + { + const int ir = ix * rho_basis->ny + iy; + zpiece[ir] = rhotot[ix * rho_basis->ny * rho_basis->nz + iy * rho_basis->nz + iz]; + // rho[ir*nczp+znow] = zpiece[ir]; + } + } + } + Pgrid.zpiece_to_all(zpiece, iz, rho_part); + } - if(GlobalV::MY_RANK==0) delete[] rhotot; - delete[] zpiece; + if (GlobalV::MY_RANK == 0) + delete[] rhotot; + delete[] zpiece; #endif - return; + return; } diff --git a/source/module_elecstate/module_charge/symmetry_rho.h b/source/module_elecstate/module_charge/symmetry_rho.h index 3e8e01a7ee..39ef442687 100644 --- a/source/module_elecstate/module_charge/symmetry_rho.h +++ b/source/module_elecstate/module_charge/symmetry_rho.h @@ -1,16 +1,15 @@ #ifndef SYMMETRY_RHO_H #define SYMMETRY_RHO_H -#include "module_elecstate/module_charge/charge.h" #include "module_basis/module_pw/pw_basis.h" -#include "module_hamilt_pw/hamilt_pwdft/parallel_grid.h" - #include "module_cell/module_symmetry/symmetry.h" +#include "module_elecstate/module_charge/charge.h" +#include "module_hamilt_pw/hamilt_pwdft/parallel_grid.h" class Symmetry_rho { - public: - Symmetry_rho(); - ~Symmetry_rho(); + public: + Symmetry_rho(); + ~Symmetry_rho(); void begin(const int& spin_now, const Charge& CHR, @@ -18,20 +17,41 @@ class Symmetry_rho Parallel_Grid& Pgrid, ModuleSymmetry::Symmetry& symm) const; + void begin(const int& spin_now, + double** rho, + std::complex** rhog, + int ngmc, + double** kin_r, + const ModulePW::PW_Basis* pw, + Parallel_Grid& Pgrid, + ModuleSymmetry::Symmetry& symm) const; + private: - // in real space: - void psymm(double *rho_part, const ModulePW::PW_Basis *pw, Parallel_Grid &Pgrid, ModuleSymmetry::Symmetry &symm) const; - // in reciprocal space: - void psymmg(std::complex* rhog_part, const ModulePW::PW_Basis *rho_basis, - Parallel_Grid &Pgrid, ModuleSymmetry::Symmetry &symm) const; + // in real space: + void psymm(double* rho_part, + const ModulePW::PW_Basis* pw, + Parallel_Grid& Pgrid, + ModuleSymmetry::Symmetry& symm) const; + // in reciprocal space: + void psymmg(std::complex* rhog_part, + const ModulePW::PW_Basis* rho_basis, + Parallel_Grid& Pgrid, + ModuleSymmetry::Symmetry& symm) const; #ifdef __MPI - void reduce_to_fullrhog(const ModulePW::PW_Basis *rho_basis, std::complex* rhogtot, - std::complex* rhogin, int* ig2isztot, const int* ig2iszin, int max_npw) const; - void rhog_piece_to_all(const ModulePW::PW_Basis *rho_basis, - std::complex* rhogtot, std::complex* rhog_part) const; + void reduce_to_fullrhog(const ModulePW::PW_Basis* rho_basis, + std::complex* rhogtot, + std::complex* rhogin, + int* ig2isztot, + const int* ig2iszin, + int max_npw) const; + void rhog_piece_to_all(const ModulePW::PW_Basis* rho_basis, + std::complex* rhogtot, + std::complex* rhog_part) const; #endif - void get_ixyz2ipw(const ModulePW::PW_Basis *rho_basis, const int* ig2isztot, - const int* fftixy2is, int* ixyz2ipw) const; //(ix, iy, iz) -> (ip, ig) + void get_ixyz2ipw(const ModulePW::PW_Basis* rho_basis, + const int* ig2isztot, + const int* fftixy2is, + int* ixyz2ipw) const; //(ix, iy, iz) -> (ip, ig) }; #endif diff --git a/source/module_esolver/esolver_ks_lcao_elec.cpp b/source/module_esolver/esolver_ks_lcao_elec.cpp index 6c93d95831..4a3b8e7e6a 100644 --- a/source/module_esolver/esolver_ks_lcao_elec.cpp +++ b/source/module_esolver/esolver_ks_lcao_elec.cpp @@ -479,8 +479,10 @@ void ESolver_KS_LCAO::others(const int istep) } else { ISC.begin(this->GK, this->pelec->charge->rho, + this->pelec->charge->rhog, this->pelec->wg, this->pelec->eferm.get_all_ef(), + this->pw_rho, this->pw_rho->nrxx, this->pw_rho->nplane, this->pw_rho->startz_current, @@ -502,7 +504,8 @@ void ESolver_KS_LCAO::others(const int istep) &GlobalC::ucell, &GlobalC::GridD, this->kv, - PARAM.inp.if_separate_k); + PARAM.inp.if_separate_k, + this->pelec->charge->ngmc); } 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 49e5a20828..f0ae136782 100644 --- a/source/module_io/istate_charge.cpp +++ b/source/module_io/istate_charge.cpp @@ -5,6 +5,7 @@ #include "module_base/global_variable.h" #include "module_base/parallel_common.h" #include "module_base/scalapack_connector.h" +#include "module_elecstate/module_charge/symmetry_rho.h" #include "module_elecstate/module_dm/cal_dm_psi.h" #include "module_elecstate/module_dm/density_matrix.h" #include "module_hamilt_lcao/module_gint/gint.h" @@ -157,8 +158,10 @@ void IState_Charge::begin(Gint_Gamma& gg, void IState_Charge::begin(Gint_k& gk, double** rho, + std::complex** rhog, const ModuleBase::matrix& wg, const std::vector& ef_all_spin, + const ModulePW::PW_Basis* rho_pw, const int rhopw_nrxx, const int rhopw_nplane, const int rhopw_startz_current, @@ -180,7 +183,8 @@ void IState_Charge::begin(Gint_k& gk, const UnitCell* ucell_in, Grid_Driver* GridD_in, const K_Vectors& kv, - const bool if_separate_k) + const bool if_separate_k, + const int ngmc) { ModuleBase::TITLE("IState_Charge", "begin"); @@ -304,6 +308,26 @@ void IState_Charge::begin(Gint_k& gk, ModuleBase::GlobalFunc::DCOPY(rho[is], rho_save[is].data(), rhopw_nrxx); // Copy data } + // Symmetrize the charge density + std::cout << " Symmetrizing charge density..." << std::endl; + Symmetry_rho srho; + for (int is = 0; is < nspin; ++is) + { + std::vector rho_save_pointers(nspin); + for (int i = 0; i < nspin; ++i) + { + rho_save_pointers[i] = rho_save[i].data(); + } + srho.begin(is, + rho_save_pointers.data(), + rhog, + ngmc, + nullptr, + rho_pw, + GlobalC::Pgrid, + GlobalC::ucell.symm); + } + std::cout << " Writting cube files..."; for (int is = 0; is < nspin; ++is) diff --git a/source/module_io/istate_charge.h b/source/module_io/istate_charge.h index 149b736e94..7f721eb605 100644 --- a/source/module_io/istate_charge.h +++ b/source/module_io/istate_charge.h @@ -57,8 +57,10 @@ class IState_Charge // For multi-k void begin(Gint_k& gk, double** rho, + std::complex** rhog, const ModuleBase::matrix& wg, const std::vector& ef_all_spin, + const ModulePW::PW_Basis* rho_pw, const int rhopw_nrxx, const int rhopw_nplane, const int rhopw_startz_current, @@ -80,7 +82,8 @@ class IState_Charge const UnitCell* ucell_in, Grid_Driver* GridD_in, const K_Vectors& kv, - const bool if_separate_k); + const bool if_separate_k, + const int ngmc); private: /** From 7faa871c015f4763b86d24ee2fa9287e29081c4e Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Fri, 26 Jul 2024 15:40:41 +0800 Subject: [PATCH 3/8] Get rid of GlobalC in istate_charge when symmetrizing charge density --- source/module_esolver/esolver_ks_lcao_elec.cpp | 1 + source/module_io/istate_charge.cpp | 13 +++++++++---- source/module_io/istate_charge.h | 4 +++- 3 files changed, 13 insertions(+), 5 deletions(-) diff --git a/source/module_esolver/esolver_ks_lcao_elec.cpp b/source/module_esolver/esolver_ks_lcao_elec.cpp index 4a3b8e7e6a..4614fb6d4b 100644 --- a/source/module_esolver/esolver_ks_lcao_elec.cpp +++ b/source/module_esolver/esolver_ks_lcao_elec.cpp @@ -505,6 +505,7 @@ void ESolver_KS_LCAO::others(const int istep) &GlobalC::GridD, this->kv, PARAM.inp.if_separate_k, + &GlobalC::Pgrid, this->pelec->charge->ngmc); } std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "getting partial charge"); diff --git a/source/module_io/istate_charge.cpp b/source/module_io/istate_charge.cpp index f0ae136782..878543f6f2 100644 --- a/source/module_io/istate_charge.cpp +++ b/source/module_io/istate_charge.cpp @@ -26,6 +26,7 @@ IState_Charge::~IState_Charge() { } +// for gamma only void IState_Charge::begin(Gint_Gamma& gg, double** rho, const ModuleBase::matrix& wg, @@ -156,6 +157,7 @@ void IState_Charge::begin(Gint_Gamma& gg, return; } +// For multi-k void IState_Charge::begin(Gint_k& gk, double** rho, std::complex** rhog, @@ -180,10 +182,11 @@ void IState_Charge::begin(Gint_k& gk, const std::string& global_out_dir, const int my_rank, std::ofstream& ofs_warning, - const UnitCell* ucell_in, + UnitCell* ucell_in, Grid_Driver* GridD_in, const K_Vectors& kv, const bool if_separate_k, + Parallel_Grid* Pgrid, const int ngmc) { ModuleBase::TITLE("IState_Charge", "begin"); @@ -308,7 +311,7 @@ void IState_Charge::begin(Gint_k& gk, ModuleBase::GlobalFunc::DCOPY(rho[is], rho_save[is].data(), rhopw_nrxx); // Copy data } - // Symmetrize the charge density + // Symmetrize the charge density, otherwise the results are incorrect if the symmetry is on std::cout << " Symmetrizing charge density..." << std::endl; Symmetry_rho srho; for (int is = 0; is < nspin; ++is) @@ -324,8 +327,8 @@ void IState_Charge::begin(Gint_k& gk, ngmc, nullptr, rho_pw, - GlobalC::Pgrid, - GlobalC::ucell.symm); + *Pgrid, + ucell_in->symm); } std::cout << " Writting cube files..."; @@ -482,6 +485,7 @@ void IState_Charge::select_bands(const int nbands_istate, } #ifdef __MPI +// for gamma only void IState_Charge::idmatrix(const int& ib, const int nspin, const double& nelec, @@ -525,6 +529,7 @@ void IState_Charge::idmatrix(const int& ib, } } +// For multi-k void IState_Charge::idmatrix(const int& ib, const int nspin, const double& nelec, diff --git a/source/module_io/istate_charge.h b/source/module_io/istate_charge.h index 7f721eb605..26e15b2579 100644 --- a/source/module_io/istate_charge.h +++ b/source/module_io/istate_charge.h @@ -5,6 +5,7 @@ #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_hamilt_pw/hamilt_pwdft/parallel_grid.h" #include "module_psi/psi.h" #include @@ -79,10 +80,11 @@ class IState_Charge const std::string& global_out_dir, const int my_rank, std::ofstream& ofs_warning, - const UnitCell* ucell_in, + UnitCell* ucell_in, Grid_Driver* GridD_in, const K_Vectors& kv, const bool if_separate_k, + Parallel_Grid* Pgrid, const int ngmc); private: From 72b89e6c54ad8b37e0794e6b13a9ca84ab562225 Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Fri, 26 Jul 2024 16:27:18 +0800 Subject: [PATCH 4/8] Fix the bug of pchg incorrect calculation under PW due to symmetry=1 --- source/module_esolver/esolver_ks_pw.cpp | 40 ++++++++++++++++++++++--- source/module_io/istate_charge.cpp | 3 +- 2 files changed, 38 insertions(+), 5 deletions(-) diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index e33414dfef..982e745553 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -795,7 +795,8 @@ void ESolver_KS_PW::after_scf(const int istep) { = new std::complex[this->pw_rho->nxyz]; double* rho_band = new double[this->pw_rho->nxyz]; - for (int ib = 0; ib < this->kspw_psi->get_nbands(); ++ib) { + for (int ib = 0; ib < this->kspw_psi->get_nbands(); ++ib) + { // Skip the loop iteration if bands_picked[ib] is 0 if (!bands_picked[ib]) { continue; @@ -816,14 +817,45 @@ void ESolver_KS_PW::after_scf(const int istep) { double w1 = static_cast(this->kv.wk[ik] / GlobalC::ucell.omega); - for (int i = 0; i < this->pw_rho->nxyz; i++) { + for (int i = 0; i < this->pw_rho->nxyz; i++) + { rho_band[i] += std::norm(wfcr[i]) * w1; } } + // Symmetrize the charge density, otherwise the results are incorrect if the symmetry is on + std::cout << " Symmetrizing band-decomposed charge density..." << std::endl; + Symmetry_rho srho; + for (int is = 0; is < GlobalV::NSPIN; is++) + { + // Create pointers to the data structures + double** rho_save_pointers = new double*[GlobalV::NSPIN]; + std::complex** rhog = new std::complex*[GlobalV::NSPIN]; + + for (int s = 0; s < GlobalV::NSPIN; s++) { + rho_save_pointers[s] = rho_band; + rhog[s] = new std::complex[this->pelec->charge->ngmc]; + } + + srho.begin(is, + rho_save_pointers, + rhog, + this->pelec->charge->ngmc, + nullptr, + this->pw_rhod, + GlobalC::Pgrid, + GlobalC::ucell.symm); + + for (int s = 0; s < GlobalV::NSPIN; s++) + { + delete[] rhog[s]; + } + delete[] rho_save_pointers; + delete[] rhog; + } + std::stringstream ssc; - ssc << GlobalV::global_out_dir << "band" << ib + 1 - << ".cube"; // band index starts from 1 + ssc << GlobalV::global_out_dir << "BAND" << ib + 1 << "_CHG.cube"; // band index starts from 1 ModuleIO::write_rho( #ifdef __MPI diff --git a/source/module_io/istate_charge.cpp b/source/module_io/istate_charge.cpp index 878543f6f2..cb33630431 100644 --- a/source/module_io/istate_charge.cpp +++ b/source/module_io/istate_charge.cpp @@ -228,6 +228,7 @@ void IState_Charge::begin(Gint_k& gk, #else ModuleBase::WARNING_QUIT("IState_Charge::begin", "The `pchg` calculation is only available for MPI now!"); #endif + // If contribution from different k-points need to be output separately if (if_separate_k) { // For multi-k, loop over all real k-points @@ -312,7 +313,7 @@ void IState_Charge::begin(Gint_k& gk, } // Symmetrize the charge density, otherwise the results are incorrect if the symmetry is on - std::cout << " Symmetrizing charge density..." << std::endl; + std::cout << " Symmetrizing band-decomposed charge density..." << std::endl; Symmetry_rho srho; for (int is = 0; is < nspin; ++is) { From b7140174e2aa7a8b081ba733c5421917882dbeef Mon Sep 17 00:00:00 2001 From: "pre-commit-ci-lite[bot]" <117423508+pre-commit-ci-lite[bot]@users.noreply.github.com> Date: Fri, 26 Jul 2024 10:37:46 +0000 Subject: [PATCH 5/8] [pre-commit.ci lite] apply automatic fixes --- source/module_elecstate/module_charge/symmetry_rho.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/source/module_elecstate/module_charge/symmetry_rho.cpp b/source/module_elecstate/module_charge/symmetry_rho.cpp index 8e9bd38be4..fada4a1749 100644 --- a/source/module_elecstate/module_charge/symmetry_rho.cpp +++ b/source/module_elecstate/module_charge/symmetry_rho.cpp @@ -18,8 +18,9 @@ void Symmetry_rho::begin(const int& spin_now, { assert(spin_now < 4); // added by zhengdy-soc - if (ModuleSymmetry::Symmetry::symm_flag != 1) + if (ModuleSymmetry::Symmetry::symm_flag != 1) { return; +} // both parallel and serial // if(symm.nrot==symm.nrotk) //pure point-group, do rho_symm in real space // { @@ -149,8 +150,9 @@ void Symmetry_rho::psymm(double* rho_part, Pgrid.zpiece_to_all(zpiece, iz, rho_part); } - if (GlobalV::MY_RANK == 0) + if (GlobalV::MY_RANK == 0) { delete[] rhotot; +} delete[] zpiece; #endif return; From b09d96c433eafee334b3236877cb04c8081bee27 Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Fri, 26 Jul 2024 18:45:05 +0800 Subject: [PATCH 6/8] Refactor symmetry_rho to replace raw pointer arrays with std::vector implementations --- .../module_charge/symmetry_rho.cpp | 44 ++++++++----------- 1 file changed, 19 insertions(+), 25 deletions(-) diff --git a/source/module_elecstate/module_charge/symmetry_rho.cpp b/source/module_elecstate/module_charge/symmetry_rho.cpp index 8e9bd38be4..24b2d716db 100644 --- a/source/module_elecstate/module_charge/symmetry_rho.cpp +++ b/source/module_elecstate/module_charge/symmetry_rho.cpp @@ -35,11 +35,11 @@ void Symmetry_rho::begin(const int& spin_now, if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) { - std::complex* kin_g = new std::complex[CHR.ngmc]; - rho_basis->real2recip(CHR.kin_r[spin_now], kin_g); - psymmg(kin_g, rho_basis, Pgrid, symm); - rho_basis->recip2real(kin_g, CHR.kin_r[spin_now]); - delete[] kin_g; + // Use std::vector to manage kin_g instead of raw pointer + std::vector> kin_g(CHR.ngmc); + rho_basis->real2recip(CHR.kin_r[spin_now], kin_g.data()); + psymmg(kin_g.data(), rho_basis, Pgrid, symm); + rho_basis->recip2real(kin_g.data(), CHR.kin_r[spin_now]); } } return; @@ -75,11 +75,11 @@ void Symmetry_rho::begin(const int& spin_now, if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) { - std::complex* kin_g = new std::complex[ngmc]; - rho_basis->real2recip(kin_r[spin_now], kin_g); - psymmg(kin_g, rho_basis, Pgrid, symm); - rho_basis->recip2real(kin_g, kin_r[spin_now]); - delete[] kin_g; + // Use std::vector to manage kin_g instead of raw pointer + std::vector> kin_g(ngmc); + rho_basis->real2recip(kin_r[spin_now], kin_g.data()); + psymmg(kin_g.data(), rho_basis, Pgrid, symm); + rho_basis->recip2real(kin_g.data(), kin_r[spin_now]); } } return; @@ -92,18 +92,18 @@ void Symmetry_rho::psymm(double* rho_part, { #ifdef __MPI // (1) reduce all rho from the first pool. - double* rhotot; + std::vector rhotot; if (GlobalV::MY_RANK == 0) { - rhotot = new double[rho_basis->nxyz]; - ModuleBase::GlobalFunc::ZEROS(rhotot, rho_basis->nxyz); + rhotot.resize(rho_basis->nxyz); + ModuleBase::GlobalFunc::ZEROS(rhotot.data(), rho_basis->nxyz); } - Pgrid.reduce_to_fullrho(rhotot, rho_part); + Pgrid.reduce_to_fullrho(rhotot.data(), rho_part); // (2) if (GlobalV::MY_RANK == 0) { - symm.rho_symmetry(rhotot, rho_basis->nx, rho_basis->ny, rho_basis->nz); + symm.rho_symmetry(rhotot.data(), rho_basis->nx, rho_basis->ny, rho_basis->nz); #else symm.rho_symmetry(rho_part, rho_basis->nx, rho_basis->ny, rho_basis->nz); #endif @@ -129,11 +129,10 @@ void Symmetry_rho::psymm(double* rho_part, // (3) const int ncxy = rho_basis->nx * rho_basis->ny; - double* zpiece = new double[ncxy]; + std::vector zpiece(ncxy); for (int iz = 0; iz < rho_basis->nz; iz++) { - // GlobalV::ofs_running << "\n iz=" << iz; - ModuleBase::GlobalFunc::ZEROS(zpiece, ncxy); + ModuleBase::GlobalFunc::ZEROS(zpiece.data(), ncxy); if (GlobalV::MY_RANK == 0) { for (int ix = 0; ix < rho_basis->nx; ix++) @@ -142,16 +141,11 @@ void Symmetry_rho::psymm(double* rho_part, { const int ir = ix * rho_basis->ny + iy; zpiece[ir] = rhotot[ix * rho_basis->ny * rho_basis->nz + iy * rho_basis->nz + iz]; - // rho[ir*nczp+znow] = zpiece[ir]; } } } - Pgrid.zpiece_to_all(zpiece, iz, rho_part); + Pgrid.zpiece_to_all(zpiece.data(), iz, rho_part); } - - if (GlobalV::MY_RANK == 0) - delete[] rhotot; - delete[] zpiece; #endif return; -} +} \ No newline at end of file From fa53e296f9ae7e90bb651b727362b10479a3055e Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Sun, 28 Jul 2024 00:16:06 +0800 Subject: [PATCH 7/8] Refactor partial charge calc in esolver_ks_pw.cpp to replace raw pointer arrays with std::vector --- source/module_esolver/esolver_ks_pw.cpp | 682 ++++++++++++------------ 1 file changed, 332 insertions(+), 350 deletions(-) diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index 46e1d2207c..707515baf7 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -43,19 +43,23 @@ #ifdef USE_PAW #include "module_cell/module_paw/paw_cell.h" #endif +#include "module_base/formatter.h" + #include #include -#include "module_base/formatter.h" -namespace ModuleESolver { +namespace ModuleESolver +{ template -ESolver_KS_PW::ESolver_KS_PW() { +ESolver_KS_PW::ESolver_KS_PW() +{ this->classname = "ESolver_KS_PW"; this->basisname = "PW"; this->device = base_device::get_device_type(this->ctx); #if ((defined __CUDA) || (defined __ROCM)) - if (this->device == base_device::GpuDevice) { + if (this->device == base_device::GpuDevice) + { hsolver::createGpuBlasHandle(); hsolver::createGpuSolverHandle(); container::kernels::createGpuBlasHandle(); @@ -65,7 +69,8 @@ ESolver_KS_PW::ESolver_KS_PW() { } template -ESolver_KS_PW::~ESolver_KS_PW() { +ESolver_KS_PW::~ESolver_KS_PW() +{ // delete HSolver and ElecState this->deallocate_hsolver(); if (this->pelec != nullptr) @@ -75,7 +80,8 @@ ESolver_KS_PW::~ESolver_KS_PW() { } // delete Hamilt this->deallocate_hamilt(); - if (this->device == base_device::GpuDevice) { + if (this->device == base_device::GpuDevice) + { #if defined(__CUDA) || defined(__ROCM) hsolver::destoryBLAShandle(); hsolver::destroyGpuSolverHandle(); @@ -84,18 +90,17 @@ ESolver_KS_PW::~ESolver_KS_PW() { #endif delete reinterpret_cast*>(this->kspw_psi); } - if (GlobalV::precision_flag == "single") { - delete reinterpret_cast, Device>*>( - this->__kspw_psi); + if (GlobalV::precision_flag == "single") + { + delete reinterpret_cast, Device>*>(this->__kspw_psi); } delete this->psi; delete this->p_wf_init; } template -void ESolver_KS_PW::Init_GlobalC(const Input_para& inp, - UnitCell& ucell, - pseudopot_cell_vnl& ppcell) { +void ESolver_KS_PW::Init_GlobalC(const Input_para& inp, UnitCell& ucell, pseudopot_cell_vnl& ppcell) +{ // GlobalC is a historically left-over namespace, it is used to store global // classes, including: pseudopot_cell_vnl: pseudopotential in cell, V // non-local UnitCell: cell information with atomic properties Grid_Driver: @@ -109,7 +114,8 @@ void ESolver_KS_PW::Init_GlobalC(const Input_para& inp, // about how to organize information stored in classes above, please feel // free to discuss with issue or pull request. - if (this->psi != nullptr) { + if (this->psi != nullptr) + { delete this->psi; } @@ -132,21 +138,21 @@ void ESolver_KS_PW::Init_GlobalC(const Input_para& inp, this->pw_wfc->npwk_max, &(this->sf)); - this->kspw_psi - = GlobalV::device_flag == "gpu" || GlobalV::precision_flag == "single" - ? new psi::Psi(this->psi[0]) - : reinterpret_cast*>(this->psi); + this->kspw_psi = GlobalV::device_flag == "gpu" || GlobalV::precision_flag == "single" + ? new psi::Psi(this->psi[0]) + : reinterpret_cast*>(this->psi); - if (GlobalV::precision_flag == "single") { - ModuleBase::Memory::record("Psi_single", - sizeof(T) * this->psi[0].size()); + if (GlobalV::precision_flag == "single") + { + ModuleBase::Memory::record("Psi_single", sizeof(T) * this->psi[0].size()); } ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT BASIS"); } template -void ESolver_KS_PW::before_all_runners(const Input_para& inp, UnitCell& ucell) { +void ESolver_KS_PW::before_all_runners(const Input_para& inp, UnitCell& ucell) +{ // 1) call before_all_runners() of ESolver_KS ESolver_KS::before_all_runners(inp, ucell); @@ -157,7 +163,8 @@ void ESolver_KS_PW::before_all_runners(const Input_para& inp, UnitCel } // 3) initialize ElecState, - if (this->pelec == nullptr) { + if (this->pelec == nullptr) + { this->pelec = new elecstate::ElecStatePW(this->pw_wfc, &(this->chr), &(this->kv), @@ -175,7 +182,8 @@ void ESolver_KS_PW::before_all_runners(const Input_para& inp, UnitCel this->pelec->omega = ucell.omega; //! 6) initialize the potential. - if (this->pelec->pot == nullptr) { + if (this->pelec->pot == nullptr) + { this->pelec->pot = new elecstate::Potential(this->pw_rhod, this->pw_rho, &ucell, @@ -205,10 +213,9 @@ void ESolver_KS_PW::before_all_runners(const Input_para& inp, UnitCel this->Init_GlobalC(inp, ucell, GlobalC::ppcell); //! 9) setup occupations - if (PARAM.inp.ocp) { - this->pelec->fixed_weights(PARAM.inp.ocp_kb, - GlobalV::NBANDS, - GlobalV::nelec); + if (PARAM.inp.ocp) + { + this->pelec->fixed_weights(PARAM.inp.ocp_kb, GlobalV::NBANDS, GlobalV::nelec); } } template @@ -234,32 +241,23 @@ void ESolver_KS_PW::deallocate_hamilt() this->p_hamilt = nullptr; } template -void ESolver_KS_PW::init_after_vc(const Input_para& inp, UnitCell& ucell) { +void ESolver_KS_PW::init_after_vc(const Input_para& inp, UnitCell& ucell) +{ ModuleBase::TITLE("ESolver_KS_PW", "init_after_vc"); ModuleBase::timer::tick("ESolver_KS_PW", "init_after_vc"); ESolver_KS::init_after_vc(inp, ucell); - if (inp.mdp.md_prec_level == 2) { - this->pw_wfc->initgrids(ucell.lat0, - ucell.latvec, - this->pw_rho->nx, - this->pw_rho->ny, - this->pw_rho->nz); + if (inp.mdp.md_prec_level == 2) + { + this->pw_wfc->initgrids(ucell.lat0, ucell.latvec, this->pw_rho->nx, this->pw_rho->ny, this->pw_rho->nz); - this->pw_wfc->initparameters(false, - inp.ecutwfc, - this->kv.get_nks(), - this->kv.kvec_d.data()); + this->pw_wfc->initparameters(false, inp.ecutwfc, this->kv.get_nks(), this->kv.kvec_d.data()); #ifdef __MPI - if (PARAM.inp.pw_seed > 0) { - MPI_Allreduce(MPI_IN_PLACE, - &this->pw_wfc->ggecut, - 1, - MPI_DOUBLE, - MPI_MAX, - MPI_COMM_WORLD); + if (PARAM.inp.pw_seed > 0) + { + MPI_Allreduce(MPI_IN_PLACE, &this->pw_wfc->ggecut, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); } // qianrui add 2021-8-13 to make different kpar parameters can get the // same results @@ -267,27 +265,25 @@ void ESolver_KS_PW::init_after_vc(const Input_para& inp, UnitCell& uc this->pw_wfc->setuptransform(); - for (int ik = 0; ik < this->kv.get_nks(); ++ik) { + for (int ik = 0; ik < this->kv.get_nks(); ++ik) + { this->kv.ngk[ik] = this->pw_wfc->npwk[ik]; } - this->pw_wfc->collect_local_pw(inp.erf_ecut, - inp.erf_height, - inp.erf_sigma); + this->pw_wfc->collect_local_pw(inp.erf_ecut, inp.erf_height, inp.erf_sigma); delete this->phsol; this->allocate_hsolver(); delete this->pelec; - this->pelec - = new elecstate::ElecStatePW(this->pw_wfc, - &(this->chr), - (K_Vectors*)(&(this->kv)), - &ucell, - &GlobalC::ppcell, - this->pw_rhod, - this->pw_rho, - this->pw_big); + this->pelec = new elecstate::ElecStatePW(this->pw_wfc, + &(this->chr), + (K_Vectors*)(&(this->kv)), + &ucell, + &GlobalC::ppcell, + this->pw_rhod, + this->pw_rho, + this->pw_big); this->pelec->charge->allocate(GlobalV::NSPIN); @@ -306,31 +302,24 @@ void ESolver_KS_PW::init_after_vc(const Input_para& inp, UnitCell& uc // temporary this->Init_GlobalC(inp, ucell, GlobalC::ppcell); - } else { + } + else + { GlobalC::ppcell.init_vnl(ucell, this->pw_rhod); - ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, - "NON-LOCAL POTENTIAL"); + ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "NON-LOCAL POTENTIAL"); - this->pw_wfc->initgrids(ucell.lat0, - ucell.latvec, - this->pw_wfc->nx, - this->pw_wfc->ny, - this->pw_wfc->nz); + this->pw_wfc->initgrids(ucell.lat0, ucell.latvec, this->pw_wfc->nx, this->pw_wfc->ny, this->pw_wfc->nz); - this->pw_wfc->initparameters(false, - PARAM.inp.ecutwfc, - this->kv.get_nks(), - this->kv.kvec_d.data()); + this->pw_wfc->initparameters(false, PARAM.inp.ecutwfc, this->kv.get_nks(), this->kv.kvec_d.data()); - this->pw_wfc->collect_local_pw(inp.erf_ecut, - inp.erf_height, - inp.erf_sigma); + this->pw_wfc->collect_local_pw(inp.erf_ecut, inp.erf_height, inp.erf_sigma); this->p_wf_init->make_table(this->kv.get_nks(), &this->sf); } #ifdef USE_PAW - if (GlobalV::use_paw) { + if (GlobalV::use_paw) + { GlobalC::paw_cell.set_libpaw_ecut(PARAM.inp.ecutwfc / 2.0, PARAM.inp.ecutwfc / 2.0); // in Hartree GlobalC::paw_cell.set_libpaw_fft(this->pw_wfc->nx, @@ -343,7 +332,8 @@ void ESolver_KS_PW::init_after_vc(const Input_para& inp, UnitCell& uc this->pw_wfc->numz); #ifdef __MPI - if (GlobalV::RANK_IN_POOL == 0) { + if (GlobalV::RANK_IN_POOL == 0) + { GlobalC::paw_cell.prepare_paw(); } #else @@ -355,10 +345,12 @@ void ESolver_KS_PW::init_after_vc(const Input_para& inp, UnitCell& uc std::vector> rhoijselect; std::vector nrhoijsel; #ifdef __MPI - if (GlobalV::RANK_IN_POOL == 0) { + if (GlobalV::RANK_IN_POOL == 0) + { GlobalC::paw_cell.get_rhoijp(rhoijp, rhoijselect, nrhoijsel); - for (int iat = 0; iat < ucell.nat; iat++) { + for (int iat = 0; iat < ucell.nat; iat++) + { GlobalC::paw_cell.set_rhoij(iat, nrhoijsel[iat], rhoijselect[iat].size(), @@ -369,7 +361,8 @@ void ESolver_KS_PW::init_after_vc(const Input_para& inp, UnitCell& uc #else GlobalC::paw_cell.get_rhoijp(rhoijp, rhoijselect, nrhoijsel); - for (int iat = 0; iat < ucell.nat; iat++) { + for (int iat = 0; iat < ucell.nat; iat++) + { GlobalC::paw_cell.set_rhoij(iat, nrhoijsel[iat], rhoijselect[iat].size(), @@ -384,13 +377,16 @@ void ESolver_KS_PW::init_after_vc(const Input_para& inp, UnitCell& uc } template -void ESolver_KS_PW::before_scf(const int istep) { +void ESolver_KS_PW::before_scf(const int istep) +{ ModuleBase::TITLE("ESolver_KS_PW", "before_scf"); - if (GlobalC::ucell.cell_parameter_updated) { + if (GlobalC::ucell.cell_parameter_updated) + { this->init_after_vc(PARAM.inp, GlobalC::ucell); } - if (GlobalC::ucell.ionic_position_updated) { + if (GlobalC::ucell.ionic_position_updated) + { this->CE.update_all_dis(GlobalC::ucell); this->CE.extrapolate_charge( #ifdef __MPI @@ -413,32 +409,33 @@ void ESolver_KS_PW::before_scf(const int istep) { // about vdw, jiyy add vdwd3 and linpz add vdwd2 //---------------------------------------------------------- auto vdw_solver = vdw::make_vdw(GlobalC::ucell, PARAM.inp); - if (vdw_solver != nullptr) { + if (vdw_solver != nullptr) + { this->pelec->f_en.evdw = vdw_solver->get_energy(); } // calculate ewald energy - if (!PARAM.inp.test_skip_ewald) { - this->pelec->f_en.ewald_energy - = H_Ewald_pw::compute_ewald(GlobalC::ucell, - this->pw_rhod, - this->sf.strucFac); + if (!PARAM.inp.test_skip_ewald) + { + this->pelec->f_en.ewald_energy = H_Ewald_pw::compute_ewald(GlobalC::ucell, this->pw_rhod, this->sf.strucFac); } //! cal_ux should be called before init_scf because //! the direction of ux is used in noncoline_rho - if (GlobalV::NSPIN == 4 && GlobalV::DOMAG) { + if (GlobalV::NSPIN == 4 && GlobalV::DOMAG) + { GlobalC::ucell.cal_ux(); } //! calculate the total local pseudopotential in real space this->pelec->init_scf(istep, this->sf.strucFac); - if (PARAM.inp.out_chg == 2) { - for (int is = 0; is < GlobalV::NSPIN; is++) { + if (PARAM.inp.out_chg == 2) + { + for (int is = 0; is < GlobalV::NSPIN; is++) + { std::stringstream ss; - ss << GlobalV::global_out_dir << "SPIN" << is + 1 - << "_CHG_INI.cube"; + ss << GlobalV::global_out_dir << "SPIN" << is + 1 << "_CHG_INI.cube"; ModuleIO::write_rho( #ifdef __MPI this->pw_big->bz, @@ -478,12 +475,9 @@ void ESolver_KS_PW::before_scf(const int istep) { //! initialized first. liuyu comment: Symmetry_rho should be located between //! init_rho and v_of_rho? Symmetry_rho srho; - for (int is = 0; is < GlobalV::NSPIN; is++) { - srho.begin(is, - *(this->pelec->charge), - this->pw_rhod, - GlobalC::Pgrid, - GlobalC::ucell.symm); + for (int is = 0; is < GlobalV::NSPIN; is++) + { + srho.begin(is, *(this->pelec->charge), this->pw_rhod, GlobalC::Pgrid, GlobalC::ucell.symm); } // liuyu move here 2023-10-09 @@ -503,54 +497,59 @@ void ESolver_KS_PW::before_scf(const int istep) { // time before scf. But for random wavefunction, we dont, because random // wavefunction is not related to atomic coordinates. What the old strategy // does is only to initialize for once... - if (((GlobalV::init_wfc == "random") && (istep == 0)) - || (GlobalV::init_wfc != "random")) { - this->p_wf_init->initialize_psi(this->psi, - this->kspw_psi, - this->p_hamilt, - GlobalV::ofs_running); + if (((GlobalV::init_wfc == "random") && (istep == 0)) || (GlobalV::init_wfc != "random")) + { + this->p_wf_init->initialize_psi(this->psi, this->kspw_psi, this->p_hamilt, GlobalV::ofs_running); } } template -void ESolver_KS_PW::others(const int istep) { +void ESolver_KS_PW::others(const int istep) +{ ModuleBase::TITLE("ESolver_KS_PW", "others"); const std::string cal_type = GlobalV::CALCULATION; - if (cal_type == "test_memory") { + if (cal_type == "test_memory") + { Cal_Test::test_memory(this->pw_rho, this->pw_wfc, this->p_chgmix->get_mixing_mode(), this->p_chgmix->get_mixing_ndim()); - } else if (cal_type == "gen_bessel") { + } + else if (cal_type == "gen_bessel") + { Numerical_Descriptor nc; nc.output_descriptor(this->psi[0], PARAM.inp.bessel_descriptor_lmax, PARAM.inp.bessel_descriptor_rcut, PARAM.inp.bessel_descriptor_tolerence, this->kv.get_nks()); - ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, - "GENERATE DESCRIPTOR FOR DEEPKS"); - } else if (cal_type == "nscf") { + ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "GENERATE DESCRIPTOR FOR DEEPKS"); + } + else if (cal_type == "nscf") + { this->nscf(); - } else { - ModuleBase::WARNING_QUIT("ESolver_KS_PW::others", - "CALCULATION type not supported"); + } + else + { + ModuleBase::WARNING_QUIT("ESolver_KS_PW::others", "CALCULATION type not supported"); } return; } template -void ESolver_KS_PW::iter_init(const int istep, const int iter) { - if (iter == 1) { +void ESolver_KS_PW::iter_init(const int istep, const int iter) +{ + if (iter == 1) + { this->p_chgmix->init_mixing(); this->p_chgmix->mixing_restart_step = GlobalV::SCF_NMAX + 1; } // for mixing restart - if (iter == this->p_chgmix->mixing_restart_step - && GlobalV::MIXING_RESTART > 0.0) { + if (iter == this->p_chgmix->mixing_restart_step && GlobalV::MIXING_RESTART > 0.0) + { this->p_chgmix->init_mixing(); } // mohan move harris functional to here, 2012-06-05 @@ -559,35 +558,34 @@ void ESolver_KS_PW::iter_init(const int istep, const int iter) { //(2) save change density as previous charge, // prepared fox mixing. - if (GlobalV::MY_STOGROUP == 0) { + if (GlobalV::MY_STOGROUP == 0) + { this->pelec->charge->save_rho_before_sum_band(); } } // Temporary, it should be replaced by hsolver later. template -void ESolver_KS_PW::hamilt2density(const int istep, - const int iter, - const double ethr) { +void ESolver_KS_PW::hamilt2density(const int istep, const int iter, const double ethr) +{ ModuleBase::timer::tick("ESolver_KS_PW", "hamilt2density"); - if (this->phsol != nullptr) { + if (this->phsol != nullptr) + { // reset energy this->pelec->f_en.eband = 0.0; this->pelec->f_en.demet = 0.0; // choose if psi should be diag in subspace // be careful that istep start from 0 and iter start from 1 // if (iter == 1) - hsolver::DiagoIterAssist::need_subspace - = ((istep == 0 || istep == 1) && iter == 1) ? false : true; + hsolver::DiagoIterAssist::need_subspace = ((istep == 0 || istep == 1) && iter == 1) ? false : true; hsolver::DiagoIterAssist::SCF_ITER = iter; hsolver::DiagoIterAssist::PW_DIAG_THR = ethr; - hsolver::DiagoIterAssist::PW_DIAG_NMAX - = GlobalV::PW_DIAG_NMAX; + hsolver::DiagoIterAssist::PW_DIAG_NMAX = GlobalV::PW_DIAG_NMAX; - this->phsol->solve(this->p_hamilt, // hamilt::Hamilt* pHamilt, - this->kspw_psi[0], // psi::Psi& psi, - this->pelec, // elecstate::ElecState* pelec, + this->phsol->solve(this->p_hamilt, // hamilt::Hamilt* pHamilt, + this->kspw_psi[0], // psi::Psi& psi, + this->pelec, // elecstate::ElecState* pelec, PARAM.inp.ks_solver, PARAM.inp.calculation, PARAM.inp.basis_type, @@ -603,19 +601,21 @@ void ESolver_KS_PW::hamilt2density(const int istep, false); - if (PARAM.inp.out_bandgap) { if (!GlobalV::TWO_EFERMI) { this->pelec->cal_bandgap(); - } else { + } + else + { this->pelec->cal_bandgap_updw(); } } - } else { - ModuleBase::WARNING_QUIT("ESolver_KS_PW", - "HSolver has not been initialed!"); + } + else + { + ModuleBase::WARNING_QUIT("ESolver_KS_PW", "HSolver has not been initialed!"); } // calculate the delta_harris energy @@ -624,12 +624,9 @@ void ESolver_KS_PW::hamilt2density(const int istep, this->pelec->cal_energies(1); Symmetry_rho srho; - for (int is = 0; is < GlobalV::NSPIN; is++) { - srho.begin(is, - *(this->pelec->charge), - this->pw_rhod, - GlobalC::Pgrid, - GlobalC::ucell.symm); + for (int is = 0; is < GlobalV::NSPIN; is++) + { + srho.begin(is, *(this->pelec->charge), this->pw_rhod, GlobalC::Pgrid, GlobalC::ucell.symm); } // compute magnetization, only for LSDA(spin==2) @@ -648,26 +645,32 @@ void ESolver_KS_PW::hamilt2density(const int istep, // Temporary, it should be rewritten with Hamilt class. template -void ESolver_KS_PW::update_pot(const int istep, const int iter) { - if (!this->conv_elec) { - if (GlobalV::NSPIN == 4) { +void ESolver_KS_PW::update_pot(const int istep, const int iter) +{ + if (!this->conv_elec) + { + if (GlobalV::NSPIN == 4) + { GlobalC::ucell.cal_ux(); } - this->pelec->pot->update_from_charge(this->pelec->charge, - &GlobalC::ucell); + this->pelec->pot->update_from_charge(this->pelec->charge, &GlobalC::ucell); this->pelec->f_en.descf = this->pelec->cal_delta_escf(); - } else { + } + else + { this->pelec->cal_converged(); } } template -void ESolver_KS_PW::iter_finish(const int iter) { +void ESolver_KS_PW::iter_finish(const int iter) +{ // liuyu 2023-10-24 // D in uspp need vloc, thus needs update when veff updated // calculate the effective coefficient matrix for non-local pseudopotential // projectors - if (GlobalV::use_uspp) { + if (GlobalV::use_uspp) + { ModuleBase::matrix veff = this->pelec->pot->get_effective_v(); GlobalC::ppcell.cal_effective_D(veff, this->pw_rhod, GlobalC::ucell); } @@ -678,30 +681,32 @@ void ESolver_KS_PW::iter_finish(const int iter) { this->pelec->cal_energies(2); bool print = false; - if (this->out_freq_elec && iter % this->out_freq_elec == 0) { + if (this->out_freq_elec && iter % this->out_freq_elec == 0) + { print = true; } - if (print == true) { - if (PARAM.inp.out_chg > 0) { - for (int is = 0; is < GlobalV::NSPIN; is++) { + if (print == true) + { + if (PARAM.inp.out_chg > 0) + { + for (int is = 0; is < GlobalV::NSPIN; is++) + { this->create_Output_Rho(is, iter, "tmp_").write(); - if (XC_Functional::get_func_type() == 3 - || XC_Functional::get_func_type() == 5) { + if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) + { this->create_Output_Kin(is, iter, "tmp_").write(); } } } // output wavefunctions - if (this->wf.out_wfc_pw == 1 || this->wf.out_wfc_pw == 2) { + if (this->wf.out_wfc_pw == 1 || this->wf.out_wfc_pw == 2) + { std::stringstream ssw; ssw << GlobalV::global_out_dir << "WAVEFUNC"; // mohan update 2011-02-21 // qianrui update 2020-10-17 - ModuleIO::write_wfc_pw(ssw.str(), - this->psi[0], - this->kv, - this->pw_wfc); + ModuleIO::write_wfc_pw(ssw.str(), this->psi[0], this->kv, this->pw_wfc); // ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running,"write wave // functions into file WAVEFUNC.dat"); } @@ -709,11 +714,13 @@ void ESolver_KS_PW::iter_finish(const int iter) { } template -void ESolver_KS_PW::after_scf(const int istep) { +void ESolver_KS_PW::after_scf(const int istep) +{ this->create_Output_Potential(istep).write(); // save charge difference into files for charge extrapolation - if (GlobalV::CALCULATION != "scf") { + if (GlobalV::CALCULATION != "scf") + { this->CE.save_files(istep, GlobalC::ucell, #ifdef __MPI @@ -723,105 +730,106 @@ void ESolver_KS_PW::after_scf(const int istep) { &this->sf); } - if (PARAM.inp.out_chg) { - for (int is = 0; is < GlobalV::NSPIN; is++) { + if (PARAM.inp.out_chg) + { + for (int is = 0; is < GlobalV::NSPIN; is++) + { this->create_Output_Rho(is, istep).write(); - if (XC_Functional::get_func_type() == 3 - || XC_Functional::get_func_type() == 5) { + if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) + { this->create_Output_Kin(is, istep).write(); } } } - if (this->wf.out_wfc_pw == 1 || this->wf.out_wfc_pw == 2) { + if (this->wf.out_wfc_pw == 1 || this->wf.out_wfc_pw == 2) + { std::stringstream ssw; ssw << GlobalV::global_out_dir << "WAVEFUNC"; ModuleIO::write_wfc_pw(ssw.str(), this->psi[0], this->kv, this->pw_wfc); } - ModuleIO::output_convergence_after_scf(this->conv_elec, - this->pelec->f_en.etot); + ModuleIO::output_convergence_after_scf(this->conv_elec, this->pelec->f_en.etot); ModuleIO::output_efermi(this->conv_elec, this->pelec->eferm.ef); - if (GlobalV::OUT_LEVEL != "m") { + if (GlobalV::OUT_LEVEL != "m") + { this->pelec->print_eigenvalue(GlobalV::ofs_running); } - if (this->device == base_device::GpuDevice) { - castmem_2d_d2h_op()( - this->psi[0].get_device(), - this->kspw_psi[0].get_device(), - this->psi[0].get_pointer() - this->psi[0].get_psi_bias(), - this->kspw_psi[0].get_pointer() - this->kspw_psi[0].get_psi_bias(), - this->psi[0].size()); + if (this->device == base_device::GpuDevice) + { + castmem_2d_d2h_op()(this->psi[0].get_device(), + this->kspw_psi[0].get_device(), + this->psi[0].get_pointer() - this->psi[0].get_psi_bias(), + this->kspw_psi[0].get_pointer() - this->kspw_psi[0].get_psi_bias(), + this->psi[0].size()); } // Get bands_to_print through public function of INPUT (returns a const // pointer to string) const std::vector bands_to_print = PARAM.inp.bands_to_print; - if (bands_to_print.size() > 0) { + if (bands_to_print.size() > 0) + { // bands_picked is a vector of 0s and 1s, where 1 means the band is // picked to output std::vector bands_picked; bands_picked.resize(this->kspw_psi->get_nbands()); - ModuleBase::GlobalFunc::ZEROS(bands_picked.data(), - this->kspw_psi->get_nbands()); + ModuleBase::GlobalFunc::ZEROS(bands_picked.data(), this->kspw_psi->get_nbands()); // Check if length of bands_to_print is valid - if (static_cast(bands_to_print.size()) - > this->kspw_psi->get_nbands()) { - ModuleBase::WARNING_QUIT( - "ESolver_KS_PW::after_scf", - "The number of bands specified by `bands_to_print` in the " - "INPUT file exceeds `nbands`!"); + if (static_cast(bands_to_print.size()) > this->kspw_psi->get_nbands()) + { + ModuleBase::WARNING_QUIT("ESolver_KS_PW::after_scf", + "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: bands_to_print) { - if (value != 0 && value != 1) { - ModuleBase::WARNING_QUIT( - "ESolver_KS_PW::after_scf", - "The elements of `bands_to_print` must be either 0 or 1. " - "Invalid values found!"); + for (int value: bands_to_print) + { + if (value != 0 && value != 1) + { + ModuleBase::WARNING_QUIT("ESolver_KS_PW::after_scf", + "The elements of `bands_to_print` must be either 0 or 1. " + "Invalid values found!"); } } // Fill bands_picked with values from bands_to_print // Remaining bands are already set to 0 - int length = std::min(static_cast(bands_to_print.size()), - this->kspw_psi->get_nbands()); - for (int i = 0; i < length; ++i) { + int length = std::min(static_cast(bands_to_print.size()), this->kspw_psi->get_nbands()); + for (int i = 0; i < length; ++i) + { // bands_to_print rely on function parse_expression // Initially designed for ocp_set, which can be double bands_picked[i] = static_cast(bands_to_print[i]); } - std::complex* wfcr - = new std::complex[this->pw_rho->nxyz]; + std::complex* wfcr = new std::complex[this->pw_rho->nxyz]; double* rho_band = new double[this->pw_rho->nxyz]; for (int ib = 0; ib < this->kspw_psi->get_nbands(); ++ib) { // Skip the loop iteration if bands_picked[ib] is 0 - if (!bands_picked[ib]) { + if (!bands_picked[ib]) + { continue; } - for (int i = 0; i < this->pw_rho->nxyz; i++) { + for (int i = 0; i < this->pw_rho->nxyz; i++) + { // Initialize rho_band to zero for each band rho_band[i] = 0.0; } - for (int ik = 0; ik < this->kv.get_nks(); ik++) { + for (int ik = 0; ik < this->kv.get_nks(); ik++) + { this->psi->fix_k(ik); - this->pw_wfc->recip_to_real(this->ctx, - &psi[0](ib, 0), - wfcr, - ik); + this->pw_wfc->recip_to_real(this->ctx, &psi[0](ib, 0), wfcr, ik); - double w1 = static_cast(this->kv.wk[ik] - / GlobalC::ucell.omega); + double w1 = static_cast(this->kv.wk[ik] / GlobalC::ucell.omega); for (int i = 0; i < this->pw_rho->nxyz; i++) { @@ -834,30 +842,27 @@ void ESolver_KS_PW::after_scf(const int istep) { Symmetry_rho srho; for (int is = 0; is < GlobalV::NSPIN; is++) { - // Create pointers to the data structures - double** rho_save_pointers = new double*[GlobalV::NSPIN]; - std::complex** rhog = new std::complex*[GlobalV::NSPIN]; - - for (int s = 0; s < GlobalV::NSPIN; s++) { - rho_save_pointers[s] = rho_band; - rhog[s] = new std::complex[this->pelec->charge->ngmc]; + // Use vector instead of raw pointers + std::vector rho_save_pointers(GlobalV::NSPIN, rho_band); + std::vector>> rhog( + GlobalV::NSPIN, + std::vector>(this->pelec->charge->ngmc)); + + // Convert vector of vectors to vector of pointers + std::vector*> rhog_pointers(GlobalV::NSPIN); + for (int s = 0; s < GlobalV::NSPIN; s++) + { + rhog_pointers[s] = rhog[s].data(); } srho.begin(is, - rho_save_pointers, - rhog, + rho_save_pointers.data(), + rhog_pointers.data(), this->pelec->charge->ngmc, nullptr, this->pw_rhod, GlobalC::Pgrid, GlobalC::ucell.symm); - - for (int s = 0; s < GlobalV::NSPIN; s++) - { - delete[] rhog[s]; - } - delete[] rho_save_pointers; - delete[] rhog; } std::stringstream ssc; @@ -888,24 +893,24 @@ void ESolver_KS_PW::after_scf(const int istep) { } template -double ESolver_KS_PW::cal_energy() { +double ESolver_KS_PW::cal_energy() +{ return this->pelec->f_en.etot; } template -void ESolver_KS_PW::cal_force(ModuleBase::matrix& force) { +void ESolver_KS_PW::cal_force(ModuleBase::matrix& force) +{ Forces ff(GlobalC::ucell.nat); - if (this->__kspw_psi != nullptr && GlobalV::precision_flag == "single") { - delete reinterpret_cast, Device>*>( - this->__kspw_psi); + if (this->__kspw_psi != nullptr && GlobalV::precision_flag == "single") + { + delete reinterpret_cast, Device>*>(this->__kspw_psi); } // Refresh __kspw_psi - this->__kspw_psi - = GlobalV::precision_flag == "single" - ? new psi::Psi, Device>(this->kspw_psi[0]) - : reinterpret_cast, Device>*>( - this->kspw_psi); + this->__kspw_psi = GlobalV::precision_flag == "single" + ? new psi::Psi, Device>(this->kspw_psi[0]) + : reinterpret_cast, Device>*>(this->kspw_psi); // Calculate forces ff.cal_force(force, @@ -919,19 +924,18 @@ void ESolver_KS_PW::cal_force(ModuleBase::matrix& force) { } template -void ESolver_KS_PW::cal_stress(ModuleBase::matrix& stress) { +void ESolver_KS_PW::cal_stress(ModuleBase::matrix& stress) +{ Stress_PW ss(this->pelec); - if (this->__kspw_psi != nullptr && GlobalV::precision_flag == "single") { - delete reinterpret_cast, Device>*>( - this->__kspw_psi); + if (this->__kspw_psi != nullptr && GlobalV::precision_flag == "single") + { + delete reinterpret_cast, Device>*>(this->__kspw_psi); } // Refresh __kspw_psi - this->__kspw_psi - = GlobalV::precision_flag == "single" - ? new psi::Psi, Device>(this->kspw_psi[0]) - : reinterpret_cast, Device>*>( - this->kspw_psi); + this->__kspw_psi = GlobalV::precision_flag == "single" + ? new psi::Psi, Device>(this->kspw_psi[0]) + : reinterpret_cast, Device>*>(this->kspw_psi); ss.cal_stress(stress, GlobalC::ucell, &GlobalC::ppcell, @@ -945,28 +949,25 @@ void ESolver_KS_PW::cal_stress(ModuleBase::matrix& stress) { // external stress double unit_transform = 0.0; - unit_transform - = ModuleBase::RYDBERG_SI / pow(ModuleBase::BOHR_RADIUS_SI, 3) * 1.0e-8; - double external_stress[3] - = {PARAM.inp.press1, PARAM.inp.press2, PARAM.inp.press3}; - for (int i = 0; i < 3; i++) { + unit_transform = ModuleBase::RYDBERG_SI / pow(ModuleBase::BOHR_RADIUS_SI, 3) * 1.0e-8; + double external_stress[3] = {PARAM.inp.press1, PARAM.inp.press2, PARAM.inp.press3}; + for (int i = 0; i < 3; i++) + { stress(i, i) -= external_stress[i] / unit_transform; } GlobalV::PRESSURE = (stress(0, 0) + stress(1, 1) + stress(2, 2)) / 3; } template -void ESolver_KS_PW::after_all_runners() { - GlobalV::ofs_running << "\n\n --------------------------------------------" - << std::endl; +void ESolver_KS_PW::after_all_runners() +{ + GlobalV::ofs_running << "\n\n --------------------------------------------" << std::endl; GlobalV::ofs_running << std::setprecision(16); - GlobalV::ofs_running << " !FINAL_ETOT_IS " - << this->pelec->f_en.etot * ModuleBase::Ry_to_eV - << " eV" << std::endl; - GlobalV::ofs_running << " --------------------------------------------\n\n" - << std::endl; + GlobalV::ofs_running << " !FINAL_ETOT_IS " << this->pelec->f_en.etot * ModuleBase::Ry_to_eV << " eV" << std::endl; + GlobalV::ofs_running << " --------------------------------------------\n\n" << std::endl; - if (PARAM.inp.out_dos != 0 || PARAM.inp.out_band[0] != 0) { + if (PARAM.inp.out_dos != 0 || PARAM.inp.out_band[0] != 0) + { GlobalV::ofs_running << "\n\n\n\n"; GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" ">>>>>>>>>>>>>>>>>>>>>>>>>" @@ -999,17 +1000,16 @@ void ESolver_KS_PW::after_all_runners() { } int nspin0 = 1; - if (GlobalV::NSPIN == 2) { + if (GlobalV::NSPIN == 2) + { nspin0 = 2; } //! print occupation in istate.info - ModuleIO::write_istate_info(this->pelec->ekb, - this->pelec->wg, - this->kv, - &(GlobalC::Pkpoints)); + ModuleIO::write_istate_info(this->pelec->ekb, this->pelec->wg, this->kv, &(GlobalC::Pkpoints)); //! compute density of states - if (PARAM.inp.out_dos) { + if (PARAM.inp.out_dos) + { ModuleIO::write_dos_pw(this->pelec->ekb, this->pelec->wg, this->kv, @@ -1017,26 +1017,26 @@ void ESolver_KS_PW::after_all_runners() { PARAM.inp.dos_scale, PARAM.inp.dos_sigma); - if (nspin0 == 1) { - GlobalV::ofs_running << " Fermi energy is " << this->pelec->eferm.ef - << " Rydberg" << std::endl; - } else if (nspin0 == 2) { - GlobalV::ofs_running << " Fermi energy (spin = 1) is " - << this->pelec->eferm.ef_up << " Rydberg" + if (nspin0 == 1) + { + GlobalV::ofs_running << " Fermi energy is " << this->pelec->eferm.ef << " Rydberg" << std::endl; + } + else if (nspin0 == 2) + { + GlobalV::ofs_running << " Fermi energy (spin = 1) is " << this->pelec->eferm.ef_up << " Rydberg" << std::endl; - GlobalV::ofs_running << " Fermi energy (spin = 2) is " - << this->pelec->eferm.ef_dw << " Rydberg" + GlobalV::ofs_running << " Fermi energy (spin = 2) is " << this->pelec->eferm.ef_dw << " Rydberg" << std::endl; } } if (PARAM.inp.out_band[0]) // pengfei 2014-10-13 { - for (int is = 0; is < nspin0; is++) { + for (int is = 0; is < nspin0; is++) + { std::stringstream ss2; ss2 << GlobalV::global_out_dir << "BANDS_" << is + 1 << ".dat"; - GlobalV::ofs_running << "\n Output bands in file: " << ss2.str() - << std::endl; + GlobalV::ofs_running << "\n Output bands in file: " << ss2.str() << std::endl; ModuleIO::nscf_band(is, ss2.str(), GlobalV::NBANDS, @@ -1048,51 +1048,38 @@ void ESolver_KS_PW::after_all_runners() { } } - if (GlobalV::BASIS_TYPE == "pw" - && winput::out_spillage) // xiaohui add 2013-09-01 + if (GlobalV::BASIS_TYPE == "pw" && winput::out_spillage) // xiaohui add 2013-09-01 { // calculate spillage value. // ! Print out overlap before spillage optimization to generate atomic // orbitals - if (winput::out_spillage <= 2) { - for (int i = 0; i < PARAM.inp.bessel_nao_rcuts.size(); i++) { - if (GlobalV::MY_RANK == 0) { - std::cout << "update value: bessel_nao_rcut <- " - << std::fixed << PARAM.inp.bessel_nao_rcuts[i] + if (winput::out_spillage <= 2) + { + for (int i = 0; i < PARAM.inp.bessel_nao_rcuts.size(); i++) + { + if (GlobalV::MY_RANK == 0) + { + std::cout << "update value: bessel_nao_rcut <- " << std::fixed << PARAM.inp.bessel_nao_rcuts[i] << " a.u." << std::endl; } Numerical_Basis numerical_basis; - numerical_basis.output_overlap(this->psi[0], - this->sf, - this->kv, - this->pw_wfc, - GlobalC::ucell, - i); + numerical_basis.output_overlap(this->psi[0], this->sf, this->kv, this->pw_wfc, GlobalC::ucell, i); } - ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, - "BASIS OVERLAP (Q and S) GENERATION."); + ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "BASIS OVERLAP (Q and S) GENERATION."); } } //! Print out wave functions in real space if (this->wf.out_wfc_r == 1) // Peize Lin add 2021.11.21 { - ModuleIO::write_psi_r_1(this->psi[0], - this->pw_wfc, - "wfc_realspace", - true, - this->kv); + ModuleIO::write_psi_r_1(this->psi[0], this->pw_wfc, "wfc_realspace", true, this->kv); } //! Use Kubo-Greenwood method to compute conductivities - if (PARAM.inp.cal_cond) { - EleCond elec_cond(&GlobalC::ucell, - &this->kv, - this->pelec, - this->pw_wfc, - this->psi, - &GlobalC::ppcell); + if (PARAM.inp.cal_cond) + { + EleCond elec_cond(&GlobalC::ucell, &this->kv, this->pelec, this->pw_wfc, this->psi, &GlobalC::ppcell); elec_cond.KG(PARAM.inp.cond_smear, PARAM.inp.cond_fwhm, PARAM.inp.cond_wcut, @@ -1104,8 +1091,10 @@ void ESolver_KS_PW::after_all_runners() { } template -void ESolver_KS_PW::hamilt2estates(const double ethr) { - if (this->phsol != nullptr) { +void ESolver_KS_PW::hamilt2estates(const double ethr) +{ + if (this->phsol != nullptr) + { hsolver::DiagoIterAssist::need_subspace = false; hsolver::DiagoIterAssist::PW_DIAG_THR = ethr; this->phsol->solve(this->p_hamilt, @@ -1125,14 +1114,16 @@ void ESolver_KS_PW::hamilt2estates(const double ethr) { hsolver::DiagoIterAssist::PW_DIAG_THR, true); - } else { - ModuleBase::WARNING_QUIT("ESolver_KS_PW", - "HSolver has not been initialed!"); + } + else + { + ModuleBase::WARNING_QUIT("ESolver_KS_PW", "HSolver has not been initialed!"); } } template -void ESolver_KS_PW::nscf() { +void ESolver_KS_PW::nscf() +{ ModuleBase::TITLE("ESolver_KS_PW", "nscf"); ModuleBase::timer::tick("ESolver_KS_PW", "nscf"); @@ -1142,10 +1133,9 @@ void ESolver_KS_PW::nscf() { //! 2) setup the parameters for diagonalization double diag_ethr = GlobalV::PW_DIAG_THR; - if (diag_ethr - 1e-2 > -1e-5) { - diag_ethr - = std::max(1e-13, - 0.1 * std::min(1e-2, PARAM.inp.scf_thr / GlobalV::nelec)); + if (diag_ethr - 1e-2 > -1e-5) + { + diag_ethr = std::max(1e-13, 0.1 * std::min(1e-2, PARAM.inp.scf_thr / GlobalV::nelec)); } GlobalV::ofs_running << " PW_DIAG_THR = " << diag_ethr << std::endl; @@ -1154,62 +1144,62 @@ void ESolver_KS_PW::nscf() { //! 3) calculate weights/Fermi energies this->pelec->calculate_weights(); - GlobalV::ofs_running << "\n End of Band Structure Calculation \n" - << std::endl; + GlobalV::ofs_running << "\n End of Band Structure Calculation \n" << std::endl; //! 4) print out band energies and weights std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "writing band energies"); const int nspin = GlobalV::NSPIN; const int nbands = GlobalV::NBANDS; - for (int ik = 0; ik < this->kv.get_nks(); ik++) { - if (nspin == 2) { - if (ik == 0) { + for (int ik = 0; ik < this->kv.get_nks(); ik++) + { + if (nspin == 2) + { + if (ik == 0) + { GlobalV::ofs_running << " spin up :" << std::endl; } - if (ik == (this->kv.get_nks() / 2)) { + if (ik == (this->kv.get_nks() / 2)) + { GlobalV::ofs_running << " spin down :" << std::endl; } } - GlobalV::ofs_running - << " k-points" << ik + 1 << "(" << this->kv.get_nkstot() - << "): " << this->kv.kvec_c[ik].x << " " << this->kv.kvec_c[ik].y - << " " << this->kv.kvec_c[ik].z << std::endl; + GlobalV::ofs_running << " k-points" << ik + 1 << "(" << this->kv.get_nkstot() << "): " << this->kv.kvec_c[ik].x + << " " << this->kv.kvec_c[ik].y << " " << this->kv.kvec_c[ik].z << std::endl; - for (int ib = 0; ib < nbands; ib++) { - GlobalV::ofs_running - << " spin" << this->kv.isk[ik] + 1 << "_final_band " << ib + 1 - << " " << this->pelec->ekb(ik, ib) * ModuleBase::Ry_to_eV << " " - << this->pelec->wg(ik, ib) * this->kv.get_nks() << std::endl; + for (int ib = 0; ib < nbands; ib++) + { + GlobalV::ofs_running << " spin" << this->kv.isk[ik] + 1 << "_final_band " << ib + 1 << " " + << this->pelec->ekb(ik, ib) * ModuleBase::Ry_to_eV << " " + << this->pelec->wg(ik, ib) * this->kv.get_nks() << std::endl; } GlobalV::ofs_running << std::endl; } std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "writing band energies"); //! 5) print out band gaps - if (PARAM.inp.out_bandgap) { + if (PARAM.inp.out_bandgap) + { std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "writing band gaps"); - if (!GlobalV::TWO_EFERMI) { + if (!GlobalV::TWO_EFERMI) + { this->pelec->cal_bandgap(); - GlobalV::ofs_running << " E_bandgap " - << this->pelec->bandgap * ModuleBase::Ry_to_eV - << " eV" << std::endl; - } else { + GlobalV::ofs_running << " E_bandgap " << this->pelec->bandgap * ModuleBase::Ry_to_eV << " eV" << std::endl; + } + else + { this->pelec->cal_bandgap_updw(); - GlobalV::ofs_running - << " E_bandgap_up " - << this->pelec->bandgap_up * ModuleBase::Ry_to_eV << " eV" - << std::endl; - GlobalV::ofs_running - << " E_bandgap_dw " - << this->pelec->bandgap_dw * ModuleBase::Ry_to_eV << " eV" - << std::endl; + GlobalV::ofs_running << " E_bandgap_up " << this->pelec->bandgap_up * ModuleBase::Ry_to_eV << " eV" + << std::endl; + GlobalV::ofs_running << " E_bandgap_dw " << this->pelec->bandgap_dw * ModuleBase::Ry_to_eV << " eV" + << std::endl; } std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "writing band gaps"); } //! 6) calculate Wannier functions - if (PARAM.inp.towannier90) { + if (PARAM.inp.towannier90) + { std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "Wannier functions calculation"); toWannier90_PW wan(PARAM.inp.out_wannier_mmn, PARAM.inp.out_wannier_amn, @@ -1219,24 +1209,16 @@ void ESolver_KS_PW::nscf() { PARAM.inp.nnkpfile, PARAM.inp.wannier_spin); - wan.calculate(this->pelec->ekb, - this->pw_wfc, - this->pw_big, - this->kv, - this->psi); + wan.calculate(this->pelec->ekb, this->pw_wfc, this->pw_big, this->kv, this->psi); std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "Wannier functions calculation"); } //! 7) calculate Berry phase polarization - if (berryphase::berry_phase_flag - && ModuleSymmetry::Symmetry::symm_flag != 1) { + if (berryphase::berry_phase_flag && ModuleSymmetry::Symmetry::symm_flag != 1) + { std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "Berry phase polarization"); berryphase bp; - bp.Macroscopic_polarization(this->pw_wfc->npwk_max, - this->psi, - this->pw_rho, - this->pw_wfc, - this->kv); + bp.Macroscopic_polarization(this->pw_wfc->npwk_max, this->psi, this->pw_rho, this->pw_wfc, this->kv); std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "Berry phase polarization"); } From dd4c86505f40abbb5d7d97ef163eca773002ecb8 Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Sun, 28 Jul 2024 00:35:30 +0800 Subject: [PATCH 8/8] Rename files: istate_charge to get_pchg, and istate_envelope to get_wf --- source/Makefile.Objects | 4 ++-- source/module_esolver/esolver_ks_lcao_elec.cpp | 4 ++-- source/module_io/CMakeLists.txt | 4 ++-- source/module_io/{istate_charge.cpp => get_pchg.cpp} | 2 +- source/module_io/{istate_charge.h => get_pchg.h} | 0 source/module_io/{istate_envelope.cpp => get_wf.cpp} | 2 +- source/module_io/{istate_envelope.h => get_wf.h} | 0 7 files changed, 8 insertions(+), 8 deletions(-) rename source/module_io/{istate_charge.cpp => get_pchg.cpp} (99%) rename source/module_io/{istate_charge.h => get_pchg.h} (100%) rename source/module_io/{istate_envelope.cpp => get_wf.cpp} (99%) rename source/module_io/{istate_envelope.h => get_wf.h} (100%) diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 8a3cb88425..d24d485b98 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -517,8 +517,8 @@ OBJS_IO_LCAO=cal_r_overlap_R.o\ write_proj_band_lcao.o\ write_istate_info.o\ nscf_fermi_surf.o\ - istate_charge.o\ - istate_envelope.o\ + get_pchg.o\ + get_wf.o\ io_dmk.o\ unk_overlap_lcao.o\ read_wfc_nao.o\ diff --git a/source/module_esolver/esolver_ks_lcao_elec.cpp b/source/module_esolver/esolver_ks_lcao_elec.cpp index 4614fb6d4b..11c63cac25 100644 --- a/source/module_esolver/esolver_ks_lcao_elec.cpp +++ b/source/module_esolver/esolver_ks_lcao_elec.cpp @@ -8,8 +8,8 @@ #include "module_cell/module_neighbor/sltk_atom_arrange.h" #include "module_cell/module_neighbor/sltk_grid_driver.h" #include "module_io/berryphase.h" -#include "module_io/istate_charge.h" -#include "module_io/istate_envelope.h" +#include "module_io/get_pchg.h" +#include "module_io/get_wf.h" #include "module_io/to_wannier90_lcao.h" #include "module_io/to_wannier90_lcao_in_pw.h" #include "module_io/write_HS_R.h" diff --git a/source/module_io/CMakeLists.txt b/source/module_io/CMakeLists.txt index d74e529eb6..9d8a66972a 100644 --- a/source/module_io/CMakeLists.txt +++ b/source/module_io/CMakeLists.txt @@ -50,8 +50,8 @@ if(ENABLE_LCAO) write_orb_info.cpp write_proj_band_lcao.cpp nscf_fermi_surf.cpp - istate_charge.cpp - istate_envelope.cpp + get_pchg.cpp + get_wf.cpp read_wfc_nao.cpp read_wfc_lcao.cpp write_wfc_nao.cpp diff --git a/source/module_io/istate_charge.cpp b/source/module_io/get_pchg.cpp similarity index 99% rename from source/module_io/istate_charge.cpp rename to source/module_io/get_pchg.cpp index cb33630431..600335112c 100644 --- a/source/module_io/istate_charge.cpp +++ b/source/module_io/get_pchg.cpp @@ -1,4 +1,4 @@ -#include "istate_charge.h" +#include "get_pchg.h" #include "module_base/blas_connector.h" #include "module_base/global_function.h" diff --git a/source/module_io/istate_charge.h b/source/module_io/get_pchg.h similarity index 100% rename from source/module_io/istate_charge.h rename to source/module_io/get_pchg.h diff --git a/source/module_io/istate_envelope.cpp b/source/module_io/get_wf.cpp similarity index 99% rename from source/module_io/istate_envelope.cpp rename to source/module_io/get_wf.cpp index 5d7d224d24..7adf301eaa 100644 --- a/source/module_io/istate_envelope.cpp +++ b/source/module_io/get_wf.cpp @@ -1,4 +1,4 @@ -#include "istate_envelope.h" +#include "get_wf.h" #include "module_base/global_function.h" #include "module_base/global_variable.h" diff --git a/source/module_io/istate_envelope.h b/source/module_io/get_wf.h similarity index 100% rename from source/module_io/istate_envelope.h rename to source/module_io/get_wf.h