From 5bb1041494a18e336c6bc2bd994c4a0ed000fe73 Mon Sep 17 00:00:00 2001 From: mohanchen Date: Sat, 11 May 2024 22:21:06 +0800 Subject: [PATCH 1/8] add explanations in esolver_ks_lcao --- source/module_esolver/esolver_ks_lcao.cpp | 64 +++++++++++++++++------ 1 file changed, 47 insertions(+), 17 deletions(-) diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 4fb606be36f..6fe75bb63a8 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -50,6 +50,8 @@ 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) { @@ -84,7 +86,6 @@ void ESolver_KS_LCAO::init(Input& inp, UnitCell& ucell) // if we are only calculating S, then there is no need // to prepare for potentials and so on - if (GlobalV::CALCULATION == "get_S") { ucell.read_pseudo(GlobalV::ofs_running); @@ -391,7 +392,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++) { @@ -698,14 +699,29 @@ void ESolver_KS_LCAO::iter_init(const int istep, const int iter) } +//------------------------------------------------------------------------------ +//! 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 +731,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 +757,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 +776,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 +794,59 @@ 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(); } +//------------------------------------------------------------------------------ +//! mohan 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 +893,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 +915,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 +924,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) From 6a10dd61035185a82e8de2f28dbcf4539b9861b0 Mon Sep 17 00:00:00 2001 From: mohanchen Date: Sat, 11 May 2024 22:35:18 +0800 Subject: [PATCH 2/8] add comments in esolver_ks_lcao --- source/module_esolver/esolver_ks_lcao.cpp | 101 ++++++++++++++++++++-- 1 file changed, 92 insertions(+), 9 deletions(-) diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 6fe75bb63a8..a5a9333f3fb 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -45,6 +45,10 @@ namespace ModuleESolver { +//------------------------------------------------------------------------------ +//! the 1st function of ESolver_KS_LCAO: constructor +//! mohan add 2024-05-11 +//------------------------------------------------------------------------------ template ESolver_KS_LCAO::ESolver_KS_LCAO() { @@ -69,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() { @@ -78,6 +86,10 @@ ESolver_KS_LCAO::~ESolver_KS_LCAO() } +//------------------------------------------------------------------------------ +//! the 3rd function of ESolver_KS_LCAO: init +//! mohan add 2024-05-11 +//------------------------------------------------------------------------------ template void ESolver_KS_LCAO::init(Input& inp, UnitCell& ucell) { @@ -242,6 +254,10 @@ void ESolver_KS_LCAO::init(Input& inp, UnitCell& ucell) } +//------------------------------------------------------------------------------ +//! 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) { @@ -288,6 +304,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() { @@ -297,6 +317,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) { @@ -339,6 +363,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) { @@ -357,6 +385,10 @@ void ESolver_KS_LCAO::cal_stress(ModuleBase::matrix& stress) } +//------------------------------------------------------------------------------ +//! the 8th function of ESolver_KS_LCAO: post_process +//! mohan add 2024-05-11 +//------------------------------------------------------------------------------ template void ESolver_KS_LCAO::post_process(void) { @@ -443,6 +475,10 @@ void ESolver_KS_LCAO::post_process(void) } +//------------------------------------------------------------------------------ +//! 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, @@ -534,6 +570,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) { @@ -700,6 +740,8 @@ 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 @@ -836,7 +878,8 @@ void ESolver_KS_LCAO::hamilt2density(int istep, int iter, double ethr) //------------------------------------------------------------------------------ -//! mohan 2024-05-11 +//! 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 @@ -941,12 +984,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) @@ -956,9 +1009,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) { @@ -967,8 +1018,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 { @@ -989,9 +1044,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) { @@ -1011,6 +1065,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) { @@ -1018,11 +1073,15 @@ 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 +//------------------------------------------------------------------------------ template void ESolver_KS_LCAO::after_scf(const int istep) { @@ -1216,6 +1275,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) { @@ -1251,6 +1314,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) { @@ -1269,6 +1336,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) { @@ -1278,6 +1349,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) { @@ -1298,6 +1373,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) { @@ -1312,6 +1391,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>; From 24fe14bf2ed39f5968f8bb6980f42e9ae938c1df Mon Sep 17 00:00:00 2001 From: mohanchen Date: Sat, 11 May 2024 22:51:58 +0800 Subject: [PATCH 3/8] add notes in esolver --- source/module_esolver/esolver_ks.cpp | 28 +++++++++++++++ source/module_esolver/esolver_ks_lcao.cpp | 42 ++++++++++++++++++++--- 2 files changed, 66 insertions(+), 4 deletions(-) diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp index e69e61c296a..89af92f1271 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -24,6 +24,10 @@ namespace ModuleESolver { +//------------------------------------------------------------------------------ +//! the 1st function of ESolver_KS: constructor +//! mohan add 2024-05-11 +//------------------------------------------------------------------------------ template ESolver_KS::ESolver_KS() { @@ -56,6 +60,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,6 +74,10 @@ ESolver_KS::~ESolver_KS() delete this->p_chgmix; } +//------------------------------------------------------------------------------ +//! the 3rd function of ESolver_KS: init +//! mohan add 2024-05-11 +//------------------------------------------------------------------------------ template void ESolver_KS::init(Input& inp, UnitCell& ucell) { @@ -287,6 +299,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 +328,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 +345,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) { @@ -373,6 +397,10 @@ void ESolver_KS::print_wfcfft(Input& inp, std::ofstream &ofs) } +//------------------------------------------------------------------------------ +//! the 7th function of ESolver_KS: run +//! mohan add 2024-05-11 +//------------------------------------------------------------------------------ template void ESolver_KS::run(const int istep, UnitCell& ucell) { diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index a5a9333f3fb..0a87279aef9 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -1081,13 +1081,31 @@ void ESolver_KS_LCAO::iter_finish(int iter) //------------------------------------------------------------------------------ //! 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 +//! 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, @@ -1099,11 +1117,13 @@ void ESolver_KS_LCAO::after_scf(const int istep) &this->sf); } + // 2) write density 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++) @@ -1116,6 +1136,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++) @@ -1124,6 +1145,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) { @@ -1148,6 +1170,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); @@ -1162,16 +1185,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); @@ -1189,10 +1218,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); @@ -1205,6 +1234,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 @@ -1223,6 +1253,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 @@ -1237,6 +1268,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(); @@ -1247,6 +1279,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) { @@ -1255,12 +1288,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); From a260077c7b0b1828fb82277436f7b73c3e89fe18 Mon Sep 17 00:00:00 2001 From: mohanchen Date: Sun, 12 May 2024 04:43:54 +0800 Subject: [PATCH 4/8] organize the functions in esolver_ks --- source/module_esolver/esolver_ks.cpp | 91 +++++++++++++++++++++++++--- 1 file changed, 84 insertions(+), 7 deletions(-) diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp index 89af92f1271..3b2153e594f 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -400,12 +400,29 @@ 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) { ModuleBase::TITLE("ESolver_KS", "run"); + // 1) run others except scf, md, relax, cell-relax if (!(GlobalV::CALCULATION == "scf" || GlobalV::CALCULATION == "md" || GlobalV::CALCULATION == "relax" @@ -417,11 +434,17 @@ void ESolver_KS::run(const int istep, UnitCell& ucell) { ModuleBase::timer::tick(this->classname, "run"); - 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. @@ -430,9 +453,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 @@ -440,10 +467,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.) @@ -484,9 +514,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); @@ -525,9 +558,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); @@ -537,6 +573,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) { @@ -544,6 +581,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( @@ -555,6 +593,7 @@ void ESolver_KS::run(const int istep, UnitCell& ucell) ); #endif //__RAPIDJSON + // 13) check convergence if (this->conv_elec) { this->niter = iter; @@ -572,9 +611,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, @@ -582,11 +623,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"); - } + }// end scf, md, relax, cell-relax + + // 16) Json again #ifdef __RAPIDJSON // add nkstot,nkstot_ibz to output json int Jnkstot = this->pelec->klist->nkstot; @@ -597,6 +641,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) { @@ -621,6 +669,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, @@ -633,6 +685,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) { @@ -645,6 +701,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() { @@ -652,6 +712,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, @@ -689,6 +753,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) { @@ -709,6 +778,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) { @@ -730,6 +803,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>; From 7f1f1979ce93b07dafe242b300bbf87ac9d3e11b Mon Sep 17 00:00:00 2001 From: mohanchen Date: Sun, 12 May 2024 05:46:53 +0800 Subject: [PATCH 5/8] add notes in esolver_ks_lcao --- source/module_esolver/esolver_ks_lcao.cpp | 67 +++++++++++++++-------- 1 file changed, 44 insertions(+), 23 deletions(-) diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 0a87279aef9..9b8830120a8 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -89,26 +89,41 @@ 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) { 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 + + // 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"); @@ -116,54 +131,57 @@ void ESolver_KS_LCAO::init(Input& inp, UnitCell& ucell) } else { + // 1) else, call init() in ESolver_KS ESolver_KS::init(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" @@ -199,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, @@ -234,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) @@ -243,6 +263,7 @@ void ESolver_KS_LCAO::init(Input& inp, UnitCell& ucell) } #endif + // 14) set occupations? // Fix this->pelec->wg by ocp_kb if (GlobalV::ocp) { @@ -1082,7 +1103,7 @@ void ESolver_KS_LCAO::iter_finish(int iter) //! 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 +//! 2) write density matrix for sparse matrix //! 3) write charge density //! 4) write density matrix //! 5) write Vxc @@ -1117,7 +1138,7 @@ void ESolver_KS_LCAO::after_scf(const int istep) &this->sf); } - // 2) write density matrix + // 2) write density matrix for sparse matrix if (this->LOC.out_dm1 == 1) { this->create_Output_DM1(istep).write(); From 9809a78525cc201fea39dd543a66a89a3876aba3 Mon Sep 17 00:00:00 2001 From: mohanchen Date: Sun, 12 May 2024 21:47:01 +0800 Subject: [PATCH 6/8] update some function names in ESolver, change init() to before_runner(), change run() to runner(), and change post_process() to after_all_runners() --- source/driver_run.cpp | 4 ++-- source/module_esolver/esolver.h | 9 ++++---- 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 | 16 +++++++------- source/module_esolver/esolver_ks.h | 4 ++-- source/module_esolver/esolver_ks_lcao.cpp | 22 +++++++++---------- 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 | 6 ++--- 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 ++-- source/module_md/md_func.cpp | 2 +- source/module_relax/relax_driver.cpp | 2 +- 22 files changed, 73 insertions(+), 72 deletions(-) diff --git a/source/driver_run.cpp b/source/driver_run.cpp index 434223166b6..de6b6e5f9dc 100644 --- a/source/driver_run.cpp +++ b/source/driver_run.cpp @@ -38,7 +38,7 @@ void Driver::driver_run(void) 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_runner(INPUT, GlobalC::ucell); #ifdef __RAPIDJSON @@ -65,7 +65,7 @@ void Driver::driver_run(void) } //! 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 142efa2a4ab..f8e8d7fae4e 100644 --- a/source/module_esolver/esolver.h +++ b/source/module_esolver/esolver.h @@ -20,10 +20,13 @@ class ESolver } //! initialize the energy solver by using input parameters and cell modules - virtual void init(Input& inp, UnitCell& cell) = 0; + virtual void before_runner(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; + + //! perform post processing calculations + virtual void after_all_runners(){}; //! deal with exx and other calculation than scf/md/relax: //! such as nscf, get_wf and get_pchg @@ -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 8352a355943..a74d9d86c2d 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_runner(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 54008eb0389..0d154c09b64 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_runner(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 22f6f32bedd..a2fca5f9d77 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_runner(Input& inp, UnitCell& cell) { - ModuleBase::TITLE("ESolver_FP", "init"); + ModuleBase::TITLE("ESolver_FP", "before_runner"); if(!GlobalV::use_paw) { diff --git a/source/module_esolver/esolver_fp.h b/source/module_esolver/esolver_fp.h index 38e492d581b..00d067bb809 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_runner(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 3b2153e594f..9391c747dff 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -75,15 +75,15 @@ ESolver_KS::~ESolver_KS() } //------------------------------------------------------------------------------ -//! the 3rd function of ESolver_KS: init +//! the 3rd function of ESolver_KS: before_runner //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -void ESolver_KS::init(Input& inp, UnitCell& ucell) +void ESolver_KS::before_runner(Input& inp, UnitCell& ucell) { - ModuleBase::TITLE("ESolver_KS", "init"); + ModuleBase::TITLE("ESolver_KS", "before_runner"); - ESolver_FP::init(inp,ucell); + ESolver_FP::before_runner(inp,ucell); //------------------Charge Mixing------------------ p_chgmix->set_mixing(GlobalV::MIXING_MODE, @@ -418,9 +418,9 @@ void ESolver_KS::print_wfcfft(Input& inp, std::ofstream &ofs) //! 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" @@ -432,7 +432,7 @@ void ESolver_KS::run(const int istep, UnitCell& ucell) } else { - ModuleBase::timer::tick(this->classname, "run"); + ModuleBase::timer::tick(this->classname, "runner"); // 2) before_scf (electronic iteration loops) this->before_scf(istep); @@ -626,7 +626,7 @@ void ESolver_KS::run(const int istep, UnitCell& ucell) // 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 diff --git a/source/module_esolver/esolver_ks.h b/source/module_esolver/esolver_ks.h index 63f79a60d87..bec25dc9bb9 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_runner(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 9b8830120a8..c760b66097b 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -105,10 +105,10 @@ ESolver_KS_LCAO::~ESolver_KS_LCAO() //! 14) set occupations? //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::init(Input& inp, UnitCell& ucell) +void ESolver_KS_LCAO::before_runner(Input& inp, UnitCell& ucell) { - ModuleBase::TITLE("ESolver_KS_LCAO", "init"); - ModuleBase::timer::tick("ESolver_KS_LCAO", "init"); + ModuleBase::TITLE("ESolver_KS_LCAO", "before_runner"); + ModuleBase::timer::tick("ESolver_KS_LCAO", "before_runner"); // 1) calculate overlap matrix S if (GlobalV::CALCULATION == "get_S") @@ -131,8 +131,8 @@ void ESolver_KS_LCAO::init(Input& inp, UnitCell& ucell) } else { - // 1) else, call init() in ESolver_KS - ESolver_KS::init(inp, ucell); + // 1) else, call before_runner() in ESolver_KS + ESolver_KS::before_runner(inp, ucell); } // end ifnot get_S @@ -270,7 +270,7 @@ void ESolver_KS_LCAO::init(Input& inp, UnitCell& ucell) this->pelec->fixed_weights(GlobalV::ocp_kb); } - ModuleBase::timer::tick("ESolver_KS_LCAO", "init"); + ModuleBase::timer::tick("ESolver_KS_LCAO", "before_runner"); return; } @@ -407,14 +407,14 @@ void ESolver_KS_LCAO::cal_stress(ModuleBase::matrix& stress) //------------------------------------------------------------------------------ -//! the 8th function of ESolver_KS_LCAO: post_process +//! 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); @@ -492,7 +492,7 @@ 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"); } diff --git a/source/module_esolver/esolver_ks_lcao.h b/source/module_esolver/esolver_ks_lcao.h index 9886081e0fc..1d27cfae6e2 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_runner(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 3e018ca01db..7cf67f7b691 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp @@ -62,9 +62,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_runner(Input& inp, UnitCell& ucell) { - ESolver_KS::init(inp, ucell); + ESolver_KS::before_runner(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 44011309f22..1f17460e95b 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_runner(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 15aca228656..60e89744592 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -220,9 +220,9 @@ 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_runner(Input& inp, UnitCell& ucell) { - ESolver_KS::init(inp, ucell); + ESolver_KS::before_runner(inp, ucell); // init HSolver if (this->phsol == nullptr) @@ -1196,7 +1196,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 41aab7892ea..e9fba701746 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_runner(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 0b5346a53f7..a059a8205b6 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_runner(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 6963e4b8547..19b858e0fa5 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_runner(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 ccd6bc3b5b8..205e9c66b4d 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_runner(Input& inp, UnitCell& ucell) { - ESolver_FP::init(inp, ucell); + ESolver_FP::before_runner(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 42722b27b16..b64f6cd1965 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_runner(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 342c513014c..b0a30969c8f 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_runner(Input& inp, UnitCell& ucell) { this->nche_sto = inp.nche_sto; this->method_sto = inp.method_sto; - ESolver_KS::init(inp, ucell); + ESolver_KS::before_runner(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 cc9b2aa8b43..2be9228bd85 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_runner(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_md/md_func.cpp b/source/module_md/md_func.cpp index 5eef0129f24..e362bb361a7 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_relax/relax_driver.cpp b/source/module_relax/relax_driver.cpp index 28f2485bc96..42e9604ead9 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 bf460b160a25567a3d16ad18a940fd66efdcf796 Mon Sep 17 00:00:00 2001 From: mohanchen Date: Sun, 12 May 2024 22:26:37 +0800 Subject: [PATCH 7/8] update esolver --- source/driver.cpp | 2 ++ source/driver_run.cpp | 13 +++++-- source/module_esolver/esolver.h | 4 +-- source/module_esolver/esolver_dp.cpp | 2 +- source/module_esolver/esolver_dp.h | 2 +- source/module_esolver/esolver_fp.cpp | 4 +-- source/module_esolver/esolver_fp.h | 2 +- source/module_esolver/esolver_ks.cpp | 36 ++++++++++++------- source/module_esolver/esolver_ks.h | 2 +- source/module_esolver/esolver_ks_lcao.cpp | 12 +++---- source/module_esolver/esolver_ks_lcao.h | 2 +- .../module_esolver/esolver_ks_lcao_tddft.cpp | 4 +-- source/module_esolver/esolver_ks_lcao_tddft.h | 2 +- source/module_esolver/esolver_ks_pw.cpp | 9 ++--- source/module_esolver/esolver_ks_pw.h | 2 +- source/module_esolver/esolver_lj.cpp | 2 +- source/module_esolver/esolver_lj.h | 2 +- source/module_esolver/esolver_of.cpp | 4 +-- source/module_esolver/esolver_of.h | 2 +- source/module_esolver/esolver_sdft_pw.cpp | 4 +-- source/module_esolver/esolver_sdft_pw.h | 2 +- 21 files changed, 68 insertions(+), 46 deletions(-) diff --git a/source/driver.cpp b/source/driver.cpp index 95d976a5002..2dedf01b8f7 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 de6b6e5f9dc..466418e47fe 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->before_runner(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,6 +70,8 @@ void Driver::driver_run(void) rl_driver.relax_driver(p_esolver); } } + // "others" in ESolver should be here. + //! 5: clean up esolver p_esolver->after_all_runners(); diff --git a/source/module_esolver/esolver.h b/source/module_esolver/esolver.h index f8e8d7fae4e..cfa4ef09bfb 100644 --- a/source/module_esolver/esolver.h +++ b/source/module_esolver/esolver.h @@ -20,7 +20,7 @@ class ESolver } //! initialize the energy solver by using input parameters and cell modules - virtual void before_runner(Input& inp, UnitCell& cell) = 0; + virtual void before_all_runners(Input& inp, UnitCell& cell) = 0; //! run energy solver virtual void runner(const int istep, UnitCell& cell) = 0; @@ -28,7 +28,7 @@ class ESolver //! perform post processing calculations virtual void after_all_runners(){}; - //! deal with exx and other calculation than scf/md/relax: + //! 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){}; diff --git a/source/module_esolver/esolver_dp.cpp b/source/module_esolver/esolver_dp.cpp index a74d9d86c2d..c902cf11709 100644 --- a/source/module_esolver/esolver_dp.cpp +++ b/source/module_esolver/esolver_dp.cpp @@ -25,7 +25,7 @@ namespace ModuleESolver { - void ESolver_DP::before_runner(Input& inp, UnitCell& ucell) + void ESolver_DP::before_all_runners(Input& inp, UnitCell& ucell) { ucell_ = &ucell; dp_potential = 0; diff --git a/source/module_esolver/esolver_dp.h b/source/module_esolver/esolver_dp.h index 0d154c09b64..d675563986f 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 before_runner(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 diff --git a/source/module_esolver/esolver_fp.cpp b/source/module_esolver/esolver_fp.cpp index a2fca5f9d77..e5a01f23e34 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::before_runner(Input& inp, UnitCell& cell) +void ESolver_FP::before_all_runners(Input& inp, UnitCell& cell) { - ModuleBase::TITLE("ESolver_FP", "before_runner"); + 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 00d067bb809..cf7cd499c88 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 before_runner(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 9391c747dff..423583c8fc1 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -27,22 +27,32 @@ 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); ///---------------------------------------------------------- @@ -75,15 +85,15 @@ ESolver_KS::~ESolver_KS() } //------------------------------------------------------------------------------ -//! the 3rd function of ESolver_KS: before_runner +//! the 3rd function of ESolver_KS: before_all_runners //! mohan add 2024-05-11 //------------------------------------------------------------------------------ template -void ESolver_KS::before_runner(Input& inp, UnitCell& ucell) +void ESolver_KS::before_all_runners(Input& inp, UnitCell& ucell) { - ModuleBase::TITLE("ESolver_KS", "before_runner"); + ModuleBase::TITLE("ESolver_KS", "before_all_runners"); - ESolver_FP::before_runner(inp,ucell); + ESolver_FP::before_all_runners(inp,ucell); //------------------Charge Mixing------------------ p_chgmix->set_mixing(GlobalV::MIXING_MODE, @@ -155,7 +165,7 @@ void ESolver_KS::before_runner(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); @@ -172,7 +182,7 @@ void ESolver_KS::before_runner(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) @@ -185,6 +195,7 @@ void ESolver_KS::before_runner(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) { @@ -200,9 +211,6 @@ void ESolver_KS::before_runner(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); @@ -215,8 +223,9 @@ void ESolver_KS::before_runner(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); } @@ -255,7 +264,7 @@ void ESolver_KS::before_runner(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); @@ -365,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; @@ -380,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] diff --git a/source/module_esolver/esolver_ks.h b/source/module_esolver/esolver_ks.h index bec25dc9bb9..cdbd392ccbd 100644 --- a/source/module_esolver/esolver_ks.h +++ b/source/module_esolver/esolver_ks.h @@ -43,7 +43,7 @@ class ESolver_KS : public ESolver_FP int out_freq_elec;// frequency for output - virtual void before_runner(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 diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index c760b66097b..eccf7642c19 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -105,10 +105,10 @@ ESolver_KS_LCAO::~ESolver_KS_LCAO() //! 14) set occupations? //------------------------------------------------------------------------------ template -void ESolver_KS_LCAO::before_runner(Input& inp, UnitCell& ucell) +void ESolver_KS_LCAO::before_all_runners(Input& inp, UnitCell& ucell) { - ModuleBase::TITLE("ESolver_KS_LCAO", "before_runner"); - ModuleBase::timer::tick("ESolver_KS_LCAO", "before_runner"); + 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") @@ -131,8 +131,8 @@ void ESolver_KS_LCAO::before_runner(Input& inp, UnitCell& ucell) } else { - // 1) else, call before_runner() in ESolver_KS - ESolver_KS::before_runner(inp, ucell); + // 1) else, call before_all_runners() in ESolver_KS + ESolver_KS::before_all_runners(inp, ucell); } // end ifnot get_S @@ -270,7 +270,7 @@ void ESolver_KS_LCAO::before_runner(Input& inp, UnitCell& ucell) this->pelec->fixed_weights(GlobalV::ocp_kb); } - ModuleBase::timer::tick("ESolver_KS_LCAO", "before_runner"); + ModuleBase::timer::tick("ESolver_KS_LCAO", "before_all_runners"); return; } diff --git a/source/module_esolver/esolver_ks_lcao.h b/source/module_esolver/esolver_ks_lcao.h index 1d27cfae6e2..ac903cbdc1f 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 before_runner(Input& inp, UnitCell& cell) override; + void before_all_runners(Input& inp, UnitCell& cell) override; void init_after_vc(Input& inp, UnitCell& cell) override; diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp index 7cf67f7b691..71bbefcb702 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp @@ -62,9 +62,9 @@ ESolver_KS_LCAO_TDDFT::~ESolver_KS_LCAO_TDDFT() } } -void ESolver_KS_LCAO_TDDFT::before_runner(Input& inp, UnitCell& ucell) +void ESolver_KS_LCAO_TDDFT::before_all_runners(Input& inp, UnitCell& ucell) { - ESolver_KS::before_runner(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 1f17460e95b..e10a9199ba6 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 before_runner(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 60e89744592..6d093545fb7 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::before_runner(Input& inp, UnitCell& ucell) +void ESolver_KS_PW::before_all_runners(Input& inp, UnitCell& ucell) { - ESolver_KS::before_runner(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, diff --git a/source/module_esolver/esolver_ks_pw.h b/source/module_esolver/esolver_ks_pw.h index e9fba701746..db59ef6de06 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 before_runner(Input& inp, UnitCell& cell) override; + void before_all_runners(Input& inp, UnitCell& cell) override; void init_after_vc(Input& inp, UnitCell& cell) override; diff --git a/source/module_esolver/esolver_lj.cpp b/source/module_esolver/esolver_lj.cpp index a059a8205b6..58fcc51ca87 100644 --- a/source/module_esolver/esolver_lj.cpp +++ b/source/module_esolver/esolver_lj.cpp @@ -7,7 +7,7 @@ namespace ModuleESolver { - void ESolver_LJ::before_runner(Input& inp, UnitCell& ucell) + void ESolver_LJ::before_all_runners(Input& inp, UnitCell& ucell) { ucell_ = &ucell; lj_potential = 0; diff --git a/source/module_esolver/esolver_lj.h b/source/module_esolver/esolver_lj.h index 19b858e0fa5..fcac58baf3b 100644 --- a/source/module_esolver/esolver_lj.h +++ b/source/module_esolver/esolver_lj.h @@ -14,7 +14,7 @@ namespace ModuleESolver classname = "ESolver_LJ"; } - void before_runner(Input& inp, UnitCell& cell) override; + void before_all_runners(Input& inp, UnitCell& cell) override; void runner(const int istep, UnitCell& cell) override; diff --git a/source/module_esolver/esolver_of.cpp b/source/module_esolver/esolver_of.cpp index 205e9c66b4d..02411dc7feb 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::before_runner(Input& inp, UnitCell& ucell) +void ESolver_OF::before_all_runners(Input& inp, UnitCell& ucell) { - ESolver_FP::before_runner(inp, ucell); + ESolver_FP::before_all_runners(inp, ucell); // save necessary parameters this->of_kinetic_ = inp.of_kinetic; diff --git a/source/module_esolver/esolver_of.h b/source/module_esolver/esolver_of.h index b64f6cd1965..31435ce9f84 100644 --- a/source/module_esolver/esolver_of.h +++ b/source/module_esolver/esolver_of.h @@ -20,7 +20,7 @@ class ESolver_OF : public ESolver_FP ESolver_OF(); ~ESolver_OF(); - virtual void before_runner(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; diff --git a/source/module_esolver/esolver_sdft_pw.cpp b/source/module_esolver/esolver_sdft_pw.cpp index b0a30969c8f..99c9dc5a4f8 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::before_runner(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::before_runner(inp, ucell); + ESolver_KS::before_all_runners(inp, ucell); this->pelec = new elecstate::ElecStatePW_SDFT(pw_wfc, &(chr), diff --git a/source/module_esolver/esolver_sdft_pw.h b/source/module_esolver/esolver_sdft_pw.h index 2be9228bd85..c5692a7303d 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 before_runner(Input& inp, UnitCell& cell) override; + void before_all_runners(Input& inp, UnitCell& cell) override; double cal_energy() override; From 4267571da2ce4c43002c0e67103779b12de86d4b Mon Sep 17 00:00:00 2001 From: mohanchen Date: Mon, 13 May 2024 09:03:35 +0800 Subject: [PATCH 8/8] fix tests in md --- source/module_esolver/test/esolver_dp_test.cpp | 8 ++++---- 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 +- 7 files changed, 10 insertions(+), 10 deletions(-) diff --git a/source/module_esolver/test/esolver_dp_test.cpp b/source/module_esolver/test/esolver_dp_test.cpp index 8f9fe878a2a..9510d9e1541 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/test/fire_test.cpp b/source/module_md/test/fire_test.cpp index 96a8e5612ef..4570b6bc715 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 8fc152c6519..02c6ad8db60 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 89a062a1b02..51a691ff5b3 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 a5d3ac97e46..9a98214688a 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 49ec2310816..7926e1ade7a 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 8852743edc2..73ad6c40b30 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);