diff --git a/source/Makefile.Objects b/source/Makefile.Objects index c58a651601..5cc976277b 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -477,7 +477,7 @@ OBJS_IO_LCAO=cal_r_overlap_R.o\ nscf_fermi_surf.o\ istate_charge.o\ istate_envelope.o\ - read_dm.o\ + io_dmk.o\ unk_overlap_lcao.o\ read_wfc_nao.o\ read_wfc_lcao.o\ @@ -485,8 +485,6 @@ OBJS_IO_LCAO=cal_r_overlap_R.o\ write_HS_sparse.o\ single_R_io.o\ write_HS_R.o\ - write_dm.o\ - output_dm.o\ write_dmr.o\ sparse_matrix.o\ output_mulliken.o\ diff --git a/source/module_elecstate/elecstate_lcao.cpp b/source/module_elecstate/elecstate_lcao.cpp index 9b2a294dd4..85a6791f30 100644 --- a/source/module_elecstate/elecstate_lcao.cpp +++ b/source/module_elecstate/elecstate_lcao.cpp @@ -10,46 +10,51 @@ #include -namespace elecstate -{ +namespace elecstate { // multi-k case template <> -void ElecStateLCAO>::psiToRho(const psi::Psi>& psi) -{ +void ElecStateLCAO>::psiToRho( + const psi::Psi>& psi) { ModuleBase::TITLE("ElecStateLCAO", "psiToRho"); ModuleBase::timer::tick("ElecStateLCAO", "psiToRho"); this->calculate_weights(); - // the calculations of dm, and dm -> rho are, technically, two separate functionalities, as we cannot - // rule out the possibility that we may have a dm from other sources, such as read from file. - // However, since we are not separating them now, I opt to add a flag to control how dm is obtained as of now - if (!GlobalV::dm_to_rho) - { + // the calculations of dm, and dm -> rho are, technically, two separate + // functionalities, as we cannot rule out the possibility that we may have a + // dm from other sources, such as read from file. However, since we are not + // separating them now, I opt to add a flag to control how dm is obtained as + // of now + if (!GlobalV::dm_to_rho) { this->calEBand(); ModuleBase::GlobalFunc::NOTE("Calculate the density matrix."); - // this part for calculating DMK in 2d-block format, not used for charge now + // this part for calculating DMK in 2d-block format, not used for charge + // now // psi::Psi> dm_k_2d(); - if (GlobalV::KS_SOLVER == "genelpa" || GlobalV::KS_SOLVER == "scalapack_gvx" || GlobalV::KS_SOLVER == "lapack" - || GlobalV::KS_SOLVER == "cusolver" || GlobalV::KS_SOLVER == "cusolvermp" + if (GlobalV::KS_SOLVER == "genelpa" + || GlobalV::KS_SOLVER == "scalapack_gvx" + || GlobalV::KS_SOLVER == "lapack" + || GlobalV::KS_SOLVER == "cusolver" + || GlobalV::KS_SOLVER == "cusolvermp" || GlobalV::KS_SOLVER == "cg_in_lcao") // Peize Lin test 2019-05-15 { // cal_dm(this->loc->ParaV, this->wg, psi, this->loc->dm_k); - elecstate::cal_dm_psi(this->DM->get_paraV_pointer(), this->wg, psi, *(this->DM)); + elecstate::cal_dm_psi(this->DM->get_paraV_pointer(), + this->wg, + psi, + *(this->DM)); this->DM->cal_DMR(); // interface for RI-related calculation, which needs loc.dm_k #ifdef __EXX - if (GlobalC::exx_info.info_global.cal_exx) - { + if (GlobalC::exx_info.info_global.cal_exx) { const K_Vectors* kv = this->DM->get_kv_pointer(); this->loc->dm_k.resize(kv->get_nks()); - for (int ik = 0; ik < kv->get_nks(); ++ik) - { + for (int ik = 0; ik < kv->get_nks(); ++ik) { this->loc->set_dm_k(ik, this->DM->get_DMK_pointer(ik)); } } @@ -57,9 +62,9 @@ void ElecStateLCAO>::psiToRho(const psi::Psicharge->rho[is], this->charge->nrxx); // mohan 2009-11-10 + for (int is = 0; is < GlobalV::NSPIN; is++) { + ModuleBase::GlobalFunc::ZEROS(this->charge->rho[is], + this->charge->nrxx); // mohan 2009-11-10 } //------------------------------------------------------------ @@ -67,15 +72,16 @@ void ElecStateLCAO>::psiToRho(const psi::Psigint_k->transfer_DM2DtoGrid(this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint + this->gint_k->transfer_DM2DtoGrid( + this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint Gint_inout inout(this->charge->rho, Gint_Tools::job_type::rho); this->gint_k->cal_gint(&inout); - if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) - { - for (int is = 0; is < GlobalV::NSPIN; is++) - { - ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[is], this->charge->nrxx); + if (XC_Functional::get_func_type() == 3 + || XC_Functional::get_func_type() == 5) { + for (int is = 0; is < GlobalV::NSPIN; is++) { + ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[is], + this->charge->nrxx); } Gint_inout inout1(this->charge->kin_r, Gint_Tools::job_type::tau); this->gint_k->cal_gint(&inout1); @@ -89,46 +95,31 @@ void ElecStateLCAO>::psiToRho(const psi::Psi -void ElecStateLCAO::psiToRho(const psi::Psi& psi) -{ +void ElecStateLCAO::psiToRho(const psi::Psi& psi) { ModuleBase::TITLE("ElecStateLCAO", "psiToRho"); ModuleBase::timer::tick("ElecStateLCAO", "psiToRho"); this->calculate_weights(); this->calEBand(); - if (GlobalV::KS_SOLVER == "genelpa" || GlobalV::KS_SOLVER == "scalapack_gvx" || GlobalV::KS_SOLVER == "lapack" - || GlobalV::KS_SOLVER == "cusolver" || GlobalV::KS_SOLVER == "cusolvermp" || GlobalV::KS_SOLVER == "cg_in_lcao") - { + if (GlobalV::KS_SOLVER == "genelpa" || GlobalV::KS_SOLVER == "scalapack_gvx" + || GlobalV::KS_SOLVER == "lapack" || GlobalV::KS_SOLVER == "cusolver" + || GlobalV::KS_SOLVER == "cusolvermp" + || GlobalV::KS_SOLVER == "cg_in_lcao") { ModuleBase::timer::tick("ElecStateLCAO", "cal_dm_2d"); // get DMK in 2d-block format // cal_dm(this->loc->ParaV, this->wg, psi, this->loc->dm_gamma); - elecstate::cal_dm_psi(this->DM->get_paraV_pointer(), this->wg, psi, *(this->DM)); + elecstate::cal_dm_psi(this->DM->get_paraV_pointer(), + this->wg, + psi, + *(this->DM)); this->DM->cal_DMR(); - if (this->loc->out_dm) // keep interface for old Output_DM until new one is ready - { - this->loc->dm_gamma.resize(GlobalV::NSPIN); - for (int is = 0; is < GlobalV::NSPIN; ++is) - { - this->loc->set_dm_gamma(is, this->DM->get_DMK_pointer(is)); - } - } - ModuleBase::timer::tick("ElecStateLCAO", "cal_dm_2d"); - for (int ik = 0; ik < psi.get_nk(); ++ik) - { - // for gamma_only case, no convertion occured, just for print. - // old 2D-to-Grid conversion has been replaced by new Gint Refactor 2023/09/25 - if (this->loc->out_dm) // keep interface for old Output_DM until new one is ready - { - this->loc->cal_dk_gamma_from_2D_pub(); - } - } } - for (int is = 0; is < GlobalV::NSPIN; is++) - { - ModuleBase::GlobalFunc::ZEROS(this->charge->rho[is], this->charge->nrxx); // mohan 2009-11-10 + for (int is = 0; is < GlobalV::NSPIN; is++) { + ModuleBase::GlobalFunc::ZEROS(this->charge->rho[is], + this->charge->nrxx); // mohan 2009-11-10 } //------------------------------------------------------------ @@ -136,17 +127,18 @@ void ElecStateLCAO::psiToRho(const psi::Psi& psi) //------------------------------------------------------------ ModuleBase::GlobalFunc::NOTE("Calculate the charge on real space grid!"); - this->gint_gamma->transfer_DM2DtoGrid(this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint + this->gint_gamma->transfer_DM2DtoGrid( + this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint Gint_inout inout(this->charge->rho, Gint_Tools::job_type::rho); this->gint_gamma->cal_gint(&inout); - if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) - { - for (int is = 0; is < GlobalV::NSPIN; is++) - { - ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[is], this->charge->nrxx); + if (XC_Functional::get_func_type() == 3 + || XC_Functional::get_func_type() == 5) { + for (int is = 0; is < GlobalV::NSPIN; is++) { + ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[is], + this->charge->nrxx); } Gint_inout inout1(this->charge->kin_r, Gint_Tools::job_type::tau); this->gint_gamma->cal_gint(&inout1); @@ -159,21 +151,21 @@ void ElecStateLCAO::psiToRho(const psi::Psi& psi) } template -void ElecStateLCAO::init_DM(const K_Vectors* kv, const Parallel_Orbitals* paraV, const int nspin) -{ +void ElecStateLCAO::init_DM(const K_Vectors* kv, + const Parallel_Orbitals* paraV, + const int nspin) { this->DM = new DensityMatrix(kv, paraV, nspin); } template <> -double ElecStateLCAO::get_spin_constrain_energy() -{ - SpinConstrain& sc = SpinConstrain::getScInstance(); +double ElecStateLCAO::get_spin_constrain_energy() { + SpinConstrain& sc + = SpinConstrain::getScInstance(); return sc.cal_escon(); } template <> -double ElecStateLCAO>::get_spin_constrain_energy() -{ +double ElecStateLCAO>::get_spin_constrain_energy() { SpinConstrain, base_device::DEVICE_CPU>& sc = SpinConstrain>::getScInstance(); return sc.cal_escon(); @@ -181,48 +173,37 @@ double ElecStateLCAO>::get_spin_constrain_energy() #ifdef __PEXSI template <> -void ElecStateLCAO::dmToRho(std::vector pexsi_DM, std::vector pexsi_EDM) -{ +void ElecStateLCAO::dmToRho(std::vector pexsi_DM, + std::vector pexsi_EDM) { ModuleBase::timer::tick("ElecStateLCAO", "dmToRho"); int nspin = GlobalV::NSPIN; - if (GlobalV::NSPIN == 4) - { + if (GlobalV::NSPIN == 4) { nspin = 1; } - // old 2D-to-Grid conversion has been replaced by new Gint Refactor 2023/09/25 - if (this->loc->out_dm) // keep interface for old Output_DM until new one is ready - { - for (int is = 0; is < nspin; ++is) - { - this->loc->set_dm_gamma(is, pexsi_DM[is]); - } - this->loc->cal_dk_gamma_from_2D_pub(); - } - this->get_DM()->pexsi_EDM = pexsi_EDM; - for (int is = 0; is < nspin; is++) - { + for (int is = 0; is < nspin; is++) { this->DM->set_DMK_pointer(is, pexsi_DM[is]); } DM->cal_DMR(); - for (int is = 0; is < GlobalV::NSPIN; is++) - { - ModuleBase::GlobalFunc::ZEROS(this->charge->rho[is], this->charge->nrxx); // mohan 2009-11-10 + for (int is = 0; is < GlobalV::NSPIN; is++) { + ModuleBase::GlobalFunc::ZEROS(this->charge->rho[is], + this->charge->nrxx); // mohan 2009-11-10 } ModuleBase::GlobalFunc::NOTE("Calculate the charge on real space grid!"); - this->gint_gamma->transfer_DM2DtoGrid(this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint + this->gint_gamma->transfer_DM2DtoGrid( + this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint Gint_inout inout(this->charge->rho, Gint_Tools::job_type::rho); this->gint_gamma->cal_gint(&inout); - if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) - { - for (int is = 0; is < GlobalV::NSPIN; is++) - { - ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[0], this->charge->nrxx); + if (XC_Functional::get_func_type() == 3 + || XC_Functional::get_func_type() == 5) { + for (int is = 0; is < GlobalV::NSPIN; is++) { + ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[0], + this->charge->nrxx); } Gint_inout inout1(this->charge->kin_r, Gint_Tools::job_type::tau); this->gint_gamma->cal_gint(&inout1); @@ -235,10 +216,11 @@ void ElecStateLCAO::dmToRho(std::vector pexsi_DM, std::vector -void ElecStateLCAO>::dmToRho(std::vector*> pexsi_DM, - std::vector*> pexsi_EDM) -{ - ModuleBase::WARNING_QUIT("ElecStateLCAO", "pexsi is not completed for multi-k case"); +void ElecStateLCAO>::dmToRho( + std::vector*> pexsi_DM, + std::vector*> pexsi_EDM) { + ModuleBase::WARNING_QUIT("ElecStateLCAO", + "pexsi is not completed for multi-k case"); } #endif diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp index e241826941..094167d61c 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -1,6 +1,6 @@ #include "esolver_ks.h" -#include +#include #ifdef __MPI #include #else @@ -23,8 +23,7 @@ #endif #include "module_io/json_output/output_info.h" -namespace ModuleESolver -{ +namespace ModuleESolver { //------------------------------------------------------------------------------ //! the 1st function of ESolver_KS: constructor @@ -33,8 +32,7 @@ namespace ModuleESolver // assumption that INPUT has been initialized, mohan 2024-05-12 //------------------------------------------------------------------------------ template -ESolver_KS::ESolver_KS() -{ +ESolver_KS::ESolver_KS() { classname = "ESolver_KS"; basisname = "PLEASE ADD BASISNAME FOR CURRENT ESOLVER."; @@ -51,8 +49,10 @@ ESolver_KS::ESolver_KS() // pw_rho = new ModuleBase::PW_Basis(); // temporary, it will be removed - pw_wfc = new ModulePW::PW_Basis_K_Big(GlobalV::device_flag, GlobalV::precision_flag); - ModulePW::PW_Basis_K_Big* tmp = static_cast(pw_wfc); + pw_wfc = new ModulePW::PW_Basis_K_Big(GlobalV::device_flag, + GlobalV::precision_flag); + ModulePW::PW_Basis_K_Big* tmp + = static_cast(pw_wfc); // should not use INPUT here, mohan 2024-05-12 tmp->setbxyz(INPUT.bx, INPUT.by, INPUT.bz); @@ -77,8 +77,7 @@ ESolver_KS::ESolver_KS() //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -ESolver_KS::~ESolver_KS() -{ +ESolver_KS::~ESolver_KS() { delete this->psi; delete this->pw_wfc; delete this->p_hamilt; @@ -91,8 +90,7 @@ ESolver_KS::~ESolver_KS() //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -void ESolver_KS::before_all_runners(Input& inp, UnitCell& ucell) -{ +void ESolver_KS::before_all_runners(Input& inp, UnitCell& ucell) { ModuleBase::TITLE("ESolver_KS", "before_all_runners"); //! 1) initialize "before_all_runniers" in ESolver_FP @@ -112,8 +110,7 @@ void ESolver_KS::before_all_runners(Input& inp, UnitCell& ucell) /// PAW Section #ifdef USE_PAW - if (GlobalV::use_paw) - { + if (GlobalV::use_paw) { int* atom_type = nullptr; double** atom_coord = nullptr; std::vector filename_list; @@ -122,16 +119,13 @@ void ESolver_KS::before_all_runners(Input& inp, UnitCell& ucell) atom_coord = new double*[ucell.nat]; filename_list.resize(ucell.ntype); - for (int ia = 0; ia < ucell.nat; ia++) - { + for (int ia = 0; ia < ucell.nat; ia++) { atom_coord[ia] = new double[3]; } int iat = 0; - for (int it = 0; it < ucell.ntype; it++) - { - for (int ia = 0; ia < ucell.atoms[it].na; ia++) - { + for (int it = 0; it < ucell.ntype; it++) { + for (int ia = 0; ia < ucell.atoms[it].na; ia++) { atom_type[iat] = it; atom_coord[iat][0] = ucell.atoms[it].taud[ia].x; atom_coord[iat][1] = ucell.atoms[it].taud[ia].y; @@ -140,30 +134,26 @@ void ESolver_KS::before_all_runners(Input& inp, UnitCell& ucell) } } - if (GlobalV::MY_RANK == 0) - { + if (GlobalV::MY_RANK == 0) { std::ifstream ifa(GlobalV::stru_file.c_str(), std::ios::in); - if (!ifa) - { - ModuleBase::WARNING_QUIT("set_libpaw_files", "can not open stru file"); + if (!ifa) { + ModuleBase::WARNING_QUIT("set_libpaw_files", + "can not open stru file"); } std::string line; - while (!ifa.eof()) - { + while (!ifa.eof()) { getline(ifa, line); if (line.find("PAW_FILES") != std::string::npos) break; } - for (int it = 0; it < ucell.ntype; it++) - { + for (int it = 0; it < ucell.ntype; it++) { ifa >> filename_list[it]; } } #ifdef __MPI - for (int it = 0; it < ucell.ntype; it++) - { + for (int it = 0; it < ucell.ntype; it++) { Parallel_Common::bcast_string(filename_list[it]); } #endif @@ -177,8 +167,7 @@ void ESolver_KS::before_all_runners(Input& inp, UnitCell& ucell) (const double**)atom_coord, filename_list); - for (int iat = 0; iat < ucell.nat; iat++) - { + for (int iat = 0; iat < ucell.nat; iat++) { delete[] atom_coord[iat]; } delete[] atom_coord; @@ -193,26 +182,30 @@ void ESolver_KS::before_all_runners(Input& inp, UnitCell& ucell) //! 4) it has been established that // xc_func is same for all elements, therefore // only the first one if used - if (GlobalV::use_paw) - { + if (GlobalV::use_paw) { XC_Functional::set_xc_type(GlobalV::DFT_FUNCTIONAL); - } - else - { + } else { XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func); } ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SETUP UNITCELL"); //! 5) ESolver depends on the Symmetry module // symmetry analysis should be performed every time the cell is changed - if (ModuleSymmetry::Symmetry::symm_flag == 1) - { - ucell.symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running); + if (ModuleSymmetry::Symmetry::symm_flag == 1) { + ucell.symm.analy_sys(ucell.lat, + ucell.st, + ucell.atoms, + GlobalV::ofs_running); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY"); } //! 6) Setup the k points according to symmetry. - this->kv.set(ucell.symm, GlobalV::global_kpoint_card, GlobalV::NSPIN, ucell.G, ucell.latvec, GlobalV::ofs_running); + this->kv.set(ucell.symm, + GlobalV::global_kpoint_card, + GlobalV::NSPIN, + ucell.G, + ucell.latvec, + GlobalV::ofs_running); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS"); @@ -221,7 +214,9 @@ void ESolver_KS::before_all_runners(Input& inp, UnitCell& ucell) //! 8) new plane wave basis, fft grids, etc. #ifdef __MPI - this->pw_wfc->initmpi(GlobalV::NPROC_IN_POOL, GlobalV::RANK_IN_POOL, POOL_WORLD); + this->pw_wfc->initmpi(GlobalV::NPROC_IN_POOL, + GlobalV::RANK_IN_POOL, + POOL_WORLD); #endif this->pw_wfc->initgrids(inp.ref_cell_factor * ucell.lat0, @@ -230,15 +225,23 @@ void ESolver_KS::before_all_runners(Input& inp, UnitCell& ucell) 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()); // the MPI allreduce should not be here, mohan 2024-05-12 #ifdef __MPI - if (inp.pw_seed > 0) - { - MPI_Allreduce(MPI_IN_PLACE, &this->pw_wfc->ggecut, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); + if (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 + // qianrui add 2021-8-13 to make different kpar parameters can get the same + // results #endif this->pw_wfc->ft.fft_mode = inp.fft_mode; @@ -246,8 +249,7 @@ void ESolver_KS::before_all_runners(Input& inp, UnitCell& ucell) this->pw_wfc->setuptransform(); //! 9) initialize the number of plane waves for each k point - 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]; } @@ -272,9 +274,9 @@ void ESolver_KS::before_all_runners(Input& inp, UnitCell& ucell) CE.Init_CE(ucell.nat); #ifdef USE_PAW - if (GlobalV::use_paw) - { - GlobalC::paw_cell.set_libpaw_ecut(inp.ecutwfc / 2.0, inp.ecutwfc / 2.0); // in Hartree + if (GlobalV::use_paw) { + GlobalC::paw_cell.set_libpaw_ecut(inp.ecutwfc / 2.0, + inp.ecutwfc / 2.0); // in Hartree GlobalC::paw_cell.set_libpaw_fft(this->pw_wfc->nx, this->pw_wfc->ny, this->pw_wfc->nz, @@ -284,8 +286,7 @@ void ESolver_KS::before_all_runners(Input& inp, UnitCell& ucell) this->pw_wfc->startz, this->pw_wfc->numz); #ifdef __MPI - if (GlobalV::RANK_IN_POOL == 0) - { + if (GlobalV::RANK_IN_POOL == 0) { GlobalC::paw_cell.prepare_paw(); } #else @@ -304,12 +305,10 @@ void ESolver_KS::before_all_runners(Input& inp, UnitCell& ucell) 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(), @@ -320,8 +319,7 @@ void ESolver_KS::before_all_runners(Input& inp, UnitCell& ucell) #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(), @@ -338,23 +336,22 @@ void ESolver_KS::before_all_runners(Input& inp, UnitCell& ucell) //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -void ESolver_KS::init_after_vc(Input& inp, UnitCell& ucell) -{ +void ESolver_KS::init_after_vc(Input& inp, UnitCell& ucell) { ModuleBase::TITLE("ESolver_KS", "init_after_vc"); ESolver_FP::init_after_vc(inp, ucell); - if (GlobalV::md_prec_level == 2) - { + if (GlobalV::md_prec_level == 2) { // initialize the real-space uniform grid for FFT and parallel // distribution of plane waves - GlobalC::Pgrid.init(this->pw_rhod->nx, - this->pw_rhod->ny, - this->pw_rhod->nz, - this->pw_rhod->nplane, - this->pw_rhod->nrxx, - pw_big->nbz, - pw_big->bz); // mohan add 2010-07-22, update 2011-05-04 + GlobalC::Pgrid.init( + this->pw_rhod->nx, + this->pw_rhod->ny, + this->pw_rhod->nz, + this->pw_rhod->nplane, + this->pw_rhod->nrxx, + pw_big->nbz, + pw_big->bz); // mohan add 2010-07-22, update 2011-05-04 // Calculate Structure factor this->sf.setup_structure_factor(&ucell, this->pw_rhod); @@ -366,8 +363,9 @@ void ESolver_KS::init_after_vc(Input& inp, UnitCell& ucell) //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -void ESolver_KS::hamilt2density(const int istep, const int iter, const double ethr) -{ +void ESolver_KS::hamilt2density(const int istep, + const int iter, + const double ethr) { ModuleBase::timer::tick(this->classname, "hamilt2density"); // Temporarily, before HSolver is constructed, it should be overrided by // LCAO, PW, SDFT and TDDFT. @@ -382,49 +380,75 @@ void ESolver_KS::hamilt2density(const int istep, const int iter, cons //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -void ESolver_KS::print_wfcfft(Input& inp, std::ofstream& ofs) -{ +void ESolver_KS::print_wfcfft(Input& inp, std::ofstream& ofs) { ofs << "\n\n\n\n"; - ofs << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl; - ofs << " | |" << std::endl; - ofs << " | Setup plane waves of wave functions: |" << std::endl; - ofs << " | Use the energy cutoff and the lattice vectors to generate the |" << std::endl; - ofs << " | dimensions of FFT grid. The number of FFT grid on each processor |" << std::endl; - ofs << " | is 'nrxx'. The number of plane wave basis in reciprocal space is |" << std::endl; - ofs << " | different for charege/potential and wave functions. We also set |" << std::endl; - ofs << " | the 'sticks' for the parallel of FFT. The number of plane wave of |" << std::endl; - ofs << " | each k-point is 'npwk[ik]' in each processor |" << std::endl; - ofs << " <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << std::endl; + ofs << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" + ">>>>" + << std::endl; + ofs << " | " + " |" + << std::endl; + ofs << " | Setup plane waves of wave functions: " + " |" + << std::endl; + ofs << " | Use the energy cutoff and the lattice vectors to generate the " + " |" + << std::endl; + ofs << " | dimensions of FFT grid. The number of FFT grid on each " + "processor |" + << std::endl; + ofs << " | is 'nrxx'. The number of plane wave basis in reciprocal space " + "is |" + << std::endl; + ofs << " | different for charege/potential and wave functions. We also set " + " |" + << std::endl; + ofs << " | the 'sticks' for the parallel of FFT. The number of plane wave " + "of |" + << std::endl; + ofs << " | each k-point is 'npwk[ik]' in each processor " + " |" + << std::endl; + ofs << " <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" + "<<<<" + << std::endl; ofs << "\n\n\n\n"; ofs << "\n SETUP PLANE WAVES FOR WAVE FUNCTIONS" << std::endl; double ecut = inp.ecutwfc; - if (std::abs(ecut - this->pw_wfc->gk_ecut * this->pw_wfc->tpiba2) > 1e-6) - { + if (std::abs(ecut - this->pw_wfc->gk_ecut * this->pw_wfc->tpiba2) > 1e-6) { ecut = this->pw_wfc->gk_ecut * this->pw_wfc->tpiba2; - ofs << "Energy cutoff for wavefunc is incompatible with nx, ny, nz and it will be reduced!" << std::endl; + ofs << "Energy cutoff for wavefunc is incompatible with nx, ny, nz and " + "it will be reduced!" + << std::endl; } - ModuleBase::GlobalFunc::OUT(ofs, "energy cutoff for wavefunc (unit:Ry)", ecut); + ModuleBase::GlobalFunc::OUT(ofs, + "energy cutoff for wavefunc (unit:Ry)", + ecut); ModuleBase::GlobalFunc::OUT(ofs, "fft grid for wave functions", this->pw_wfc->nx, this->pw_wfc->ny, this->pw_wfc->nz); - ModuleBase::GlobalFunc::OUT(ofs, "number of plane waves", this->pw_wfc->npwtot); + ModuleBase::GlobalFunc::OUT(ofs, + "number of plane waves", + this->pw_wfc->npwtot); ModuleBase::GlobalFunc::OUT(ofs, "number of sticks", this->pw_wfc->nstot); ofs << "\n PARALLEL PW FOR WAVE FUNCTIONS" << std::endl; - ofs << " " << std::setw(8) << "PROC" << std::setw(15) << "COLUMNS(POT)" << std::setw(15) << "PW" << std::endl; + ofs << " " << std::setw(8) << "PROC" << std::setw(15) << "COLUMNS(POT)" + << std::setw(15) << "PW" << std::endl; - for (int i = 0; i < GlobalV::NPROC_IN_POOL; ++i) - { - ofs << " " << std::setw(8) << i + 1 << std::setw(15) << this->pw_wfc->nst_per[i] << std::setw(15) + for (int i = 0; i < GlobalV::NPROC_IN_POOL; ++i) { + ofs << " " << std::setw(8) << i + 1 << std::setw(15) + << this->pw_wfc->nst_per[i] << std::setw(15) << this->pw_wfc->npw_per[i] << std::endl; } ofs << " --------------- sum -------------------" << std::endl; - ofs << " " << std::setw(8) << GlobalV::NPROC_IN_POOL << std::setw(15) << this->pw_wfc->nstot << std::setw(15) - << this->pw_wfc->npwtot << std::endl; + ofs << " " << std::setw(8) << GlobalV::NPROC_IN_POOL << std::setw(15) + << this->pw_wfc->nstot << std::setw(15) << this->pw_wfc->npwtot + << std::endl; ModuleBase::GlobalFunc::DONE(ofs, "INIT PLANEWAVE"); } @@ -448,8 +472,7 @@ void ESolver_KS::print_wfcfft(Input& inp, std::ofstream& ofs) //! 16) Json again //------------------------------------------------------------------------------ template -void ESolver_KS::runner(const int istep, UnitCell& ucell) -{ +void ESolver_KS::runner(const int istep, UnitCell& ucell) { ModuleBase::TITLE("ESolver_KS", "runner"); ModuleBase::timer::tick(this->classname, "runner"); @@ -458,8 +481,7 @@ void ESolver_KS::runner(const int istep, UnitCell& ucell) this->before_scf(istep); // 3) write charge density - if (GlobalV::dm_to_rho) - { + if (GlobalV::dm_to_rho) { ModuleBase::timer::tick(this->classname, "runner"); return; // nothing further is needed } @@ -472,8 +494,7 @@ void ESolver_KS::runner(const int istep, UnitCell& ucell) // 4) SCF iterations std::cout << " * * * * * *\n << Start SCF iteration." << std::endl; - for (int iter = 1; iter <= this->maxniter; ++iter) - { + for (int iter = 1; iter <= this->maxniter; ++iter) { // 5) write head this->write_head(GlobalV::ofs_running, istep, iter); @@ -491,27 +512,31 @@ void ESolver_KS::runner(const int istep, UnitCell& ucell) this->hamilt2density(istep, iter, diag_ethr); // 8) for MPI: STOGROUP? need to rewrite - // It may be changed when more clever parallel algorithm is put forward. - // When parallel algorithm for bands are adopted. Density will only be treated in the first group. - //(Different ranks should have abtained the same, but small differences always exist in practice.) - // Maybe in the future, density and wavefunctions should use different parallel algorithms, in which - // they do not occupy all processors, for example wavefunctions uses 20 processors while density uses 10. - if (GlobalV::MY_STOGROUP == 0) - { + // It may be changed when more clever parallel algorithm is + //put forward. + // When parallel algorithm for bands are adopted. Density will only be + // treated in the first group. + //(Different ranks should have abtained the same, but small differences + //always exist in practice.) + // Maybe in the future, density and wavefunctions should use different + // parallel algorithms, in which they do not occupy all processors, for + // example wavefunctions uses 20 processors while density uses 10. + if (GlobalV::MY_STOGROUP == 0) { // double drho = this->estate.caldr2(); // EState should be used after it is constructed. drho = p_chgmix->get_drho(pelec->charge, GlobalV::nelec); double hsolver_error = 0.0; - if (firstscf) - { + if (firstscf) { firstscf = false; hsolver_error = this->phsol->cal_hsolerror(); // The error of HSolver is larger than drho, // so a more precise HSolver should be excuconv_elected. - if (hsolver_error > drho) - { - diag_ethr = this->phsol->reset_diagethr(GlobalV::ofs_running, hsolver_error, drho); + if (hsolver_error > drho) { + diag_ethr + = this->phsol->reset_diagethr(GlobalV::ofs_running, + hsolver_error, + drho); this->hamilt2density(istep, iter, diag_ethr); drho = p_chgmix->get_drho(pelec->charge, GlobalV::nelec); hsolver_error = this->phsol->cal_hsolerror(); @@ -519,48 +544,51 @@ void ESolver_KS::runner(const int istep, UnitCell& ucell) } // mixing will restart at this->p_chgmix->mixing_restart steps if (drho <= GlobalV::MIXING_RESTART && GlobalV::MIXING_RESTART > 0.0 - && this->p_chgmix->mixing_restart_step > iter) - { + && this->p_chgmix->mixing_restart_step > iter) { this->p_chgmix->mixing_restart_step = iter + 1; } - // drho will be 0 at this->p_chgmix->mixing_restart step, which is not ground state - bool not_restart_step = !(iter == this->p_chgmix->mixing_restart_step && GlobalV::MIXING_RESTART > 0.0); + // drho will be 0 at this->p_chgmix->mixing_restart step, which is + // not ground state + bool not_restart_step + = !(iter == this->p_chgmix->mixing_restart_step + && GlobalV::MIXING_RESTART > 0.0); // SCF will continue if U is not converged for uramping calculation bool is_U_converged = true; // to avoid unnecessary dependence on dft+u, refactor is needed #ifdef __LCAO - if (GlobalV::dft_plus_u) - { + if (GlobalV::dft_plus_u) { is_U_converged = GlobalC::dftu.u_converged(); } #endif - this->conv_elec = (drho < this->scf_thr && not_restart_step && is_U_converged); + this->conv_elec + = (drho < this->scf_thr && not_restart_step && is_U_converged); - // If drho < hsolver_error in the first iter or drho < scf_thr, we do not change rho. - if (drho < hsolver_error || this->conv_elec) - { - if (drho < hsolver_error) - { - GlobalV::ofs_warning << " drho < hsolver_error, keep charge density unchanged." << std::endl; + // If drho < hsolver_error in the first iter or drho < scf_thr, we + // do not change rho. + if (drho < hsolver_error || this->conv_elec) { + if (drho < hsolver_error) { + GlobalV::ofs_warning << " drho < hsolver_error, keep " + "charge density unchanged." + << std::endl; } - } - else - { + } else { //----------charge mixing--------------- - // mixing will restart after this->p_chgmix->mixing_restart steps - if (GlobalV::MIXING_RESTART > 0 && iter == this->p_chgmix->mixing_restart_step - 1) - { + // mixing will restart after this->p_chgmix->mixing_restart + // steps + if (GlobalV::MIXING_RESTART > 0 + && iter == this->p_chgmix->mixing_restart_step - 1 + && drho <= GlobalV::MIXING_RESTART) { // do not mix charge density + } else { + p_chgmix->mix_rho( + pelec->charge); // update chr->rho by mixing } - else - { - p_chgmix->mix_rho(pelec->charge); // update chr->rho by mixing - } - if (GlobalV::SCF_THR_TYPE == 2) - { - pelec->charge->renormalize_rho(); // renormalize rho in R-space would induce a error in K-space + if (GlobalV::SCF_THR_TYPE == 2) { + pelec->charge + ->renormalize_rho(); // renormalize rho in R-space would + // induce a error in K-space } //----------charge mixing done----------- } @@ -568,7 +596,11 @@ void ESolver_KS::runner(const int istep, UnitCell& ucell) #ifdef __MPI MPI_Bcast(&drho, 1, MPI_DOUBLE, 0, PARAPW_WORLD); MPI_Bcast(&this->conv_elec, 1, MPI_DOUBLE, 0, PARAPW_WORLD); - MPI_Bcast(pelec->charge->rho[0], this->pw_rhod->nrxx, MPI_DOUBLE, 0, PARAPW_WORLD); + MPI_Bcast(pelec->charge->rho[0], + this->pw_rhod->nrxx, + MPI_DOUBLE, + 0, + PARAPW_WORLD); #endif // 9) update potential @@ -582,15 +614,16 @@ void ESolver_KS::runner(const int istep, UnitCell& ucell) double duration = (double)(MPI_Wtime() - iterstart); #else double duration - = (std::chrono::duration_cast(std::chrono::system_clock::now() - iterstart)) + = (std::chrono::duration_cast( + std::chrono::system_clock::now() - iterstart)) .count() / static_cast(1e6); #endif // 11) get mtaGGA related parameters double dkin = 0.0; // for meta-GGA - 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) { dkin = p_chgmix->get_dkin(pelec->charge, GlobalV::nelec); } this->print_iter(iter, drho, dkin, duration, diag_ethr); @@ -601,25 +634,25 @@ void ESolver_KS::runner(const int istep, UnitCell& ucell) Json::add_output_scf_mag(GlobalC::ucell.magnet.tot_magnetization, GlobalC::ucell.magnet.abs_magnetization, this->pelec->f_en.etot * ModuleBase::Ry_to_eV, - this->pelec->f_en.etot_delta * ModuleBase::Ry_to_eV, + this->pelec->f_en.etot_delta + * ModuleBase::Ry_to_eV, drho, duration); #endif //__RAPIDJSON // 13) check convergence - if (this->conv_elec) - { + if (this->conv_elec) { this->niter = iter; bool stop = this->do_after_converge(iter); - if (stop) - { + if (stop) { break; } } // notice for restart - if (GlobalV::MIXING_RESTART > 0 && iter == this->p_chgmix->mixing_restart_step - 1 && iter != GlobalV::SCF_NMAX) - { + if (GlobalV::MIXING_RESTART > 0 + && iter == this->p_chgmix->mixing_restart_step - 1 + && iter != GlobalV::SCF_NMAX) { std::cout << " SCF restart after this step!" << std::endl; } } // end scf iterations @@ -627,7 +660,9 @@ void ESolver_KS::runner(const int istep, UnitCell& ucell) #ifdef __RAPIDJSON // 14) add Json of efermi energy converge - Json::add_output_efermi_converge(this->pelec->eferm.ef * ModuleBase::Ry_to_eV, this->conv_elec); + Json::add_output_efermi_converge(this->pelec->eferm.ef + * ModuleBase::Ry_to_eV, + this->conv_elec); #endif //__RAPIDJSON // 15) after scf @@ -649,12 +684,10 @@ void ESolver_KS::runner(const int istep, UnitCell& ucell) //! mohan add 2024-05-12 //------------------------------------------------------------------------------ template -void ESolver_KS::print_head(void) -{ +void ESolver_KS::print_head() { std::cout << " " << std::setw(7) << "ITER"; - if (GlobalV::NSPIN == 2) - { + if (GlobalV::NSPIN == 2) { std::cout << std::setw(10) << "TMAG"; std::cout << std::setw(10) << "AMAG"; } @@ -663,8 +696,8 @@ void ESolver_KS::print_head(void) std::cout << std::setw(15) << "EDIFF(eV)"; std::cout << std::setw(11) << "DRHO"; - 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) { std::cout << std::setw(11) << "DKIN"; } @@ -680,9 +713,14 @@ void ESolver_KS::print_iter(const int iter, const double drho, const double dkin, const double duration, - const double ethr) -{ - this->pelec->print_etot(this->conv_elec, iter, drho, dkin, duration, INPUT.printe, ethr); + const double ethr) { + this->pelec->print_etot(this->conv_elec, + iter, + drho, + dkin, + duration, + INPUT.printe, + ethr); } //------------------------------------------------------------------------------ @@ -690,10 +728,13 @@ void ESolver_KS::print_iter(const int iter, //! mohan add 2024-05-12 //------------------------------------------------------------------------------ template -void ESolver_KS::write_head(std::ofstream& ofs_running, const int istep, const int iter) -{ - ofs_running << "\n " << this->basisname << " ALGORITHM --------------- ION=" << std::setw(4) << istep + 1 - << " ELEC=" << std::setw(4) << iter << "--------------------------------\n"; +void ESolver_KS::write_head(std::ofstream& ofs_running, + const int istep, + const int iter) { + ofs_running << "\n " << this->basisname + << " ALGORITHM --------------- ION=" << std::setw(4) + << istep + 1 << " ELEC=" << std::setw(4) << iter + << "--------------------------------\n"; } //------------------------------------------------------------------------------ @@ -701,8 +742,7 @@ void ESolver_KS::write_head(std::ofstream& ofs_running, const int ist //! mohan add 2024-05-12 //------------------------------------------------------------------------------ template -int ESolver_KS::get_niter() -{ +int ESolver_KS::get_niter() { return this->niter; } @@ -711,8 +751,7 @@ int ESolver_KS::get_niter() //! tqzhao add 2024-05-15 //------------------------------------------------------------------------------ template -int ESolver_KS::get_maxniter() -{ +int ESolver_KS::get_maxniter() { return this->maxniter; } @@ -721,8 +760,7 @@ int ESolver_KS::get_maxniter() //! tqzhao add 2024-05-15 //------------------------------------------------------------------------------ template -bool ESolver_KS::get_conv_elec() -{ +bool ESolver_KS::get_conv_elec() { return this->conv_elec; } @@ -731,12 +769,13 @@ bool ESolver_KS::get_conv_elec() //! mohan add 2024-05-12 //------------------------------------------------------------------------------ template -ModuleIO::Output_Rho ESolver_KS::create_Output_Rho(int is, int iter, const std::string& prefix) -{ +ModuleIO::Output_Rho + ESolver_KS::create_Output_Rho(int is, + int iter, + const std::string& prefix) { const int precision = 3; std::string tag = "CHG"; - if (GlobalV::dm_to_rho) - { + if (GlobalV::dm_to_rho) { return ModuleIO::Output_Rho(this->pw_big, this->pw_rhod, is, @@ -769,8 +808,10 @@ ModuleIO::Output_Rho ESolver_KS::create_Output_Rho(int is, int iter, //! mohan add 2024-05-12 //------------------------------------------------------------------------------ template -ModuleIO::Output_Rho ESolver_KS::create_Output_Kin(int is, int iter, const std::string& prefix) -{ +ModuleIO::Output_Rho + ESolver_KS::create_Output_Kin(int is, + int iter, + const std::string& prefix) { const int precision = 11; std::string tag = "TAU"; return ModuleIO::Output_Rho(this->pw_big, @@ -792,8 +833,9 @@ ModuleIO::Output_Rho ESolver_KS::create_Output_Kin(int is, int iter, //! mohan add 2024-05-12 //------------------------------------------------------------------------------ template -ModuleIO::Output_Potential ESolver_KS::create_Output_Potential(int iter, const std::string& prefix) -{ +ModuleIO::Output_Potential + ESolver_KS::create_Output_Potential(int iter, + const std::string& prefix) { const int precision = 3; std::string tag = "POT"; return ModuleIO::Output_Potential(this->pw_big, diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 682eb2c3d6..0a8ac41f9e 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -44,34 +44,35 @@ //--------------------------------------------------- #include "module_hamilt_lcao/module_deltaspin/spin_constrain.h" +#include "module_io/io_dmk.h" #include "module_io/write_dmr.h" #include "module_io/write_wfc_nao.h" -namespace ModuleESolver -{ +namespace ModuleESolver { //------------------------------------------------------------------------------ //! the 1st function of ESolver_KS_LCAO: constructor //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -ESolver_KS_LCAO::ESolver_KS_LCAO() -{ +ESolver_KS_LCAO::ESolver_KS_LCAO() { this->classname = "ESolver_KS_LCAO"; this->basisname = "LCAO"; // the following EXX code should be removed to other places, mohan 2024/05/11 #ifdef __EXX - if (GlobalC::exx_info.info_ri.real_number) - { - this->exx_lri_double = std::make_shared>(GlobalC::exx_info.info_ri); - this->exd = std::make_shared>(this->exx_lri_double); + if (GlobalC::exx_info.info_ri.real_number) { + this->exx_lri_double + = std::make_shared>(GlobalC::exx_info.info_ri); + this->exd = std::make_shared>( + this->exx_lri_double); this->LM.Hexxd = &this->exd->get_Hexxs(); - } - else - { - this->exx_lri_complex = std::make_shared>>(GlobalC::exx_info.info_ri); - this->exc = std::make_shared>>(this->exx_lri_complex); + } else { + this->exx_lri_complex = std::make_shared>>( + GlobalC::exx_info.info_ri); + this->exc + = std::make_shared>>( + this->exx_lri_complex); this->LM.Hexxc = &this->exc->get_Hexxs(); } #endif @@ -82,9 +83,7 @@ ESolver_KS_LCAO::ESolver_KS_LCAO() //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -ESolver_KS_LCAO::~ESolver_KS_LCAO() -{ -} +ESolver_KS_LCAO::~ESolver_KS_LCAO() {} //------------------------------------------------------------------------------ //! the 3rd function of ESolver_KS_LCAO: init @@ -105,50 +104,53 @@ ESolver_KS_LCAO::~ESolver_KS_LCAO() //! 14) set occupations? //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::before_all_runners(Input& inp, UnitCell& ucell) -{ +void ESolver_KS_LCAO::before_all_runners(Input& inp, UnitCell& ucell) { ModuleBase::TITLE("ESolver_KS_LCAO", "before_all_runners"); ModuleBase::timer::tick("ESolver_KS_LCAO", "before_all_runners"); // 1) calculate overlap matrix S - if (GlobalV::CALCULATION == "get_S") - { + if (GlobalV::CALCULATION == "get_S") { // 1.1) read pseudopotentials ucell.read_pseudo(GlobalV::ofs_running); // 1.2) symmetrize things - if (ModuleSymmetry::Symmetry::symm_flag == 1) - { - ucell.symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running); + if (ModuleSymmetry::Symmetry::symm_flag == 1) { + ucell.symm.analy_sys(ucell.lat, + ucell.st, + ucell.atoms, + GlobalV::ofs_running); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY"); } // 1.3) Setup k-points according to symmetry. - this->kv - .set(ucell.symm, GlobalV::global_kpoint_card, GlobalV::NSPIN, ucell.G, ucell.latvec, GlobalV::ofs_running); + this->kv.set(ucell.symm, + GlobalV::global_kpoint_card, + GlobalV::NSPIN, + ucell.G, + ucell.latvec, + GlobalV::ofs_running); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS"); Print_Info::setup_parameters(ucell, this->kv); - } - else - { + } else { // 1) else, call before_all_runners() in ESolver_KS ESolver_KS::before_all_runners(inp, ucell); } // end ifnot get_S // 2) init ElecState - // autoset nbands in ElecState, it should before basis_init (for Psi 2d divid) - if (this->pelec == nullptr) - { + // autoset nbands in ElecState, it should before basis_init (for Psi 2d + // divid) + if (this->pelec == nullptr) { // TK stands for double and complex? - this->pelec = new elecstate::ElecStateLCAO(&(this->chr), // use which parameter? - &(this->kv), - this->kv.get_nks(), - &(this->LOC), // use which parameter? - &(this->GG), // mohan add 2024-04-01 - &(this->GK), // mohan add 2024-04-01 - this->pw_rho, - this->pw_big); + this->pelec = new elecstate::ElecStateLCAO( + &(this->chr), // use which parameter? + &(this->kv), + this->kv.get_nks(), + &(this->LOC), // use which parameter? + &(this->GG), // mohan add 2024-04-01 + &(this->GK), // mohan add 2024-04-01 + this->pw_rho, + this->pw_big); } // 3) init LCAO basis @@ -166,8 +168,7 @@ void ESolver_KS_LCAO::before_all_runners(Input& inp, UnitCell& ucell) ->init_DM(&this->kv, &(this->orb_con.ParaV), GlobalV::NSPIN); // this function should be removed outside of the function - if (GlobalV::CALCULATION == "get_S") - { + if (GlobalV::CALCULATION == "get_S") { ModuleBase::timer::tick("ESolver_KS_LCAO", "init"); return; } @@ -175,36 +176,33 @@ void ESolver_KS_LCAO::before_all_runners(Input& inp, UnitCell& ucell) // 6) initialize Hamilt in LCAO // * allocate H and S matrices according to computational resources // * set the 'trace' between local H/S and global H/S - this->LM.divide_HS_in_frag(GlobalV::GAMMA_ONLY_LOCAL, orb_con.ParaV, this->kv.get_nks()); + this->LM.divide_HS_in_frag(GlobalV::GAMMA_ONLY_LOCAL, + orb_con.ParaV, + this->kv.get_nks()); #ifdef __EXX // 7) initialize exx // PLEASE simplify the Exx_Global interface - if (GlobalV::CALCULATION == "scf" || GlobalV::CALCULATION == "relax" || GlobalV::CALCULATION == "cell-relax" - || GlobalV::CALCULATION == "md") - { - if (GlobalC::exx_info.info_global.cal_exx) - { + if (GlobalV::CALCULATION == "scf" || GlobalV::CALCULATION == "relax" + || GlobalV::CALCULATION == "cell-relax" + || GlobalV::CALCULATION == "md") { + if (GlobalC::exx_info.info_global.cal_exx) { /* In the special "two-level" calculation case, - first scf iteration only calculate the functional without exact exchange. - but in "nscf" calculation, there is no need of "two-level" method. */ - if (ucell.atoms[0].ncpp.xc_func == "HF" || ucell.atoms[0].ncpp.xc_func == "PBE0" - || ucell.atoms[0].ncpp.xc_func == "HSE") - { + first scf iteration only calculate the functional without exact + exchange. but in "nscf" calculation, there is no need of "two-level" + method. */ + if (ucell.atoms[0].ncpp.xc_func == "HF" + || ucell.atoms[0].ncpp.xc_func == "PBE0" + || ucell.atoms[0].ncpp.xc_func == "HSE") { XC_Functional::set_xc_type("pbe"); - } - else if (ucell.atoms[0].ncpp.xc_func == "SCAN0") - { + } else if (ucell.atoms[0].ncpp.xc_func == "SCAN0") { XC_Functional::set_xc_type("scan"); } // GlobalC::exx_lcao.init(); - if (GlobalC::exx_info.info_ri.real_number) - { + if (GlobalC::exx_info.info_ri.real_number) { this->exx_lri_double->init(MPI_COMM_WORLD, this->kv); - } - else - { + } else { this->exx_lri_complex->init(MPI_COMM_WORLD, this->kv); } } @@ -212,8 +210,7 @@ void ESolver_KS_LCAO::before_all_runners(Input& inp, UnitCell& ucell) #endif // 8) initialize DFT+U - if (GlobalV::dft_plus_u) - { + if (GlobalV::dft_plus_u) { GlobalC::dftu.init(ucell, this->LM, this->kv.get_nks()); } @@ -221,8 +218,7 @@ void ESolver_KS_LCAO::before_all_runners(Input& inp, UnitCell& ucell) GlobalC::ppcell.init_vloc(GlobalC::ppcell.vloc, this->pw_rho); // 10) initialize the HSolver - if (this->phsol == nullptr) - { + if (this->phsol == nullptr) { this->phsol = new hsolver::HSolverLCAO(&(this->orb_con.ParaV)); this->phsol->method = GlobalV::KS_SOLVER; } @@ -232,8 +228,7 @@ void ESolver_KS_LCAO::before_all_runners(Input& inp, UnitCell& ucell) this->pelec->omega = GlobalC::ucell.omega; // 12) 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, &GlobalC::ucell, @@ -245,17 +240,17 @@ void ESolver_KS_LCAO::before_all_runners(Input& inp, UnitCell& ucell) #ifdef __DEEPKS // 13) initialize deepks - if (GlobalV::deepks_scf) - { + if (GlobalV::deepks_scf) { // load the DeePKS model from deep neural network GlobalC::ld.load_model(INPUT.deepks_model); } #endif // 14) set occupations - if (GlobalV::ocp) - { - this->pelec->fixed_weights(GlobalV::ocp_kb, GlobalV::NBANDS, GlobalV::nelec); + if (GlobalV::ocp) { + this->pelec->fixed_weights(GlobalV::ocp_kb, + GlobalV::NBANDS, + GlobalV::nelec); } ModuleBase::timer::tick("ESolver_KS_LCAO", "before_all_runners"); @@ -267,25 +262,25 @@ void ESolver_KS_LCAO::before_all_runners(Input& inp, UnitCell& ucell) //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::init_after_vc(Input& inp, UnitCell& ucell) -{ +void ESolver_KS_LCAO::init_after_vc(Input& inp, UnitCell& ucell) { ModuleBase::TITLE("ESolver_KS_LCAO", "init_after_vc"); ModuleBase::timer::tick("ESolver_KS_LCAO", "init_after_vc"); ESolver_KS::init_after_vc(inp, ucell); - if (GlobalV::md_prec_level == 2) - { + if (GlobalV::md_prec_level == 2) { delete this->pelec; - this->pelec = new elecstate::ElecStateLCAO(&(this->chr), - &(this->kv), - this->kv.get_nks(), - &(this->LOC), - &(this->GG), // mohan add 2024-04-01 - &(this->GK), // mohan add 2024-04-01 - this->pw_rho, - this->pw_big); - - dynamic_cast*>(this->pelec)->init_DM(&this->kv, this->LM.ParaV, GlobalV::NSPIN); + this->pelec = new elecstate::ElecStateLCAO( + &(this->chr), + &(this->kv), + this->kv.get_nks(), + &(this->LOC), + &(this->GG), // mohan add 2024-04-01 + &(this->GK), // mohan add 2024-04-01 + this->pw_rho, + this->pw_big); + + dynamic_cast*>(this->pelec) + ->init_DM(&this->kv, this->LM.ParaV, GlobalV::NSPIN); GlobalC::ppcell.init_vloc(GlobalC::ppcell.vloc, this->pw_rho); @@ -293,15 +288,15 @@ void ESolver_KS_LCAO::init_after_vc(Input& inp, UnitCell& ucell) this->pelec->omega = GlobalC::ucell.omega; // Initialize the potential. - if (this->pelec->pot == nullptr) - { - this->pelec->pot = new elecstate::Potential(this->pw_rhod, - this->pw_rho, - &GlobalC::ucell, - &(GlobalC::ppcell.vloc), - &(this->sf), - &(this->pelec->f_en.etxc), - &(this->pelec->f_en.vtxc)); + if (this->pelec->pot == nullptr) { + this->pelec->pot + = new elecstate::Potential(this->pw_rhod, + this->pw_rho, + &GlobalC::ucell, + &(GlobalC::ppcell.vloc), + &(this->sf), + &(this->pelec->f_en.etxc), + &(this->pelec->f_en.vtxc)); } } @@ -314,8 +309,7 @@ void ESolver_KS_LCAO::init_after_vc(Input& inp, UnitCell& ucell) //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -double ESolver_KS_LCAO::cal_energy() -{ +double ESolver_KS_LCAO::cal_energy() { ModuleBase::TITLE("ESolver_KS_LCAO", "cal_energy"); return this->pelec->f_en.etot; @@ -326,8 +320,7 @@ double ESolver_KS_LCAO::cal_energy() //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::cal_force(ModuleBase::matrix& force) -{ +void ESolver_KS_LCAO::cal_force(ModuleBase::matrix& force) { ModuleBase::TITLE("ESolver_KS_LCAO", "cal_force"); ModuleBase::timer::tick("ESolver_KS_LCAO", "cal_force"); @@ -369,13 +362,11 @@ void ESolver_KS_LCAO::cal_force(ModuleBase::matrix& force) //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::cal_stress(ModuleBase::matrix& stress) -{ +void ESolver_KS_LCAO::cal_stress(ModuleBase::matrix& stress) { ModuleBase::TITLE("ESolver_KS_LCAO", "cal_stress"); ModuleBase::timer::tick("ESolver_KS_LCAO", "cal_stress"); - if (!this->have_force) - { + if (!this->have_force) { ModuleBase::matrix fcs; this->cal_force(fcs); } @@ -390,45 +381,68 @@ void ESolver_KS_LCAO::cal_stress(ModuleBase::matrix& stress) //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::after_all_runners() -{ +void ESolver_KS_LCAO::after_all_runners() { ModuleBase::TITLE("ESolver_KS_LCAO", "after_all_runners"); ModuleBase::timer::tick("ESolver_KS_LCAO", "after_all_runners"); - GlobalV::ofs_running << "\n\n --------------------------------------------" << std::endl; + 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; - - if (INPUT.out_dos != 0 || INPUT.out_band[0] != 0 || INPUT.out_proj_band != 0) - { + 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 (INPUT.out_dos != 0 || INPUT.out_band[0] != 0 + || INPUT.out_proj_band != 0) { GlobalV::ofs_running << "\n\n\n\n"; - GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl; - GlobalV::ofs_running << " | |" << std::endl; - GlobalV::ofs_running << " | Post-processing of data: |" << std::endl; - GlobalV::ofs_running << " | DOS (density of states) and bands will be output here. |" << std::endl; - GlobalV::ofs_running << " | If atomic orbitals are used, Mulliken charge analysis can be done. |" << std::endl; - GlobalV::ofs_running << " | Also the .bxsf file containing fermi surface information can be |" << std::endl; - GlobalV::ofs_running << " | done here. |" << std::endl; - GlobalV::ofs_running << " | |" << std::endl; - GlobalV::ofs_running << " <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << std::endl; + GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" + ">>>>>>>>>>>>>>>>>>>>>>>>>" + << std::endl; + GlobalV::ofs_running << " | " + " |" + << std::endl; + GlobalV::ofs_running << " | Post-processing of data: " + " |" + << std::endl; + GlobalV::ofs_running << " | DOS (density of states) and bands will be " + "output here. |" + << std::endl; + GlobalV::ofs_running << " | If atomic orbitals are used, Mulliken " + "charge analysis can be done. |" + << std::endl; + GlobalV::ofs_running << " | Also the .bxsf file containing fermi " + "surface information can be |" + << std::endl; + GlobalV::ofs_running << " | done here. " + " |" + << std::endl; + GlobalV::ofs_running << " | " + " |" + << std::endl; + GlobalV::ofs_running << " <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" + "<<<<<<<<<<<<<<<<<<<<<<<<<" + << std::endl; GlobalV::ofs_running << "\n\n\n\n"; } // qianrui modify 2020-10-18 - if (GlobalV::CALCULATION == "scf" || GlobalV::CALCULATION == "md" || GlobalV::CALCULATION == "relax") - { - ModuleIO::write_istate_info(this->pelec->ekb, this->pelec->wg, this->kv, &(GlobalC::Pkpoints)); + if (GlobalV::CALCULATION == "scf" || GlobalV::CALCULATION == "md" + || GlobalV::CALCULATION == "relax") { + ModuleIO::write_istate_info(this->pelec->ekb, + this->pelec->wg, + this->kv, + &(GlobalC::Pkpoints)); } const int nspin0 = (GlobalV::NSPIN == 2) ? 2 : 1; - if (INPUT.out_band[0]) - { - for (int is = 0; is < nspin0; is++) - { + if (INPUT.out_band[0]) { + 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, @@ -442,11 +456,15 @@ void ESolver_KS_LCAO::after_all_runners() if (INPUT.out_proj_band) // Projeced band structure added by jiyy-2022-4-20 { - ModuleIO::write_proj_band_lcao(this->psi, this->LM, this->pelec, this->kv, GlobalC::ucell, this->p_hamilt); + ModuleIO::write_proj_band_lcao(this->psi, + this->LM, + this->pelec, + this->kv, + GlobalC::ucell, + this->p_hamilt); } - if (INPUT.out_dos) - { + if (INPUT.out_dos) { ModuleIO::out_dos_nao(this->psi, this->LM, this->orb_con.ParaV, @@ -463,8 +481,7 @@ void ESolver_KS_LCAO::after_all_runners() this->p_hamilt); } - if (INPUT.out_mat_xc) - { + if (INPUT.out_mat_xc) { ModuleIO::write_Vxc(GlobalV::NSPIN, GlobalV::NLOCAL, GlobalV::DRANK, @@ -491,23 +508,20 @@ void ESolver_KS_LCAO::after_all_runners() //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::init_basis_lcao(ORB_control& orb_con, Input& inp, UnitCell& ucell) -{ +void ESolver_KS_LCAO::init_basis_lcao(ORB_control& orb_con, + Input& inp, + UnitCell& ucell) { ModuleBase::TITLE("ESolver_KS_LCAO", "init_basis_lcao"); // autoset NB2D first - if (GlobalV::NB2D == 0) - { - if (GlobalV::NLOCAL > 0) - { + if (GlobalV::NB2D == 0) { + if (GlobalV::NLOCAL > 0) { GlobalV::NB2D = (GlobalV::NSPIN == 4) ? 2 : 1; } - if (GlobalV::NLOCAL > 500) - { + if (GlobalV::NLOCAL > 500) { GlobalV::NB2D = 32; } - if (GlobalV::NLOCAL > 1000) - { + if (GlobalV::NLOCAL > 1000) { GlobalV::NB2D = 64; } } @@ -529,15 +543,24 @@ void ESolver_KS_LCAO::init_basis_lcao(ORB_control& orb_con, Input& inp, // * construct the interpolation tables. two_center_bundle_.build_orb(ucell.ntype, ucell.orbital_fn); - two_center_bundle_.build_alpha(GlobalV::deepks_setorb, &ucell.descriptor_file); + two_center_bundle_.build_alpha(GlobalV::deepks_setorb, + &ucell.descriptor_file); two_center_bundle_.build_orb_onsite(ucell.ntype, GlobalV::onsite_radius); - // currently deepks only use one descriptor file, so cast bool to int is fine + // currently deepks only use one descriptor file, so cast bool to int is + // fine // TODO Due to the omnipresence of GlobalC::ORB, we still have to rely // on the old interface for now. - two_center_bundle_.to_LCAO_Orbitals(GlobalC::ORB, inp.lcao_ecut, inp.lcao_dk, inp.lcao_dr, inp.lcao_rmax); + two_center_bundle_.to_LCAO_Orbitals(GlobalC::ORB, + inp.lcao_ecut, + inp.lcao_dk, + inp.lcao_dr, + inp.lcao_rmax); - ucell.infoNL.setupNonlocal(ucell.ntype, ucell.atoms, GlobalV::ofs_running, GlobalC::ORB); + ucell.infoNL.setupNonlocal(ucell.ntype, + ucell.atoms, + GlobalV::ofs_running, + GlobalC::ORB); two_center_bundle_.build_beta(ucell.ntype, ucell.infoNL.Beta); @@ -549,13 +572,18 @@ void ESolver_KS_LCAO::init_basis_lcao(ORB_control& orb_con, Input& inp, #ifdef USE_NEW_TWO_CENTER two_center_bundle_.tabulate(); #else - two_center_bundle_.tabulate(inp.lcao_ecut, inp.lcao_dk, inp.lcao_dr, inp.lcao_rmax); + two_center_bundle_.tabulate(inp.lcao_ecut, + inp.lcao_dk, + inp.lcao_dr, + inp.lcao_rmax); #endif - if (this->orb_con.setup_2d) - { - this->orb_con.setup_2d_division(GlobalV::ofs_running, GlobalV::ofs_warning); - this->orb_con.ParaV.set_atomic_trace(GlobalC::ucell.get_iat2iwt(), GlobalC::ucell.nat, GlobalV::NLOCAL); + if (this->orb_con.setup_2d) { + this->orb_con.setup_2d_division(GlobalV::ofs_running, + GlobalV::ofs_warning); + this->orb_con.ParaV.set_atomic_trace(GlobalC::ucell.get_iat2iwt(), + GlobalC::ucell.nat, + GlobalV::NLOCAL); } return; @@ -566,45 +594,39 @@ void ESolver_KS_LCAO::init_basis_lcao(ORB_control& orb_con, Input& inp, //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::iter_init(const int istep, const int iter) -{ +void ESolver_KS_LCAO::iter_init(const int istep, const int iter) { ModuleBase::TITLE("ESolver_KS_LCAO", "iter_init"); - if (iter == 1) - { + if (iter == 1) { this->p_chgmix->init_mixing(); // init mixing this->p_chgmix->mixing_restart_step = GlobalV::SCF_NMAX + 1; this->p_chgmix->mixing_restart_count = 0; // this output will be removed once the feeature is stable - if (GlobalC::dftu.uramping > 0.01) - { + if (GlobalC::dftu.uramping > 0.01) { std::cout << " U-Ramping! Current U = "; - for (int i = 0; i < GlobalC::dftu.U0.size(); i++) - { + for (int i = 0; i < GlobalC::dftu.U0.size(); i++) { std::cout << GlobalC::dftu.U[i] * ModuleBase::Ry_to_eV << " "; } std::cout << " eV " << std::endl; } } // 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(); this->p_chgmix->mixing_restart_count++; - if (GlobalV::dft_plus_u) - { - GlobalC::dftu.uramping_update(); // update U by uramping if uramping > 0.01 - if (GlobalC::dftu.uramping > 0.01) - { + if (GlobalV::dft_plus_u) { + GlobalC::dftu + .uramping_update(); // update U by uramping if uramping > 0.01 + if (GlobalC::dftu.uramping > 0.01) { std::cout << " U-Ramping! Current U = "; - for (int i = 0; i < GlobalC::dftu.U0.size(); i++) - { - std::cout << GlobalC::dftu.U[i] * ModuleBase::Ry_to_eV << " "; + for (int i = 0; i < GlobalC::dftu.U0.size(); i++) { + std::cout << GlobalC::dftu.U[i] * ModuleBase::Ry_to_eV + << " "; } std::cout << " eV " << std::endl; } - if (GlobalC::dftu.uramping > 0.01 && !GlobalC::dftu.u_converged()) - { + if (GlobalC::dftu.uramping > 0.01 && !GlobalC::dftu.u_converged()) { this->p_chgmix->mixing_restart_step = GlobalV::SCF_NMAX + 1; } } @@ -612,7 +634,8 @@ void ESolver_KS_LCAO::iter_init(const int istep, const int iter) { // allocate memory for dmr_mdata const elecstate::DensityMatrix* dm - = dynamic_cast*>(this->pelec)->get_DM(); + = dynamic_cast*>(this->pelec) + ->get_DM(); int nnr_tmp = dm->get_DMR_pointer(1)->get_nnr(); this->p_chgmix->allocate_mixing_dmr(nnr_tmp); } @@ -624,10 +647,8 @@ void ESolver_KS_LCAO::iter_init(const int istep, const int iter) // mohan move it outside 2011-01-13 // first need to calculate the weight according to // electrons number. - if (istep == 0 && this->wf.init_wfc == "file") - { - if (iter == 1) - { + if (istep == 0 && this->wf.init_wfc == "file") { + if (iter == 1) { std::cout << " WAVEFUN -> CHARGE " << std::endl; // calculate the density matrix using read in wave functions @@ -659,13 +680,13 @@ void ESolver_KS_LCAO::iter_init(const int istep, const int iter) // rho1 and rho2 are the same rho. // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - if (GlobalV::NSPIN == 4) - { + if (GlobalV::NSPIN == 4) { GlobalC::ucell.cal_ux(); } //! update the potentials by using new electron charge density - this->pelec->pot->update_from_charge(this->pelec->charge, &GlobalC::ucell); + this->pelec->pot->update_from_charge(this->pelec->charge, + &GlobalC::ucell); //! compute the correction energy for metals this->pelec->f_en.descf = this->pelec->cal_delta_escf(); @@ -674,24 +695,28 @@ void ESolver_KS_LCAO::iter_init(const int istep, const int iter) #ifdef __EXX // calculate exact-exchange - if (GlobalC::exx_info.info_ri.real_number) - { - this->exd->exx_eachiterinit(*dynamic_cast*>(this->pelec)->get_DM(), iter); - } - else - { - this->exc->exx_eachiterinit(*dynamic_cast*>(this->pelec)->get_DM(), iter); + if (GlobalC::exx_info.info_ri.real_number) { + this->exd->exx_eachiterinit( + *dynamic_cast*>(this->pelec) + ->get_DM(), + iter); + } else { + this->exc->exx_eachiterinit( + *dynamic_cast*>(this->pelec) + ->get_DM(), + iter); } #endif - if (GlobalV::dft_plus_u) - { - if (istep != 0 || iter != 1) - { - GlobalC::dftu.set_dmr(dynamic_cast*>(this->pelec)->get_DM()); + if (GlobalV::dft_plus_u) { + if (istep != 0 || iter != 1) { + GlobalC::dftu.set_dmr( + dynamic_cast*>(this->pelec) + ->get_DM()); } // Calculate U and J if Yukawa potential is used - GlobalC::dftu.cal_slater_UJ(this->pelec->charge->rho, this->pw_rho->nrxx); + GlobalC::dftu.cal_slater_UJ(this->pelec->charge->rho, + this->pw_rho->nrxx); } #ifdef __DEEPKS @@ -699,27 +724,25 @@ void ESolver_KS_LCAO::iter_init(const int istep, const int iter) GlobalC::ld.set_hr_cal(true); // HR in HamiltLCAO should be recalculate - if (GlobalV::deepks_scf) - { + if (GlobalV::deepks_scf) { this->p_hamilt->refresh(); } #endif - if (GlobalV::VL_IN_H) - { + if (GlobalV::VL_IN_H) { // update Gint_K - if (!GlobalV::GAMMA_ONLY_LOCAL) - { + if (!GlobalV::GAMMA_ONLY_LOCAL) { this->GK.renew(); } // update real space Hamiltonian this->p_hamilt->refresh(); } - // run the inner lambda loop to contrain atomic moments with the DeltaSpin method - if (GlobalV::sc_mag_switch && iter > GlobalV::sc_scf_nmin) - { - SpinConstrain& sc = SpinConstrain::getScInstance(); + // run the inner lambda loop to contrain atomic moments with the DeltaSpin + // method + if (GlobalV::sc_mag_switch && iter > GlobalV::sc_scf_nmin) { + SpinConstrain& sc + = SpinConstrain::getScInstance(); sc.run_lambda_loop(iter - 1); } } @@ -741,75 +764,69 @@ void ESolver_KS_LCAO::iter_init(const int istep, const int iter) //! 12) calculate delta energy //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::hamilt2density(int istep, int iter, double ethr) -{ +void ESolver_KS_LCAO::hamilt2density(int istep, int iter, double ethr) { ModuleBase::TITLE("ESolver_KS_LCAO", "hamilt2density"); // 1) save input rho this->pelec->charge->save_rho_before_sum_band(); // 2) save density matrix DMR for mixing - if (GlobalV::MIXING_RESTART > 0 && GlobalV::MIXING_DMR && this->p_chgmix->mixing_restart_count > 0) - { - elecstate::DensityMatrix* dm = dynamic_cast*>(this->pelec)->get_DM(); + if (GlobalV::MIXING_RESTART > 0 && GlobalV::MIXING_DMR + && this->p_chgmix->mixing_restart_count > 0) { + elecstate::DensityMatrix* dm + = dynamic_cast*>(this->pelec) + ->get_DM(); dm->save_DMR(); } // 3) solve the Hamiltonian and output band gap - if (this->phsol != nullptr) - { + if (this->phsol != nullptr) { // reset energy this->pelec->f_en.eband = 0.0; this->pelec->f_en.demet = 0.0; - this->phsol->solve(this->p_hamilt, this->psi[0], this->pelec, GlobalV::KS_SOLVER); + this->phsol->solve(this->p_hamilt, + this->psi[0], + this->pelec, + GlobalV::KS_SOLVER); - if (GlobalV::out_bandgap) - { - if (!GlobalV::TWO_EFERMI) - { + if (GlobalV::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!"); } // 4) print bands for each k-point and each band - for (int ik = 0; ik < this->kv.get_nks(); ++ik) - { + for (int ik = 0; ik < this->kv.get_nks(); ++ik) { this->pelec->print_band(ik, INPUT.printe, iter); } // 5) what's the exd used for? #ifdef __EXX - if (GlobalC::exx_info.info_ri.real_number) - { + if (GlobalC::exx_info.info_ri.real_number) { this->exd->exx_hamilt2density(*this->pelec, this->orb_con.ParaV, iter); - } - else - { + } else { this->exc->exx_hamilt2density(*this->pelec, this->orb_con.ParaV, iter); } #endif - // 6) calculate the local occupation number matrix and energy correction in DFT+U - if (GlobalV::dft_plus_u) - { + // 6) calculate the local occupation number matrix and energy correction in + // DFT+U + if (GlobalV::dft_plus_u) { // only old DFT+U method should calculated energy correction in esolver, // new DFT+U method will calculate energy in calculating Hamiltonian - if (GlobalV::dft_plus_u == 2) - { - if (GlobalC::dftu.omc != 2) - { + if (GlobalV::dft_plus_u == 2) { + if (GlobalC::dftu.omc != 2) { const std::vector>& tmp_dm - = dynamic_cast*>(this->pelec)->get_DM()->get_DMK_vector(); + = dynamic_cast*>(this->pelec) + ->get_DM() + ->get_DMK_vector(); this->dftu_cal_occup_m(iter, tmp_dm); } GlobalC::dftu.cal_energy_correction(istep); @@ -819,19 +836,20 @@ void ESolver_KS_LCAO::hamilt2density(int istep, int iter, double ethr) // (7) for deepks, calculate delta_e #ifdef __DEEPKS - if (GlobalV::deepks_scf) - { + if (GlobalV::deepks_scf) { const std::vector>& dm - = dynamic_cast*>(this->pelec)->get_DM()->get_DMK_vector(); + = dynamic_cast*>(this->pelec) + ->get_DM() + ->get_DMK_vector(); this->dpks_cal_e_delta_band(dm); } #endif // 8) for delta spin - if (GlobalV::sc_mag_switch) - { - SpinConstrain& sc = SpinConstrain::getScInstance(); + if (GlobalV::sc_mag_switch) { + SpinConstrain& sc + = SpinConstrain::getScInstance(); sc.cal_MW(iter, &(this->LM)); } @@ -840,9 +858,12 @@ void ESolver_KS_LCAO::hamilt2density(int istep, int iter, double ethr) // 10) symmetrize the charge density Symmetry_rho srho; - for (int is = 0; is < GlobalV::NSPIN; is++) - { - srho.begin(is, *(this->pelec->charge), this->pw_rho, GlobalC::Pgrid, GlobalC::ucell.symm); + for (int is = 0; is < GlobalV::NSPIN; is++) { + srho.begin(is, + *(this->pelec->charge), + this->pw_rho, + GlobalC::Pgrid, + GlobalC::ucell.symm); } // 11) compute magnetization, only for spin==2 @@ -863,31 +884,26 @@ void ESolver_KS_LCAO::hamilt2density(int istep, int iter, double ethr) //! 3) print potential //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::update_pot(const int istep, const int iter) -{ +void ESolver_KS_LCAO::update_pot(const int istep, const int iter) { ModuleBase::TITLE("ESolver_KS_LCAO", "update_pot"); // 1) print Hamiltonian and Overlap matrix - if (this->conv_elec || iter == GlobalV::SCF_NMAX) - { - if (!GlobalV::GAMMA_ONLY_LOCAL && hsolver::HSolverLCAO::out_mat_hs[0]) - { + if (this->conv_elec || iter == GlobalV::SCF_NMAX) { + if (!GlobalV::GAMMA_ONLY_LOCAL + && hsolver::HSolverLCAO::out_mat_hs[0]) { this->GK.renew(true); } - for (int ik = 0; ik < this->kv.get_nks(); ++ik) - { - if (hsolver::HSolverLCAO::out_mat_hs[0]) - { + for (int ik = 0; ik < this->kv.get_nks(); ++ik) { + if (hsolver::HSolverLCAO::out_mat_hs[0]) { this->p_hamilt->updateHk(ik); } bool bit = false; // LiuXh, 2017-03-21 - // if set bit = true, there would be error in soc-multi-core calculation, noted by zhengdy-soc - if (this->psi != nullptr && (istep % GlobalV::out_interval == 0)) - { + // if set bit = true, there would be error in soc-multi-core + // calculation, noted by zhengdy-soc + if (this->psi != nullptr && (istep % GlobalV::out_interval == 0)) { hamilt::MatrixBlock h_mat, s_mat; this->p_hamilt->matrix(h_mat, s_mat); - if (hsolver::HSolverLCAO::out_mat_hs[0]) - { + if (hsolver::HSolverLCAO::out_mat_hs[0]) { ModuleIO::save_mat(istep, h_mat.p, GlobalV::NLOCAL, @@ -916,9 +932,9 @@ void ESolver_KS_LCAO::update_pot(const int istep, const int iter) } // 2) print wavefunctions - if (elecstate::ElecStateLCAO::out_wfc_lcao && (this->conv_elec || iter == GlobalV::SCF_NMAX) - && (istep % GlobalV::out_interval == 0)) - { + if (elecstate::ElecStateLCAO::out_wfc_lcao + && (this->conv_elec || iter == GlobalV::SCF_NMAX) + && (istep % GlobalV::out_interval == 0)) { ModuleIO::write_wfc_nao(elecstate::ElecStateLCAO::out_wfc_lcao, this->psi[0], this->pelec->ekb, @@ -929,25 +945,21 @@ void ESolver_KS_LCAO::update_pot(const int istep, const int iter) } // 3) print potential - if (this->conv_elec || iter == GlobalV::SCF_NMAX) - { + if (this->conv_elec || iter == GlobalV::SCF_NMAX) { if (GlobalV::out_pot < 0) // mohan add 2011-10-10 { GlobalV::out_pot = -2; } } - if (!this->conv_elec) - { - if (GlobalV::NSPIN == 4) - { + 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(); } } @@ -963,47 +975,58 @@ void ESolver_KS_LCAO::update_pot(const int istep, const int iter) //! 6) calculate the total energy? //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::iter_finish(int iter) -{ +void ESolver_KS_LCAO::iter_finish(int iter) { ModuleBase::TITLE("ESolver_KS_LCAO", "iter_finish"); - // 1) mix density matrix if mixing_restart + mixing_dmr + not first mixing_restart at every iter - if (GlobalV::MIXING_RESTART > 0 && this->p_chgmix->mixing_restart_count > 0 && GlobalV::MIXING_DMR) - { - elecstate::DensityMatrix* dm = dynamic_cast*>(this->pelec)->get_DM(); + // 1) mix density matrix if mixing_restart + mixing_dmr + not first + // mixing_restart at every iter + if (GlobalV::MIXING_RESTART > 0 && this->p_chgmix->mixing_restart_count > 0 + && GlobalV::MIXING_DMR) { + elecstate::DensityMatrix* dm + = dynamic_cast*>(this->pelec) + ->get_DM(); this->p_chgmix->mix_dmr(dm); } // 2) save charge density // Peize Lin add 2020.04.04 - if (GlobalC::restart.info_save.save_charge) - { - for (int is = 0; is < GlobalV::NSPIN; ++is) - { - GlobalC::restart.save_disk("charge", is, this->pelec->charge->nrxx, this->pelec->charge->rho[is]); + if (GlobalC::restart.info_save.save_charge) { + for (int is = 0; is < GlobalV::NSPIN; ++is) { + GlobalC::restart.save_disk("charge", + is, + this->pelec->charge->nrxx, + this->pelec->charge->rho[is]); } } #ifdef __EXX // 3) save exx matrix - int two_level_step = GlobalC::exx_info.info_ri.real_number ? this->exd->two_level_step : this->exc->two_level_step; + int two_level_step = GlobalC::exx_info.info_ri.real_number + ? this->exd->two_level_step + : this->exc->two_level_step; if (GlobalC::restart.info_save.save_H && two_level_step > 0 - && (!GlobalC::exx_info.info_global.separate_loop || iter == 1)) // to avoid saving the same value repeatedly + && (!GlobalC::exx_info.info_global.separate_loop + || iter == 1)) // to avoid saving the same value repeatedly { std::vector Hexxk_save(this->orb_con.ParaV.get_local_size()); - for (int ik = 0; ik < this->kv.get_nks(); ++ik) - { + for (int ik = 0; ik < this->kv.get_nks(); ++ik) { ModuleBase::GlobalFunc::ZEROS(Hexxk_save.data(), Hexxk_save.size()); - hamilt::OperatorEXX> opexx_save(&this->LM, nullptr, &Hexxk_save, this->kv); + hamilt::OperatorEXX> opexx_save( + &this->LM, + nullptr, + &Hexxk_save, + this->kv); opexx_save.contributeHk(ik); - GlobalC::restart.save_disk("Hexx", ik, this->orb_con.ParaV.get_local_size(), Hexxk_save.data()); + GlobalC::restart.save_disk("Hexx", + ik, + this->orb_con.ParaV.get_local_size(), + Hexxk_save.data()); } - if (GlobalV::MY_RANK == 0) - { + if (GlobalV::MY_RANK == 0) { GlobalC::restart.save_disk("Eexx", 0, 1, &this->pelec->f_en.exx); } } @@ -1011,29 +1034,26 @@ void ESolver_KS_LCAO::iter_finish(int iter) // 4) output charge density and density matrix 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) - { - for (int is = 0; is < GlobalV::NSPIN; is++) - { + if (print) { + for (int is = 0; is < GlobalV::NSPIN; is++) { this->create_Output_Rho(is, iter, "tmp_").write(); - this->create_Output_DM(is, iter).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(); } } } // 5) cal_MW? - // escon: energy of spin constraint depends on Mi, so cal_energies should be called after cal_MW - if (GlobalV::sc_mag_switch) - { - SpinConstrain& sc = SpinConstrain::getScInstance(); + // escon: energy of spin constraint depends on Mi, so cal_energies should be + // called after cal_MW + if (GlobalV::sc_mag_switch) { + SpinConstrain& sc + = SpinConstrain::getScInstance(); sc.cal_MW(iter, &(this->LM)); } @@ -1064,13 +1084,11 @@ void ESolver_KS_LCAO::iter_finish(int iter) //! 18) write quasi-orbitals //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::after_scf(const int istep) -{ +void ESolver_KS_LCAO::after_scf(const int istep) { ModuleBase::TITLE("ESolver_KS_LCAO", "after_scf"); // 1) write charge difference into files for charge extrapolation - if (GlobalV::CALCULATION != "scf") - { + if (GlobalV::CALCULATION != "scf") { this->CE.save_files(istep, GlobalC::ucell, #ifdef __MPI @@ -1081,46 +1099,53 @@ void ESolver_KS_LCAO::after_scf(const int istep) } // 2) write density matrix for sparse matrix - ModuleIO::write_dmr(dynamic_cast*>(this->pelec)->get_DM()->get_DMR_vector(), - this->orb_con.ParaV, - INPUT.out_dm1, - false, - GlobalV::out_app_flag, - istep); + ModuleIO::write_dmr( + dynamic_cast*>(this->pelec) + ->get_DM() + ->get_DMR_vector(), + this->orb_con.ParaV, + INPUT.out_dm1, + false, + GlobalV::out_app_flag, + istep); // 3) write charge density - if (GlobalV::out_chg) - { - for (int is = 0; is < GlobalV::NSPIN; is++) - { + if (GlobalV::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(); } } } // 4) write density matrix - if (this->LOC.out_dm) - { - for (int is = 0; is < GlobalV::NSPIN; is++) - { - this->create_Output_DM(is, istep).write(); + if (INPUT.out_dm) { + std::vector efermis(GlobalV::NSPIN == 2 ? 2 : 1); + for (int ispin = 0; ispin < efermis.size(); ispin++) { + efermis[ispin] = this->pelec->eferm.get_efval(ispin); } + const int precision = 3; + ModuleIO::write_dmk( + dynamic_cast*>(this->pelec) + ->get_DM() + ->get_DMK_vector(), + precision, + efermis, + &(GlobalC::ucell), + this->orb_con.ParaV); } #ifdef __EXX // 5) write Exx matrix if (GlobalC::exx_info.info_global.cal_exx) // Peize Lin add if 2022.11.14 { - const std::string file_name_exx = GlobalV::global_out_dir + "HexxR" + std::to_string(GlobalV::MY_RANK); - if (GlobalC::exx_info.info_ri.real_number) - { + const std::string file_name_exx = GlobalV::global_out_dir + "HexxR" + + std::to_string(GlobalV::MY_RANK); + if (GlobalC::exx_info.info_ri.real_number) { this->exd->write_Hexxs_csr(file_name_exx, GlobalC::ucell); - } - else - { + } else { this->exc->write_Hexxs_csr(file_name_exx, GlobalC::ucell); } } @@ -1130,116 +1155,124 @@ void ESolver_KS_LCAO::after_scf(const int istep) this->create_Output_Potential(istep).write(); // 7) write convergence - 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); // 8) write fermi energy ModuleIO::output_efermi(this->conv_elec, this->pelec->eferm.ef); // 9) write eigenvalues - if (GlobalV::OUT_LEVEL != "m") - { + if (GlobalV::OUT_LEVEL != "m") { this->pelec->print_eigenvalue(GlobalV::ofs_running); } // 10) write deepks information #ifdef __DEEPKS - std::shared_ptr ld_shared_ptr(&GlobalC::ld, [](LCAO_Deepks*) {}); + std::shared_ptr ld_shared_ptr(&GlobalC::ld, + [](LCAO_Deepks*) {}); LCAO_Deepks_Interface LDI = LCAO_Deepks_Interface(ld_shared_ptr); ModuleBase::timer::tick("ESolver_KS_LCAO", "out_deepks_labels"); - LDI.out_deepks_labels(this->pelec->f_en.etot, - this->pelec->klist->get_nks(), - GlobalC::ucell.nat, - this->pelec->ekb, - this->pelec->klist->kvec_d, - GlobalC::ucell, - GlobalC::ORB, - GlobalC::GridD, - &(this->orb_con.ParaV), - *(this->psi), - dynamic_cast*>(this->pelec)->get_DM()); + LDI.out_deepks_labels( + this->pelec->f_en.etot, + this->pelec->klist->get_nks(), + GlobalC::ucell.nat, + this->pelec->ekb, + this->pelec->klist->kvec_d, + GlobalC::ucell, + GlobalC::ORB, + GlobalC::GridD, + &(this->orb_con.ParaV), + *(this->psi), + dynamic_cast*>(this->pelec) + ->get_DM()); ModuleBase::timer::tick("ESolver_KS_LCAO", "out_deepks_labels"); #endif #ifdef __EXX // 11) write rpa information - if (INPUT.rpa) - { - // ModuleRPA::DFT_RPA_interface rpa_interface(GlobalC::exx_info.info_global); + if (INPUT.rpa) { + // ModuleRPA::DFT_RPA_interface + // rpa_interface(GlobalC::exx_info.info_global); // rpa_interface.rpa_exx_lcao().info.files_abfs = GlobalV::rpa_orbitals; RPA_LRI rpa_lri_double(GlobalC::exx_info.info_ri); - rpa_lri_double.cal_postSCF_exx(*dynamic_cast*>(this->pelec)->get_DM(), - MPI_COMM_WORLD, - this->kv); + rpa_lri_double.cal_postSCF_exx( + *dynamic_cast*>(this->pelec) + ->get_DM(), + MPI_COMM_WORLD, + this->kv); rpa_lri_double.init(MPI_COMM_WORLD, this->kv); - rpa_lri_double.out_for_RPA(this->orb_con.ParaV, *(this->psi), this->pelec); + rpa_lri_double.out_for_RPA(this->orb_con.ParaV, + *(this->psi), + this->pelec); } #endif // 12) write HR in npz format - if (GlobalV::out_hr_npz) - { + if (GlobalV::out_hr_npz) { this->p_hamilt->updateHk(0); // first k point, up spin hamilt::HamiltLCAO, double>* p_ham_lcao - = dynamic_cast, double>*>(this->p_hamilt); + = dynamic_cast, double>*>( + this->p_hamilt); std::string zipname = "output_HR0.npz"; this->output_mat_npz(zipname, *(p_ham_lcao->getHR())); - if (GlobalV::NSPIN == 2) - { - this->p_hamilt->updateHk(this->kv.get_nks() / 2); // the other half of k points, down spin + if (GlobalV::NSPIN == 2) { + this->p_hamilt->updateHk( + this->kv.get_nks() + / 2); // the other half of k points, down spin hamilt::HamiltLCAO, double>* p_ham_lcao - = dynamic_cast, double>*>(this->p_hamilt); + = dynamic_cast< + hamilt::HamiltLCAO, double>*>( + this->p_hamilt); zipname = "output_HR1.npz"; this->output_mat_npz(zipname, *(p_ham_lcao->getHR())); } } // 13) write dm in npz format - if (GlobalV::out_dm_npz) - { + if (GlobalV::out_dm_npz) { const elecstate::DensityMatrix* dm - = dynamic_cast*>(this->pelec)->get_DM(); + = dynamic_cast*>(this->pelec) + ->get_DM(); std::string zipname = "output_DM0.npz"; this->output_mat_npz(zipname, *(dm->get_DMR_pointer(1))); - if (GlobalV::NSPIN == 2) - { + if (GlobalV::NSPIN == 2) { zipname = "output_DM1.npz"; this->output_mat_npz(zipname, *(dm->get_DMR_pointer(2))); } } // 14) write md related - if (!md_skip_out(GlobalV::CALCULATION, istep, GlobalV::out_interval)) - { + if (!md_skip_out(GlobalV::CALCULATION, istep, GlobalV::out_interval)) { this->create_Output_Mat_Sparse(istep).write(); // mulliken charge analysis - if (GlobalV::out_mul) - { + if (GlobalV::out_mul) { this->cal_mag(istep, true); } } // 15) write spin constrian MW? // spin constrain calculations, added by Tianqi Zhao. - if (GlobalV::sc_mag_switch) - { - SpinConstrain& sc = SpinConstrain::getScInstance(); + if (GlobalV::sc_mag_switch) { + SpinConstrain& sc + = SpinConstrain::getScInstance(); sc.cal_MW(istep, &(this->LM), true); sc.print_Mag_Force(); } // 16) delete grid - if (!GlobalV::CAL_FORCE && !GlobalV::CAL_STRESS) - { + if (!GlobalV::CAL_FORCE && !GlobalV::CAL_STRESS) { RA.delete_grid(); } // 17) write quasi-orbitals, added by Yike Huang. - if (GlobalV::qo_switch) - { - toQO tqo(GlobalV::qo_basis, GlobalV::qo_strategy, GlobalV::qo_thr, GlobalV::qo_screening_coeff); + if (GlobalV::qo_switch) { + toQO tqo(GlobalV::qo_basis, + GlobalV::qo_strategy, + GlobalV::qo_thr, + GlobalV::qo_screening_coeff); tqo.initialize(GlobalV::global_out_dir, GlobalV::global_pseudo_dir, GlobalV::global_orbital_dir, @@ -1257,31 +1290,30 @@ void ESolver_KS_LCAO::after_scf(const int istep) //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -bool ESolver_KS_LCAO::do_after_converge(int& iter) -{ +bool ESolver_KS_LCAO::do_after_converge(int& iter) { ModuleBase::TITLE("ESolver_KS_LCAO", "do_after_converge"); #ifdef __EXX - if (GlobalC::exx_info.info_ri.real_number) - { - return this->exd->exx_after_converge(*this->p_hamilt, - this->LM, - *dynamic_cast*>(this->pelec)->get_DM(), - this->kv, - iter); - } - else - { - return this->exc->exx_after_converge(*this->p_hamilt, - this->LM, - *dynamic_cast*>(this->pelec)->get_DM(), - this->kv, - iter); + if (GlobalC::exx_info.info_ri.real_number) { + return this->exd->exx_after_converge( + *this->p_hamilt, + this->LM, + *dynamic_cast*>(this->pelec) + ->get_DM(), + this->kv, + iter); + } else { + return this->exc->exx_after_converge( + *this->p_hamilt, + this->LM, + *dynamic_cast*>(this->pelec) + ->get_DM(), + this->kv, + iter); } #endif // __EXX - if (GlobalV::dft_plus_u) - { + if (GlobalV::dft_plus_u) { // use the converged occupation matrix for next MD/Relax SCF calculation GlobalC::dftu.initialed_locale = true; } @@ -1293,22 +1325,6 @@ bool ESolver_KS_LCAO::do_after_converge(int& iter) //! the 16th function of ESolver_KS_LCAO: create_Output_DM //! mohan add 2024-05-11 //------------------------------------------------------------------------------ -template -ModuleIO::Output_DM ESolver_KS_LCAO::create_Output_DM(int is, int iter) -{ - const int precision = 3; - - return ModuleIO::Output_DM(this->GridT, - is, - iter, - precision, - this->LOC.out_dm, - this->LOC.DM, - this->pelec->eferm.get_efval(is), - &(GlobalC::ucell), - GlobalV::global_out_dir, - GlobalV::GAMMA_ONLY_LOCAL); -} //------------------------------------------------------------------------------ //! the 17th function of ESolver_KS_LCAO: create_Output_DM1 @@ -1320,21 +1336,22 @@ ModuleIO::Output_DM ESolver_KS_LCAO::create_Output_DM(int is, int iter) //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -ModuleIO::Output_Mat_Sparse ESolver_KS_LCAO::create_Output_Mat_Sparse(int istep) -{ - return ModuleIO::Output_Mat_Sparse(hsolver::HSolverLCAO::out_mat_hsR, - hsolver::HSolverLCAO::out_mat_dh, - hsolver::HSolverLCAO::out_mat_t, - INPUT.out_mat_r, - istep, - this->pelec->pot->get_effective_v(), - this->orb_con.ParaV, - this->GK, // mohan add 2024-04-01 - two_center_bundle_, - this->LM, - GlobalC::GridD, // mohan add 2024-04-06 - this->kv, - this->p_hamilt); +ModuleIO::Output_Mat_Sparse + ESolver_KS_LCAO::create_Output_Mat_Sparse(int istep) { + return ModuleIO::Output_Mat_Sparse( + hsolver::HSolverLCAO::out_mat_hsR, + hsolver::HSolverLCAO::out_mat_dh, + hsolver::HSolverLCAO::out_mat_t, + INPUT.out_mat_r, + istep, + this->pelec->pot->get_effective_v(), + this->orb_con.ParaV, + this->GK, // mohan add 2024-04-01 + two_center_bundle_, + this->LM, + GlobalC::GridD, // mohan add 2024-04-06 + this->kv, + this->p_hamilt); } //------------------------------------------------------------------------------ @@ -1342,12 +1359,11 @@ ModuleIO::Output_Mat_Sparse ESolver_KS_LCAO::create_Output_Mat_Spars //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -bool ESolver_KS_LCAO::md_skip_out(std::string calculation, int istep, int interval) -{ - if (calculation == "md") - { - if (istep % interval != 0) - { +bool ESolver_KS_LCAO::md_skip_out(std::string calculation, + int istep, + int interval) { + if (calculation == "md") { + if (istep % interval != 0) { return true; } } @@ -1355,8 +1371,7 @@ bool ESolver_KS_LCAO::md_skip_out(std::string calculation, int istep, in } template -void ESolver_KS_LCAO::cal_mag(const int istep, const bool print) -{ +void ESolver_KS_LCAO::cal_mag(const int istep, const bool print) { auto cell_index = CellIndex(GlobalC::ucell.get_atomLabels(), GlobalC::ucell.get_atomCounts(), GlobalC::ucell.get_lnchiCounts(), @@ -1366,10 +1381,12 @@ void ESolver_KS_LCAO::cal_mag(const int istep, const bool print) &(this->orb_con.ParaV), GlobalV::NSPIN, this->kv.get_nks()); - auto out_dmk = ModuleIO::Output_DMK(dynamic_cast*>(this->pelec)->get_DM(), - &(this->orb_con.ParaV), - GlobalV::NSPIN, - this->kv.get_nks()); + auto out_dmk = ModuleIO::Output_DMK( + dynamic_cast*>(this->pelec) + ->get_DM(), + &(this->orb_con.ParaV), + GlobalV::NSPIN, + this->kv.get_nks()); auto mulp = ModuleIO::Output_Mulliken(&(out_sk), &(out_dmk), &(this->orb_con.ParaV), @@ -1379,8 +1396,7 @@ void ESolver_KS_LCAO::cal_mag(const int istep, const bool print) auto atom_chg = mulp.get_atom_chg(); /// used in updating mag info in STRU file GlobalC::ucell.atom_mulliken = mulp.get_atom_mulliken(atom_chg); - if (print && GlobalV::MY_RANK == 0) - { + if (print && GlobalV::MY_RANK == 0) { /// write the Orbital file cell_index.write_orb_info(GlobalV::global_out_dir); /// write mulliken.txt diff --git a/source/module_esolver/esolver_ks_lcao.h b/source/module_esolver/esolver_ks_lcao.h index 2ecdb67de9..d7aa653f3f 100644 --- a/source/module_esolver/esolver_ks_lcao.h +++ b/source/module_esolver/esolver_ks_lcao.h @@ -12,16 +12,13 @@ #include "module_ri/Mix_DMk_2D.h" #endif #include "module_basis/module_nao/two_center_bundle.h" -#include "module_io/output_dm.h" #include "module_io/output_mat_sparse.h" #include -namespace ModuleESolver -{ +namespace ModuleESolver { template -class ESolver_KS_LCAO : public ESolver_KS -{ +class ESolver_KS_LCAO : public ESolver_KS { public: ESolver_KS_LCAO(); ~ESolver_KS_LCAO(); @@ -49,7 +46,9 @@ class ESolver_KS_LCAO : public ESolver_KS virtual void iter_init(const int istep, const int iter) override; - virtual void hamilt2density(const int istep, const int iter, const double ethr) override; + virtual void hamilt2density(const int istep, + const int iter, + const double ethr) override; virtual void update_pot(const int istep, const int iter) override; @@ -99,14 +98,13 @@ class ESolver_KS_LCAO : public ESolver_KS void beforesolver(const int istep); //---------------------------------------------------------------------- - /// @brief create ModuleIO::Output_DM object to output density matrix - ModuleIO::Output_DM create_Output_DM(int is, int iter); - - /// @brief create ModuleIO::Output_Mat_Sparse object to output sparse density matrix of H, S, T, r + /// @brief create ModuleIO::Output_Mat_Sparse object to output sparse + /// density matrix of H, S, T, r ModuleIO::Output_Mat_Sparse create_Output_Mat_Sparse(int istep); void read_mat_npz(std::string& zipname, hamilt::HContainer& hR); - void output_mat_npz(std::string& zipname, const hamilt::HContainer& hR); + void output_mat_npz(std::string& zipname, + const hamilt::HContainer& hR); /// @brief check if skip the corresponding output in md calculation bool md_skip_out(std::string calculation, int istep, int interval); @@ -120,12 +118,14 @@ class ESolver_KS_LCAO : public ESolver_KS private: // tmp interfaces before sub-modules are refactored - void dftu_cal_occup_m(const int& iter, const std::vector>& dm) const; + void dftu_cal_occup_m(const int& iter, + const std::vector>& dm) const; #ifdef __DEEPKS void dpks_cal_e_delta_band(const std::vector>& dm) const; - void dpks_cal_projected_DM(const elecstate::DensityMatrix* dm) const; + void dpks_cal_projected_DM( + const elecstate::DensityMatrix* dm) const; #endif }; } // namespace ModuleESolver diff --git a/source/module_esolver/esolver_ks_lcao_elec.cpp b/source/module_esolver/esolver_ks_lcao_elec.cpp index 41426c335c..ff54dbfd7d 100644 --- a/source/module_esolver/esolver_ks_lcao_elec.cpp +++ b/source/module_esolver/esolver_ks_lcao_elec.cpp @@ -24,25 +24,23 @@ #include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.h" #include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h" #include "module_hamilt_lcao/module_deltaspin/spin_constrain.h" -#include "module_io/dm_io.h" #include "module_io/rho_io.h" #include "module_io/write_pot.h" -namespace ModuleESolver -{ +namespace ModuleESolver { template -void ESolver_KS_LCAO::set_matrix_grid(Record_adj& ra) -{ +void ESolver_KS_LCAO::set_matrix_grid(Record_adj& ra) { ModuleBase::TITLE("ESolver_KS_LCAO", "set_matrix_grid"); ModuleBase::timer::tick("ESolver_KS_LCAO", "set_matrix_grid"); // (1) Find adjacent atoms for each atom. - GlobalV::SEARCH_RADIUS = atom_arrange::set_sr_NL(GlobalV::ofs_running, - GlobalV::OUT_LEVEL, - GlobalC::ORB.get_rcutmax_Phi(), - GlobalC::ucell.infoNL.get_rcutmax_Beta(), - GlobalV::GAMMA_ONLY_LOCAL); + GlobalV::SEARCH_RADIUS + = atom_arrange::set_sr_NL(GlobalV::ofs_running, + GlobalV::OUT_LEVEL, + GlobalC::ORB.get_rcutmax_Phi(), + GlobalC::ucell.infoNL.get_rcutmax_Beta(), + GlobalV::GAMMA_ONLY_LOCAL); atom_arrange::search(GlobalV::SEARCH_PBC, GlobalV::ofs_running, @@ -51,7 +49,8 @@ void ESolver_KS_LCAO::set_matrix_grid(Record_adj& ra) GlobalV::SEARCH_RADIUS, GlobalV::test_atom_input); - // ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running,"SEARCH ADJACENT ATOMS"); + // ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running,"SEARCH ADJACENT + // ATOMS"); // (3) Periodic condition search for each grid. double dr_uniform = 0.001; @@ -60,7 +59,12 @@ void ESolver_KS_LCAO::set_matrix_grid(Record_adj& ra) std::vector> dpsi_u; std::vector> d2psi_u; - Gint_Tools::init_orb(dr_uniform, rcuts, GlobalC::ucell, psi_u, dpsi_u, d2psi_u); + Gint_Tools::init_orb(dr_uniform, + rcuts, + GlobalC::ucell, + psi_u, + dpsi_u, + d2psi_u); this->GridT.set_pbc_grid(this->pw_rho->nx, this->pw_rho->ny, @@ -95,8 +99,7 @@ void ESolver_KS_LCAO::set_matrix_grid(Record_adj& ra) // If k point is used here, allocate HlocR after atom_arrange. Parallel_Orbitals* pv = this->LM.ParaV; ra.for_2d(*pv, GlobalV::GAMMA_ONLY_LOCAL); - if (!GlobalV::GAMMA_ONLY_LOCAL) - { + if (!GlobalV::GAMMA_ONLY_LOCAL) { // need to first calculae lgd. // using GridT.init. this->GridT.cal_nnrg(pv); @@ -107,8 +110,7 @@ void ESolver_KS_LCAO::set_matrix_grid(Record_adj& ra) } template -void ESolver_KS_LCAO::beforesolver(const int istep) -{ +void ESolver_KS_LCAO::beforesolver(const int istep) { ModuleBase::TITLE("ESolver_KS_LCAO", "beforesolver"); ModuleBase::timer::tick("ESolver_KS_LCAO", "beforesolver"); @@ -123,22 +125,20 @@ void ESolver_KS_LCAO::beforesolver(const int istep) // the force. // init psi - if (this->psi == nullptr) - { + if (this->psi == nullptr) { int nsk = 0; int ncol = 0; - if (GlobalV::GAMMA_ONLY_LOCAL) - { + if (GlobalV::GAMMA_ONLY_LOCAL) { nsk = GlobalV::NSPIN; ncol = this->orb_con.ParaV.ncol_bands; - if (GlobalV::KS_SOLVER == "genelpa" || GlobalV::KS_SOLVER == "lapack_gvx" || GlobalV::KS_SOLVER == "pexsi" - || GlobalV::KS_SOLVER == "cusolver" || GlobalV::KS_SOLVER == "cusolvermp") - { + if (GlobalV::KS_SOLVER == "genelpa" + || GlobalV::KS_SOLVER == "lapack_gvx" + || GlobalV::KS_SOLVER == "pexsi" + || GlobalV::KS_SOLVER == "cusolver" + || GlobalV::KS_SOLVER == "cusolvermp") { ncol = this->orb_con.ParaV.ncol; } - } - else - { + } else { nsk = this->kv.get_nks(); #ifdef __MPI ncol = this->orb_con.ParaV.ncol_bands; @@ -146,21 +146,26 @@ void ESolver_KS_LCAO::beforesolver(const int istep) ncol = GlobalV::NBANDS; #endif } - this->psi = new psi::Psi(nsk, ncol, this->orb_con.ParaV.nrow, nullptr); + this->psi + = new psi::Psi(nsk, ncol, this->orb_con.ParaV.nrow, nullptr); } // prepare grid in Gint - LCAO_domain::grid_prepare(this->GridT, this->GG, this->GK, *this->pw_rho, *this->pw_big); + LCAO_domain::grid_prepare(this->GridT, + this->GG, + this->GK, + *this->pw_rho, + *this->pw_big); // init Hamiltonian - if (this->p_hamilt != nullptr) - { + if (this->p_hamilt != nullptr) { delete this->p_hamilt; this->p_hamilt = nullptr; } - if (this->p_hamilt == nullptr) - { - elecstate::DensityMatrix* DM = dynamic_cast*>(this->pelec)->get_DM(); + if (this->p_hamilt == nullptr) { + elecstate::DensityMatrix* DM + = dynamic_cast*>(this->pelec) + ->get_DM(); this->p_hamilt = new hamilt::HamiltLCAO( GlobalV::GAMMA_ONLY_LOCAL ? &(this->GG) : nullptr, GlobalV::GAMMA_ONLY_LOCAL ? nullptr : &(this->GK), @@ -171,19 +176,23 @@ void ESolver_KS_LCAO::beforesolver(const int istep) two_center_bundle_, #ifdef __EXX DM, - GlobalC::exx_info.info_ri.real_number ? &this->exd->two_level_step : &this->exc->two_level_step); + GlobalC::exx_info.info_ri.real_number ? &this->exd->two_level_step + : &this->exc->two_level_step); #else DM); #endif } // init density kernel and wave functions. - this->LOC.allocate_dm_wfc(this->GridT, this->pelec, this->psi, this->kv, istep); + this->LOC.allocate_dm_wfc(this->GridT, + this->pelec, + this->psi, + this->kv, + istep); #ifdef __DEEPKS // for each ionic step, the overlap must be rebuilt // since it depends on ionic positions - if (GlobalV::deepks_setorb) - { + if (GlobalV::deepks_setorb) { const Parallel_Orbitals* pv = this->LM.ParaV; // build and save at beginning GlobalC::ld.build_psialpha(GlobalV::CAL_FORCE, @@ -192,15 +201,17 @@ void ESolver_KS_LCAO::beforesolver(const int istep) GlobalC::GridD, *(two_center_bundle_.overlap_orb_alpha)); - if (GlobalV::deepks_out_unittest) - { - GlobalC::ld.check_psialpha(GlobalV::CAL_FORCE, GlobalC::ucell, GlobalC::ORB, GlobalC::GridD); + if (GlobalV::deepks_out_unittest) { + GlobalC::ld.check_psialpha(GlobalV::CAL_FORCE, + GlobalC::ucell, + GlobalC::ORB, + GlobalC::GridD); } } #endif - if (GlobalV::sc_mag_switch) - { - SpinConstrain& sc = SpinConstrain::getScInstance(); + if (GlobalV::sc_mag_switch) { + SpinConstrain& sc + = SpinConstrain::getScInstance(); sc.init_sc(GlobalV::sc_thr, GlobalV::nsc, GlobalV::nsc_min, @@ -224,25 +235,21 @@ void ESolver_KS_LCAO::beforesolver(const int istep) // 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(); } ModuleBase::timer::tick("ESolver_KS_LCAO", "beforesolver"); } template -void ESolver_KS_LCAO::before_scf(const int istep) -{ +void ESolver_KS_LCAO::before_scf(const int istep) { ModuleBase::TITLE("ESolver_KS_LCAO", "before_scf"); ModuleBase::timer::tick("ESolver_KS_LCAO", "before_scf"); - if (GlobalC::ucell.cell_parameter_updated) - { + if (GlobalC::ucell.cell_parameter_updated) { this->init_after_vc(INPUT, 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 @@ -257,8 +264,7 @@ void ESolver_KS_LCAO::before_scf(const int istep) // about vdw, jiyy add vdwd3 and linpz add vdwd2 //---------------------------------------------------------- auto vdw_solver = vdw::make_vdw(GlobalC::ucell, INPUT); - if (vdw_solver != nullptr) - { + if (vdw_solver != nullptr) { this->pelec->f_en.evdw = vdw_solver->get_energy(); } @@ -266,23 +272,19 @@ void ESolver_KS_LCAO::before_scf(const int istep) // Peize Lin add 2016-12-03 #ifdef __EXX // set xc type before the first cal of xc in pelec->init_scf - if (GlobalC::exx_info.info_ri.real_number) - { + if (GlobalC::exx_info.info_ri.real_number) { this->exd->exx_beforescf(this->kv, *this->p_chgmix); - } - else - { + } else { this->exc->exx_beforescf(this->kv, *this->p_chgmix); } #endif // __EXX this->pelec->init_scf(istep, this->sf.strucFac); - if (GlobalV::out_chg == 2) - { - for (int is = 0; is < GlobalV::NSPIN; is++) - { + if (GlobalV::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, // bz first, then nbz @@ -322,16 +324,16 @@ void ESolver_KS_LCAO::before_scf(const int istep) // DMR should be same size with Hamiltonian(R) dynamic_cast*>(this->pelec) ->get_DM() - ->init_DMR(*(dynamic_cast*>(this->p_hamilt)->getHR())); + ->init_DMR(*(dynamic_cast*>(this->p_hamilt) + ->getHR())); - if (GlobalV::dm_to_rho) - { + if (GlobalV::dm_to_rho) { std::string zipname = "output_DM0.npz"; elecstate::DensityMatrix* dm - = dynamic_cast*>(this->pelec)->get_DM(); + = dynamic_cast*>(this->pelec) + ->get_DM(); this->read_mat_npz(zipname, *(dm->get_DMR_pointer(1))); - if (GlobalV::NSPIN == 2) - { + if (GlobalV::NSPIN == 2) { zipname = "output_DM1.npz"; this->read_mat_npz(zipname, *(dm->get_DMR_pointer(2))); } @@ -339,8 +341,7 @@ void ESolver_KS_LCAO::before_scf(const int istep) this->pelec->psiToRho(*this->psi); this->create_Output_Rho(0, istep).write(); - if (GlobalV::NSPIN == 2) - { + if (GlobalV::NSPIN == 2) { this->create_Output_Rho(1, istep).write(); } @@ -350,16 +351,21 @@ void ESolver_KS_LCAO::before_scf(const int istep) // the electron charge density should be symmetrized, // here is the initialization Symmetry_rho srho; - for (int is = 0; is < GlobalV::NSPIN; is++) - { - srho.begin(is, *(this->pelec->charge), this->pw_rho, GlobalC::Pgrid, GlobalC::ucell.symm); + for (int is = 0; is < GlobalV::NSPIN; is++) { + srho.begin(is, + *(this->pelec->charge), + this->pw_rho, + GlobalC::Pgrid, + GlobalC::ucell.symm); } // 1. calculate ewald energy. // mohan update 2021-02-25 - if (!GlobalV::test_skip_ewald) - { - this->pelec->f_en.ewald_energy = H_Ewald_pw::compute_ewald(GlobalC::ucell, this->pw_rho, this->sf.strucFac); + if (!GlobalV::test_skip_ewald) { + this->pelec->f_en.ewald_energy + = H_Ewald_pw::compute_ewald(GlobalC::ucell, + this->pw_rho, + this->sf.strucFac); } this->p_hamilt->non_first_scf = istep; @@ -369,15 +375,13 @@ void ESolver_KS_LCAO::before_scf(const int istep) } template -void ESolver_KS_LCAO::others(const int istep) -{ +void ESolver_KS_LCAO::others(const int istep) { ModuleBase::TITLE("ESolver_KS_LCAO", "others"); ModuleBase::timer::tick("ESolver_KS_LCAO", "others"); const std::string cal_type = GlobalV::CALCULATION; - if (cal_type == "get_S") - { + if (cal_type == "get_S") { std::cout << "\n * * * * * *" << std::endl; std::cout << " << Start writing the overlap matrix." << std::endl; this->get_S(); @@ -386,22 +390,19 @@ void ESolver_KS_LCAO::others(const int istep) ModuleBase::QUIT(); - // return; // use 'return' will cause segmentation fault. by mohan 2024-06-09 - } - else if (cal_type == "test_memory") - { + // return; // use 'return' will cause segmentation fault. by mohan + // 2024-06-09 + } else 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()); return; - } - else if (cal_type == "test_neighbour") - { + } else if (cal_type == "test_neighbour") { // test_search_neighbor(); - if (GlobalV::SEARCH_RADIUS < 0) - { - std::cout << " SEARCH_RADIUS : " << GlobalV::SEARCH_RADIUS << std::endl; + if (GlobalV::SEARCH_RADIUS < 0) { + std::cout << " SEARCH_RADIUS : " << GlobalV::SEARCH_RADIUS + << std::endl; std::cout << " please make sure search_radius > 0" << std::endl; } @@ -419,12 +420,9 @@ void ESolver_KS_LCAO::others(const int istep) // pelec should be initialized before these calculations this->pelec->init_scf(istep, this->sf.strucFac); // self consistent calculations for electronic ground state - if (GlobalV::CALCULATION == "nscf") - { + if (GlobalV::CALCULATION == "nscf") { this->nscf(); - } - else if (cal_type == "get_pchg") - { + } else if (cal_type == "get_pchg") { IState_Charge ISC(this->psi, &(this->orb_con.ParaV)); ISC.begin(this->GG, this->pelec->charge->rho, @@ -450,12 +448,9 @@ void ESolver_KS_LCAO::others(const int istep) GlobalV::ofs_warning, &GlobalC::ucell, &GlobalC::GridD); - } - else if (cal_type == "get_wf") - { + } else if (cal_type == "get_wf") { IState_Envelope IEP(this->pelec); - if (GlobalV::GAMMA_ONLY_LOCAL) - { + if (GlobalV::GAMMA_ONLY_LOCAL) { IEP.begin(this->psi, this->pw_rho, this->pw_wfc, @@ -472,9 +467,7 @@ void ESolver_KS_LCAO::others(const int istep) GlobalV::NSPIN, GlobalV::NLOCAL, GlobalV::global_out_dir); - } - else - { + } else { IEP.begin(this->psi, this->pw_rho, this->pw_wfc, @@ -492,10 +485,9 @@ void ESolver_KS_LCAO::others(const int istep) GlobalV::NLOCAL, GlobalV::global_out_dir); } - } - else - { - ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::others", "CALCULATION type not supported"); + } else { + ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::others", + "CALCULATION type not supported"); } ModuleBase::timer::tick("ESolver_KS_LCAO", "others"); @@ -503,22 +495,22 @@ void ESolver_KS_LCAO::others(const int istep) } template <> -void ESolver_KS_LCAO::get_S(void) -{ +void ESolver_KS_LCAO::get_S(void) { ModuleBase::TITLE("ESolver_KS_LCAO", "get_S"); - ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::get_S", "not implemented for"); + ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::get_S", + "not implemented for"); } template <> -void ESolver_KS_LCAO, double>::get_S(void) -{ +void ESolver_KS_LCAO, double>::get_S(void) { ModuleBase::TITLE("ESolver_KS_LCAO", "get_S"); // (1) Find adjacent atoms for each atom. - GlobalV::SEARCH_RADIUS = atom_arrange::set_sr_NL(GlobalV::ofs_running, - GlobalV::OUT_LEVEL, - GlobalC::ORB.get_rcutmax_Phi(), - GlobalC::ucell.infoNL.get_rcutmax_Beta(), - GlobalV::GAMMA_ONLY_LOCAL); + GlobalV::SEARCH_RADIUS + = atom_arrange::set_sr_NL(GlobalV::ofs_running, + GlobalV::OUT_LEVEL, + GlobalC::ORB.get_rcutmax_Phi(), + GlobalC::ucell.infoNL.get_rcutmax_Beta(), + GlobalV::GAMMA_ONLY_LOCAL); atom_arrange::search(GlobalV::SEARCH_PBC, GlobalV::ofs_running, @@ -531,12 +523,14 @@ void ESolver_KS_LCAO, double>::get_S(void) this->LM.ParaV = &this->orb_con.ParaV; - if (this->p_hamilt == nullptr) - { - this->p_hamilt = new hamilt::HamiltLCAO, double>(&this->LM, - this->kv, - *(two_center_bundle_.overlap_orb)); - dynamic_cast, double>*>(this->p_hamilt->ops)->contributeHR(); + if (this->p_hamilt == nullptr) { + this->p_hamilt = new hamilt::HamiltLCAO, double>( + &this->LM, + this->kv, + *(two_center_bundle_.overlap_orb)); + dynamic_cast, double>*>( + this->p_hamilt->ops) + ->contributeHR(); } // mohan add 2024-06-09 @@ -544,21 +538,25 @@ void ESolver_KS_LCAO, double>::get_S(void) std::cout << " The file is saved in " << fn << std::endl; - ModuleIO::output_SR(orb_con.ParaV, this->LM, GlobalC::GridD, this->p_hamilt, fn); + ModuleIO::output_SR(orb_con.ParaV, + this->LM, + GlobalC::GridD, + this->p_hamilt, + fn); return; } template <> -void ESolver_KS_LCAO, std::complex>::get_S(void) -{ +void ESolver_KS_LCAO, std::complex>::get_S(void) { ModuleBase::TITLE("ESolver_KS_LCAO", "get_S"); // (1) Find adjacent atoms for each atom. - GlobalV::SEARCH_RADIUS = atom_arrange::set_sr_NL(GlobalV::ofs_running, - GlobalV::OUT_LEVEL, - GlobalC::ORB.get_rcutmax_Phi(), - GlobalC::ucell.infoNL.get_rcutmax_Beta(), - GlobalV::GAMMA_ONLY_LOCAL); + GlobalV::SEARCH_RADIUS + = atom_arrange::set_sr_NL(GlobalV::ofs_running, + GlobalV::OUT_LEVEL, + GlobalC::ORB.get_rcutmax_Phi(), + GlobalC::ucell.infoNL.get_rcutmax_Beta(), + GlobalV::GAMMA_ONLY_LOCAL); atom_arrange::search(GlobalV::SEARCH_PBC, GlobalV::ofs_running, @@ -569,13 +567,15 @@ void ESolver_KS_LCAO, std::complex>::get_S(void) this->RA.for_2d(this->orb_con.ParaV, GlobalV::GAMMA_ONLY_LOCAL); this->LM.ParaV = &this->orb_con.ParaV; - if (this->p_hamilt == nullptr) - { - this->p_hamilt - = new hamilt::HamiltLCAO, std::complex>(&this->LM, - this->kv, - *(two_center_bundle_.overlap_orb)); - dynamic_cast, std::complex>*>(this->p_hamilt->ops) + if (this->p_hamilt == nullptr) { + this->p_hamilt = new hamilt::HamiltLCAO, + std::complex>( + &this->LM, + this->kv, + *(two_center_bundle_.overlap_orb)); + dynamic_cast< + hamilt::OperatorLCAO, std::complex>*>( + this->p_hamilt->ops) ->contributeHR(); } @@ -584,14 +584,17 @@ void ESolver_KS_LCAO, std::complex>::get_S(void) std::cout << " The file is saved in " << fn << std::endl; - ModuleIO::output_SR(orb_con.ParaV, this->LM, GlobalC::GridD, this->p_hamilt, fn); + ModuleIO::output_SR(orb_con.ParaV, + this->LM, + GlobalC::GridD, + this->p_hamilt, + fn); return; } template -void ESolver_KS_LCAO::nscf(void) -{ +void ESolver_KS_LCAO::nscf(void) { ModuleBase::TITLE("ESolver_KS_LCAO", "nscf"); std::cout << " NON-SELF CONSISTENT CALCULATIONS" << std::endl; @@ -601,15 +604,12 @@ void ESolver_KS_LCAO::nscf(void) #ifdef __EXX #ifdef __MPI // Peize Lin add 2018-08-14 - if (GlobalC::exx_info.info_global.cal_exx) - { - const std::string file_name_exx = GlobalV::global_out_dir + "HexxR" + std::to_string(GlobalV::MY_RANK); - if (GlobalC::exx_info.info_ri.real_number) - { + if (GlobalC::exx_info.info_global.cal_exx) { + const std::string file_name_exx = GlobalV::global_out_dir + "HexxR" + + std::to_string(GlobalV::MY_RANK); + if (GlobalC::exx_info.info_ri.real_number) { this->exd->read_Hexxs_csr(file_name_exx, GlobalC::ucell); - } - else - { + } else { this->exc->read_Hexxs_csr(file_name_exx, GlobalC::ucell); } } @@ -621,72 +621,73 @@ void ESolver_KS_LCAO::nscf(void) // then when the istep is a variable of scf or nscf, // istep becomes istep-1, this should be fixed in future int istep = 0; - if (this->phsol != nullptr) - { - this->phsol->solve(this->p_hamilt, this->psi[0], this->pelec, GlobalV::KS_SOLVER, true); - } - else - { - ModuleBase::WARNING_QUIT("ESolver_KS_PW", "HSolver has not been initialed!"); + if (this->phsol != nullptr) { + this->phsol->solve(this->p_hamilt, + this->psi[0], + this->pelec, + GlobalV::KS_SOLVER, + true); + } else { + ModuleBase::WARNING_QUIT("ESolver_KS_PW", + "HSolver has not been initialed!"); } time_t time_finish = std::time(NULL); ModuleBase::GlobalFunc::OUT_TIME("cal_bands", time_start, time_finish); GlobalV::ofs_running << " end of band structure calculation " << std::endl; - GlobalV::ofs_running << " band eigenvalue in this processor (eV) :" << std::endl; + GlobalV::ofs_running << " band eigenvalue in this processor (eV) :" + << std::endl; 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_state " << 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_state " << 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; } - if (GlobalV::out_bandgap) - { - if (!GlobalV::TWO_EFERMI) - { + if (GlobalV::out_bandgap) { + 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; } } // add by jingan in 2018.11.7 - if (GlobalV::CALCULATION == "nscf" && INPUT.towannier90) - { + if (GlobalV::CALCULATION == "nscf" && INPUT.towannier90) { #ifdef __LCAO - if (INPUT.wannier_method == 1) - { + if (INPUT.wannier_method == 1) { toWannier90_LCAO_IN_PW myWannier(INPUT.out_wannier_mmn, INPUT.out_wannier_amn, INPUT.out_wannier_unk, @@ -702,9 +703,7 @@ void ESolver_KS_LCAO::nscf(void) this->kv, this->psi, &(this->orb_con.ParaV)); - } - else if (INPUT.wannier_method == 2) - { + } else if (INPUT.wannier_method == 2) { toWannier90_LCAO myWannier(INPUT.out_wannier_mmn, INPUT.out_wannier_amn, INPUT.out_wannier_unk, @@ -713,27 +712,35 @@ void ESolver_KS_LCAO::nscf(void) INPUT.nnkpfile, INPUT.wannier_spin); - myWannier.calculate(this->pelec->ekb, this->kv, *(this->psi), &(this->orb_con.ParaV)); + myWannier.calculate(this->pelec->ekb, + this->kv, + *(this->psi), + &(this->orb_con.ParaV)); } #endif } // add by jingan - if (berryphase::berry_phase_flag && ModuleSymmetry::Symmetry::symm_flag != 1) - { + if (berryphase::berry_phase_flag + && ModuleSymmetry::Symmetry::symm_flag != 1) { berryphase bp(this->LOC); - bp.lcao_init( - this->kv, - this->GridT); // additional step before calling macroscopic_polarization (why capitalize the function name?) - bp.Macroscopic_polarization(this->pw_wfc->npwk_max, this->psi, this->pw_rho, this->pw_wfc, this->kv); + bp.lcao_init(this->kv, + this->GridT); // additional step before calling + // macroscopic_polarization (why capitalize + // the function name?) + bp.Macroscopic_polarization(this->pw_wfc->npwk_max, + this->psi, + this->pw_rho, + this->pw_wfc, + this->kv); } // below is for DeePKS NSCF calculation #ifdef __DEEPKS - if (GlobalV::deepks_out_labels || GlobalV::deepks_scf) - { + if (GlobalV::deepks_out_labels || GlobalV::deepks_scf) { const elecstate::DensityMatrix* dm - = dynamic_cast*>(this->pelec)->get_DM(); + = dynamic_cast*>(this->pelec) + ->get_DM(); this->dpks_cal_projected_DM(dm); GlobalC::ld.cal_descriptor(GlobalC::ucell.nat); // final descriptor GlobalC::ld.cal_gedm(GlobalC::ucell.nat); @@ -743,12 +750,15 @@ void ESolver_KS_LCAO::nscf(void) this->create_Output_Mat_Sparse(0).write(); // mulliken charge analysis - if (GlobalV::out_mul) - { - elecstate::ElecStateLCAO* pelec_lcao = dynamic_cast*>(this->pelec); + if (GlobalV::out_mul) { + elecstate::ElecStateLCAO* pelec_lcao + = dynamic_cast*>(this->pelec); this->pelec->calculate_weights(); this->pelec->calEBand(); - elecstate::cal_dm_psi(&(this->orb_con.ParaV), pelec_lcao->wg, *(this->psi), *(pelec_lcao->get_DM())); + elecstate::cal_dm_psi(&(this->orb_con.ParaV), + pelec_lcao->wg, + *(this->psi), + *(pelec_lcao->get_DM())); this->cal_mag(istep, true); } diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp index e05512f81c..628e76c374 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp @@ -2,7 +2,6 @@ #include "module_io/cal_r_overlap_R.h" #include "module_io/dipole_io.h" -#include "module_io/dm_io.h" #include "module_io/rho_io.h" #include "module_io/td_current_io.h" #include "module_io/write_HS.h" @@ -32,38 +31,30 @@ //--------------------------------------------------- -namespace ModuleESolver -{ +namespace ModuleESolver { -ESolver_KS_LCAO_TDDFT::ESolver_KS_LCAO_TDDFT() -{ +ESolver_KS_LCAO_TDDFT::ESolver_KS_LCAO_TDDFT() { classname = "ESolver_KS_LCAO_TDDFT"; basisname = "LCAO"; } -ESolver_KS_LCAO_TDDFT::~ESolver_KS_LCAO_TDDFT() -{ +ESolver_KS_LCAO_TDDFT::~ESolver_KS_LCAO_TDDFT() { delete psi_laststep; - if (Hk_laststep != nullptr) - { - for (int ik = 0; ik < kv.get_nks(); ++ik) - { + if (Hk_laststep != nullptr) { + for (int ik = 0; ik < kv.get_nks(); ++ik) { delete[] Hk_laststep[ik]; } delete[] Hk_laststep; } - if (Sk_laststep != nullptr) - { - for (int ik = 0; ik < kv.get_nks(); ++ik) - { + if (Sk_laststep != nullptr) { + for (int ik = 0; ik < kv.get_nks(); ++ik) { delete[] Sk_laststep[ik]; } delete[] Sk_laststep; } } -void ESolver_KS_LCAO_TDDFT::before_all_runners(Input& inp, UnitCell& ucell) -{ +void ESolver_KS_LCAO_TDDFT::before_all_runners(Input& inp, UnitCell& ucell) { // 1) run "before_all_runners" in ESolver_KS ESolver_KS::before_all_runners(inp, ucell); @@ -71,22 +62,24 @@ void ESolver_KS_LCAO_TDDFT::before_all_runners(Input& inp, UnitCell& ucell) GlobalC::ppcell.init_vloc(GlobalC::ppcell.vloc, pw_rho); // 3) initialize the electronic states for TDDFT - if (this->pelec == nullptr) - { - this->pelec = new elecstate::ElecStateLCAO_TDDFT(&this->chr, - &kv, - kv.get_nks(), - &this->LOC, - &this->GK, // mohan add 2024-04-01 - this->pw_rho, - pw_big); + if (this->pelec == nullptr) { + this->pelec = new elecstate::ElecStateLCAO_TDDFT( + &this->chr, + &kv, + kv.get_nks(), + &this->LOC, + &this->GK, // mohan add 2024-04-01 + this->pw_rho, + pw_big); } // 4) read the local orbitals and construct the interpolation tables. this->init_basis_lcao(this->orb_con, inp, ucell); // 5) allocate H and S matrices according to computational resources - this->LM.divide_HS_in_frag(GlobalV::GAMMA_ONLY_LOCAL, orb_con.ParaV, kv.get_nks()); + this->LM.divide_HS_in_frag(GlobalV::GAMMA_ONLY_LOCAL, + orb_con.ParaV, + kv.get_nks()); // this part will be updated soon // pass Hamilt-pointer to Operator @@ -97,9 +90,9 @@ void ESolver_KS_LCAO_TDDFT::before_all_runners(Input& inp, UnitCell& ucell) ->init_DM(&kv, this->LM.ParaV, GlobalV::NSPIN); // 7) initialize Hsolver - if (this->phsol == nullptr) - { - this->phsol = new hsolver::HSolverLCAO>(this->LM.ParaV); + if (this->phsol == nullptr) { + this->phsol + = new hsolver::HSolverLCAO>(this->LM.ParaV); this->phsol->method = GlobalV::KS_SOLVER; } @@ -120,14 +113,13 @@ void ESolver_KS_LCAO_TDDFT::before_all_runners(Input& inp, UnitCell& ucell) this->pelec_td = dynamic_cast(this->pelec); } -void ESolver_KS_LCAO_TDDFT::hamilt2density(const int istep, const int iter, const double ethr) -{ +void ESolver_KS_LCAO_TDDFT::hamilt2density(const int istep, + const int iter, + const double ethr) { pelec->charge->save_rho_before_sum_band(); - if (wf.init_wfc == "file") - { - if (istep >= 1) - { + if (wf.init_wfc == "file") { + if (istep >= 1) { module_tddft::Evolve_elec::solve_psi(istep, GlobalV::NBANDS, GlobalV::NLOCAL, @@ -144,9 +136,7 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density(const int istep, const int iter, cons this->pelec_td->psiToRho_td(this->psi[0]); } this->pelec_td->psiToRho_td(this->psi[0]); - } - else if (istep >= 2) - { + } else if (istep >= 2) { module_tddft::Evolve_elec::solve_psi(istep, GlobalV::NBANDS, GlobalV::NLOCAL, @@ -161,49 +151,46 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density(const int istep, const int iter, cons INPUT.propagator, kv.get_nks()); this->pelec_td->psiToRho_td(this->psi[0]); - } - else if (this->phsol != nullptr) - { + } else if (this->phsol != nullptr) { // reset energy this->pelec->f_en.eband = 0.0; this->pelec->f_en.demet = 0.0; - if (this->psi != nullptr) - { - this->phsol->solve(this->p_hamilt, this->psi[0], this->pelec_td, GlobalV::KS_SOLVER); + if (this->psi != nullptr) { + this->phsol->solve(this->p_hamilt, + this->psi[0], + this->pelec_td, + GlobalV::KS_SOLVER); } - } - else - { - ModuleBase::WARNING_QUIT("ESolver_KS_LCAO", "HSolver has not been initialed!"); + } else { + ModuleBase::WARNING_QUIT("ESolver_KS_LCAO", + "HSolver has not been initialed!"); } // print occupation of each band - if (iter == 1 && istep <= 2) - { + if (iter == 1 && istep <= 2) { GlobalV::ofs_running - << "------------------------------------------------------------------------------------------------" + << "---------------------------------------------------------------" + "---------------------------------" << std::endl; GlobalV::ofs_running << "occupation : " << std::endl; GlobalV::ofs_running << "ik iband occ " << std::endl; GlobalV::ofs_running << std::setprecision(6); GlobalV::ofs_running << std::setiosflags(std::ios::showpoint); - for (int ik = 0; ik < kv.get_nks(); ik++) - { - for (int ib = 0; ib < GlobalV::NBANDS; ib++) - { + for (int ik = 0; ik < kv.get_nks(); ik++) { + for (int ib = 0; ib < GlobalV::NBANDS; ib++) { std::setprecision(6); - GlobalV::ofs_running << ik + 1 << " " << ib + 1 << " " << this->pelec_td->wg(ik, ib) - << std::endl; + GlobalV::ofs_running << ik + 1 << " " << ib + 1 << " " + << this->pelec_td->wg(ik, ib) << std::endl; } } GlobalV::ofs_running << std::endl; GlobalV::ofs_running - << "------------------------------------------------------------------------------------------------" + << "---------------------------------------------------------------" + "---------------------------------" << std::endl; } - for (int ik = 0; ik < kv.get_nks(); ++ik) - { + for (int ik = 0; ik < kv.get_nks(); ++ik) { this->pelec_td->print_band(ik, INPUT.printe, iter); } @@ -211,12 +198,14 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density(const int istep, const int iter, cons this->pelec->cal_energies(1); // symmetrize the charge density only for ground state - if (istep <= 1) - { + if (istep <= 1) { Symmetry_rho srho; - for (int is = 0; is < GlobalV::NSPIN; is++) - { - srho.begin(is, *(pelec->charge), pw_rho, GlobalC::Pgrid, GlobalC::ucell.symm); + for (int is = 0; is < GlobalV::NSPIN; is++) { + srho.begin(is, + *(pelec->charge), + pw_rho, + GlobalC::Pgrid, + GlobalC::ucell.symm); } } @@ -230,34 +219,29 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density(const int istep, const int iter, cons this->pelec->f_en.deband = this->pelec->cal_delta_eband(); } -void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter) -{ +void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter) { // print Hamiltonian and Overlap matrix - if (this->conv_elec) - { - if (!GlobalV::GAMMA_ONLY_LOCAL) - { + if (this->conv_elec) { + if (!GlobalV::GAMMA_ONLY_LOCAL) { this->GK.renew(true); } - for (int ik = 0; ik < kv.get_nks(); ++ik) - { - if (hsolver::HSolverLCAO>::out_mat_hs[0]) - { + for (int ik = 0; ik < kv.get_nks(); ++ik) { + if (hsolver::HSolverLCAO>::out_mat_hs[0]) { this->p_hamilt->updateHk(ik); } bool bit = false; // LiuXh, 2017-03-21 - // if set bit = true, there would be error in soc-multi-core calculation, noted by zhengdy-soc - if (this->psi != nullptr && (istep % GlobalV::out_interval == 0)) - { + // if set bit = true, there would be error in soc-multi-core + // calculation, noted by zhengdy-soc + if (this->psi != nullptr && (istep % GlobalV::out_interval == 0)) { hamilt::MatrixBlock> h_mat, s_mat; this->p_hamilt->matrix(h_mat, s_mat); - if (hsolver::HSolverLCAO>::out_mat_hs[0]) - { + if (hsolver::HSolverLCAO>::out_mat_hs[0]) { ModuleIO::save_mat(istep, h_mat.p, GlobalV::NLOCAL, bit, - hsolver::HSolverLCAO>::out_mat_hs[1], + hsolver::HSolverLCAO< + std::complex>::out_mat_hs[1], 1, GlobalV::out_app_flag, "H", @@ -269,7 +253,8 @@ void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter) s_mat.p, GlobalV::NLOCAL, bit, - hsolver::HSolverLCAO>::out_mat_hs[1], + hsolver::HSolverLCAO< + std::complex>::out_mat_hs[1], 1, GlobalV::out_app_flag, "S", @@ -281,30 +266,28 @@ void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter) } } - if (elecstate::ElecStateLCAO>::out_wfc_lcao && (this->conv_elec || iter == GlobalV::SCF_NMAX) - && (istep % GlobalV::out_interval == 0)) - { - ModuleIO::write_wfc_nao(elecstate::ElecStateLCAO>::out_wfc_lcao, - this->psi[0], - this->pelec->ekb, - this->pelec->wg, - this->pelec->klist->kvec_c, - this->orb_con.ParaV, - istep); + if (elecstate::ElecStateLCAO>::out_wfc_lcao + && (this->conv_elec || iter == GlobalV::SCF_NMAX) + && (istep % GlobalV::out_interval == 0)) { + ModuleIO::write_wfc_nao( + elecstate::ElecStateLCAO>::out_wfc_lcao, + this->psi[0], + this->pelec->ekb, + this->pelec->wg, + this->pelec->klist->kvec_c, + this->orb_con.ParaV, + istep); } // Calculate new potential according to new Charge Density - if (!this->conv_elec) - { - if (GlobalV::NSPIN == 4) - { + 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(); } @@ -315,52 +298,51 @@ void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter) const int nlocal = GlobalV::NLOCAL; // store wfc and Hk laststep - if (istep >= (wf.init_wfc == "file" ? 0 : 1) && this->conv_elec) - { - if (this->psi_laststep == nullptr) - { + if (istep >= (wf.init_wfc == "file" ? 0 : 1) && this->conv_elec) { + if (this->psi_laststep == nullptr) { #ifdef __MPI - this->psi_laststep = new psi::Psi>(kv.get_nks(), ncol_nbands, nrow, nullptr); + this->psi_laststep + = new psi::Psi>(kv.get_nks(), + ncol_nbands, + nrow, + nullptr); #else - this->psi_laststep = new psi::Psi>(kv.get_nks(), nbands, nlocal, nullptr); + this->psi_laststep + = new psi::Psi>(kv.get_nks(), + nbands, + nlocal, + nullptr); #endif } - if (td_htype == 1) - { - if (this->Hk_laststep == nullptr) - { + if (td_htype == 1) { + if (this->Hk_laststep == nullptr) { this->Hk_laststep = new std::complex*[kv.get_nks()]; - for (int ik = 0; ik < kv.get_nks(); ++ik) - { + for (int ik = 0; ik < kv.get_nks(); ++ik) { this->Hk_laststep[ik] = new std::complex[nloc]; ModuleBase::GlobalFunc::ZEROS(Hk_laststep[ik], nloc); } } - if (this->Sk_laststep == nullptr) - { + if (this->Sk_laststep == nullptr) { this->Sk_laststep = new std::complex*[kv.get_nks()]; - for (int ik = 0; ik < kv.get_nks(); ++ik) - { + for (int ik = 0; ik < kv.get_nks(); ++ik) { this->Sk_laststep[ik] = new std::complex[nloc]; ModuleBase::GlobalFunc::ZEROS(Sk_laststep[ik], nloc); } } } - for (int ik = 0; ik < kv.get_nks(); ++ik) - { + for (int ik = 0; ik < kv.get_nks(); ++ik) { this->psi->fix_k(ik); this->psi_laststep->fix_k(ik); int size0 = psi->get_nbands() * psi->get_nbasis(); - for (int index = 0; index < size0; ++index) - { - psi_laststep[0].get_pointer()[index] = psi[0].get_pointer()[index]; + for (int index = 0; index < size0; ++index) { + psi_laststep[0].get_pointer()[index] + = psi[0].get_pointer()[index]; } // store Hamiltonian - if (td_htype == 1) - { + if (td_htype == 1) { this->p_hamilt->updateHk(ik); hamilt::MatrixBlock> h_mat, s_mat; this->p_hamilt->matrix(h_mat, s_mat); @@ -370,53 +352,57 @@ void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter) } // calculate energy density matrix for tddft - if (istep >= (wf.init_wfc == "file" ? 0 : 2) && module_tddft::Evolve_elec::td_edm == 0) - { + if (istep >= (wf.init_wfc == "file" ? 0 : 2) + && module_tddft::Evolve_elec::td_edm == 0) { this->cal_edm_tddft(); } } // print "eigen value" for tddft - if (this->conv_elec) - { + if (this->conv_elec) { GlobalV::ofs_running - << "------------------------------------------------------------------------------------------------" + << "---------------------------------------------------------------" + "---------------------------------" << std::endl; GlobalV::ofs_running << "Eii : " << std::endl; GlobalV::ofs_running << "ik iband Eii (eV)" << std::endl; GlobalV::ofs_running << std::setprecision(6); GlobalV::ofs_running << std::setiosflags(std::ios::showpoint); - for (int ik = 0; ik < kv.get_nks(); ik++) - { - for (int ib = 0; ib < GlobalV::NBANDS; ib++) - { - GlobalV::ofs_running << ik + 1 << " " << ib + 1 << " " - << this->pelec_td->ekb(ik, ib) * ModuleBase::Ry_to_eV << std::endl; + for (int ik = 0; ik < kv.get_nks(); ik++) { + for (int ib = 0; ib < GlobalV::NBANDS; ib++) { + GlobalV::ofs_running + << ik + 1 << " " << ib + 1 << " " + << this->pelec_td->ekb(ik, ib) * ModuleBase::Ry_to_eV + << std::endl; } } GlobalV::ofs_running << std::endl; GlobalV::ofs_running - << "------------------------------------------------------------------------------------------------" + << "---------------------------------------------------------------" + "---------------------------------" << std::endl; } } -void ESolver_KS_LCAO_TDDFT::after_scf(const int istep) -{ - for (int is = 0; is < GlobalV::NSPIN; is++) - { - if (module_tddft::Evolve_elec::out_dipole == 1) - { +void ESolver_KS_LCAO_TDDFT::after_scf(const int istep) { + for (int is = 0; is < GlobalV::NSPIN; is++) { + if (module_tddft::Evolve_elec::out_dipole == 1) { std::stringstream ss_dipole; - ss_dipole << GlobalV::global_out_dir << "SPIN" << is + 1 << "_DIPOLE"; - ModuleIO::write_dipole(pelec->charge->rho_save[is], pelec->charge->rhopw, is, istep, ss_dipole.str()); + ss_dipole << GlobalV::global_out_dir << "SPIN" << is + 1 + << "_DIPOLE"; + ModuleIO::write_dipole(pelec->charge->rho_save[is], + pelec->charge->rhopw, + is, + istep, + ss_dipole.str()); } } - if (module_tddft::Evolve_elec::out_current == 1) - { + if (module_tddft::Evolve_elec::out_current == 1) { elecstate::DensityMatrix, double>* tmp_DM - = dynamic_cast>*>(this->pelec)->get_DM(); + = dynamic_cast>*>( + this->pelec) + ->get_DM(); ModuleIO::write_current(istep, this->psi, pelec, @@ -429,25 +415,35 @@ void ESolver_KS_LCAO_TDDFT::after_scf(const int istep) ESolver_KS_LCAO, double>::after_scf(istep); } -// use the original formula (Hamiltonian matrix) to calculate energy density matrix -void ESolver_KS_LCAO_TDDFT::cal_edm_tddft() -{ +// use the original formula (Hamiltonian matrix) to calculate energy density +// matrix +void ESolver_KS_LCAO_TDDFT::cal_edm_tddft() { // mohan add 2024-03-27 const int nlocal = GlobalV::NLOCAL; assert(nlocal >= 0); // this->LOC.edm_k_tddft.resize(kv.get_nks()); - dynamic_cast>*>(this->pelec)->get_DM()->EDMK.resize(kv.get_nks()); - for (int ik = 0; ik < kv.get_nks(); ++ik) - { + dynamic_cast>*>(this->pelec) + ->get_DM() + ->EDMK.resize(kv.get_nks()); + for (int ik = 0; ik < kv.get_nks(); ++ik) { std::complex* tmp_dmk - = dynamic_cast>*>(this->pelec)->get_DM()->get_DMK_pointer(ik); + = dynamic_cast>*>( + this->pelec) + ->get_DM() + ->get_DMK_pointer(ik); ModuleBase::ComplexMatrix& tmp_edmk - = dynamic_cast>*>(this->pelec)->get_DM()->EDMK[ik]; + = dynamic_cast>*>( + this->pelec) + ->get_DM() + ->EDMK[ik]; const Parallel_Orbitals* tmp_pv - = dynamic_cast>*>(this->pelec)->get_DM()->get_paraV_pointer(); + = dynamic_cast>*>( + this->pelec) + ->get_DM() + ->get_paraV_pointer(); #ifdef __MPI @@ -486,7 +482,14 @@ void ESolver_KS_LCAO_TDDFT::cal_edm_tddft() int info = 0; const int one_int = 1; - pzgetrf_(&nlocal, &nlocal, Sinv, &one_int, &one_int, this->orb_con.ParaV.desc, ipiv.data(), &info); + pzgetrf_(&nlocal, + &nlocal, + Sinv, + &one_int, + &one_int, + this->orb_con.ParaV.desc, + ipiv.data(), + &info); int lwork = -1; int liwork = -1; @@ -645,10 +648,8 @@ void ESolver_KS_LCAO_TDDFT::cal_edm_tddft() p_hamilt->matrix(h_mat, s_mat); // cout<<"hmat "<* + // I just use ModuleBase::ComplexMatrix temporarily, and will change it + // to complex* ModuleBase::ComplexMatrix tmp_dmk_base(nlocal, nlocal); - for (int i = 0; i < nlocal; i++) - { - for (int j = 0; j < nlocal; j++) - { + for (int i = 0; i < nlocal; i++) { + for (int j = 0; j < nlocal; j++) { tmp_dmk_base(i, j) = tmp_dmk[i * nlocal + j]; } } - tmp_edmk = 0.5 * (Sinv * Htmp * tmp_dmk_base + tmp_dmk_base * Htmp * Sinv); + tmp_edmk + = 0.5 * (Sinv * Htmp * tmp_dmk_base + tmp_dmk_base * Htmp * Sinv); delete[] work; #endif } diff --git a/source/module_io/CMakeLists.txt b/source/module_io/CMakeLists.txt index 47a4e30316..1a5337ff6f 100644 --- a/source/module_io/CMakeLists.txt +++ b/source/module_io/CMakeLists.txt @@ -53,14 +53,12 @@ if(ENABLE_LCAO) nscf_fermi_surf.cpp istate_charge.cpp istate_envelope.cpp - read_dm.cpp read_wfc_nao.cpp read_wfc_lcao.cpp write_wfc_nao.cpp - write_dm.cpp + io_dmk.cpp write_dmr.cpp dos_nao.cpp - output_dm.cpp sparse_matrix.cpp file_reader.cpp csr_reader.cpp diff --git a/source/module_io/dm_io.h b/source/module_io/dm_io.h deleted file mode 100644 index 346763b2a6..0000000000 --- a/source/module_io/dm_io.h +++ /dev/null @@ -1,44 +0,0 @@ -#ifndef DM_IO_H -#define DM_IO_H - -#include -#include "module_cell/unitcell.h" - -namespace ModuleIO -{ - -void read_dm( -#ifdef __MPI - const int nnrg, - const int* trace_lo, -#endif - const bool gamma_only_local, - const int nlocal, - const int nspin, - const int &is, - const std::string &fn, - double*** DM, - double** DM_R, - double& ef, - const UnitCell* ucell); - -void write_dm( -#ifdef __MPI - const int* trace_lo, -#endif - const int &is, - const int &iter, - const std::string &fn, - int precision, - int out_dm, - double*** DM, - const double& ef, - const UnitCell* ucell, - const int my_rank, - const int nspin, - const int nlocal); - -} - -#endif - diff --git a/source/module_io/io_dmk.cpp b/source/module_io/io_dmk.cpp new file mode 100644 index 0000000000..23776ac3d4 --- /dev/null +++ b/source/module_io/io_dmk.cpp @@ -0,0 +1,368 @@ +#include "module_io/io_dmk.h" + +#include "module_base/parallel_common.h" +#include "module_base/scalapack_connector.h" +#include "module_base/timer.h" + +/* +The format of the DMK file is as follows: +''' + + + + + + ... + ... +Direct + + +... + + +... + + + (fermi energy) + + + +... +''' + + +Example: +''' +sc + 5.29177 + 1 0 0 + 0 1 0 + 0 0 1 + H + 2 +Direct + 0 0 0.933859999999186 + 0 0 0.0661400000008143 + + 1 + -0.0883978533958687 (fermi energy) + 10 10 + + 5.773e-01 3.902e-02 1.661e-02 4.797e-17 -2.255e-17 5.773e-01 3.902e-02 +-1.661e-02 -1.461e-17 -4.414e-17 + ... + ''' + */ + +std::string ModuleIO::dmk_gen_fname(const bool gamma_only, + const int ispin, + const int ik) { + if (gamma_only) { + return "SPIN" + std::to_string(ispin + 1) + "_DM"; + } else { + // this case is not implemented now. + ModuleBase::WARNING_QUIT("dmk_gen_fname", + "Not implemented for non-gamma_only case."); + } +} + +void ModuleIO::dmk_write_ucell(std::ofstream& ofs, const UnitCell* ucell) { + // write the UnitCell information + ofs << ucell->latName << std::endl; + ofs << " " << ucell->lat0 * ModuleBase::BOHR_TO_A << std::endl; + ofs << " " << ucell->latvec.e11 << " " << ucell->latvec.e12 << " " + << ucell->latvec.e13 << std::endl; + ofs << " " << ucell->latvec.e21 << " " << ucell->latvec.e22 << " " + << ucell->latvec.e23 << std::endl; + ofs << " " << ucell->latvec.e31 << " " << ucell->latvec.e32 << " " + << ucell->latvec.e33 << std::endl; + for (int it = 0; it < ucell->ntype; it++) { + ofs << " " << ucell->atoms[it].label; + } + ofs << std::endl; + for (int it = 0; it < ucell->ntype; it++) { + ofs << " " << ucell->atoms[it].na; + } + ofs << std::endl; + ofs << "Direct" << std::endl; + for (int it = 0; it < ucell->ntype; it++) { + Atom* atom = &ucell->atoms[it]; + ofs << std::setprecision(15); + for (int ia = 0; ia < ucell->atoms[it].na; ia++) { + ofs << " " << atom->taud[ia].x << " " << atom->taud[ia].y << " " + << atom->taud[ia].z << std::endl; + } + } +} + +void ModuleIO::dmk_read_ucell(std::ifstream& ifs) { + std::string tmp; + for (int i = 0; i < 6; i++) { + std::getline(ifs, tmp); // latName + lat0 + latvec + atom label + } + std::getline(ifs, tmp); // atom number of each type + + std::istringstream iss(tmp); + int natom = 0; + int total_natom = 0; + while (iss >> natom) { + total_natom += natom; + } + for (int i = 0; i < total_natom + 1; i++) { + std::getline(ifs, tmp); // Direct + atom coordinates + } +} + +void ModuleIO::dmk_readData(std::ifstream& ifs, double& data) { ifs >> data; } + +void ModuleIO::dmk_readData(std::ifstream& ifs, std::complex& data) { + double real, imag; + ifs >> real; + ifs >> imag; + data = std::complex(real, imag); +} + +template +bool ModuleIO::read_dmk(const int nspin, + const int nk, + const Parallel_2D& pv, + const std::string& dmk_dir, + std::vector>& dmk) { + ModuleBase::TITLE("ModuleIO", "read_dmk"); + ModuleBase::timer::tick("ModuleIO", "read_dmk"); + + int my_rank = 0; +#ifdef __MPI + MPI_Comm_rank(pv.comm_2D, &my_rank); +#endif + + int nlocal = pv.get_global_row_size(); + bool gamma_only = std::is_same::value; + std::vector> dmk_global; + + // write a lambda function to check the consistency of the data + auto check_consistency = [&](const std::string& fn, + const std::string& name, + const std::string& value, + const int& target) { + if (std::stoi(value) != target) { + ModuleBase::WARNING("ModuleIO::read_dmk", + name + " is not consistent in file < " + fn + + " >."); + std::cout << name << " = " << target << ", " << name + << " in file = " << value << std::endl; + return false; + } + return true; + }; + + bool read_success = true; + std::string tmp; + if (my_rank == 0) { + dmk_global.resize(nspin * nk, std::vector(nlocal * nlocal)); + + for (int ispin = 0; ispin < nspin; ispin++) { + for (int ik = 0; ik < nk; ik++) { + std::string fn = dmk_dir + dmk_gen_fname(gamma_only, ispin, ik); + std::ifstream ifs(fn.c_str()); + + if (!ifs) { + ModuleBase::WARNING("ModuleIO::read_dmk", + "Can't open DENSITY MATRIX File < " + fn + + " >."); + read_success = false; + break; + } + + // read the UnitCell + dmk_read_ucell(ifs); + + ifs >> tmp; // nspin + if (!check_consistency(fn, "ispin", tmp, ispin + 1)) { + read_success = false; + ifs.close(); + break; + } + ifs >> tmp; + ifs >> tmp; + ifs >> tmp; // fermi energy + ifs >> tmp; // nlocal + if (!check_consistency(fn, "nlocal", tmp, nlocal)) { + read_success = false; + ifs.close(); + break; + } + ifs >> tmp; // nlocal + if (!check_consistency(fn, "nlocal", tmp, nlocal)) { + read_success = false; + ifs.close(); + break; + } + + // read the DMK data + for (int i = 0; i < nlocal; ++i) { + for (int j = 0; j < nlocal; ++j) { + dmk_readData( + ifs, + dmk_global[ik + nk * ispin][i * nlocal + j]); + } + } + ifs.close(); + } // ik + if (!read_success) { + break; + } + } // ispin + } // rank0 + +#ifdef __MPI + MPI_Bcast(&read_success, 1, MPI_C_BOOL, 0, pv.comm_2D); +#endif + + if (read_success) { +#ifdef __MPI + // seperate dmk data to each processor with 2D block distribution + dmk.resize(nspin * nk, + std::vector(pv.get_row_size() * pv.get_col_size())); + Parallel_2D pv_glb; + pv_glb.set(nlocal, nlocal, nlocal, pv.comm_2D, pv.blacs_ctxt); + for (int ik = 0; ik < nspin * nk; ik++) { + Cpxgemr2d(nlocal, + nlocal, + dmk_global[ik].data(), + 1, + 1, + pv_glb.desc, + dmk[ik].data(), + 1, + 1, + const_cast(pv.desc), + pv_glb.blacs_ctxt); + } +#else + dmk = dmk_global; +#endif + } + ModuleBase::timer::tick("ModuleIO", "read_dmk"); + return read_success; +} + +template +void ModuleIO::write_dmk(const std::vector>& dmk, + const int precision, + const std::vector& efs, + const UnitCell* ucell, + const Parallel_2D& pv) { + ModuleBase::TITLE("ModuleIO", "write_dmk"); + ModuleBase::timer::tick("ModuleIO", "write_dmk"); + + int my_rank = 0; +#ifdef __MPI + MPI_Comm_rank(pv.comm_2D, &my_rank); +#endif + + bool gamma_only = std::is_same::value; + int nlocal = pv.get_global_row_size(); + int nspin = efs.size(); + int nk = dmk.size() / nspin; + if (nk * nspin != dmk.size()) { + ModuleBase::WARNING_QUIT( + "write_dmk", + "The size of dmk is not consistent with nspin and nk."); + } + Parallel_2D pv_glb; + + // when nspin == 2, assume the order of K in dmk is K1_up, K2_up, ..., + // K1_down, K2_down, ... + for (int ispin = 0; ispin < nspin; ispin++) { + for (int ik = 0; ik < nk; ik++) { + // gather dmk[ik] to dmk_global + std::vector dmk_global(my_rank == 0 ? nlocal * nlocal : 0); +#ifdef __MPI + pv_glb.set(nlocal, nlocal, nlocal, pv.comm_2D, pv.blacs_ctxt); + Cpxgemr2d(nlocal, + nlocal, + const_cast(dmk[ik + nk * ispin].data()), + 1, + 1, + const_cast(pv.desc), + dmk_global.data(), + 1, + 1, + pv_glb.desc, + pv_glb.blacs_ctxt); +#else + dmk_global = dmk[ik + nk * ispin]; +#endif + + if (my_rank == 0) { + std::string fn = GlobalV::global_out_dir + + dmk_gen_fname(gamma_only, ispin, ik); + std::ofstream ofs(fn.c_str()); + + if (!ofs) { + ModuleBase::WARNING("ModuleIO::write_dmk", + "Can't create DENSITY MATRIX File < " + + fn + " >."); + continue; + } + + // write the UnitCell information + dmk_write_ucell(ofs, ucell); + + ofs << "\n " << dmk.size(); // nspin + ofs << "\n " << std::fixed << std::setprecision(5) << efs[ispin] + << " (fermi energy)"; + ofs << "\n " << nlocal << " " << nlocal << std::endl; + + ofs << std::setprecision(precision); + ofs << std::scientific; + for (int i = 0; i < nlocal; ++i) { + for (int j = 0; j < nlocal; ++j) { + if (j % 8 == 0) { + ofs << "\n"; + } + if (std::is_same::value) { + ofs << " " << dmk_global[i * nlocal + j]; + } else if (std::is_same, + T>::value) { + ofs << " (" << std::real(dmk_global[i * nlocal + j]) + << "," << std::imag(dmk_global[i * nlocal + j]) + << ")"; + } + } + } + ofs.close(); + } // rank0 + } // ik + } // ispin + + ModuleBase::timer::tick("ModuleIO", "write_dmk"); +} + +template bool ModuleIO::read_dmk(const int nspin, + const int nk, + const Parallel_2D& pv, + const std::string& dmk_dir, + std::vector>& dmk); + +template bool ModuleIO::read_dmk>( + const int nspin, + const int nk, + const Parallel_2D& pv, + const std::string& dmk_dir, + std::vector>>& dmk); + +template void + ModuleIO::write_dmk(const std::vector>& dmk, + const int precision, + const std::vector& efs, + const UnitCell* ucell, + const Parallel_2D& pv); + +template void ModuleIO::write_dmk>( + const std::vector>>& dmk, + const int precision, + const std::vector& efs, + const UnitCell* ucell, + const Parallel_2D& pv); \ No newline at end of file diff --git a/source/module_io/io_dmk.h b/source/module_io/io_dmk.h new file mode 100644 index 0000000000..bbd35946d8 --- /dev/null +++ b/source/module_io/io_dmk.h @@ -0,0 +1,89 @@ +#ifndef DM_IO_H +#define DM_IO_H + +#include "module_basis/module_ao/parallel_2d.h" +#include "module_cell/unitcell.h" + +#include +#include + +namespace ModuleIO { + +/** + * @brief Generates the filename for the DMK file based on the given parameters. + * + * @param gamma_only Whether the calculation is gamma_only. + * @param ispin The index of the spin component. + * @param ik The index of the k-point. + * @return The generated filename. + */ +std::string dmk_gen_fname(const bool gamma_only, const int ispin, const int ik); + +/** + * @brief Writes the unit cell information to a DMK file. + * + * @param ofs The output file stream. + * @param ucell A pointer to the UnitCell object. + */ +void dmk_write_ucell(std::ofstream& ofs, const UnitCell* ucell); + +/** + * @brief Reads the unit cell information lines in a DMK file. + * + * @param ifs The input file stream. + */ +void dmk_read_ucell(std::ifstream& ifs); + +/** + * @brief Read one double from a file. + */ +void dmk_readData(std::ifstream& ifs, double& data); + +/** + * @brief Read one complex double from a file. + */ +void dmk_readData(std::ifstream& ifs, std::complex& data); + +/** + * @brief Reads the DMK data from a file. + * + * @tparam T The type of the DMK data. + * @param nspin The number of spin components. + * @param nk The number of k-points. + * @param pv The Parallel_2D object. Will get the global size and local size + * from it, and seperate the data into different processors accordingly. + * @param dmk_dir The directory path of the DMK file. + * @param dmk A vector to store the DMK data. If use MPI parallel, the data will + * be seperated into different processors based on the Parallel_2D object. + * @return True if the DMK data is successfully read, false otherwise. + */ +template +bool read_dmk(const int nspin, + const int nk, + const Parallel_2D& pv, + const std::string& dmk_dir, + std::vector>& dmk); + +/** + * @brief Writes the DMK data to a file. + * + * @tparam T The type of the DMK data. + * @param dmk A vector containing the DMK data. The first dimension is nspin*nk, + * and the second dimension is nlocal*nlocal. DMK is parallel in 2d-block type + * if using MPI. + * @param precision The precision of the output of DMK. + * @param efs A vector containing the Fermi energies, and should have the same + * size as the number of SPIN. + * @param ucell A pointer to the UnitCell object. + * @param pv The Parallel_2D object. The 2d-block parallel information of DMK. + */ +template +void write_dmk(const std::vector>& dmk, + const int precision, + const std::vector& efs, + const UnitCell* ucell, + const Parallel_2D& pv); + +} // namespace ModuleIO + +#endif // IO_DMK_H diff --git a/source/module_io/output_dm.cpp b/source/module_io/output_dm.cpp deleted file mode 100644 index 2b08c85a7a..0000000000 --- a/source/module_io/output_dm.cpp +++ /dev/null @@ -1,55 +0,0 @@ -#include "module_io/output_dm.h" -#include "module_io/dm_io.h" - -namespace ModuleIO -{ - -Output_DM::Output_DM(const Grid_Technique& GridT, - int is, - int iter, - int precision, - int out_dm, - double*** DM, - const double& ef, - const UnitCell* ucell, - const std::string& directory, - bool gamma_only_local) - : _GridT(GridT), - _is(is), - _iter(iter), - _precision(precision), - _out_dm(out_dm), - _DM(DM), - _ef(ef), - _ucell(ucell), - _directory(directory), - _gamma_only_local(gamma_only_local) -{ - if (gamma_only_local) - { - this->_fn = this->_directory + "/SPIN" + std::to_string(this->_is + 1) + "_DM"; - } - else - { - this->_fn = this->_directory + "/SPIN" + std::to_string(this->_is + 1) + "_DM_R"; - } -} -void Output_DM::write() -{ - ModuleIO::write_dm( -#ifdef __MPI - _GridT.trace_lo.data(), -#endif - _is, - _iter, - _fn, - _precision, - _out_dm, - _DM, - _ef, - _ucell, - GlobalV::MY_RANK, - GlobalV::NSPIN, - GlobalV::NLOCAL); -} -} // namespace ModuleIO diff --git a/source/module_io/output_dm.h b/source/module_io/output_dm.h deleted file mode 100644 index b1415f4c75..0000000000 --- a/source/module_io/output_dm.h +++ /dev/null @@ -1,43 +0,0 @@ -#ifndef OUTPUT_DM_H -#define OUTPUT_DM_H - -#include "module_cell/unitcell.h" -#include "module_hamilt_lcao/module_gint/grid_technique.h" - -#include - -namespace ModuleIO -{ - -/// @brief the output interface to write the density matrix -class Output_DM -{ - public: - Output_DM(const Grid_Technique& GridT, - int is, - int iter, - int precision, - int out_dm, - double*** DM, - const double& ef, - const UnitCell* ucell, - const std::string& directory, - bool gamma_only_local); - void write(); - - private: - const Grid_Technique& _GridT; - int _is; - int _iter; - std::string _fn; - int _precision; - int _out_dm; - double*** _DM; - const double& _ef; - const UnitCell* _ucell; - const std::string& _directory; - bool _gamma_only_local; -}; -} // namespace ModuleIO - -#endif // OUTPUT_DM_H \ No newline at end of file diff --git a/source/module_io/read_dm.cpp b/source/module_io/read_dm.cpp deleted file mode 100644 index 2332866cb3..0000000000 --- a/source/module_io/read_dm.cpp +++ /dev/null @@ -1,173 +0,0 @@ -#include "module_base/parallel_common.h" -#include "module_base/timer.h" -#include "module_io/dm_io.h" - -void ModuleIO::read_dm( -#ifdef __MPI - const int nnrg, - const int* trace_lo, -#endif - const bool gamma_only_local, - const int nlocal, - const int nspin, - const int& is, - const std::string& fn, - double*** DM, - double** DM_R, - double& ef, - const UnitCell* ucell) -{ - ModuleBase::TITLE("ModuleIO", "read_dm"); - ModuleBase::timer::tick("ModuleIO", "read_dm"); - - GlobalV::ofs_running << "\n processor 0 is reading density matrix from file < " << fn << " > " << std::endl; - // xiaohui modify 2015-03-25 - // bool quit_mesia = false; - bool quit_abacus = false; - - std::ifstream ifs; - if (GlobalV::MY_RANK == 0) - { - ifs.open(fn.c_str()); - if (!ifs) - { - // xiaohui modify 2015-03-25 - // quit_mesia = true; - quit_abacus = true; - } - else - { - // if the number is not match, - // quit the program or not. - bool quit = false; - - std::string name; - ifs >> name; - - // check lattice constant, unit is Angstrom - ModuleBase::CHECK_DOUBLE(ifs, ucell->lat0 * ModuleBase::BOHR_TO_A, quit); - ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e11, quit); - ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e12, quit); - ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e13, quit); - ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e21, quit); - ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e22, quit); - ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e23, quit); - ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e31, quit); - ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e32, quit); - ModuleBase::CHECK_DOUBLE(ifs, ucell->latvec.e33, quit); - - for (int it = 0; it < ucell->ntype; it++) - { - ModuleBase::CHECK_STRING(ifs, ucell->atoms[it].label, quit); - } - - for (int it = 0; it < ucell->ntype; it++) - { - ModuleBase::CHECK_DOUBLE(ifs, ucell->atoms[it].na, quit); - } - - std::string coordinate; - ifs >> coordinate; - - for (int it = 0; it < ucell->ntype; it++) - { - for (int ia = 0; ia < ucell->atoms[it].na; ia++) - { - ModuleBase::CHECK_DOUBLE(ifs, ucell->atoms[it].taud[ia].x, quit); - ModuleBase::CHECK_DOUBLE(ifs, ucell->atoms[it].taud[ia].y, quit); - ModuleBase::CHECK_DOUBLE(ifs, ucell->atoms[it].taud[ia].z, quit); - } - } - - ModuleBase::CHECK_INT(ifs, nspin); - ModuleBase::GlobalFunc::READ_VALUE(ifs, ef); - ModuleBase::CHECK_INT(ifs, nlocal); - ModuleBase::CHECK_INT(ifs, nlocal); - } // If file exist, read in data. - } // Finish reading the first part of density matrix. - -#ifndef __MPI - GlobalV::ofs_running << " Read SPIN = " << is + 1 << " density matrix now." << std::endl; - - if (gamma_only_local) - { - for (int i = 0; i < nlocal; ++i) - { - for (int j = 0; j < nlocal; ++j) - { - ifs >> DM[is][i][j]; - } - } - } - else - { -#ifdef __MPI - ModuleBase::WARNING_QUIT("ModuleIO::read_dm", "The nnrg should not be update"); - ModuleBase::CHECK_INT(ifs, nnrg); - - for (int i = 0; i < nnrg; ++i) - { - ifs >> DM_R[is][i]; - } -#endif - } -#else - - // distribution of necessary data - // xiaohui modify 2015-03-25 - // Parallel_Common::bcast_bool(quit_mesia); - Parallel_Common::bcast_bool(quit_abacus); - // xiaohui modify 2015-03-25 - // if(quit_mesia) - if (quit_abacus) - { - ModuleBase::WARNING_QUIT("ModuleIO::read_dm", "Can not find the density matrix file."); - } - - Parallel_Common::bcast_double(ef); - - if (gamma_only_local) - { - - double* tmp = new double[nlocal]; - for (int i = 0; i < nlocal; ++i) - { - // GlobalV::ofs_running << " i=" << i << std::endl; - ModuleBase::GlobalFunc::ZEROS(tmp, nlocal); - if (GlobalV::MY_RANK == 0) - { - for (int j = 0; j < nlocal; ++j) - { - ifs >> tmp[j]; - } - } - Parallel_Common::bcast_double(tmp, nlocal); - - const int mu = trace_lo[i]; - if (mu >= 0) - { - for (int j = 0; j < nlocal; ++j) - { - const int nu = trace_lo[j]; - if (nu >= 0) - { - DM[is][mu][nu] = tmp[j]; - } - } - } - } // i - delete[] tmp; - } - else - { - ModuleBase::WARNING_QUIT("ModuleIO::read_dm", "not ready to readin DM_R"); - } -#endif - if (GlobalV::MY_RANK == 0) - ifs.close(); - - GlobalV::ofs_running << " Finish reading density matrix." << std::endl; - - ModuleBase::timer::tick("ModuleIO", "read_dm"); - return; -} diff --git a/source/module_io/test_serial/CMakeLists.txt b/source/module_io/test_serial/CMakeLists.txt index 1e6aebe635..8a643d6308 100644 --- a/source/module_io/test_serial/CMakeLists.txt +++ b/source/module_io/test_serial/CMakeLists.txt @@ -12,9 +12,9 @@ AddTest( ) AddTest( - TARGET io_dm_io + TARGET io_dmk_io LIBS ${math_libs} base device cell_info - SOURCES dm_io_test.cpp ../read_dm.cpp ../write_dm.cpp ../output.cpp + SOURCES io_dmk_test.cpp ../io_dmk.cpp ../output.cpp ../../module_basis/module_ao/parallel_2d.cpp ) AddTest( diff --git a/source/module_io/test_serial/dm_io_test.cpp b/source/module_io/test_serial/dm_io_test.cpp deleted file mode 100644 index 63ddf38a53..0000000000 --- a/source/module_io/test_serial/dm_io_test.cpp +++ /dev/null @@ -1,123 +0,0 @@ -#include "gtest/gtest.h" -#include "gmock/gmock.h" -#include "module_io/dm_io.h" -#include "module_base/global_variable.h" -#include "prepare_unitcell.h" - -#ifdef __LCAO -InfoNonlocal::InfoNonlocal(){} -InfoNonlocal::~InfoNonlocal(){} -LCAO_Orbitals::LCAO_Orbitals(){} -LCAO_Orbitals::~LCAO_Orbitals(){} -#endif -Magnetism::Magnetism() -{ - this->tot_magnetization = 0.0; - this->abs_magnetization = 0.0; - this->start_magnetization = nullptr; -} -Magnetism::~Magnetism() -{ - delete[] this->start_magnetization; -} - -/************************************************ - * unit test of read_dm and write_dm - ***********************************************/ - -/** - * - Tested Functions: - * - read_dm() - * - the function to read density matrix from file - * - the serial version without MPI - * - write_dm() - * - the function to write density matrix to file - * - the serial version without MPI - */ - -class DMIOTest : public ::testing::Test -{ -protected: - int nspin = 1; - int lgd = 26; - int nnrg = 26*26; - double*** DM; - double** DM_R; - UnitCell* ucell; - void SetUp() - { - DM = new double**[nspin]; - DM_R = new double*[nspin]; - ucell = new UnitCell; - for(int is=0; is(ifs)),std::istreambuf_iterator()); - EXPECT_THAT(str, testing::HasSubstr("0.570336288802337 (fermi energy)")); - ifs.close(); - //remove("SPIN1_DM"); -} diff --git a/source/module_io/test_serial/io_dmk_test.cpp b/source/module_io/test_serial/io_dmk_test.cpp new file mode 100644 index 0000000000..9d4baca6e5 --- /dev/null +++ b/source/module_io/test_serial/io_dmk_test.cpp @@ -0,0 +1,139 @@ +#include "module_io/io_dmk.h" + +#include "module_base/global_variable.h" +#include "prepare_unitcell.h" + +#include "gmock/gmock.h" +#include "gtest/gtest.h" + +#ifdef __LCAO +InfoNonlocal::InfoNonlocal() {} +InfoNonlocal::~InfoNonlocal() {} +LCAO_Orbitals::LCAO_Orbitals() {} +LCAO_Orbitals::~LCAO_Orbitals() {} +#endif +Magnetism::Magnetism() { + this->tot_magnetization = 0.0; + this->abs_magnetization = 0.0; + this->start_magnetization = nullptr; +} +Magnetism::~Magnetism() { delete[] this->start_magnetization; } + +/************************************************ + * unit test of read_dmk and write_dmk + ***********************************************/ + +/** + * - Tested Functions: + * - read_dmk() + * - the function to read density matrix K from file + * - the serial version without MPI + * - write_dmk() + * - the function to write density matrix K to file + * - the serial version without MPI + */ + +TEST(DMKTest, GenFileName) { + std::string fname = ModuleIO::dmk_gen_fname(true, 0, 0); + EXPECT_EQ(fname, "SPIN1_DM"); + fname = ModuleIO::dmk_gen_fname(true, 1, 1); + EXPECT_EQ(fname, "SPIN2_DM"); + + // do not support non-gamma-only case now + std::string output; + testing::internal::CaptureStdout(); + EXPECT_EXIT(ModuleIO::dmk_gen_fname(false, 2, 0), + ::testing::ExitedWithCode(0), + ""); + output = testing::internal::GetCapturedStdout(); +}; + +class DMKIOTest : public ::testing::Test { + protected: + int nspin = 2; + int nk = 1; + int nlocal = 20; + std::vector> dmk; + Parallel_2D pv; + std::vector efs; + + void SetUp() { + dmk.resize(nspin * nk, std::vector(nlocal * nlocal, 0.0)); + for (int i = 0; i < nspin * nk; i++) { + for (int j = 0; j < nlocal * nlocal; j++) { + dmk[i][j] = 1.0 * i + 0.1 * j; + } + } + + efs.resize(nspin, 0.0); + for (int i = 0; i < nspin; i++) { + efs[i] = 0.1 * i; + } + + pv.nrow = nlocal; + pv.ncol = nlocal; + + GlobalV::global_out_dir = "./"; + } +}; + +TEST_F(DMKIOTest, WriteDMK) { + UnitCell* ucell; + UcellTestPrepare utp = UcellTestLib["Si"]; + ucell = utp.SetUcellInfo(); + ModuleIO::write_dmk(dmk, 3, efs, ucell, pv); + std::ifstream ifs; + + std::string fn = "SPIN1_DM"; + ifs.open(fn); + std::string str((std::istreambuf_iterator(ifs)), + std::istreambuf_iterator()); + EXPECT_THAT(str, testing::HasSubstr("0.00000 (fermi energy)")); + EXPECT_THAT(str, testing::HasSubstr("20 20")); + EXPECT_THAT( + str, + testing::HasSubstr("0.000e+00 1.000e-01 2.000e-01 3.000e-01 4.000e-01 " + "5.000e-01 6.000e-01 7.000e-01\n")); + EXPECT_THAT( + str, + testing::HasSubstr("8.000e-01 9.000e-01 1.000e+00 1.100e+00 1.200e+00 " + "1.300e+00 1.400e+00 1.500e+00\n")); + EXPECT_THAT( + str, + testing::HasSubstr("1.600e+00 1.700e+00 1.800e+00 1.900e+00\n")); + ifs.close(); + + fn = "SPIN2_DM"; + ifs.open(fn); + str = std::string((std::istreambuf_iterator(ifs)), + std::istreambuf_iterator()); + EXPECT_THAT(str, testing::HasSubstr("0.10000 (fermi energy)")); + EXPECT_THAT(str, testing::HasSubstr("20 20")); + EXPECT_THAT( + str, + testing::HasSubstr("1.000e+00 1.100e+00 1.200e+00 1.300e+00 1.400e+00 " + "1.500e+00 1.600e+00 1.700e+00\n")); + EXPECT_THAT( + str, + testing::HasSubstr("1.800e+00 1.900e+00 2.000e+00 2.100e+00 2.200e+00 " + "2.300e+00 2.400e+00 2.500e+00\n")); + EXPECT_THAT( + str, + testing::HasSubstr("2.600e+00 2.700e+00 2.800e+00 2.900e+00\n")); + ifs.close(); + + delete ucell; + // remove the generated files + remove("SPIN1_DM"); + remove("SPIN2_DM"); +}; + +TEST_F(DMKIOTest, ReadDMK) { + pv.nrow = 26; + pv.ncol = 26; + EXPECT_TRUE(ModuleIO::read_dmk(1, 1, pv, "./support/", dmk)); + EXPECT_EQ(dmk.size(), 1); + EXPECT_EQ(dmk[0].size(), 26 * 26); + EXPECT_NEAR(dmk[0][0], 3.904e-01, 1e-6); + EXPECT_NEAR(dmk[0][25 * 26 + 25], 3.445e-02, 1e-6); +} diff --git a/source/module_io/write_dm.cpp b/source/module_io/write_dm.cpp deleted file mode 100644 index ebfd7532b6..0000000000 --- a/source/module_io/write_dm.cpp +++ /dev/null @@ -1,195 +0,0 @@ -#include "module_base/parallel_reduce.h" -#include "module_base/timer.h" -#include "module_io/dm_io.h" - -//------------------------------------------------- -// NOTE for ModuleIO::write_dm -// I will give an example here, suppose we have a 4*4 -// density matrix (symmetry) which is -// 1.1 2.3 3.6 4.2 -// 2.3 5.2 7.1 8.9 -// 3.6 7.1 1.3 3.2 -// 4.2 8.9 3.2 2.4 -// we use two processors, each one has 3 orbitals -// processor 1 has orbital index 1, 2, 4: -// ('no' means no value on this processor) -// 1.1 2.3 no 4.2 -// 2.3 5.2 no 8.9 -// no no no no -// 4.2 8.9 no 2.4 -// processor 2 has orbital index 1, 3, 4; -// 1.1 no 3.6 4.2 -// no no no no -// 3.6 no 1.3 3.2 -// 4.2 no 3.2 2.4 -// now we want to reduce them and print out, -// we plan to reduce them one row by one row, -// then for the first row, we need to set the -// temparary array to 4 (GlobalV::NLOCAL in code), -// then we reduce first row, it will become -// 2.2 2.3 3.6 8.4, -// we notice that the first element and fourth -// element is doubled, that's because the density -// may have overlap, so we need to first count -// for each element, how many times it is duplicated -// on other processors, this is why there is -// a 'count' integer array in the code. -// UPDATED BY MOHAN 2014-05-18 -void ModuleIO::write_dm( -#ifdef __MPI - const int* trace_lo, -#endif - const int &is, - const int &iter, - const std::string &fn, - int precision, - int out_dm, - double*** DM, - const double& ef, - const UnitCell* ucell, - const int my_rank, - const int nspin, - const int nlocal) -{ - ModuleBase::TITLE("ModuleIO","write_dm"); - - if (out_dm==0) - { - return; - } - - ModuleBase::timer::tick("ModuleIO","write_dm"); - - time_t start, end; - std::ofstream ofs; - - if(my_rank==0) - { - start = time(NULL); - - ofs.open(fn.c_str()); - if (!ofs) - { - ModuleBase::WARNING("ModuleIO::write_dm","Can't create DENSITY MATRIX File!"); - } - - //GlobalV::ofs_running << "\n Output charge file." << std::endl; - - ofs << ucell->latName << std::endl;//1 - ofs << " " << ucell->lat0 * ModuleBase::BOHR_TO_A << std::endl; - ofs << " " << ucell->latvec.e11 << " " << ucell->latvec.e12 << " " << ucell->latvec.e13 << std::endl; - ofs << " " << ucell->latvec.e21 << " " << ucell->latvec.e22 << " " << ucell->latvec.e23 << std::endl; - ofs << " " << ucell->latvec.e31 << " " << ucell->latvec.e32 << " " << ucell->latvec.e33 << std::endl; - for(int it=0; itntype; it++) - { - ofs << " " << ucell->atoms[it].label; - } - ofs << std::endl; - for(int it=0; itntype; it++) - { - ofs << " " << ucell->atoms[it].na; - } - ofs << std::endl; - ofs << "Direct" << std::endl; - - for(int it=0; itntype; it++) - { - Atom* atom = &ucell->atoms[it]; - ofs << std::setprecision(15); - for(int ia=0; iaatoms[it].na; ia++) - { - ofs << " " << atom->taud[ia].x - << " " << atom->taud[ia].y - << " " << atom->taud[ia].z << std::endl; - } - } - - ofs << "\n " << nspin; - ofs << "\n " << ef << " (fermi energy)"; - - ofs << "\n " << nlocal << " " << nlocal << std::endl; - - ofs << std::setprecision(precision); - ofs << std::scientific; - - } - - //ofs << "\n " << GlobalV::GAMMA_ONLY_LOCAL << " (GAMMA ONLY LOCAL)" << std::endl; -#ifndef __MPI - - for(int i=0; i tmp(nlocal,0); - std::vector count(nlocal,0); - for (int i=0; i= 0) - { - for (int j=0; j= 0) - { - count[j]=1; - } - } - } - Parallel_Reduce::reduce_all(count.data(), nlocal); - - // reduce the density matrix for 'i' line. - tmp.assign(tmp.size(),0); - if (mu >= 0) - { - for (int j=0; j=0) - { - tmp[j] = DM[is][mu][nu]; - //GlobalV::ofs_running << " dmi=" << i << " j=" << j << " " << DM[is][mu][nu] << std::endl; - } - } - } - Parallel_Reduce::reduce_all(tmp.data(), nlocal); - - if(my_rank==0) - { - for (int j=0; j0) - { - ofs << " " << tmp[j]/(double)count[j]; - } - else - { - ofs << " 0"; - } - } - } - } - std::vector().swap(tmp); - std::vector().swap(count); -#endif - if(my_rank==0) - { - end = time(NULL); - ModuleBase::GlobalFunc::OUT_TIME("write_dm",start,end); - ofs.close(); - } - ModuleBase::timer::tick("ModuleIO","write_dm"); - - return; -} diff --git a/tests/integrate/314_NO_GO_dm_out/SPIN1_DM.ref b/tests/integrate/314_NO_GO_dm_out/SPIN1_DM.ref new file mode 100644 index 0000000000..cdc3ba9105 --- /dev/null +++ b/tests/integrate/314_NO_GO_dm_out/SPIN1_DM.ref @@ -0,0 +1,35 @@ +sc + 5.29177 + 1 0 0 + 0 1 0 + 0 0 1 + H + 2 +Direct + 0 0 0.933859999999186 + 0 0 0.0661400000008143 + + 1 + -0.0883978533958687 (fermi energy) + 10 10 + + 5.773e-01 3.902e-02 1.661e-02 4.797e-17 -2.255e-17 5.773e-01 3.902e-02 -1.661e-02 + -1.461e-17 -4.414e-17 + 3.902e-02 2.637e-03 1.122e-03 3.242e-18 -1.524e-18 3.902e-02 2.637e-03 -1.122e-03 + -9.878e-19 -2.984e-18 + 1.661e-02 1.122e-03 4.777e-04 1.380e-18 -6.487e-19 1.661e-02 1.122e-03 -4.777e-04 + -4.204e-19 -1.270e-18 + 4.797e-17 3.242e-18 1.380e-18 3.986e-33 -1.874e-33 4.797e-17 3.242e-18 -1.380e-18 + -1.214e-33 -3.668e-33 + -2.255e-17 -1.524e-18 -6.487e-19 -1.874e-33 8.809e-34 -2.255e-17 -1.524e-18 6.487e-19 + 5.709e-34 1.724e-33 + 5.773e-01 3.902e-02 1.661e-02 4.797e-17 -2.255e-17 5.773e-01 3.902e-02 -1.661e-02 + -1.461e-17 -4.414e-17 + 3.902e-02 2.637e-03 1.122e-03 3.242e-18 -1.524e-18 3.902e-02 2.637e-03 -1.122e-03 + -9.878e-19 -2.984e-18 + -1.661e-02 -1.122e-03 -4.777e-04 -1.380e-18 6.487e-19 -1.661e-02 -1.122e-03 4.777e-04 + 4.204e-19 1.270e-18 + -1.461e-17 -9.878e-19 -4.204e-19 -1.214e-33 5.709e-34 -1.461e-17 -9.878e-19 4.204e-19 + 3.700e-34 1.118e-33 + -4.414e-17 -2.984e-18 -1.270e-18 -3.668e-33 1.724e-33 -4.414e-17 -2.984e-18 1.270e-18 + 1.118e-33 3.376e-33 \ No newline at end of file diff --git a/tests/integrate/314_NO_GO_dm_out/result.ref b/tests/integrate/314_NO_GO_dm_out/result.ref index 16627e485b..670d9daff7 100644 --- a/tests/integrate/314_NO_GO_dm_out/result.ref +++ b/tests/integrate/314_NO_GO_dm_out/result.ref @@ -1,3 +1,4 @@ etotref -31.58774083965451 etotperatomref -15.7938704198 +DM_different 0 totaltimeref 0.18486 diff --git a/tests/integrate/tools/catch_properties.sh b/tests/integrate/tools/catch_properties.sh index 8f5aa446e7..56ac8fb216 100755 --- a/tests/integrate/tools/catch_properties.sh +++ b/tests/integrate/tools/catch_properties.sh @@ -21,46 +21,53 @@ sum_file(){ #echo $answer #exit 0 +get_input_key_value(){ + key=$1 + inputf=$2 + value=$(awk -v key=$key '{if($1==key) a=$2} END {print a}' $inputf) + echo $value +} + file=$1 #echo $1 calculation=`grep calculation INPUT | awk '{print $2}' | sed s/[[:space:]]//g` running_path=`echo "OUT.autotest/running_$calculation"".log"` natom=`grep -En '(^|[[:space:]])TOTAL ATOM NUMBER($|[[:space:]])' $running_path | tail -1 | awk '{print $6}'` -has_force=`grep -En '(^|[[:space:]])cal_force($|[[:space:]])' INPUT | awk '{print $2}'` -has_stress=`grep -En '(^|[[:space:]])cal_stress($|[[:space:]])' INPUT | awk '{print $2}'` -has_dftu=`grep -En '(^|[[:space:]])dft_plus_u($|[[:space:]])' INPUT | awk '{print $2}'` -has_band=`grep -En '(^|[[:space:]])out_band($|[[:space:]])' INPUT | awk '{print $2}'` -has_dos=`grep -En '(^|[[:space:]])out_dos($|[[:space:]])' INPUT | awk '{print $2}'` -has_cond=`grep -En '(^|[[:space:]])cal_cond($|[[:space:]])' INPUT | awk '{print $2}'` -has_hs=`grep -En '(^|[[:space:]])out_mat_hs($|[[:space:]])' INPUT | awk '{print $2}'` -has_hs2=`grep -En '(^|[[:space:]])out_mat_hs2($|[[:space:]])' INPUT | awk '{print $2}'` -has_xc=`grep -En '(^|[[:space:]])out_mat_xc($|[[:space:]])' INPUT | awk '{print $2}'` -has_r=`grep -En '(^|[[:space:]])out_mat_r($|[[:space:]])' INPUT | awk '{print $2}'` -deepks_out_labels=`grep deepks_out_labels INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -deepks_bandgap=`grep deepks_bandgap INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -has_lowf=`grep out_wfc_lcao INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -out_app_flag=`grep out_app_flag INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -has_wfc_r=`grep out_wfc_r INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -has_wfc_pw=`grep out_wfc_pw INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -out_dm=`grep "out_dm " INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -out_mul=`grep out_mul INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -gamma_only=`grep gamma_only INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -imp_sol=`grep imp_sol INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -run_rpa=`grep rpa INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -out_pot=`grep out_pot INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -out_dm1=`grep out_dm1 INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -get_s=`grep calculation INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -out_pband=`grep out_proj_band INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -toW90=`grep towannier90 INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -has_mat_r=`grep out_mat_r INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -has_mat_t=`grep out_mat_t INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -has_mat_dh=`grep out_mat_dh INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -has_scan=`grep dft_functional INPUT | awk '{print $2}' | sed s/[[:space:]]//g` -out_chg=`grep out_chg INPUT | awk '{print $2}' | sed s/[[:space:]]//g` +has_force=$(get_input_key_value "cal_force" "INPUT") +has_stress=$(get_input_key_value "cal_stress" "INPUT") +has_dftu=$(get_input_key_value "dft_plus_u" "INPUT") +has_band=$(get_input_key_value "out_band" "INPUT") +has_dos=$(get_input_key_value "out_dos" "INPUT") +has_cond=$(get_input_key_value "cal_cond" "INPUT") +has_hs=$(get_input_key_value "out_mat_hs" "INPUT") +has_hs2=$(get_input_key_value "out_mat_hs2" "INPUT") +has_xc=$(get_input_key_value "out_mat_xc" "INPUT") +has_r=$(get_input_key_value "out_mat_r" "INPUT") +deepks_out_labels=$(get_input_key_value "deepks_out_labels" "INPUT") +deepks_bandgap=$(get_input_key_value "deepks_bandgap" "INPUT") +has_lowf=$(get_input_key_value "out_wfc_lcao" "INPUT") +out_app_flag=$(get_input_key_value "out_app_flag" "INPUT") +has_wfc_r=$(get_input_key_value "out_wfc_r" "INPUT") +has_wfc_pw=$(get_input_key_value "out_wfc_pw" "INPUT") +out_dm=$(get_input_key_value "out_dm" "INPUT") +out_mul=$(get_input_key_value "out_mul" "INPUT") +gamma_only=$(get_input_key_value "gamma_only" "INPUT") +imp_sol=$(get_input_key_value "imp_sol" "INPUT") +run_rpa=$(get_input_key_value "rpa" "INPUT") +out_pot=$(get_input_key_value "out_pot" "INPUT") +out_dm1=$(get_input_key_value "out_dm1" "INPUT") +get_s=$(get_input_key_value "calculation" "INPUT") +out_pband=$(get_input_key_value "out_proj_band" "INPUT") +toW90=$(get_input_key_value "towannier90" "INPUT") +has_mat_r=$(get_input_key_value "out_mat_r" "INPUT") +has_mat_t=$(get_input_key_value "out_mat_t" "INPUT") +has_mat_dh=$(get_input_key_value "out_mat_dh" "INPUT") +has_scan=$(get_input_key_value "dft_functional" "INPUT") +out_chg=$(get_input_key_value "out_chg" "INPUT") #echo $running_path -base=`grep -En '(^|[[:space:]])basis_type($|[[:space:]])' INPUT | awk '{print $2}' | sed s/[[:space:]]//g` +base=$(get_input_key_value "basis_type" "INPUT") word="driver_line" -symmetry=`grep -w "symmetry" INPUT | awk '{print $2}' | sed s/[[:space:]]//g` +symmetry=$(get_input_key_value "symmetry" "INPUT") test -e $1 && rm $1 #-------------------------------------------- # if NOT non-self-consistent calculations @@ -341,29 +348,14 @@ if ! test -z "$has_lowf" && [ $has_lowf == 1 ]; then fi if ! test -z "$out_dm" && [ $out_dm == 1 ]; then - dmfile=`ls OUT.autotest/ | grep "^SPIN1_DM"` + dmfile=OUT.autotest/SPIN1_DM + dmref=SPIN1_DM.ref if test -z "$dmfile"; then echo "Can't find DM files" exit 1 else - for dm in $dmfile; - do - if ! test -f OUT.autotest/$dm; then - echo "Irregular DM file found" - exit 1 - else - total_dm=$(awk 'BEGIN {sum=0.0;startline=999} - { - if(NR==7){startline=$1+14;} - else if(NR>=startline) - { - for(i=1;i<=NF;i++){sum+=sqrt($i*$i)} - } - } - END {printf"%.6f",sum}' OUT.autotest/$dm) - echo "$dm $total_dm" >>$1 - fi - done + python3 ../tools/CompareFile.py $dmref $dmfile 5 + echo "DM_different $?" >>$1 fi fi