diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 1aea12100a..91e4f34bea 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -634,6 +634,7 @@ OBJS_LCAO=evolve_elec.o\ spar_hsr.o\ spar_st.o\ spar_u.o\ + LCAO_set.o\ LCAO_set_fs.o\ LCAO_set_st.o\ LCAO_nl_mu.o\ diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index d92bf2ab3b..62f00adcb6 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -17,10 +17,16 @@ #include "source_io/ctrl_runner_lcao.h" // use ctrl_runner_lcao() #include "source_io/ctrl_iter_lcao.h" // use ctrl_iter_lcao() #include "source_io/ctrl_scf_lcao.h" // use ctrl_scf_lcao() -#include "source_psi/setup_psi.h" // mohan add 20251019 -#include "source_io/read_wfc_nao.h" #include "source_io/print_info.h" #include "source_lcao/rho_tau_lcao.h" // mohan add 20251024 +#include "source_lcao/LCAO_set.h" // mohan add 20251111 + + +// tmp +#include "source_psi/setup_psi.h" // use Setup_Psi +#include "source_io/read_wfc_nao.h" // use read_wfc_nao +#include "source_estate/elecstate_tools.h" // use fixed_weights + namespace ModuleESolver { @@ -74,62 +80,17 @@ void ESolver_KS_LCAO::before_all_runners(UnitCell& ucell, const Input_pa return; } - // 5) init electronic wave function psi - Setup_Psi::allocate_psi(this->psi, this->kv, this->pv, inp); - - //! read psi from file - if (inp.init_wfc == "file" && inp.esolver_type != "tddft") - { - if (!ModuleIO::read_wfc_nao(PARAM.globalv.global_readin_dir, - this->pv, *this->psi, this->pelec->ekb, this->pelec->wg, this->kv.ik2iktot, - this->kv.get_nkstot(), inp.nspin)) - { - ModuleBase::WARNING_QUIT("ESolver_KS_LCAO", "read electronic wave functions failed"); - } - } - + LCAO_domain::set_psi_occ_dm_chg(this->kv, this->psi, this->pv, this->pelec, + this->dmat, this->chr, inp); - // 7) init DMK, but DMR is constructed in before_scf() - this->dmat.allocate_dm(&this->kv, &this->pv, inp.nspin); + LCAO_domain::set_pot(ucell, this->kv, this->sf, *this->pw_rho, *this->pw_rhod, + this->pelec, this->orb_, this->pv, this->locpp, this->dftu, + this->solvent, this->exx_nao, this->deepks, inp); - // 8) init exact exchange calculations - this->exx_nao.before_runner(ucell, this->kv, this->orb_, this->pv, inp); - - // 9) initialize DFT+U - if (inp.dft_plus_u) - { - this->dftu.init(ucell, &this->pv, this->kv.get_nks(), &orb_); - } - - // 10) init local pseudopotentials - this->locpp.init_vloc(ucell, this->pw_rho); - ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL"); - - // 11) init charge density - this->chr.allocate(inp.nspin); - - // 12) init potentials - if (this->pelec->pot == nullptr) - { - this->pelec->pot = new elecstate::Potential(this->pw_rhod, this->pw_rho, - &ucell, &(this->locpp.vloc), &(this->sf), &(this->solvent), - &(this->pelec->f_en.etxc), &(this->pelec->f_en.vtxc)); - } - - // 13) init deepks - this->deepks.before_runner(ucell, this->kv.get_nks(), this->orb_, this->pv, inp); - - // 14) set occupations, tddft does not need to set occupations in the first scf - if (inp.ocp && inp.esolver_type != "tddft") - { - elecstate::fixed_weights(inp.ocp_kb, inp.nbands, inp.nelec, - this->pelec->klist, this->pelec->wg, this->pelec->skip_weights); - } - - // 15) if kpar is not divisible by nks, print a warning + //! if kpar is not divisible by nks, print a warning ModuleIO::print_kpar(this->kv.get_nks(), PARAM.globalv.kpar_lcao); - // 16) init rdmft, added by jghan + //! init rdmft, added by jghan if (inp.rdmft == true) { rdmft_solver.init(this->pv, ucell, @@ -310,8 +271,6 @@ void ESolver_KS_LCAO::after_all_runners(UnitCell& ucell) ESolver_KS::after_all_runners(ucell); - const int nspin0 = (PARAM.inp.nspin == 2) ? 2 : 1; - auto* hamilt_lcao = dynamic_cast*>(this->p_hamilt); if(!hamilt_lcao) { diff --git a/source/source_lcao/CMakeLists.txt b/source/source_lcao/CMakeLists.txt index 5199625c41..7b090ab275 100644 --- a/source/source_lcao/CMakeLists.txt +++ b/source/source_lcao/CMakeLists.txt @@ -33,6 +33,7 @@ if(ENABLE_LCAO) spar_hsr.cpp spar_st.cpp spar_u.cpp + LCAO_set.cpp LCAO_set_fs.cpp LCAO_set_st.cpp LCAO_nl_mu.cpp diff --git a/source/source_lcao/LCAO_set.cpp b/source/source_lcao/LCAO_set.cpp new file mode 100644 index 0000000000..2c87852fde --- /dev/null +++ b/source/source_lcao/LCAO_set.cpp @@ -0,0 +1,147 @@ +#include "source_lcao/LCAO_set.h" +#include "source_io/module_parameter/parameter.h" +#include "source_psi/setup_psi.h" // use Setup_Psi +#include "source_io/read_wfc_nao.h" // use read_wfc_nao +#include "source_estate/elecstate_tools.h" // use fixed_weights + +template +void LCAO_domain::set_psi_occ_dm_chg( + const K_Vectors &kv, // k-points + psi::Psi* &psi, // coefficients of NAO basis + const Parallel_Orbitals &pv, // parallel scheme of NAO basis + elecstate::ElecState* pelec, // eigen values and weights + LCAO_domain::Setup_DM &dmat, // density matrix + Charge &chr, // charge density + const Input_para &inp) // input parameters +{ + + //! 1) init electronic wave function psi + Setup_Psi::allocate_psi(psi, kv, pv, inp); + + //! 2) read psi from file + if (inp.init_wfc == "file" && inp.esolver_type != "tddft") + { + if (!ModuleIO::read_wfc_nao(PARAM.globalv.global_readin_dir, + pv, *psi, pelec->ekb, pelec->wg, kv.ik2iktot, + kv.get_nkstot(), inp.nspin)) + { + ModuleBase::WARNING_QUIT("set_psi_occ_dm_chg", "read electronic wave functions failed"); + } + } + + //! 3) set occupations, tddft does not need to set occupations in the first scf + if (inp.ocp && inp.esolver_type != "tddft") + { + elecstate::fixed_weights(inp.ocp_kb, inp.nbands, inp.nelec, + &kv, pelec->wg, pelec->skip_weights); + } + + //! 4) init DMK, but DMR is constructed in before_scf() + dmat.allocate_dm(&kv, &pv, inp.nspin); + + //! 5) init charge density + chr.allocate(inp.nspin); + + ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "CHARGE"); + + return; +} + + +template +void LCAO_domain::set_pot( + UnitCell &ucell, // not const because of dftu + K_Vectors &kv, // not const due to exx + Structure_Factor& sf, // will be modified in potential + const ModulePW::PW_Basis &pw_rho, + const ModulePW::PW_Basis &pw_rhod, + elecstate::ElecState* pelec, + const LCAO_Orbitals& orb, + Parallel_Orbitals &pv, // not const due to deepks + pseudopot_cell_vl &locpp, + Plus_U &dftu, + surchem& solvent, + Exx_NAO &exx_nao, + Setup_DeePKS &deepks, + const Input_para &inp) +{ + //! 1) init local pseudopotentials + locpp.init_vloc(ucell, &pw_rho); + + //! 2) init potentials + if (pelec->pot == nullptr) + { + // where is the pot deleted? + pelec->pot = new elecstate::Potential(&pw_rhod, &pw_rho, + &ucell, &locpp.vloc, &sf, &solvent, + &(pelec->f_en.etxc), &(pelec->f_en.vtxc)); + } + + //! 3) initialize DFT+U + if (inp.dft_plus_u) + { + dftu.init(ucell, &pv, kv.get_nks(), &orb); + } + + //! 4) init exact exchange calculations + exx_nao.before_runner(ucell, kv, orb, pv, inp); + + //! 5) init deepks + deepks.before_runner(ucell, kv.get_nks(), orb, pv, inp); + + ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "POTENTIALS"); + + return; +} + + + +template void LCAO_domain::set_psi_occ_dm_chg( + const K_Vectors &kv, // k-points + psi::Psi* &psi, // coefficients of NAO basis + const Parallel_Orbitals &pv, // parallel scheme of NAO basis + elecstate::ElecState* pelec, // eigen values and weights + LCAO_domain::Setup_DM &dmat, // density matrix + Charge &chr, // charge density + const Input_para &inp); + +template void LCAO_domain::set_psi_occ_dm_chg>( + const K_Vectors &kv, // k-points + psi::Psi>* &psi, // coefficients of NAO basis + const Parallel_Orbitals &pv, // parallel scheme of NAO basis + elecstate::ElecState* pelec, // eigen values and weights + LCAO_domain::Setup_DM> &dmat, // density matrix + Charge &chr, // charge density + const Input_para &inp); + +template void LCAO_domain::set_pot( + UnitCell &ucell, + K_Vectors &kv, + Structure_Factor& sf, + const ModulePW::PW_Basis &pw_rho, + const ModulePW::PW_Basis &pw_rhod, + elecstate::ElecState* pelec, + const LCAO_Orbitals& orb, + Parallel_Orbitals &pv, + pseudopot_cell_vl &locpp, + Plus_U &dftu, + surchem& solvent, + Exx_NAO &exx_nao, + Setup_DeePKS &deepks, + const Input_para &inp); + +template void LCAO_domain::set_pot>( + UnitCell &ucell, + K_Vectors &kv, + Structure_Factor& sf, + const ModulePW::PW_Basis &pw_rho, + const ModulePW::PW_Basis &pw_rhod, + elecstate::ElecState* pelec, + const LCAO_Orbitals& orb, + Parallel_Orbitals &pv, + pseudopot_cell_vl &locpp, + Plus_U &dftu, + surchem& solvent, + Exx_NAO> &exx_nao, + Setup_DeePKS> &deepks, + const Input_para &inp); diff --git a/source/source_lcao/LCAO_set.h b/source/source_lcao/LCAO_set.h new file mode 100644 index 0000000000..f5df2df054 --- /dev/null +++ b/source/source_lcao/LCAO_set.h @@ -0,0 +1,58 @@ +#ifndef LCAO_SET_H +#define LCAO_SET_H + +#include "source_cell/unitcell.h" +#include "source_cell/klist.h" +#include "source_psi/psi.h" +#include "source_estate/elecstate.h" +#include "source_lcao/setup_dm.h" +#include "source_pw/module_pwdft/structure_factor.h" +#include "source_basis/module_pw/pw_basis.h" +#include "source_hamilt/module_surchem/surchem.h" +#include "source_pw/module_pwdft/VL_in_pw.h" +#include "source_lcao/module_deepks/LCAO_deepks.h" +#include "source_lcao/module_dftu/dftu.h" +#include "source_lcao/setup_exx.h" +#include "source_lcao/setup_deepks.h" + +namespace LCAO_domain +{ + +/** + * @brief set up wave functions, occupation numbers, + * density matrix and charge density + */ +template +void set_psi_occ_dm_chg( + const K_Vectors &kv, // k-points + psi::Psi* &psi, // coefficients of NAO basis + const Parallel_Orbitals &pv, // parallel scheme of NAO basis + elecstate::ElecState* pelec, // eigen values and weights + LCAO_domain::Setup_DM &dmat, // density matrix + Charge &chr, // charge density + const Input_para& inp); // input parameters + +/** + * @brief set up potentials, including local pseudopotentials, + * +U potential, solvent potential, exx potential and deepks potential + */ +template +void set_pot( + UnitCell &ucell, + K_Vectors &kv, + Structure_Factor& sf, + const ModulePW::PW_Basis &pw_rho, + const ModulePW::PW_Basis &pw_rhod, + elecstate::ElecState* pelec, + const LCAO_Orbitals& orb, + Parallel_Orbitals &pv, + pseudopot_cell_vl &locpp, + Plus_U &dftu, + surchem& solvent, + Exx_NAO &exx_nao, + Setup_DeePKS &deepks, + const Input_para &inp); + +} // end namespace + +#endif