From c56f2776ed69a385b53e1a0dcfca1baccea38273 Mon Sep 17 00:00:00 2001 From: Mohan Chen Date: Tue, 14 May 2024 09:49:26 +0800 Subject: [PATCH 1/2] Update ESolver (#4157) * add explanations in esolver_ks_lcao * add comments in esolver_ks_lcao * add notes in esolver * organize the functions in esolver_ks * add notes in esolver_ks_lcao * update some function names in ESolver, change init() to before_runner(), change run() to runner(), and change post_process() to after_all_runners() * update esolver * fix tests in md --- source/driver.cpp | 2 + source/driver_run.cpp | 15 +- source/module_esolver/esolver.h | 11 +- source/module_esolver/esolver_dp.cpp | 12 +- source/module_esolver/esolver_dp.h | 6 +- source/module_esolver/esolver_fp.cpp | 4 +- source/module_esolver/esolver_fp.h | 2 +- source/module_esolver/esolver_ks.cpp | 161 ++++++++-- source/module_esolver/esolver_ks.h | 4 +- source/module_esolver/esolver_ks_lcao.cpp | 286 ++++++++++++++---- source/module_esolver/esolver_ks_lcao.h | 4 +- .../module_esolver/esolver_ks_lcao_tddft.cpp | 4 +- source/module_esolver/esolver_ks_lcao_tddft.h | 2 +- source/module_esolver/esolver_ks_pw.cpp | 11 +- source/module_esolver/esolver_ks_pw.h | 4 +- source/module_esolver/esolver_lj.cpp | 8 +- source/module_esolver/esolver_lj.h | 6 +- source/module_esolver/esolver_of.cpp | 12 +- source/module_esolver/esolver_of.h | 6 +- source/module_esolver/esolver_sdft_pw.cpp | 6 +- source/module_esolver/esolver_sdft_pw.h | 4 +- .../module_esolver/test/esolver_dp_test.cpp | 8 +- source/module_md/md_func.cpp | 2 +- source/module_md/test/fire_test.cpp | 2 +- source/module_md/test/langevin_test.cpp | 2 +- source/module_md/test/lj_pot_test.cpp | 2 +- source/module_md/test/msst_test.cpp | 2 +- source/module_md/test/nhchain_test.cpp | 2 +- source/module_md/test/verlet_test.cpp | 2 +- source/module_relax/relax_driver.cpp | 2 +- 30 files changed, 445 insertions(+), 149 deletions(-) diff --git a/source/driver.cpp b/source/driver.cpp index 95d976a500..2dedf01b8f 100644 --- a/source/driver.cpp +++ b/source/driver.cpp @@ -29,6 +29,8 @@ void Driver::init() ModuleBase::timer::start(); // (1) read the input parameters. + // INPUT should be initalized here and then pass to atomic world, mohan 2024-05-12 + // INPUT should not be GlobalC, mohan 2024-05-12 this->reading(); // (2) welcome to the atomic world! diff --git a/source/driver_run.cpp b/source/driver_run.cpp index 434223166b..466418e47f 100644 --- a/source/driver_run.cpp +++ b/source/driver_run.cpp @@ -29,18 +29,23 @@ void Driver::driver_run(void) ModuleESolver::init_esolver(p_esolver); //! 2: setup cell and atom information + + // this warning should not be here, mohan 2024-05-22 #ifndef __LCAO if(GlobalV::BASIS_TYPE == "lcao_in_pw" || GlobalV::BASIS_TYPE == "lcao") { ModuleBase::WARNING_QUIT("driver","to use LCAO basis, compile with __LCAO"); } #endif + + // the life of ucell should begin here, mohan 2024-05-12 + // delete ucell as a GlobalC in near future GlobalC::ucell.setup_cell(GlobalV::stru_file, GlobalV::ofs_running); //! 3: initialize Esolver and fill json-structure - p_esolver->init(INPUT, GlobalC::ucell); - + p_esolver->before_all_runners(INPUT, GlobalC::ucell); + // this Json part should be moved to before_all_runners, mohan 2024-05-12 #ifdef __RAPIDJSON Json::gen_stru_wrapper(&GlobalC::ucell); #endif @@ -52,6 +57,8 @@ void Driver::driver_run(void) } else //! scf; cell relaxation; nscf; etc { + // mixed-precision should not be like this, mohan 2024-05-12, + // DEVICE should not depend on psi if (GlobalV::precision_flag == "single") { Relax_Driver rl_driver; @@ -63,9 +70,11 @@ void Driver::driver_run(void) rl_driver.relax_driver(p_esolver); } } + // "others" in ESolver should be here. + //! 5: clean up esolver - p_esolver->post_process(); + p_esolver->after_all_runners(); ModuleESolver::clean_esolver(p_esolver); ModuleBase::timer::tick("Driver", "driver_line"); diff --git a/source/module_esolver/esolver.h b/source/module_esolver/esolver.h index 142efa2a4a..cfa4ef09bf 100644 --- a/source/module_esolver/esolver.h +++ b/source/module_esolver/esolver.h @@ -20,12 +20,15 @@ class ESolver } //! initialize the energy solver by using input parameters and cell modules - virtual void init(Input& inp, UnitCell& cell) = 0; + virtual void before_all_runners(Input& inp, UnitCell& cell) = 0; //! run energy solver - virtual void run(int istep, UnitCell& cell) = 0; + virtual void runner(const int istep, UnitCell& cell) = 0; - //! deal with exx and other calculation than scf/md/relax: + //! perform post processing calculations + virtual void after_all_runners(){}; + + //! deal with exx and other calculation than scf/md/relax/cell-relax: //! such as nscf, get_wf and get_pchg virtual void others(const int istep){}; @@ -38,8 +41,6 @@ class ESolver //! calcualte stress of given cell virtual void cal_stress(ModuleBase::matrix& stress) = 0; - //! perform post processing calculations - virtual void post_process(){}; // Print current classname. void printname(); diff --git a/source/module_esolver/esolver_dp.cpp b/source/module_esolver/esolver_dp.cpp index 8352a35594..c902cf1170 100644 --- a/source/module_esolver/esolver_dp.cpp +++ b/source/module_esolver/esolver_dp.cpp @@ -25,7 +25,7 @@ namespace ModuleESolver { - void ESolver_DP::init(Input& inp, UnitCell& ucell) + void ESolver_DP::before_all_runners(Input& inp, UnitCell& ucell) { ucell_ = &ucell; dp_potential = 0; @@ -59,10 +59,10 @@ namespace ModuleESolver assert(ucell.nat == iat); } - void ESolver_DP::run(const int istep, UnitCell& ucell) + void ESolver_DP::runner(const int istep, UnitCell& ucell) { - ModuleBase::TITLE("ESolver_DP", "Run"); - ModuleBase::timer::tick("ESolver_DP", "Run"); + ModuleBase::TITLE("ESolver_DP", "runner"); + ModuleBase::timer::tick("ESolver_DP", "runner"); cell[0] = ucell.latvec.e11 * ucell.lat0_angstrom; cell[1] = ucell.latvec.e12 * ucell.lat0_angstrom; @@ -119,7 +119,7 @@ namespace ModuleESolver #else ModuleBase::WARNING_QUIT("ESolver_DP", "Please recompile with -D__DPMD"); #endif - ModuleBase::timer::tick("ESolver_DP", "Run"); + ModuleBase::timer::tick("ESolver_DP", "runner"); } double ESolver_DP::cal_energy() @@ -148,7 +148,7 @@ namespace ModuleESolver ModuleIO::print_stress("TOTAL-STRESS", stress, true, false); } - void ESolver_DP::post_process(void) + void ESolver_DP::after_all_runners(void) { GlobalV::ofs_running << "\n\n --------------------------------------------" << std::endl; GlobalV::ofs_running << std::setprecision(16); diff --git a/source/module_esolver/esolver_dp.h b/source/module_esolver/esolver_dp.h index 54008eb038..d675563986 100644 --- a/source/module_esolver/esolver_dp.h +++ b/source/module_esolver/esolver_dp.h @@ -36,7 +36,7 @@ class ESolver_DP : public ESolver * @param inp input parameters * @param cell unitcell information */ - void init(Input& inp, UnitCell& cell) override; + void before_all_runners(Input& inp, UnitCell& cell) override; /** * @brief Run the DP solver for a given ion/md step and unit cell @@ -44,7 +44,7 @@ class ESolver_DP : public ESolver * @param istep the current ion/md step * @param cell unitcell information */ - void run(const int istep, UnitCell& cell) override; + void runner(const int istep, UnitCell& cell) override; /** * @brief get the total energy without ion kinetic energy @@ -73,7 +73,7 @@ class ESolver_DP : public ESolver * * This function prints the final total energy of the DP model in eV to the output file along with some formatting. */ - void post_process() override; + void after_all_runners() override; private: /** diff --git a/source/module_esolver/esolver_fp.cpp b/source/module_esolver/esolver_fp.cpp index 22f6f32bed..e5a01f23e3 100644 --- a/source/module_esolver/esolver_fp.cpp +++ b/source/module_esolver/esolver_fp.cpp @@ -40,9 +40,9 @@ ESolver_FP::~ESolver_FP() } -void ESolver_FP::init(Input& inp, UnitCell& cell) +void ESolver_FP::before_all_runners(Input& inp, UnitCell& cell) { - ModuleBase::TITLE("ESolver_FP", "init"); + ModuleBase::TITLE("ESolver_FP", "before_all_runners"); if(!GlobalV::use_paw) { diff --git a/source/module_esolver/esolver_fp.h b/source/module_esolver/esolver_fp.h index 38e492d581..cf7cd499c8 100644 --- a/source/module_esolver/esolver_fp.h +++ b/source/module_esolver/esolver_fp.h @@ -40,7 +40,7 @@ namespace ModuleESolver virtual ~ESolver_FP(); //! Initialize of the first-principels energy solver - virtual void init(Input& inp, UnitCell& cell) override; + virtual void before_all_runners(Input& inp, UnitCell& cell) override; virtual void init_after_vc(Input& inp, UnitCell& cell); // liuyu add 2023-03-09 diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp index e69e61c296..423583c8fc 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -24,21 +24,35 @@ namespace ModuleESolver { +//------------------------------------------------------------------------------ +//! the 1st function of ESolver_KS: constructor +//! mohan add 2024-05-11 +// in future, the initialize of ESolver_KS should not be based on the +// assumption that INPUT has been initialized, mohan 2024-05-12 +//------------------------------------------------------------------------------ template ESolver_KS::ESolver_KS() { classname = "ESolver_KS"; basisname = "PLEASE ADD BASISNAME FOR CURRENT ESOLVER."; + + // should not use GlobalV here, mohan 2024-05-12 scf_thr = GlobalV::SCF_THR; drho = 0.0; + + // should not use GlobalV here, mohan 2024-05-12 maxniter = GlobalV::SCF_NMAX; niter = maxniter; + + // should not use GlobalV here, mohan 2024-05-12 out_freq_elec = GlobalV::OUT_FREQ_ELEC; // 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); + + // should not use INPUT here, mohan 2024-05-12 tmp->setbxyz(INPUT.bx,INPUT.by,INPUT.bz); ///---------------------------------------------------------- @@ -56,6 +70,10 @@ ESolver_KS::ESolver_KS() this->wf.out_wfc_r = INPUT.out_wfc_r; } +//------------------------------------------------------------------------------ +//! the 2nd function of ESolver_KS: deconstructor +//! mohan add 2024-05-11 +//------------------------------------------------------------------------------ template ESolver_KS::~ESolver_KS() { @@ -66,12 +84,16 @@ ESolver_KS::~ESolver_KS() delete this->p_chgmix; } +//------------------------------------------------------------------------------ +//! the 3rd function of ESolver_KS: before_all_runners +//! mohan add 2024-05-11 +//------------------------------------------------------------------------------ template -void ESolver_KS::init(Input& inp, UnitCell& ucell) +void ESolver_KS::before_all_runners(Input& inp, UnitCell& ucell) { - ModuleBase::TITLE("ESolver_KS", "init"); + ModuleBase::TITLE("ESolver_KS", "before_all_runners"); - ESolver_FP::init(inp,ucell); + ESolver_FP::before_all_runners(inp,ucell); //------------------Charge Mixing------------------ p_chgmix->set_mixing(GlobalV::MIXING_MODE, @@ -143,7 +165,7 @@ void ESolver_KS::init(Input& inp, UnitCell& ucell) } #endif - GlobalC::paw_cell.init_paw_cell(INPUT.ecutwfc, INPUT.cell_factor, + GlobalC::paw_cell.init_paw_cell(inp.ecutwfc, inp.cell_factor, ucell.omega,ucell.nat,ucell.ntype, atom_type,(const double **) atom_coord, filename_list); @@ -160,7 +182,7 @@ void ESolver_KS::init(Input& inp, UnitCell& ucell) ucell.cal_nelec(GlobalV::nelec); - /* it has been established that that + /* it has been established that xc_func is same for all elements, therefore only the first one if used*/ if(GlobalV::use_paw) @@ -173,6 +195,7 @@ void ESolver_KS::init(Input& inp, UnitCell& ucell) } ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SETUP UNITCELL"); + // ESolver depends on the Symmetry module // symmetry analysis should be performed every time the cell is changed if (ModuleSymmetry::Symmetry::symm_flag == 1) { @@ -188,9 +211,6 @@ void ESolver_KS::init(Input& inp, UnitCell& ucell) // mohan add 2021-01-30 Print_Info::setup_parameters(ucell, this->kv); - //if(GlobalV::BASIS_TYPE=="pw" || GlobalV::CALCULATION=="get_wf") - //{ - //Envelope function is calculated as lcao_in_pw //new plane wave basis #ifdef __MPI this->pw_wfc->initmpi(GlobalV::NPROC_IN_POOL, GlobalV::RANK_IN_POOL, POOL_WORLD); @@ -203,8 +223,9 @@ void ESolver_KS::init(Input& inp, UnitCell& ucell) this->pw_rho->nz); this->pw_wfc->initparameters(false, inp.ecutwfc, this->kv.nks, this->kv.kvec_d.data()); + // the MPI allreduce should not be here, mohan 2024-05-12 #ifdef __MPI - if (INPUT.pw_seed > 0) + if (inp.pw_seed > 0) { MPI_Allreduce(MPI_IN_PLACE, &this->pw_wfc->ggecut, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); } @@ -243,7 +264,7 @@ void ESolver_KS::init(Input& inp, UnitCell& ucell) #ifdef USE_PAW if(GlobalV::use_paw) { - GlobalC::paw_cell.set_libpaw_ecut(INPUT.ecutwfc/2.0,INPUT.ecutwfc/2.0); //in Hartree + 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, this->pw_wfc->nx,this->pw_wfc->ny,this->pw_wfc->nz, this->pw_wfc->startz,this->pw_wfc->numz); @@ -287,6 +308,10 @@ void ESolver_KS::init(Input& inp, UnitCell& ucell) } +//------------------------------------------------------------------------------ +//! the 4th function of ESolver_KS: init_after_vc +//! mohan add 2024-05-11 +//------------------------------------------------------------------------------ template void ESolver_KS::init_after_vc(Input& inp, UnitCell& ucell) { @@ -312,6 +337,10 @@ void ESolver_KS::init_after_vc(Input& inp, UnitCell& ucell) } +//------------------------------------------------------------------------------ +//! the 5th function of ESolver_KS: hamilt2density +//! mohan add 2024-05-11 +//------------------------------------------------------------------------------ template void ESolver_KS::hamilt2density(const int istep, const int iter, const double ethr) { @@ -325,6 +354,10 @@ void ESolver_KS::hamilt2density(const int istep, const int iter, cons } +//------------------------------------------------------------------------------ +//! the 6th function of ESolver_KS: print_wfcfft +//! mohan add 2024-05-11 +//------------------------------------------------------------------------------ template void ESolver_KS::print_wfcfft(Input& inp, std::ofstream &ofs) { @@ -341,7 +374,8 @@ void ESolver_KS::print_wfcfft(Input& inp, std::ofstream &ofs) ofs << " <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << std::endl; ofs << "\n\n\n\n"; ofs << "\n SETUP PLANE WAVES FOR WAVE FUNCTIONS" << std::endl; - double ecut = INPUT.ecutwfc; + + double ecut = inp.ecutwfc; 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; @@ -356,7 +390,7 @@ void ESolver_KS::print_wfcfft(Input& inp, std::ofstream &ofs) 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; - for (int i = 0; i < GlobalV::NPROC_IN_POOL ; ++i) + for (int i = 0; i < GlobalV::NPROC_IN_POOL; ++i) { ofs <<" "<pw_wfc->nst_per[i] @@ -373,11 +407,32 @@ void ESolver_KS::print_wfcfft(Input& inp, std::ofstream &ofs) } +//------------------------------------------------------------------------------ +//! the 7th function of ESolver_KS: run +//! mohan add 2024-05-11 +//! 1) run others except scf, md, relax, cell-relax +//! 2) before_scf (electronic iteration loops) +//! 3) print head +//! 4) SCF iterations +//! 5) write head +//! 6) initialization of SCF iterations +//! 7) use Hamiltonian to obtain charge density +//! 8) for MPI: STOGROUP? need to rewrite +//! 9) update potential +//! 10) finish scf iterations +//! 11) get mtaGGA related parameters +//! 12) Json, need to be moved to somewhere else +//! 13) check convergence +//! 14) add Json of efermi energy converge +//! 15) after scf +//! 16) Json again +//------------------------------------------------------------------------------ template -void ESolver_KS::run(const int istep, UnitCell& ucell) +void ESolver_KS::runner(const int istep, UnitCell& ucell) { - ModuleBase::TITLE("ESolver_KS", "run"); + ModuleBase::TITLE("ESolver_KS", "runner"); + // 1) run others except scf, md, relax, cell-relax if (!(GlobalV::CALCULATION == "scf" || GlobalV::CALCULATION == "md" || GlobalV::CALCULATION == "relax" @@ -387,13 +442,19 @@ void ESolver_KS::run(const int istep, UnitCell& ucell) } else { - ModuleBase::timer::tick(this->classname, "run"); + ModuleBase::timer::tick(this->classname, "runner"); - this->before_scf(istep); //Something else to do before the iter loop - if(GlobalV::dm_to_rho) return; //nothing further is needed + // 2) before_scf (electronic iteration loops) + this->before_scf(istep); + + if(GlobalV::dm_to_rho) + { + return; //nothing further is needed + } ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT SCF"); + // 3) print head if(this->maxniter > 0) { this->print_head(); //print the headline on the screen. @@ -402,9 +463,13 @@ void ESolver_KS::run(const int istep, UnitCell& ucell) bool firstscf = true; this->conv_elec = false; this->niter = this->maxniter; + + // 4) SCF iterations for (int iter = 1; iter <= this->maxniter; ++iter) { + // 5) write head this->write_head(GlobalV::ofs_running, istep, iter); + #ifdef __MPI auto iterstart = MPI_Wtime(); #else @@ -412,10 +477,13 @@ void ESolver_KS::run(const int istep, UnitCell& ucell) #endif double diag_ethr = this->phsol->set_diagethr(istep, iter, drho); + // 6) initialization of SCF iterations this->iter_init(istep, iter); + // 7) use Hamiltonian to obtain charge density 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.) @@ -456,9 +524,12 @@ void ESolver_KS::run(const int istep, UnitCell& ucell) bool is_U_converged = true; // to avoid unnecessary dependence on dft+u, refactor is needed #ifdef __LCAO - if (GlobalV::dft_plus_u) is_U_converged = GlobalC::dftu.u_converged(); + 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); @@ -497,9 +568,12 @@ void ESolver_KS::run(const int istep, UnitCell& ucell) MPI_Bcast(pelec->charge->rho[0], this->pw_rhod->nrxx, MPI_DOUBLE, 0, PARAPW_WORLD); #endif + // 9) update potential // Hamilt should be used after it is constructed. // this->phamilt->update(conv_elec); this->update_pot(istep, iter); + + // 10) finish scf iterations this->iter_finish(iter); #ifdef __MPI double duration = (double)(MPI_Wtime() - iterstart); @@ -509,6 +583,7 @@ void ESolver_KS::run(const int istep, UnitCell& ucell) - 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) { @@ -516,6 +591,7 @@ void ESolver_KS::run(const int istep, UnitCell& ucell) } this->print_iter(iter, drho, dkin, duration, diag_ethr); + // 12) Json, need to be moved to somewhere else #ifdef __RAPIDJSON //add Json of scf mag Json::add_output_scf_mag( @@ -527,6 +603,7 @@ void ESolver_KS::run(const int istep, UnitCell& ucell) ); #endif //__RAPIDJSON + // 13) check convergence if (this->conv_elec) { this->niter = iter; @@ -544,9 +621,11 @@ void ESolver_KS::run(const int istep, UnitCell& ucell) { std::cout<<" SCF restart after this step!"<pelec->eferm.ef * ModuleBase::Ry_to_eV, this->pelec->f_en.etot * ModuleBase::Ry_to_eV, @@ -554,11 +633,14 @@ void ESolver_KS::run(const int istep, UnitCell& ucell) ); #endif //__RAPIDJSON + // 15) after scf this->after_scf(istep); - ModuleBase::timer::tick(this->classname, "run"); - } + ModuleBase::timer::tick(this->classname, "runner"); + }// end scf, md, relax, cell-relax + + // 16) Json again #ifdef __RAPIDJSON // add nkstot,nkstot_ibz to output json int Jnkstot = this->pelec->klist->nkstot; @@ -569,6 +651,10 @@ void ESolver_KS::run(const int istep, UnitCell& ucell) }; +//------------------------------------------------------------------------------ +//! the 8th function of ESolver_KS: print_head +//! mohan add 2024-05-12 +//------------------------------------------------------------------------------ template void ESolver_KS::print_head(void) { @@ -593,6 +679,10 @@ void ESolver_KS::print_head(void) } +//------------------------------------------------------------------------------ +//! the 8th function of ESolver_KS: print_iter +//! mohan add 2024-05-12 +//------------------------------------------------------------------------------ template void ESolver_KS::print_iter( const int iter, @@ -605,6 +695,10 @@ void ESolver_KS::print_iter( } +//------------------------------------------------------------------------------ +//! the 9th function of ESolver_KS: write_head +//! mohan add 2024-05-12 +//------------------------------------------------------------------------------ template void ESolver_KS::write_head(std::ofstream& ofs_running, const int istep, const int iter) { @@ -617,6 +711,10 @@ void ESolver_KS::write_head(std::ofstream& ofs_running, const int ist } +//------------------------------------------------------------------------------ +//! the 10th function of ESolver_KS: getnieter +//! mohan add 2024-05-12 +//------------------------------------------------------------------------------ template int ESolver_KS::getniter() { @@ -624,6 +722,10 @@ int ESolver_KS::getniter() } +//------------------------------------------------------------------------------ +//! the 11th function of ESolver_KS: create_Output_Rho +//! mohan add 2024-05-12 +//------------------------------------------------------------------------------ template ModuleIO::Output_Rho ESolver_KS::create_Output_Rho( int is, @@ -661,6 +763,11 @@ ModuleIO::Output_Rho ESolver_KS::create_Output_Rho( prefix); } + +//------------------------------------------------------------------------------ +//! the 12th function of ESolver_KS: create_Output_Kin +//! mohan add 2024-05-12 +//------------------------------------------------------------------------------ template ModuleIO::Output_Rho ESolver_KS::create_Output_Kin(int is, int iter, const std::string& prefix) { @@ -681,6 +788,10 @@ ModuleIO::Output_Rho ESolver_KS::create_Output_Kin(int is, int iter, } +//------------------------------------------------------------------------------ +//! the 13th function of ESolver_KS: create_Output_Potential +//! mohan add 2024-05-12 +//------------------------------------------------------------------------------ template ModuleIO::Output_Potential ESolver_KS::create_Output_Potential(int iter, const std::string& prefix) { @@ -702,6 +813,10 @@ ModuleIO::Output_Potential ESolver_KS::create_Output_Potential(int it } +//------------------------------------------------------------------------------ +//! the 14th-18th functions of ESolver_KS +//! mohan add 2024-05-12 +//------------------------------------------------------------------------------ //! This is for mixed-precision pw/LCAO basis sets. template class ESolver_KS, psi::DEVICE_CPU>; template class ESolver_KS, psi::DEVICE_CPU>; diff --git a/source/module_esolver/esolver_ks.h b/source/module_esolver/esolver_ks.h index 63f79a60d8..cdbd392ccb 100644 --- a/source/module_esolver/esolver_ks.h +++ b/source/module_esolver/esolver_ks.h @@ -43,11 +43,11 @@ class ESolver_KS : public ESolver_FP int out_freq_elec;// frequency for output - virtual void init(Input& inp, UnitCell& cell) override; + virtual void before_all_runners(Input& inp, UnitCell& cell) override; virtual void init_after_vc(Input& inp, UnitCell& cell) override; // liuyu add 2023-03-09 - virtual void run(const int istep, UnitCell& cell) override; + virtual void runner(const int istep, UnitCell& cell) override; // calculate electron density from a specific Hamiltonian virtual void hamilt2density(const int istep, const int iter, const double ethr); diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 4fb606be36..eccf7642c1 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -45,11 +45,17 @@ namespace ModuleESolver { +//------------------------------------------------------------------------------ +//! the 1st function of ESolver_KS_LCAO: constructor +//! mohan add 2024-05-11 +//------------------------------------------------------------------------------ template 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) { @@ -67,6 +73,10 @@ ESolver_KS_LCAO::ESolver_KS_LCAO() } +//------------------------------------------------------------------------------ +//! the 2nd function of ESolver_KS_LCAO: deconstructor +//! mohan add 2024-05-11 +//------------------------------------------------------------------------------ template ESolver_KS_LCAO::~ESolver_KS_LCAO() { @@ -76,26 +86,44 @@ ESolver_KS_LCAO::~ESolver_KS_LCAO() } +//------------------------------------------------------------------------------ +//! the 3rd function of ESolver_KS_LCAO: init +//! mohan add 2024-05-11 +//! 1) calculate overlap matrix S or initialize +//! 2) init ElecState +//! 3) init LCAO basis +//! 4) redundant ParaV and LM pointers +//! 5) initialize Density Matrix +//! 6) initialize Hamilt in LCAO +//! 7) initialize exx +//! 8) Quxin added for DFT+U +//! 9) ppcell +//! 10) init HSolver +//! 11) inititlize the charge density. +//! 12) initialize the potential. +//! 13) initialize deepks +//! 14) set occupations? +//------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::init(Input& inp, UnitCell& ucell) +void ESolver_KS_LCAO::before_all_runners(Input& inp, UnitCell& ucell) { - ModuleBase::TITLE("ESolver_KS_LCAO", "init"); - ModuleBase::timer::tick("ESolver_KS_LCAO", "init"); - - // if we are only calculating S, then there is no need - // to prepare for potentials and so on - + 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") { + // 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); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY"); } - // Setup the k points according to symmetry. + // 1.3) Setup k-points according to symmetry. this->kv.set(ucell.symm, GlobalV::global_kpoint_card, GlobalV::NSPIN, ucell.G, ucell.latvec); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS"); @@ -103,54 +131,57 @@ void ESolver_KS_LCAO::init(Input& inp, UnitCell& ucell) } else { - ESolver_KS::init(inp, ucell); + // 1) else, call before_all_runners() in ESolver_KS + ESolver_KS::before_all_runners(inp, ucell); } // end ifnot get_S - // init ElecState + + // 2) init ElecState // 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), + &(this->chr), // use which parameter? &(this->kv), this->kv.nks, - &(this->LOC), + &(this->LOC), // use which parameter? &(this->GG), // mohan add 2024-04-01 &(this->GK), // mohan add 2024-04-01 - &(this->LOWF), - this->pw_rho, + &(this->LOWF), // use which parameter? + this->pw_rho, this->pw_big); } - //------------------init Basis_lcao---------------------- - // Init Basis should be put outside of Ensolver. - // * reading the localized orbitals/projectors - // * construct the interpolation tables. + // 3) init LCAO basis + // reading the localized orbitals/projectors + // construct the interpolation tables. this->init_basis_lcao(this->orb_con, inp, ucell); //------------------init Basis_lcao---------------------- - //! pass Hamilt-pointer to Operator + // 4) redundant ParaV and LM pointers this->gen_h.LM = &this->LM; //! pass basis-pointer to EState and Psi this->LOC.ParaV = this->LOWF.ParaV = this->LM.ParaV = &(this->orb_con.ParaV); - //! initialize DensityMatrix + // 5) initialize Density Matrix dynamic_cast*>(this->pelec)->init_DM(&this->kv, this->LM.ParaV, GlobalV::NSPIN); + if (GlobalV::CALCULATION == "get_S") { ModuleBase::timer::tick("ESolver_KS_LCAO", "init"); return; } - //------------------init Hamilt_lcao---------------------- + // 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.nks); - //------------------init Hamilt_lcao---------------------- #ifdef __EXX + // 7) initialize exx // PLEASE simplify the Exx_Global interface if (GlobalV::CALCULATION == "scf" || GlobalV::CALCULATION == "relax" @@ -186,29 +217,30 @@ void ESolver_KS_LCAO::init(Input& inp, UnitCell& ucell) } #endif - // Quxin added for DFT+U + // 8) Quxin added for DFT+U if (GlobalV::dft_plus_u) { GlobalC::dftu.init(ucell, this->LM, this->kv.nks); } + // 9) ppcell // output is GlobalC::ppcell.vloc 3D local pseudopotentials // without structure factors // this function belongs to cell LOOP GlobalC::ppcell.init_vloc(GlobalC::ppcell.vloc, this->pw_rho); - // init HSolver + // 10) init HSolver if (this->phsol == nullptr) { this->phsol = new hsolver::HSolverLCAO(this->LOWF.ParaV); this->phsol->method = GlobalV::KS_SOLVER; } - // Inititlize the charge density. + // 11) inititlize the charge density. this->pelec->charge->allocate(GlobalV::NSPIN); this->pelec->omega = GlobalC::ucell.omega; - // Initialize the potential. + // 12) initialize the potential. if (this->pelec->pot == nullptr) { this->pelec->pot = new elecstate::Potential(this->pw_rhod, @@ -221,6 +253,7 @@ void ESolver_KS_LCAO::init(Input& inp, UnitCell& ucell) } #ifdef __DEEPKS + // 13) initialize deepks // wenfei 2021-12-19 // if we are performing DeePKS calculations, we need to load a model if (GlobalV::deepks_scf) @@ -230,17 +263,22 @@ void ESolver_KS_LCAO::init(Input& inp, UnitCell& ucell) } #endif + // 14) set occupations? // Fix this->pelec->wg by ocp_kb if (GlobalV::ocp) { this->pelec->fixed_weights(GlobalV::ocp_kb); } - ModuleBase::timer::tick("ESolver_KS_LCAO", "init"); + ModuleBase::timer::tick("ESolver_KS_LCAO", "before_all_runners"); return; } +//------------------------------------------------------------------------------ +//! the 4th function of ESolver_KS_LCAO: init_after_vc +//! mohan add 2024-05-11 +//------------------------------------------------------------------------------ template void ESolver_KS_LCAO::init_after_vc(Input& inp, UnitCell& ucell) { @@ -287,6 +325,10 @@ void ESolver_KS_LCAO::init_after_vc(Input& inp, UnitCell& ucell) } +//------------------------------------------------------------------------------ +//! the 5th function of ESolver_KS_LCAO: cal_energy +//! mohan add 2024-05-11 +//------------------------------------------------------------------------------ template double ESolver_KS_LCAO::cal_energy() { @@ -296,6 +338,10 @@ double ESolver_KS_LCAO::cal_energy() } +//------------------------------------------------------------------------------ +//! the 6th function of ESolver_KS_LCAO: cal_force +//! mohan add 2024-05-11 +//------------------------------------------------------------------------------ template void ESolver_KS_LCAO::cal_force(ModuleBase::matrix& force) { @@ -338,6 +384,10 @@ void ESolver_KS_LCAO::cal_force(ModuleBase::matrix& force) } +//------------------------------------------------------------------------------ +//! the 7th function of ESolver_KS_LCAO: cal_stress +//! mohan add 2024-05-11 +//------------------------------------------------------------------------------ template void ESolver_KS_LCAO::cal_stress(ModuleBase::matrix& stress) { @@ -356,11 +406,15 @@ void ESolver_KS_LCAO::cal_stress(ModuleBase::matrix& stress) } +//------------------------------------------------------------------------------ +//! the 8th function of ESolver_KS_LCAO: after_all_runners +//! mohan add 2024-05-11 +//------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::post_process(void) +void ESolver_KS_LCAO::after_all_runners(void) { - ModuleBase::TITLE("ESolver_KS_LCAO", "post_process"); - ModuleBase::timer::tick("ESolver_KS_LCAO", "post_process"); + 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 << std::setprecision(16); @@ -391,7 +445,7 @@ void ESolver_KS_LCAO::post_process(void) const int nspin0 = (GlobalV::NSPIN == 2) ? 2 : 1; - if (INPUT.out_band[0]) // pengfei 2014-10-13 + if (INPUT.out_band[0]) { for (int is = 0; is < nspin0; is++) { @@ -438,10 +492,14 @@ void ESolver_KS_LCAO::post_process(void) GlobalV::NBANDS, this->p_hamilt); } - ModuleBase::timer::tick("ESolver_KS_LCAO", "post_process"); + ModuleBase::timer::tick("ESolver_KS_LCAO", "after_all_runners"); } +//------------------------------------------------------------------------------ +//! the 9th function of ESolver_KS_LCAO: init_basis_lcao +//! mohan add 2024-05-11 +//------------------------------------------------------------------------------ template void ESolver_KS_LCAO::init_basis_lcao( ORB_control& orb_con, @@ -533,6 +591,10 @@ void ESolver_KS_LCAO::init_basis_lcao( } +//------------------------------------------------------------------------------ +//! the 10th function of ESolver_KS_LCAO: iter_init +//! mohan add 2024-05-11 +//------------------------------------------------------------------------------ template void ESolver_KS_LCAO::iter_init(const int istep, const int iter) { @@ -698,14 +760,31 @@ void ESolver_KS_LCAO::iter_init(const int istep, const int iter) } +//------------------------------------------------------------------------------ +//! the 11th function of ESolver_KS_LCAO: hamilt2density +//! mohan add 2024-05-11 +//! 1) save input rho +//! 2) save density matrix DMR for mixing +//! 3) solve the Hamiltonian and output band gap +//! 4) print bands for each k-point and each band +//! 5) EXX: +//! 6) DFT+U: compute local occupation number matrix and energy correction +//! 7) DeePKS: compute delta_e +//! 8) DeltaSpin: +//! 9) use new charge density to calculate energy +//! 10) symmetrize the charge density +//! 11) compute magnetization, only for spin==2 +//! 12) calculate delta energy +//------------------------------------------------------------------------------ template void ESolver_KS_LCAO::hamilt2density(int istep, int iter, double ethr) { ModuleBase::TITLE("ESolver_KS_LCAO", "hamilt2density"); - // save input rho + // 1) save input rho this->pelec->charge->save_rho_before_sum_band(); - // save density matrix for mixing + + // 2) save density matrix DMR for mixing if (GlobalV::MIXING_RESTART > 0 && GlobalV::MIXING_DMR && this->p_chgmix->mixing_restart_count > 0) @@ -715,7 +794,7 @@ void ESolver_KS_LCAO::hamilt2density(int istep, int iter, double ethr) dm->save_DMR(); } - // using HSolverLCAO::solve() + // 3) solve the Hamiltonian and output band gap if (this->phsol != nullptr) { // reset energy @@ -741,12 +820,14 @@ void ESolver_KS_LCAO::hamilt2density(int istep, int iter, double ethr) ModuleBase::WARNING_QUIT("ESolver_KS_PW", "HSolver has not been initialed!"); } - // print ekb for each k point and each band + + // 4) print bands for each k-point and each band for (int ik = 0; ik < this->kv.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) { @@ -758,11 +839,11 @@ void ESolver_KS_LCAO::hamilt2density(int istep, int iter, double ethr) } #endif - // if DFT+U calculation is needed, this function will calculate - // the local occupation number matrix and energy correction + // 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 + // 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) @@ -776,49 +857,60 @@ void ESolver_KS_LCAO::hamilt2density(int istep, int iter, double ethr) GlobalC::dftu.output(); } + // (7) for deepks, calculate delta_e #ifdef __DEEPKS if (GlobalV::deepks_scf) { const Parallel_Orbitals* pv = this->LOWF.ParaV; + const std::vector>& dm = 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(); sc.cal_MW(iter, &(this->LM)); } - // (4) mohan add 2010-06-24 - // using new charge density. + // 9) use new charge density to calculate energy this->pelec->cal_energies(1); - // (5) symmetrize the charge density + // 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); } - // (6) compute magnetization, only for spin==2 + // 11) compute magnetization, only for spin==2 GlobalC::ucell.magnet.compute_magnetization(this->pelec->charge->nrxx, this->pelec->charge->nxyz, this->pelec->charge->rho, this->pelec->nelec_spin.data()); - // (7) calculate delta energy + // 12) calculate delta energy this->pelec->f_en.deband = this->pelec->cal_delta_eband(); } +//------------------------------------------------------------------------------ +//! the 12th function of ESolver_KS_LCAO: update_pot +//! mohan add 2024-05-11 +//! 1) print Hamiltonian and Overlap matrix (why related to update_pot()?) +//! 2) print wavefunctions (why related to update_pot()?) +//! 3) print potential +//------------------------------------------------------------------------------ template void ESolver_KS_LCAO::update_pot(const int istep, const int iter) { ModuleBase::TITLE("ESolver_KS_LCAO", "update_pot"); - // print Hamiltonian and Overlap matrix + // 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]) @@ -865,7 +957,8 @@ void ESolver_KS_LCAO::update_pot(const int istep, const int iter) } } } - // print wavefunctions + + // 2) print wavefunctions if (this->conv_elec) { if (elecstate::ElecStateLCAO::out_wfc_lcao) @@ -886,8 +979,8 @@ void ESolver_KS_LCAO::update_pot(const int istep, const int iter) elecstate::ElecStateLCAO::out_wfc_flag = 0; } } - // (9) Calculate new potential according to new Charge Density. + // 3) print potential if (this->conv_elec || iter == GlobalV::SCF_NMAX) { if (GlobalV::out_pot < 0) // mohan add 2011-10-10 @@ -895,6 +988,7 @@ void ESolver_KS_LCAO::update_pot(const int istep, const int iter) GlobalV::out_pot = -2; } } + if (!this->conv_elec) { if (GlobalV::NSPIN == 4) @@ -911,12 +1005,22 @@ void ESolver_KS_LCAO::update_pot(const int istep, const int iter) } +//------------------------------------------------------------------------------ +//! the 13th function of ESolver_KS_LCAO: iter_finish +//! mohan add 2024-05-11 +//! 1) mix density matrix +//! 2) output charge density +//! 3) output exx matrix +//! 4) output charge density and density matrix +//! 5) cal_MW? (why put it here?) +//! 6) calculate the total energy? +//------------------------------------------------------------------------------ template void ESolver_KS_LCAO::iter_finish(int iter) { ModuleBase::TITLE("ESolver_KS_LCAO", "iter_finish"); - // mix density matrix if mixing_restart + mixing_dmr + not first mixing_restart at every iter + // 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) @@ -926,9 +1030,7 @@ void ESolver_KS_LCAO::iter_finish(int iter) this->p_chgmix->mix_dmr(dm); } - //----------------------------------- - // save charge density - //----------------------------------- + // 2) save charge density // Peize Lin add 2020.04.04 if (GlobalC::restart.info_save.save_charge) { @@ -937,8 +1039,12 @@ void ESolver_KS_LCAO::iter_finish(int iter) 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; + 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 { @@ -959,9 +1065,8 @@ void ESolver_KS_LCAO::iter_finish(int iter) } } #endif - //----------------------------------- - // output charge density for tmp - //----------------------------------- + + // 4) output charge density and density matrix bool print = false; if (this->out_freq_elec && iter % this->out_freq_elec == 0) { @@ -981,6 +1086,7 @@ void ESolver_KS_LCAO::iter_finish(int iter) } } + // 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) { @@ -988,17 +1094,39 @@ void ESolver_KS_LCAO::iter_finish(int iter) sc.cal_MW(iter, &(this->LM)); } - // (11) calculate the total energy. + // 6) calculate the total energy. this->pelec->cal_energies(2); } +//------------------------------------------------------------------------------ +//! the 14th function of ESolver_KS_LCAO: after_scf +//! mohan add 2024-05-11 +//! 1) write charge difference into files for charge extrapolation +//! 2) write density matrix for sparse matrix +//! 3) write charge density +//! 4) write density matrix +//! 5) write Vxc +//! 6) write Exx matrix +//! 7) write potential +//! 8) write convergence +//! 9) write fermi energy +//! 10) write eigenvalues +//! 11) write deepks information +//! 12) write rpa information +//! 13) write HR in npz format +//! 14) write dm in npz format +//! 15) write md related +//! 16) write spin constrian MW? +//! 17) delete grid +//! 18) write quasi-orbitals +//------------------------------------------------------------------------------ template void ESolver_KS_LCAO::after_scf(const int istep) { ModuleBase::TITLE("ESolver_KS_LCAO", "after_scf"); - // save charge difference into files for charge extrapolation + // 1) write charge difference into files for charge extrapolation if (GlobalV::CALCULATION != "scf") { this->CE.save_files(istep, @@ -1010,11 +1138,13 @@ void ESolver_KS_LCAO::after_scf(const int istep) &this->sf); } + // 2) write density matrix for sparse matrix if (this->LOC.out_dm1 == 1) { this->create_Output_DM1(istep).write(); } + // 3) write charge density if (GlobalV::out_chg) { for (int is = 0; is < GlobalV::NSPIN; is++) @@ -1027,6 +1157,7 @@ void ESolver_KS_LCAO::after_scf(const int istep) } } + // 4) write density matrix if (this->LOC.out_dm) { for (int is = 0; is < GlobalV::NSPIN; is++) @@ -1035,6 +1166,7 @@ void ESolver_KS_LCAO::after_scf(const int istep) } } + // 5) write Vxc bool out_exc = true; // tmp, add parameter! if (GlobalV::out_mat_xc) { @@ -1059,6 +1191,7 @@ void ESolver_KS_LCAO::after_scf(const int istep) } #ifdef __EXX + // 6) 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); @@ -1073,16 +1206,22 @@ void ESolver_KS_LCAO::after_scf(const int istep) } #endif + // 7) write potential this->create_Output_Potential(istep).write(); + // 8) write convergence ModuleIO::output_convergence_after_scf(this->conv_elec, this->pelec->f_en.etot); + + // 9) write fermi energy ModuleIO::output_efermi(this->conv_elec, this->pelec->eferm.ef); + // 10) write eigenvalues if (GlobalV::OUT_LEVEL != "m") { this->pelec->print_eigenvalue(GlobalV::ofs_running); } + // 11) write deepks information #ifdef __DEEPKS std::shared_ptr ld_shared_ptr(&GlobalC::ld,[](LCAO_Deepks*){}); LCAO_Deepks_Interface LDI = LCAO_Deepks_Interface(ld_shared_ptr); @@ -1100,10 +1239,10 @@ void ESolver_KS_LCAO::after_scf(const int istep) dynamic_cast*>(this->pelec)->get_DM()); ModuleBase::timer::tick("ESolver_KS_LCAO", "out_deepks_labels"); - #endif - // 3. some outputs + #ifdef __EXX + // 12) write rpa information if (INPUT.rpa) { // ModuleRPA::DFT_RPA_interface rpa_interface(GlobalC::exx_info.info_global); @@ -1116,6 +1255,7 @@ void ESolver_KS_LCAO::after_scf(const int istep) } #endif + // 13) write HR in npz format if(GlobalV::out_hr_npz) { this->p_hamilt->updateHk(0); // first k point, up spin @@ -1134,6 +1274,7 @@ void ESolver_KS_LCAO::after_scf(const int istep) } } + // 14) write dm in npz format if(GlobalV::out_dm_npz) { const elecstate::DensityMatrix* dm @@ -1148,6 +1289,7 @@ void ESolver_KS_LCAO::after_scf(const int istep) } } + // 15) write md related if (!md_skip_out(GlobalV::CALCULATION, istep, GlobalV::out_interval)) { this->create_Output_Mat_Sparse(istep).write(); @@ -1158,6 +1300,7 @@ void ESolver_KS_LCAO::after_scf(const int istep) } // qifeng add 2019/9/10, jiyy modify 2023/2/27, liuyu move here 2023-04-18 } + // 16) write spin constrian MW? // spin constrain calculations, added by Tianqi Zhao. if (GlobalV::sc_mag_switch) { @@ -1166,12 +1309,13 @@ void ESolver_KS_LCAO::after_scf(const int istep) sc.print_Mag_Force(); } + // 17) delete grid if (!GlobalV::CAL_FORCE && !GlobalV::CAL_STRESS) { RA.delete_grid(); } - // quasi-orbitals, added by Yike Huang. + // 18) 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); @@ -1186,6 +1330,10 @@ void ESolver_KS_LCAO::after_scf(const int istep) } +//------------------------------------------------------------------------------ +//! the 15th function of ESolver_KS_LCAO: do_after_converge +//! mohan add 2024-05-11 +//------------------------------------------------------------------------------ template bool ESolver_KS_LCAO::do_after_converge(int& iter) { @@ -1221,6 +1369,10 @@ 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) { @@ -1239,6 +1391,10 @@ ModuleIO::Output_DM ESolver_KS_LCAO::create_Output_DM(int is, int iter) } +//------------------------------------------------------------------------------ +//! the 17th function of ESolver_KS_LCAO: create_Output_DM1 +//! mohan add 2024-05-11 +//------------------------------------------------------------------------------ template ModuleIO::Output_DM1 ESolver_KS_LCAO::create_Output_DM1(int istep) { @@ -1248,6 +1404,10 @@ ModuleIO::Output_DM1 ESolver_KS_LCAO::create_Output_DM1(int istep) } +//------------------------------------------------------------------------------ +//! the 18th function of ESolver_KS_LCAO: create_Output_Mat_Sparse +//! mohan add 2024-05-11 +//------------------------------------------------------------------------------ template ModuleIO::Output_Mat_Sparse ESolver_KS_LCAO::create_Output_Mat_Sparse(int istep) { @@ -1268,6 +1428,10 @@ ModuleIO::Output_Mat_Sparse ESolver_KS_LCAO::create_Output_Mat_Spars } +//------------------------------------------------------------------------------ +//! the 19th function of ESolver_KS_LCAO: md_skip_out +//! mohan add 2024-05-11 +//------------------------------------------------------------------------------ template bool ESolver_KS_LCAO::md_skip_out(std::string calculation, int istep, int interval) { @@ -1282,6 +1446,10 @@ bool ESolver_KS_LCAO::md_skip_out(std::string calculation, int istep, in } +//------------------------------------------------------------------------------ +//! the 20th,21th,22th functions of ESolver_KS_LCAO +//! mohan add 2024-05-11 +//------------------------------------------------------------------------------ template class ESolver_KS_LCAO; template class ESolver_KS_LCAO, double>; template class ESolver_KS_LCAO, std::complex>; diff --git a/source/module_esolver/esolver_ks_lcao.h b/source/module_esolver/esolver_ks_lcao.h index 9886081e0f..ac903cbdc1 100644 --- a/source/module_esolver/esolver_ks_lcao.h +++ b/source/module_esolver/esolver_ks_lcao.h @@ -29,7 +29,7 @@ namespace ModuleESolver ESolver_KS_LCAO(); ~ESolver_KS_LCAO(); - void init(Input& inp, UnitCell& cell) override; + void before_all_runners(Input& inp, UnitCell& cell) override; void init_after_vc(Input& inp, UnitCell& cell) override; @@ -39,7 +39,7 @@ namespace ModuleESolver void cal_stress(ModuleBase::matrix &stress) override; - void post_process() override; + void after_all_runners() override; void nscf() override; diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp index 8c2d379ca0..3626a5ecf1 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp @@ -61,9 +61,9 @@ ESolver_KS_LCAO_TDDFT::~ESolver_KS_LCAO_TDDFT() } } -void ESolver_KS_LCAO_TDDFT::init(Input& inp, UnitCell& ucell) +void ESolver_KS_LCAO_TDDFT::before_all_runners(Input& inp, UnitCell& ucell) { - ESolver_KS::init(inp, ucell); + ESolver_KS::before_all_runners(inp, ucell); // Initialize the FFT. // this function belongs to cell LOOP diff --git a/source/module_esolver/esolver_ks_lcao_tddft.h b/source/module_esolver/esolver_ks_lcao_tddft.h index 44011309f2..e10a9199ba 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.h +++ b/source/module_esolver/esolver_ks_lcao_tddft.h @@ -21,7 +21,7 @@ class ESolver_KS_LCAO_TDDFT : public ESolver_KS_LCAO, doubl ~ESolver_KS_LCAO_TDDFT(); - void init(Input& inp, UnitCell& cell) override; + void before_all_runners(Input& inp, UnitCell& cell) override; psi::Psi>* psi_laststep = nullptr; diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index 15aca22865..6d093545fb 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -220,17 +220,18 @@ void ESolver_KS_PW::Init_GlobalC(Input& inp, UnitCell& cell) template -void ESolver_KS_PW::init(Input& inp, UnitCell& ucell) +void ESolver_KS_PW::before_all_runners(Input& inp, UnitCell& ucell) { - ESolver_KS::init(inp, ucell); + // 1) call before_all_runners() of ESolver_KS + ESolver_KS::before_all_runners(inp, ucell); - // init HSolver + // 2) initialize HSolver if (this->phsol == nullptr) { this->phsol = new hsolver::HSolverPW(this->pw_wfc, &this->wf); } - // init ElecState, + // 3) initialize ElecState, if (this->pelec == nullptr) { this->pelec = new elecstate::ElecStatePW(this->pw_wfc, @@ -1196,7 +1197,7 @@ void ESolver_KS_PW::cal_stress(ModuleBase::matrix& stress) template -void ESolver_KS_PW::post_process(void) +void ESolver_KS_PW::after_all_runners(void) { GlobalV::ofs_running << "\n\n --------------------------------------------" << std::endl; diff --git a/source/module_esolver/esolver_ks_pw.h b/source/module_esolver/esolver_ks_pw.h index 41aab7892e..db59ef6de0 100644 --- a/source/module_esolver/esolver_ks_pw.h +++ b/source/module_esolver/esolver_ks_pw.h @@ -28,7 +28,7 @@ class ESolver_KS_PW : public ESolver_KS ~ESolver_KS_PW(); - void init(Input& inp, UnitCell& cell) override; + void before_all_runners(Input& inp, UnitCell& cell) override; void init_after_vc(Input& inp, UnitCell& cell) override; @@ -44,7 +44,7 @@ class ESolver_KS_PW : public ESolver_KS virtual void nscf() override; - void post_process() override; + void after_all_runners() override; /** * @brief calculate Onsager coefficients Lmn(\omega) and conductivities with Kubo-Greenwood formula diff --git a/source/module_esolver/esolver_lj.cpp b/source/module_esolver/esolver_lj.cpp index 0b5346a53f..58fcc51ca8 100644 --- a/source/module_esolver/esolver_lj.cpp +++ b/source/module_esolver/esolver_lj.cpp @@ -7,7 +7,7 @@ namespace ModuleESolver { - void ESolver_LJ::init(Input& inp, UnitCell& ucell) + void ESolver_LJ::before_all_runners(Input& inp, UnitCell& ucell) { ucell_ = &ucell; lj_potential = 0; @@ -24,7 +24,7 @@ namespace ModuleESolver lj_sigma *= ModuleBase::ANGSTROM_AU; } - void ESolver_LJ::run(const int istep, UnitCell& ucell) + void ESolver_LJ::runner(const int istep, UnitCell& ucell) { Grid_Driver grid_neigh(GlobalV::test_deconstructor, GlobalV::test_grid_driver, GlobalV::test_grid); atom_arrange::search( @@ -35,7 +35,7 @@ namespace ModuleESolver GlobalV::SEARCH_RADIUS, GlobalV::test_atom_input); - double distance; + double distance=0.0; int index = 0; // Important! potential, force, virial must be zero per step @@ -119,7 +119,7 @@ namespace ModuleESolver ModuleIO::print_stress("TOTAL-STRESS", stress, true, false); } - void ESolver_LJ::post_process(void) + void ESolver_LJ::after_all_runners(void) { GlobalV::ofs_running << "\n\n --------------------------------------------" << std::endl; GlobalV::ofs_running << std::setprecision(16); diff --git a/source/module_esolver/esolver_lj.h b/source/module_esolver/esolver_lj.h index 6963e4b854..fcac58baf3 100644 --- a/source/module_esolver/esolver_lj.h +++ b/source/module_esolver/esolver_lj.h @@ -14,9 +14,9 @@ namespace ModuleESolver classname = "ESolver_LJ"; } - void init(Input& inp, UnitCell& cell) override; + void before_all_runners(Input& inp, UnitCell& cell) override; - void run(const int istep, UnitCell& cell) override; + void runner(const int istep, UnitCell& cell) override; double cal_energy() override; @@ -24,7 +24,7 @@ namespace ModuleESolver void cal_stress(ModuleBase::matrix& stress) override; - void post_process() override; + void after_all_runners() override; private: diff --git a/source/module_esolver/esolver_of.cpp b/source/module_esolver/esolver_of.cpp index ccd6bc3b5b..02411dc7fe 100644 --- a/source/module_esolver/esolver_of.cpp +++ b/source/module_esolver/esolver_of.cpp @@ -57,9 +57,9 @@ ESolver_OF::~ESolver_OF() delete this->opt_cg_mag_; } -void ESolver_OF::init(Input& inp, UnitCell& ucell) +void ESolver_OF::before_all_runners(Input& inp, UnitCell& ucell) { - ESolver_FP::init(inp, ucell); + ESolver_FP::before_all_runners(inp, ucell); // save necessary parameters this->of_kinetic_ = inp.of_kinetic; @@ -217,9 +217,9 @@ void ESolver_OF::init_after_vc(Input& inp, UnitCell& ucell) } } -void ESolver_OF::run(int istep, UnitCell& ucell) +void ESolver_OF::runner(int istep, UnitCell& ucell) { - ModuleBase::timer::tick("ESolver_OF", "run"); + ModuleBase::timer::tick("ESolver_OF", "runner"); // get Ewald energy, initial rho and phi if necessary this->before_opt(istep, ucell); this->iter_ = 0; @@ -251,7 +251,7 @@ void ESolver_OF::run(int istep, UnitCell& ucell) this->after_opt(istep, ucell); - ModuleBase::timer::tick("ESolver_OF", "run"); + ModuleBase::timer::tick("ESolver_OF", "runner"); } /** @@ -587,7 +587,7 @@ void ESolver_OF::after_opt(const int istep, UnitCell& ucell) /** * @brief Output the FINAL_ETOT */ -void ESolver_OF::post_process() +void ESolver_OF::after_all_runners() { GlobalV::ofs_running << "\n\n --------------------------------------------" << std::endl; diff --git a/source/module_esolver/esolver_of.h b/source/module_esolver/esolver_of.h index 42722b27b1..31435ce9f8 100644 --- a/source/module_esolver/esolver_of.h +++ b/source/module_esolver/esolver_of.h @@ -20,13 +20,13 @@ class ESolver_OF : public ESolver_FP ESolver_OF(); ~ESolver_OF(); - virtual void init(Input& inp, UnitCell& ucell) override; + virtual void before_all_runners(Input& inp, UnitCell& ucell) override; virtual void init_after_vc(Input& inp, UnitCell& ucell) override; - virtual void run(int istep, UnitCell& ucell) override; + virtual void runner(const int istep, UnitCell& ucell) override; - virtual void post_process() override; + virtual void after_all_runners() override; virtual double cal_energy() override; diff --git a/source/module_esolver/esolver_sdft_pw.cpp b/source/module_esolver/esolver_sdft_pw.cpp index 342c513014..99c9dc5a4f 100644 --- a/source/module_esolver/esolver_sdft_pw.cpp +++ b/source/module_esolver/esolver_sdft_pw.cpp @@ -36,12 +36,12 @@ ESolver_SDFT_PW::~ESolver_SDFT_PW() { } -void ESolver_SDFT_PW::init(Input& inp, UnitCell& ucell) +void ESolver_SDFT_PW::before_all_runners(Input& inp, UnitCell& ucell) { this->nche_sto = inp.nche_sto; this->method_sto = inp.method_sto; - ESolver_KS::init(inp, ucell); + ESolver_KS::before_all_runners(inp, ucell); this->pelec = new elecstate::ElecStatePW_SDFT(pw_wfc, &(chr), @@ -260,7 +260,7 @@ void ESolver_SDFT_PW::cal_stress(ModuleBase::matrix& stress) } -void ESolver_SDFT_PW::post_process(void) +void ESolver_SDFT_PW::after_all_runners(void) { GlobalV::ofs_running << "\n\n --------------------------------------------" << std::endl; diff --git a/source/module_esolver/esolver_sdft_pw.h b/source/module_esolver/esolver_sdft_pw.h index cc9b2aa8b4..c5692a7303 100644 --- a/source/module_esolver/esolver_sdft_pw.h +++ b/source/module_esolver/esolver_sdft_pw.h @@ -16,7 +16,7 @@ class ESolver_SDFT_PW : public ESolver_KS_PW> ESolver_SDFT_PW(); ~ESolver_SDFT_PW(); - void init(Input& inp, UnitCell& cell) override; + void before_all_runners(Input& inp, UnitCell& cell) override; double cal_energy() override; @@ -41,7 +41,7 @@ class ESolver_SDFT_PW : public ESolver_KS_PW> virtual void after_scf(const int istep) override; - virtual void post_process() override; + virtual void after_all_runners() override; public: /** diff --git a/source/module_esolver/test/esolver_dp_test.cpp b/source/module_esolver/test/esolver_dp_test.cpp index 8f9fe878a2..9510d9e154 100644 --- a/source/module_esolver/test/esolver_dp_test.cpp +++ b/source/module_esolver/test/esolver_dp_test.cpp @@ -39,7 +39,7 @@ class ESolverDPTest : public ::testing::Test { // Initialize variables before each test esolver = new ModuleESolver::ESolver_DP("./support/case_1.pb"); - esolver->init(inp, ucell); + esolver->before_all_runners(inp, ucell); } void TearDown() override @@ -85,7 +85,7 @@ TEST_F(ESolverDPTest, InitCase2) esolver->dp_type[0] = 0; esolver->dp_type[1] = 0; esolver->dp_file = "./support/case_2.pb"; - esolver->init(inp, ucell); + esolver->before_all_runners(inp, ucell); // Check the initialized variables EXPECT_EQ(esolver->dp_type[0], 0); @@ -100,7 +100,7 @@ TEST_F(ESolverDPTest, RunWarningQuit) int istep = 0; testing::internal::CaptureStdout(); - EXPECT_EXIT(esolver->run(istep, ucell), ::testing::ExitedWithCode(0), ""); + EXPECT_EXIT(esolver->runner(istep, ucell), ::testing::ExitedWithCode(0), ""); std::string output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output, testing::HasSubstr("Please recompile with -D__DPMD")); } @@ -171,7 +171,7 @@ TEST_F(ESolverDPTest, Postprocess) // Check the results GlobalV::ofs_running.open("log"); - esolver->post_process(); + esolver->after_all_runners(); GlobalV::ofs_running.close(); std::string expected_output = "\n\n --------------------------------------------\n !FINAL_ETOT_IS 133.3358404 eV\n " diff --git a/source/module_md/md_func.cpp b/source/module_md/md_func.cpp index 5eef0129f2..e362bb361a 100644 --- a/source/module_md/md_func.cpp +++ b/source/module_md/md_func.cpp @@ -241,7 +241,7 @@ void force_virial(ModuleESolver::ESolver* p_esolver, ModuleBase::TITLE("MD_func", "force_virial"); ModuleBase::timer::tick("MD_func", "force_virial"); - p_esolver->run(istep, unit_in); + p_esolver->runner(istep, unit_in); potential = p_esolver->cal_energy(); diff --git a/source/module_md/test/fire_test.cpp b/source/module_md/test/fire_test.cpp index 96a8e5612e..4570b6bc71 100644 --- a/source/module_md/test/fire_test.cpp +++ b/source/module_md/test/fire_test.cpp @@ -46,7 +46,7 @@ class FIREtest : public testing::Test Setcell::parameters(); ModuleESolver::ESolver* p_esolver = new ModuleESolver::ESolver_LJ(); - p_esolver->init(INPUT, ucell); + p_esolver->before_all_runners(INPUT, ucell); mdrun = new FIRE(INPUT.mdp, ucell); mdrun->setup(p_esolver, GlobalV::global_readin_dir); diff --git a/source/module_md/test/langevin_test.cpp b/source/module_md/test/langevin_test.cpp index 8fc152c651..02c6ad8db6 100644 --- a/source/module_md/test/langevin_test.cpp +++ b/source/module_md/test/langevin_test.cpp @@ -46,7 +46,7 @@ class Langevin_test : public testing::Test Setcell::parameters(); ModuleESolver::ESolver* p_esolver = new ModuleESolver::ESolver_LJ(); - p_esolver->init(INPUT, ucell); + p_esolver->before_all_runners(INPUT, ucell); mdrun = new Langevin(INPUT.mdp, ucell); mdrun->setup(p_esolver, GlobalV::global_readin_dir); diff --git a/source/module_md/test/lj_pot_test.cpp b/source/module_md/test/lj_pot_test.cpp index 89a062a1b0..51a691ff5b 100644 --- a/source/module_md/test/lj_pot_test.cpp +++ b/source/module_md/test/lj_pot_test.cpp @@ -35,7 +35,7 @@ class LJ_pot_test : public testing::Test Setcell::parameters(); ModuleESolver::ESolver* p_esolver = new ModuleESolver::ESolver_LJ(); - p_esolver->init(INPUT, ucell); + p_esolver->before_all_runners(INPUT, ucell); MD_func::force_virial(p_esolver, 0, ucell, potential, force, true, stress); } diff --git a/source/module_md/test/msst_test.cpp b/source/module_md/test/msst_test.cpp index a5d3ac97e4..9a98214688 100644 --- a/source/module_md/test/msst_test.cpp +++ b/source/module_md/test/msst_test.cpp @@ -46,7 +46,7 @@ class MSST_test : public testing::Test Setcell::parameters(); ModuleESolver::ESolver* p_esolver = new ModuleESolver::ESolver_LJ(); - p_esolver->init(INPUT, ucell); + p_esolver->before_all_runners(INPUT, ucell); mdrun = new MSST(INPUT.mdp, ucell); mdrun->setup(p_esolver, GlobalV::global_readin_dir); diff --git a/source/module_md/test/nhchain_test.cpp b/source/module_md/test/nhchain_test.cpp index 49ec231081..7926e1ade7 100644 --- a/source/module_md/test/nhchain_test.cpp +++ b/source/module_md/test/nhchain_test.cpp @@ -46,7 +46,7 @@ class NHC_test : public testing::Test Setcell::parameters(); ModuleESolver::ESolver* p_esolver = new ModuleESolver::ESolver_LJ(); - p_esolver->init(INPUT, ucell); + p_esolver->before_all_runners(INPUT, ucell); INPUT.mdp.md_type = "npt"; INPUT.mdp.md_pmode = "tri"; diff --git a/source/module_md/test/verlet_test.cpp b/source/module_md/test/verlet_test.cpp index 8852743edc..73ad6c40b3 100644 --- a/source/module_md/test/verlet_test.cpp +++ b/source/module_md/test/verlet_test.cpp @@ -46,7 +46,7 @@ class Verlet_test : public testing::Test Setcell::parameters(); ModuleESolver::ESolver* p_esolver = new ModuleESolver::ESolver_LJ(); - p_esolver->init(INPUT, ucell); + p_esolver->before_all_runners(INPUT, ucell); mdrun = new Verlet(INPUT.mdp, ucell); mdrun->setup(p_esolver, GlobalV::global_readin_dir); diff --git a/source/module_relax/relax_driver.cpp b/source/module_relax/relax_driver.cpp index 28f2485bc9..42e9604ead 100644 --- a/source/module_relax/relax_driver.cpp +++ b/source/module_relax/relax_driver.cpp @@ -48,7 +48,7 @@ void Relax_Driver::relax_driver(ModuleESolver::ESolver *p_esolve #endif //__RAPIDJSON // mohan added eiter to count for the electron iteration number, 2021-01-28 - p_esolver->run(istep - 1, GlobalC::ucell); + p_esolver->runner(istep - 1, GlobalC::ucell); time_t eend = time(NULL); time_t fstart = time(NULL); From c4db90485cd1591a18467d1fe283b5481c341c4a Mon Sep 17 00:00:00 2001 From: Denghui Lu Date: Tue, 14 May 2024 15:24:20 +0800 Subject: [PATCH 2/2] fix intel compiler warnings (#4166) --- source/module_hsolver/diago_bpcg.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/source/module_hsolver/diago_bpcg.cpp b/source/module_hsolver/diago_bpcg.cpp index a1efc98a8a..0017de4e96 100644 --- a/source/module_hsolver/diago_bpcg.cpp +++ b/source/module_hsolver/diago_bpcg.cpp @@ -1,3 +1,5 @@ +#include + #include "module_hsolver/diago_bpcg.h" #include @@ -264,7 +266,7 @@ void DiagoBPCG::diag( setmem_complex_op()(this->grad_old.template data(), 0, this->n_basis * this->n_band); - setmem_var_op()(this->beta.template data(), 1E+40, this->n_band); + setmem_var_op()(this->beta.template data(), std::numeric_limits::infinity(), this->n_band); int ntry = 0; int max_iter = current_scf_iter > 1 ?