diff --git a/source/source_esolver/esolver_gets.h b/source/source_esolver/esolver_gets.h index 564fd55035..7a7fb1d34b 100644 --- a/source/source_esolver/esolver_gets.h +++ b/source/source_esolver/esolver_gets.h @@ -10,7 +10,7 @@ namespace ModuleESolver { -class ESolver_GetS : public ESolver_KS> +class ESolver_GetS : public ESolver_KS { public: ESolver_GetS(); diff --git a/source/source_esolver/esolver_ks.cpp b/source/source_esolver/esolver_ks.cpp index fc99b8a572..cc94510a66 100644 --- a/source/source_esolver/esolver_ks.cpp +++ b/source/source_esolver/esolver_ks.cpp @@ -15,35 +15,29 @@ #include "source_io/module_output/output_log.h" // use write_head #include "source_estate/elecstate_print.h" // print_etot #include "source_io/module_output/print_info.h" // print_parameters -#include "source_psi/setup_psi.h" // mohan add 20251009 #include "source_lcao/module_dftu/dftu.h" // mohan add 2025-11-07 namespace ModuleESolver { -template -ESolver_KS::ESolver_KS(){} +ESolver_KS::ESolver_KS() {} -template -ESolver_KS::~ESolver_KS() +ESolver_KS::~ESolver_KS() { - //**************************************************** - // do not add any codes in this deconstructor funcion - //**************************************************** - Setup_Psi::deallocate_psi(this->psi); - + //**************************************************** + // do not add any codes in this deconstructor funcion + //**************************************************** delete this->p_hamilt; delete this->p_chgmix; this->ppcell.release_memory(); - + // mohan add 2025-10-18, should be put int clean() function pw::teardown_pwwfc(this->pw_wfc); } -template -void ESolver_KS::before_all_runners(UnitCell& ucell, const Input_para& inp) +void ESolver_KS::before_all_runners(UnitCell& ucell, const Input_para& inp) { ModuleBase::TITLE("ESolver_KS", "before_all_runners"); @@ -81,12 +75,10 @@ void ESolver_KS::before_all_runners(UnitCell& ucell, const Input_para } -template -void ESolver_KS::hamilt2rho_single(UnitCell& ucell, const int istep, const int iter, const double ethr) +void ESolver_KS::hamilt2rho_single(UnitCell& ucell, const int istep, const int iter, const double ethr) {} -template -void ESolver_KS::hamilt2rho(UnitCell& ucell, const int istep, const int iter, const double ethr) +void ESolver_KS::hamilt2rho(UnitCell& ucell, const int istep, const int iter, const double ethr) { // 1) use Hamiltonian to obtain charge density this->hamilt2rho_single(ucell, istep, iter, diag_ethr); @@ -126,8 +118,7 @@ void ESolver_KS::hamilt2rho(UnitCell& ucell, const int istep, const i } } -template -void ESolver_KS::runner(UnitCell& ucell, const int istep) +void ESolver_KS::runner(UnitCell& ucell, const int istep) { ModuleBase::TITLE("ESolver_KS", "runner"); ModuleBase::timer::tick(this->classname, "runner"); @@ -142,14 +133,14 @@ void ESolver_KS::runner(UnitCell& ucell, const int istep) this->diag_ethr = PARAM.inp.pw_diag_thr; this->scf_nmax_flag = false; // mohan add 2025-09-21 for (int iter = 1; iter <= this->maxniter; ++iter) - { - if(iter == this->maxniter) - { - this->scf_nmax_flag=true; - } + { + if(iter == this->maxniter) + { + this->scf_nmax_flag=true; + } - // 3) initialization of SCF iterations - this->iter_init(ucell, istep, iter); + // 3) initialization of SCF iterations + this->iter_init(ucell, istep, iter); // 4) use Hamiltonian to obtain charge density this->hamilt2rho(ucell, istep, iter, diag_ethr); @@ -169,22 +160,20 @@ void ESolver_KS::runner(UnitCell& ucell, const int istep) } } // end scf iterations - // 7) after scf + // 7) after scf this->after_scf(ucell, istep, conv_esolver); ModuleBase::timer::tick(this->classname, "runner"); return; }; -template -void ESolver_KS::before_scf(UnitCell& ucell, const int istep) +void ESolver_KS::before_scf(UnitCell& ucell, const int istep) { ModuleBase::TITLE("ESolver_KS", "before_scf"); ESolver_FP::before_scf(ucell, istep); } -template -void ESolver_KS::iter_init(UnitCell& ucell, const int istep, const int iter) +void ESolver_KS::iter_init(UnitCell& ucell, const int istep, const int iter) { if(PARAM.inp.esolver_type != "tddft") { @@ -210,8 +199,7 @@ void ESolver_KS::iter_init(UnitCell& ucell, const int istep, const in this->chr.save_rho_before_sum_band(); } -template -void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& iter, bool &conv_esolver) +void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& iter, bool &conv_esolver) { // 1.1) print out band gap @@ -227,25 +215,25 @@ void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& i // 1.2) print out eigenvalues and occupations if (PARAM.inp.out_band[0]) { - if (iter % PARAM.inp.out_freq_elec == 0 || iter == PARAM.inp.scf_nmax || conv_esolver) - { - ModuleIO::write_eig_iter(this->pelec->ekb,this->pelec->wg,*this->pelec->klist); - } + if (iter % PARAM.inp.out_freq_elec == 0 || iter == PARAM.inp.scf_nmax || conv_esolver) + { + ModuleIO::write_eig_iter(this->pelec->ekb,this->pelec->wg,*this->pelec->klist); + } } // 2.1) compute magnetization, only for spin==2 ucell.magnet.compute_mag(ucell.omega, this->chr.nrxx, this->chr.nxyz, this->chr.rho, this->pelec->nelec_spin.data()); - // 2.2) charge mixing + // 2.2) charge mixing // SCF will continue if U is not converged for uramping calculation - bool converged_u = true; - // to avoid unnecessary dependence on dft+u, refactor is needed + bool converged_u = true; + // to avoid unnecessary dependence on dft+u, refactor is needed #ifdef __LCAO - if (PARAM.inp.dft_plus_u) - { - converged_u = this->dftu.u_converged(); - } + if (PARAM.inp.dft_plus_u) + { + converged_u = this->dftu.u_converged(); + } #endif module_charge::chgmixing_ks(iter, ucell, this->pelec, this->chr, this->p_chgmix, @@ -296,8 +284,7 @@ void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& i } //! Something to do after SCF iterations when SCF is converged or comes to the max iter step. -template -void ESolver_KS::after_scf(UnitCell& ucell, const int istep, const bool conv_esolver) +void ESolver_KS::after_scf(UnitCell& ucell, const int istep, const bool conv_esolver) { ModuleBase::TITLE("ESolver_KS", "after_scf"); @@ -321,29 +308,10 @@ void ESolver_KS::after_scf(UnitCell& ucell, const int istep, const bo } -template -void ESolver_KS::after_all_runners(UnitCell& ucell) +void ESolver_KS::after_all_runners(UnitCell& ucell) { // 1) write Etot information ESolver_FP::after_all_runners(ucell); } -//------------------------------------------------------------------------------ -//! the 16th-20th functions of ESolver_KS -//! mohan add 2024-05-12 -//------------------------------------------------------------------------------ -//! This is for mixed-precision pw/LCAO basis sets. -template class ESolver_KS, base_device::DEVICE_CPU>; -template class ESolver_KS, base_device::DEVICE_CPU>; - -//! This is for GPU codes. -#if ((defined __CUDA) || (defined __ROCM)) -template class ESolver_KS, base_device::DEVICE_GPU>; -template class ESolver_KS, base_device::DEVICE_GPU>; -#endif - -//! This is for LCAO basis set. -#ifdef __LCAO -template class ESolver_KS; -#endif } // namespace ModuleESolver diff --git a/source/source_esolver/esolver_ks.h b/source/source_esolver/esolver_ks.h index 1913aa4101..b6affc7b0c 100644 --- a/source/source_esolver/esolver_ks.h +++ b/source/source_esolver/esolver_ks.h @@ -5,7 +5,6 @@ #include "source_basis/module_pw/pw_basis_k.h" // use plane wave #include "source_cell/klist.h" // use k-points in Brillouin zone #include "source_estate/module_charge/charge_mixing.h" // use charge mixing -#include "source_psi/psi.h" // use electronic wave functions #include "source_hamilt/hamilt.h" // use Hamiltonian #include "source_hamilt/hamilt_base.h" // use Hamiltonian base class #include "source_lcao/module_dftu/dftu.h" // mohan add 20251107 @@ -13,7 +12,6 @@ namespace ModuleESolver { -template class ESolver_KS : public ESolver_FP { public: @@ -60,9 +58,6 @@ class ESolver_KS : public ESolver_FP //! nonlocal pseudopotentials pseudopot_cell_vnl ppcell; - //! Electronic wavefunctions - psi::Psi* psi = nullptr; - //! DFT+U method, mohan add 2025-11-07 Plus_U dftu; diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index bad1cec7b6..0558942e91 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -22,6 +22,7 @@ #include "source_io/module_output/print_info.h" #include "source_lcao/rho_tau_lcao.h" // mohan add 20251024 #include "source_lcao/LCAO_set.h" // mohan add 20251111 +#include "source_psi/setup_psi.h" // use Setup_Psi for deallocate_psi namespace ModuleESolver { @@ -40,6 +41,7 @@ ESolver_KS_LCAO::~ESolver_KS_LCAO() //**************************************************** // do not add any codes in this deconstructor funcion //**************************************************** + Setup_Psi::deallocate_psi(this->psi); } template @@ -49,7 +51,7 @@ void ESolver_KS_LCAO::before_all_runners(UnitCell& ucell, const Input_pa ModuleBase::timer::tick("ESolver_KS_LCAO", "before_all_runners"); // 1) before_all_runners in ESolver_KS - ESolver_KS::before_all_runners(ucell, inp); + ESolver_KS::before_all_runners(ucell, inp); // 2) autoset nbands in ElecState before init_basis (for Psi 2d division) if (this->pelec == nullptr) @@ -105,7 +107,7 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) ModuleBase::timer::tick("ESolver_KS_LCAO", "before_scf"); //! 1) call before_scf() of ESolver_KS. - ESolver_KS::before_scf(ucell, istep); + ESolver_KS::before_scf(ucell, istep); //! 2) find search radius double search_radius = atom_arrange::set_sr_NL(GlobalV::ofs_running, @@ -269,7 +271,7 @@ void ESolver_KS_LCAO::after_all_runners(UnitCell& ucell) ModuleBase::TITLE("ESolver_KS_LCAO", "after_all_runners"); ModuleBase::timer::tick("ESolver_KS_LCAO", "after_all_runners"); - ESolver_KS::after_all_runners(ucell); + ESolver_KS::after_all_runners(ucell); auto* hamilt_lcao = dynamic_cast*>(this->p_hamilt); if(!hamilt_lcao) @@ -301,7 +303,7 @@ void ESolver_KS_LCAO::iter_init(UnitCell& ucell, const int istep, const ModuleBase::TITLE("ESolver_KS_LCAO", "iter_init"); // call iter_init() of ESolver_KS - ESolver_KS::iter_init(ucell, istep, iter); + ESolver_KS::iter_init(ucell, istep, iter); module_charge::chgmixing_ks_lcao(iter, this->p_chgmix, this->dftu, this->dmat.dm->get_DMR_pointer(1)->get_nnr(), PARAM.inp); @@ -436,7 +438,7 @@ void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& // eig and occ are printed, magnetization is calculated, // charge mixing is performed, potential is updated, // HF and kS energies are computed, meta-GGA, Jason and restart - ESolver_KS::iter_finish(ucell, istep, iter, conv_esolver); + ESolver_KS::iter_finish(ucell, istep, iter, conv_esolver); // mix density matrix if mixing_restart + mixing_dmr + not first // mixing_restart at every iter except the last iter @@ -474,7 +476,7 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep, const } //! 1) call after_scf() of ESolver_KS - ESolver_KS::after_scf(ucell, istep, conv_esolver); + ESolver_KS::after_scf(ucell, istep, conv_esolver); //! 2) output of lcao every few ionic steps ModuleIO::ctrl_scf_lcao(ucell, diff --git a/source/source_esolver/esolver_ks_lcao.h b/source/source_esolver/esolver_ks_lcao.h index 4191306788..143f7089ba 100644 --- a/source/source_esolver/esolver_ks_lcao.h +++ b/source/source_esolver/esolver_ks_lcao.h @@ -28,7 +28,7 @@ namespace ModuleESolver { template -class ESolver_KS_LCAO : public ESolver_KS +class ESolver_KS_LCAO : public ESolver_KS { public: ESolver_KS_LCAO(); @@ -57,6 +57,9 @@ class ESolver_KS_LCAO : public ESolver_KS virtual void others(UnitCell& ucell, const int istep) override; + //! Electronic wave functions (moved from base class) + psi::Psi* psi = nullptr; + //! Store information about Adjacent Atoms Record_adj RA; diff --git a/source/source_esolver/esolver_ks_lcao_tddft.cpp b/source/source_esolver/esolver_ks_lcao_tddft.cpp index 05dc8c9233..2e463acdcd 100644 --- a/source/source_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/source_esolver/esolver_ks_lcao_tddft.cpp @@ -40,7 +40,11 @@ ESolver_KS_LCAO_TDDFT::~ESolver_KS_LCAO_TDDFT() //************************************************* // Do not add any code in this destructor function //************************************************* - delete psi_laststep; + if (psi_laststep != nullptr) + { + delete psi_laststep; + psi_laststep = nullptr; + } if (td_p != nullptr) { diff --git a/source/source_esolver/esolver_ks_pw.cpp b/source/source_esolver/esolver_ks_pw.cpp index b2976733bf..74507e09c7 100644 --- a/source/source_esolver/esolver_ks_pw.cpp +++ b/source/source_esolver/esolver_ks_pw.cpp @@ -84,7 +84,7 @@ void ESolver_KS_PW::deallocate_hamilt() template void ESolver_KS_PW::before_all_runners(UnitCell& ucell, const Input_para& inp) { - ESolver_KS::before_all_runners(ucell, inp); + ESolver_KS::before_all_runners(ucell, inp); //! setup and allocation for pelec, potentials, etc. elecstate::setup_estate_pw(ucell, this->kv, this->sf, this->pelec, this->chr, @@ -105,7 +105,7 @@ void ESolver_KS_PW::before_scf(UnitCell& ucell, const int istep) ModuleBase::TITLE("ESolver_KS_PW", "before_scf"); ModuleBase::timer::tick("ESolver_KS_PW", "before_scf"); - ESolver_KS::before_scf(ucell, istep); + ESolver_KS::before_scf(ucell, istep); //! Init variables (once the cell has changed) pw::update_cell_pw(ucell, this->ppcell, this->kv, this->pw_wfc, PARAM.inp); @@ -142,7 +142,7 @@ void ESolver_KS_PW::before_scf(UnitCell& ucell, const int istep) template void ESolver_KS_PW::iter_init(UnitCell& ucell, const int istep, const int iter) { - ESolver_KS::iter_init(ucell, istep, iter); + ESolver_KS::iter_init(ucell, istep, iter); module_charge::chgmixing_ks_pw(iter, this->p_chgmix, this->dftu, PARAM.inp); @@ -212,7 +212,7 @@ void ESolver_KS_PW::iter_finish(UnitCell& ucell, const int istep, int this->pelec->f_en.deband = this->pelec->cal_delta_eband(ucell); // Call iter_finish() of ESolver_KS - ESolver_KS::iter_finish(ucell, istep, iter, conv_esolver); + ESolver_KS::iter_finish(ucell, istep, iter, conv_esolver); // D in USPP needs vloc, thus needs update when veff updated // calculate the effective coefficient matrix for non-local @@ -240,16 +240,13 @@ void ESolver_KS_PW::after_scf(UnitCell& ucell, const int istep, const ModuleBase::TITLE("ESolver_KS_PW", "after_scf"); ModuleBase::timer::tick("ESolver_KS_PW", "after_scf"); - // Since ESolver_KS::psi is hidden by ESolver_KS_PW::psi, - // we need to copy the data from ESolver_KS::psi to ESolver_KS_PW::psi. - // sunliang 2025-04-10 + // Calculate kinetic energy density tau for ELF if needed if (PARAM.inp.out_elf[0] > 0) { - this->ESolver_KS::psi = new psi::Psi(this->stp.psi_cpu[0]); - this->pelec->cal_tau(*(this->psi)); + this->pelec->cal_tau(*(this->stp.psi_cpu)); } - ESolver_KS::after_scf(ucell, istep, conv_esolver); + ESolver_KS::after_scf(ucell, istep, conv_esolver); // Output quantities ModuleIO::ctrl_scf_pw(istep, ucell, this->pelec, this->chr, this->kv, this->pw_wfc, @@ -303,7 +300,7 @@ void ESolver_KS_PW::cal_stress(UnitCell& ucell, ModuleBase::matrix& s template void ESolver_KS_PW::after_all_runners(UnitCell& ucell) { - ESolver_KS::after_all_runners(ucell); + ESolver_KS::after_all_runners(ucell); ModuleIO::ctrl_runner_pw(ucell, this->pelec, this->pw_wfc, this->pw_rho, this->pw_rhod, this->chr, this->kv, this->stp, diff --git a/source/source_esolver/esolver_ks_pw.h b/source/source_esolver/esolver_ks_pw.h index 01e1027d79..6a6be52b73 100644 --- a/source/source_esolver/esolver_ks_pw.h +++ b/source/source_esolver/esolver_ks_pw.h @@ -13,7 +13,7 @@ namespace ModuleESolver { template -class ESolver_KS_PW : public ESolver_KS +class ESolver_KS_PW : public ESolver_KS { private: using Real = typename GetTypeReal::type; diff --git a/source/source_esolver/esolver_sdft_pw.cpp b/source/source_esolver/esolver_sdft_pw.cpp index 597825cd6d..f7f9a29983 100644 --- a/source/source_esolver/esolver_sdft_pw.cpp +++ b/source/source_esolver/esolver_sdft_pw.cpp @@ -119,7 +119,7 @@ template void ESolver_SDFT_PW::iter_finish(UnitCell& ucell, const int istep, int& iter, bool& conv_esolver) { // call iter_finish() of ESolver_KS - ESolver_KS::iter_finish(ucell, istep, iter, conv_esolver); + ESolver_KS::iter_finish(ucell, istep, iter, conv_esolver); } template diff --git a/source/source_esolver/pw_others.cpp b/source/source_esolver/pw_others.cpp index fc42df14bd..49f7465b46 100644 --- a/source/source_esolver/pw_others.cpp +++ b/source/source_esolver/pw_others.cpp @@ -32,7 +32,7 @@ void ESolver_KS_PW::others(UnitCell& ucell, const int istep) { Numerical_Descriptor nc; nc.output_descriptor(ucell, - this->psi[0], + *(this->stp.psi_cpu), PARAM.inp.bessel_descriptor_lmax, PARAM.inp.bessel_descriptor_rcut, PARAM.inp.bessel_descriptor_tolerence,