diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 537ad1a4dd..d43b91ee3d 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -581,7 +581,9 @@ OBJS_IO=input_conv.o\ write_libxc_r.o\ output_log.o\ output_mat_sparse.o\ - ctrl_output_lcao.o\ + ctrl_scf_lcao.o\ + ctrl_runner_lcao.o\ + ctrl_iter_lcao.o\ ctrl_output_fp.o\ ctrl_output_pw.o\ para_json.o\ @@ -704,6 +706,7 @@ OBJS_SRCPW=H_Ewald_pw.o\ VL_in_pw.o\ VNL_in_pw.o\ VNL_grad_pw.o\ + chgmixing.o\ charge.o\ charge_init.o\ charge_mpi.o\ diff --git a/source/source_esolver/esolver_ks.cpp b/source/source_esolver/esolver_ks.cpp index b477c89656..a00abc02e6 100644 --- a/source/source_esolver/esolver_ks.cpp +++ b/source/source_esolver/esolver_ks.cpp @@ -1,7 +1,5 @@ #include "esolver_ks.h" - -// To setup plane wave for electronic wave functions -#include "pw_setup.h" +#include "pw_setup.h" // setup plane wave #include "source_base/timer.h" #include "source_base/global_variable.h" @@ -152,10 +150,6 @@ void ESolver_KS::before_all_runners(UnitCell& ucell, const Input_para this->sf.setup_structure_factor(&ucell, Pgrid, this->pw_rhod); } -//------------------------------------------------------------------------------ -//! the 5th function of ESolver_KS: hamilt2rho_single -//! mohan add 2024-05-11 -//------------------------------------------------------------------------------ template void ESolver_KS::hamilt2rho_single(UnitCell& ucell, const int istep, const int iter, const double ethr) { @@ -248,24 +242,16 @@ void ESolver_KS::runner(UnitCell& ucell, const int istep) this->scf_nmax_flag=true; } - //---------------------------------------------------------------- // 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); - //---------------------------------------------------------------- // 5) finish scf iterations - //---------------------------------------------------------------- this->iter_finish(ucell, istep, iter, conv_esolver); - //---------------------------------------------------------------- // 6) check convergence - //---------------------------------------------------------------- if (conv_esolver || this->oscillate_esolver) { this->niter = iter; @@ -277,9 +263,7 @@ void ESolver_KS::runner(UnitCell& ucell, const int istep) } } // end scf iterations - //---------------------------------------------------------------- // 7) after scf - //---------------------------------------------------------------- this->after_scf(ucell, istep, conv_esolver); ModuleBase::timer::tick(this->classname, "runner"); @@ -598,9 +582,6 @@ void ESolver_KS::after_scf(UnitCell& ucell, const int istep, const bo this->kv); } } - - - } template diff --git a/source/source_esolver/esolver_ks.h b/source/source_esolver/esolver_ks.h index 974351a4a7..d54794dd89 100644 --- a/source/source_esolver/esolver_ks.h +++ b/source/source_esolver/esolver_ks.h @@ -1,21 +1,12 @@ #ifndef ESOLVER_KS_H #define ESOLVER_KS_H -#include -//#include - -// for first-principles esolver -#include "esolver_fp.h" -// for plane wave basis set -#include "source_basis/module_pw/pw_basis_k.h" -// for k-points in Brillouin zone -#include "source_cell/klist.h" -// for charge mixing -#include "source_estate/module_charge/charge_mixing.h" -// for electronic wave functions -#include "source_psi/psi.h" -// for Hamiltonian -#include "source_hamilt/hamilt.h" +#include "esolver_fp.h" // first-principles esolver +#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 namespace ModuleESolver { diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index b7ef11ed55..456aafee7a 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -1,8 +1,6 @@ #include "esolver_ks_lcao.h" -#include "source_io/write_proj_band_lcao.h" // projcted band structure - -#include "source_base/formatter.h" +//#include "source_base/formatter.h" #include "source_base/global_variable.h" #include "source_base/tool_title.h" #include "source_estate/elecstate_tools.h" @@ -10,28 +8,24 @@ #include "source_estate/module_dm/cal_dm_psi.h" #include "source_lcao/module_deltaspin/spin_constrain.h" #include "source_lcao/module_dftu/dftu.h" -#include "source_io/berryphase.h" -#include "source_io/cal_ldos.h" +//#include "source_io/berryphase.h" #include "source_io/cube_io.h" //#include "source_io/io_npz.h" -#include "source_io/output_dmk.h" +//#include "source_io/output_dmk.h" #include "source_io/output_log.h" -#include "source_io/output_mat_sparse.h" -#include "source_io/output_mulliken.h" -#include "source_io/output_sk.h" +//#include "source_io/output_mat_sparse.h" +//#include "source_io/output_mulliken.h" +//#include "source_io/output_sk.h" #include "source_io/read_wfc_nao.h" -#include "source_io/to_qo.h" -#include "source_io/to_wannier90_lcao.h" -#include "source_io/to_wannier90_lcao_in_pw.h" -#include "source_io/write_HS.h" +//#include "source_io/to_qo.h" +//#include "source_io/to_wannier90_lcao.h" +//#include "source_io/to_wannier90_lcao_in_pw.h" +//#include "source_io/write_HS.h" #include "source_io/write_elecstat_pot.h" #include "source_io/module_parameter/parameter.h" // be careful of hpp, there may be multiple definitions of functions, 20250302, mohan #include "source_lcao/hs_matrix_k.hpp" -#include "source_io/write_eband_terms.hpp" -#include "source_io/write_vxc.hpp" -#include "source_io/write_vxc_r.hpp" #include "source_base/global_function.h" #include "source_cell/module_neighbor/sltk_grid_driver.h" @@ -63,10 +57,12 @@ // test RDMFT #include "source_lcao/module_rdmft/rdmft.h" - #include "source_lcao/module_gint/temp_gint/gint_info.h" -#include +#include "source_estate/module_charge/chgmixing.h" // use charge mixing, mohan add 20251006 +#include "source_io/ctrl_runner_lcao.h" // use ctrl_runner_lcao() +#include "source_io/ctrl_iter_lcao.h" // use ctrl_iter_lcao() + namespace ModuleESolver { @@ -112,32 +108,19 @@ void ESolver_KS_LCAO::before_all_runners(UnitCell& ucell, const Input_pa ESolver_KS::before_all_runners(ucell, inp); // 2) init ElecState - // autoset nbands in ElecState before basis_init (for Psi 2d division) + // autoset nbands in ElecState before init_basis (for Psi 2d division) if (this->pelec == nullptr) { // TK stands for double and std::complex? - this->pelec = new elecstate::ElecStateLCAO(&(this->chr), // use which parameter? - &(this->kv), - this->kv.get_nks(), - &(this->GG), - &(this->GK), - this->pw_rho, - this->pw_big); - } - - // 3) init LCAO basis - // reading the localized orbitals/projectors - // construct the interpolation tables. - LCAO_domain::init_basis_lcao(this->pv, - inp.onsite_radius, - inp.lcao_ecut, - inp.lcao_dk, - inp.lcao_dr, - inp.lcao_rmax, - ucell, - two_center_bundle_, - orb_); + this->pelec = new elecstate::ElecStateLCAO(&(this->chr), &(this->kv), + this->kv.get_nks(), &(this->GG), &(this->GK), this->pw_rho, this->pw_big); + } + // 3) read the LCAO orbitals/projectors and construct the interpolation tables. + LCAO_domain::init_basis_lcao(this->pv, inp.onsite_radius, inp.lcao_ecut, + inp.lcao_dk, inp.lcao_dr, inp.lcao_rmax, ucell, two_center_bundle_, orb_); + + // 4) setup EXX calculations if (PARAM.inp.calculation == "gen_opt_abfs") { #ifdef __EXX @@ -149,7 +132,7 @@ void ESolver_KS_LCAO::before_all_runners(UnitCell& ucell, const Input_pa return; } - // 4) initialize electronic wave function psi + // 5) initialize electronic wave function psi if (this->psi == nullptr) { int nsk = 0; @@ -177,7 +160,7 @@ void ESolver_KS_LCAO::before_all_runners(UnitCell& ucell, const Input_pa this->psi = new psi::Psi(nsk, ncol, this->pv.nrow, this->kv.ngk, true); } - // 5) read psi from file + // 6) read psi from file if (inp.init_wfc == "file" && inp.esolver_type != "tddft") { if (!ModuleIO::read_wfc_nao(PARAM.globalv.global_readin_dir, @@ -192,12 +175,11 @@ void ESolver_KS_LCAO::before_all_runners(UnitCell& ucell, const Input_pa } } - // 6) initialize the density matrix - // DensityMatrix is allocated here, DMK is also initialized here - // DMR is not initialized here, it will be constructed in each before_scf + // 7) initialize the density matrix + // DMK are allocated here, but DMR is constructed in before_scf() dynamic_cast*>(this->pelec)->init_DM(&this->kv, &(this->pv), inp.nspin); - // 7) initialize exact exchange calculations + // 8) initialize exact exchange calculations #ifdef __EXX if (inp.calculation == "scf" || inp.calculation == "relax" || inp.calculation == "cell-relax" || inp.calculation == "md") @@ -224,35 +206,30 @@ void ESolver_KS_LCAO::before_all_runners(UnitCell& ucell, const Input_pa } #endif - // 8) initialize DFT+U + // 9) initialize DFT+U if (inp.dft_plus_u) { auto* dftu = ModuleDFTU::DFTU::get_instance(); dftu->init(ucell, &this->pv, this->kv.get_nks(), &orb_); } - // 9) initialize local pseudopotentials + // 10) initialize local pseudopotentials this->locpp.init_vloc(ucell, this->pw_rho); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL"); - // 10) inititlize the charge density + // 11) inititlize the charge density this->chr.allocate(inp.nspin); this->pelec->omega = ucell.omega; - // 11) initialize the potential + // 12) initialize the potential 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)); + 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)); } - // 12) initialize deepks + // 13) initialize deepks #ifdef __MLALGO LCAO_domain::DeePKS_init(ucell, pv, this->kv.get_nks(), orb_, this->ld, GlobalV::ofs_running); if (inp.deepks_scf) @@ -260,31 +237,21 @@ void ESolver_KS_LCAO::before_all_runners(UnitCell& ucell, const Input_pa // load the DeePKS model from deep neural network DeePKS_domain::load_model(inp.deepks_model, ld.model_deepks); // read pdm from file for NSCF or SCF-restart, do it only once in whole calculation - DeePKS_domain::read_pdm((inp.init_chg == "file"), - inp.deepks_equiv, - ld.init_pdm, - ucell.nat, - orb_.Alpha[0].getTotal_nchi() * ucell.nat, - ld.lmaxd, - ld.inl2l, - *orb_.Alpha, - ld.pdm); + DeePKS_domain::read_pdm((inp.init_chg == "file"), inp.deepks_equiv, + ld.init_pdm, ucell.nat, orb_.Alpha[0].getTotal_nchi() * ucell.nat, + ld.lmaxd, ld.inl2l, *orb_.Alpha, ld.pdm); } #endif - // 13) set occupations + // 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); + elecstate::fixed_weights(inp.ocp_kb, inp.nbands, inp.nelec, + this->pelec->klist, this->pelec->wg, this->pelec->skip_weights); } - // 14) if kpar is not divisible by nks, print a warning + // 15) if kpar is not divisible by nks, print a warning if (PARAM.globalv.kpar_lcao > 1) { if (this->kv.get_nks() % PARAM.globalv.kpar_lcao != 0) @@ -305,20 +272,12 @@ void ESolver_KS_LCAO::before_all_runners(UnitCell& ucell, const Input_pa } } - // 15) initialize rdmft, added by jghan + // 16) initialize rdmft, added by jghan if (inp.rdmft == true) { - rdmft_solver.init(this->GG, - this->GK, - this->pv, - ucell, - this->gd, - this->kv, - *(this->pelec), - this->orb_, - two_center_bundle_, - inp.dft_functional, - inp.rdmft_power_alpha); + rdmft_solver.init(this->GG, this->GK, this->pv, ucell, + this->gd, this->kv, *(this->pelec), this->orb_, + two_center_bundle_, inp.dft_functional, inp.rdmft_power_alpha); } ModuleBase::timer::tick("ESolver_KS_LCAO", "before_all_runners"); @@ -414,103 +373,30 @@ void ESolver_KS_LCAO::after_all_runners(UnitCell& ucell) const int nspin0 = (PARAM.inp.nspin == 2) ? 2 : 1; - // 1) write projected band structure - if (PARAM.inp.out_proj_band) - { - ModuleIO::write_proj_band_lcao(this->psi, this->pv, this->pelec, this->kv, ucell, this->p_hamilt); - } - - // 2) out ldos - if (PARAM.inp.out_ldos[0]) - { - ModuleIO::Cal_ldos::cal_ldos_lcao(reinterpret_cast*>(this->pelec), - this->psi[0], - this->Pgrid, - ucell); - } - - // 3) print out exchange-correlation potential - if (PARAM.inp.out_mat_xc) - { - ModuleIO::write_Vxc(PARAM.inp.nspin, - PARAM.globalv.nlocal, - GlobalV::DRANK, - &this->pv, - *this->psi, - ucell, - this->sf, - this->solvent, - *this->pw_rho, - *this->pw_rhod, - this->locpp.vloc, - this->chr, - this->GG, - this->GK, - this->kv, - orb_.cutoffs(), - this->pelec->wg, - this->gd -#ifdef __EXX - , - this->exd ? &this->exd->get_Hexxs() : nullptr, - this->exc ? &this->exc->get_Hexxs() : nullptr -#endif - ); - } - - if (PARAM.inp.out_mat_xc2) - { - ModuleIO::write_Vxc_R(PARAM.inp.nspin, - &this->pv, - ucell, - this->sf, - this->solvent, - *this->pw_rho, - *this->pw_rhod, - this->locpp.vloc, - this->chr, - this->GG, - this->GK, - this->kv, - orb_.cutoffs(), - this->gd -#ifdef __EXX - , - this->exd ? &this->exd->get_Hexxs() : nullptr, - this->exc ? &this->exc->get_Hexxs() : nullptr -#endif - ); - } - - // write eband terms - if (PARAM.inp.out_eband_terms) - { - ModuleIO::write_eband_terms(PARAM.inp.nspin, - PARAM.globalv.nlocal, - GlobalV::DRANK, - &this->pv, - *this->psi, - ucell, - this->sf, - this->solvent, - *this->pw_rho, - *this->pw_rhod, - this->locpp.vloc, - this->chr, - this->GG, - this->GK, - this->kv, - this->pelec->wg, - this->gd, - orb_.cutoffs(), - this->two_center_bundle_ + auto* estate = dynamic_cast*>(this->pelec); + auto* hamilt_lcao = dynamic_cast*>(this->p_hamilt); + + if(!estate) + { + ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::after_all_runners","pelec does not exist"); + } + + if(!hamilt_lcao) + { + ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::after_all_runners","p_hamilt does not exist"); + } + + ModuleIO::ctrl_runner_lcao(ucell, + PARAM.inp, this->kv, estate, this->pv, this->Pgrid, + this->gd, this->psi, this->chr, hamilt_lcao, + this->two_center_bundle_, this->GG, this->GK, + this->orb_, this->pw_rho, this->pw_rhod, + this->sf, this->locpp.vloc, #ifdef __EXX - , - this->exd ? &this->exd->get_Hexxs() : nullptr, - this->exc ? &this->exc->get_Hexxs() : nullptr + this->exd, + this->exc, #endif - ); - } + this->solvent); ModuleBase::timer::tick("ESolver_KS_LCAO", "after_all_runners"); } @@ -523,54 +409,10 @@ void ESolver_KS_LCAO::iter_init(UnitCell& ucell, const int istep, const // call iter_init() of ESolver_KS ESolver_KS::iter_init(ucell, istep, iter); - if (iter == 1) - { - this->p_chgmix->mix_reset(); // init mixing - this->p_chgmix->mixing_restart_step = PARAM.inp.scf_nmax + 1; - this->p_chgmix->mixing_restart_count = 0; - // this output will be removed once the feeature is stable - if (GlobalC::dftu.uramping > 0.01) - { - std::cout << " U-Ramping! Current U = "; - for (int i = 0; i < GlobalC::dftu.U0.size(); i++) - { - std::cout << GlobalC::dftu.U[i] * ModuleBase::Ry_to_eV << " "; - } - std::cout << " eV " << std::endl; - } - } + elecstate::DensityMatrix* dm + = dynamic_cast*>(this->pelec)->get_DM(); - // for mixing restart - if (iter == this->p_chgmix->mixing_restart_step && PARAM.inp.mixing_restart > 0.0) - { - this->p_chgmix->init_mixing(); - this->p_chgmix->mixing_restart_count++; - if (PARAM.inp.dft_plus_u) - { - GlobalC::dftu.uramping_update(); // update U by uramping if uramping > 0.01 - if (GlobalC::dftu.uramping > 0.01) - { - std::cout << " U-Ramping! Current U = "; - for (int i = 0; i < GlobalC::dftu.U0.size(); i++) - { - std::cout << GlobalC::dftu.U[i] * ModuleBase::Ry_to_eV << " "; - } - std::cout << " eV " << std::endl; - } - if (GlobalC::dftu.uramping > 0.01 && !GlobalC::dftu.u_converged()) - { - this->p_chgmix->mixing_restart_step = PARAM.inp.scf_nmax + 1; - } - } - if (PARAM.inp.mixing_dmr) // for mixing_dmr - { - // allocate memory for dmr_mdata - const elecstate::DensityMatrix* dm - = dynamic_cast*>(this->pelec)->get_DM(); - int nnr_tmp = dm->get_DMR_pointer(1)->get_nnr(); - this->p_chgmix->allocate_mixing_dmr(nnr_tmp); - } - } + module_charge::chgmixing_ks_lcao(iter, this->p_chgmix, dm->get_DMR_pointer(1)->get_nnr(), PARAM.inp); // mohan update 2012-06-05 this->pelec->f_en.deband_harris = this->pelec->cal_delta_eband(ucell); @@ -612,28 +454,6 @@ void ESolver_KS_LCAO::iter_init(UnitCell& ucell, const int istep, const this->pelec->psiToRho(*this->psi); this->pelec->skip_weights = false; - // calculate the local potential(rho) again. - // the grid integration will do in later grid integration. - - // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - // a puzzle remains here. - // if I don't renew potential, - // The scf_thr is very small. - // OneElectron, Hartree and - // Exc energy are all correct - // except the band energy. - // - // solved by mohan 2010-09-10 - // there are there rho here: - // rho1: formed by read in orbitals. - // rho2: atomic rho, used to construct H - // rho3: generated by after diagonalize - // here converged because rho3 and rho1 - // are very close. - // so be careful here, make sure - // rho1 and rho2 are the same rho. - // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - elecstate::cal_ux(ucell); //! update the potentials by using new electron charge density @@ -650,19 +470,11 @@ void ESolver_KS_LCAO::iter_init(UnitCell& ucell, const int istep, const { if (GlobalC::exx_info.info_ri.real_number) { - this->exd->exx_eachiterinit(istep, - ucell, - *dynamic_cast*>(this->pelec)->get_DM(), - this->kv, - iter); + this->exd->exx_eachiterinit(istep, ucell, *dm, this->kv, iter); } else { - this->exc->exx_eachiterinit(istep, - ucell, - *dynamic_cast*>(this->pelec)->get_DM(), - this->kv, - iter); + this->exc->exx_eachiterinit(istep, ucell, *dm, this->kv, iter); } } #endif @@ -671,7 +483,7 @@ void ESolver_KS_LCAO::iter_init(UnitCell& ucell, const int istep, const { if (istep != 0 || iter != 1) { - GlobalC::dftu.set_dmr(dynamic_cast*>(this->pelec)->get_DM()); + GlobalC::dftu.set_dmr(dm); } // Calculate U and J if Yukawa potential is used GlobalC::dftu.cal_slater_UJ(ucell, this->chr.rho, this->pw_rho->nrxx); @@ -702,7 +514,6 @@ void ESolver_KS_LCAO::iter_init(UnitCell& ucell, const int istep, const // save density matrix DMR for mixing if (PARAM.inp.mixing_restart > 0 && PARAM.inp.mixing_dmr && this->p_chgmix->mixing_restart_count > 0) { - elecstate::DensityMatrix* dm = dynamic_cast*>(this->pelec)->get_DM(); dm->save_DMR(); } } @@ -792,8 +603,23 @@ void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& { ModuleBase::TITLE("ESolver_KS_LCAO", "iter_finish"); - // 1) calculate the local occupation number matrix and energy correction - // in DFT+U + auto* estate = dynamic_cast*>(this->pelec); + auto* hamilt_lcao = dynamic_cast*>(this->p_hamilt); + + if(!estate) + { + ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::iter_finish","pelec does not exist"); + } + + if(!hamilt_lcao) + { + ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::iter_finish","p_hamilt does not exist"); + } + + const std::vector>& dm_vec = estate->get_DM()->get_DMK_vector(); + + + // 1) calculate the local occupation number matrix and energy correction in DFT+U if (PARAM.inp.dft_plus_u) { // only old DFT+U method should calculated energy correction in esolver, @@ -802,14 +628,12 @@ void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& { if (GlobalC::dftu.omc != 2) { - const std::vector>& tmp_dm - = dynamic_cast*>(this->pelec)->get_DM()->get_DMK_vector(); ModuleDFTU::dftu_cal_occup_m(iter, ucell, - tmp_dm, + dm_vec, this->kv, this->p_chgmix->get_mixing_beta(), - this->p_hamilt); + hamilt_lcao); } GlobalC::dftu.cal_energy_correction(ucell, istep); } @@ -820,13 +644,10 @@ void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& #ifdef __MLALGO if (PARAM.inp.deepks_scf) { - const std::vector>& dm - = dynamic_cast*>(this->pelec)->get_DM()->get_DMK_vector(); - - ld.dpks_cal_e_delta_band(dm, this->kv.get_nks()); - DeePKS_domain::update_dmr(this->kv.kvec_d, dm, ucell, orb_, this->pv, this->gd, ld.dm_r); - this->pelec->f_en.edeepks_scf = ld.E_delta - ld.e_delta_band; - this->pelec->f_en.edeepks_delta = ld.E_delta; + ld.dpks_cal_e_delta_band(dm_vec, this->kv.get_nks()); + DeePKS_domain::update_dmr(this->kv.kvec_d, dm_vec, ucell, orb_, this->pv, this->gd, ld.dm_r); + estate->f_en.edeepks_scf = ld.E_delta - ld.e_delta_band; + estate->f_en.edeepks_delta = ld.E_delta; } #endif @@ -837,93 +658,42 @@ void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& sc.cal_mi_lcao(iter); } - // 4) call iter_finish() of ESolver_KS + // call iter_finish() of ESolver_KS, where band gap is printed, + // 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); - // 5) mix density matrix if mixing_restart + mixing_dmr + not first + // mix density matrix if mixing_restart + mixing_dmr + not first // mixing_restart at every iter except the last iter if(iter != PARAM.inp.scf_nmax && !conv_esolver) { if (PARAM.inp.mixing_restart > 0 && this->p_chgmix->mixing_restart_count > 0 && PARAM.inp.mixing_dmr) { - elecstate::DensityMatrix* dm = dynamic_cast*>(this->pelec)->get_DM(); + elecstate::DensityMatrix* dm = estate->get_DM(); this->p_chgmix->mix_dmr(dm); } } - // 6) save charge density - // Peize Lin add 2020.04.04 - if (GlobalC::restart.info_save.save_charge) - { - for (int is = 0; is < PARAM.inp.nspin; ++is) - { - GlobalC::restart.save_disk("charge", is, this->chr.nrxx, this->chr.rho[is]); - } - } - -#ifdef __EXX - // 7) save exx matrix - if (PARAM.inp.calculation != "nscf") - { - if (GlobalC::exx_info.info_global.cal_exx) - { - GlobalC::exx_info.info_ri.real_number ? this->exd->exx_iter_finish(this->kv, - ucell, - *this->p_hamilt, - *this->pelec, - *this->p_chgmix, - this->scf_ene_thr, - iter, - istep, - conv_esolver) - : this->exc->exx_iter_finish(this->kv, - ucell, - *this->p_hamilt, - *this->pelec, - *this->p_chgmix, - this->scf_ene_thr, - iter, - istep, - conv_esolver); - } - } -#endif - // 8) use the converged occupation matrix for next MD/Relax SCF calculation + // use the converged occupation matrix for next MD/Relax SCF calculation if (PARAM.inp.dft_plus_u && conv_esolver) { GlobalC::dftu.initialed_locale = true; } - // 9) for deepks, output labels during electronic steps (after conv_esolver is renewed) + // control the output related to the finished iteration + ModuleIO::ctrl_iter_lcao(ucell, PARAM.inp, this->kv, estate, + this->pv, this->gd, this->psi, this->chr, this->p_chgmix, + hamilt_lcao, this->orb_, #ifdef __MLALGO - if (PARAM.inp.deepks_out_labels >0 && PARAM.inp.deepks_out_freq_elec) - { - if (iter % PARAM.inp.deepks_out_freq_elec == 0 ) - { - hamilt::HamiltLCAO* p_ham_deepks = dynamic_cast*>(this->p_hamilt); - std::shared_ptr> ld_shared_ptr(&ld, [](LCAO_Deepks*) {}); - LCAO_Deepks_Interface deepks_interface(ld_shared_ptr); - - deepks_interface.out_deepks_labels(this->pelec->f_en.etot, - this->kv.get_nks(), - ucell.nat, - PARAM.globalv.nlocal, - this->pelec->ekb, - this->kv.kvec_d, - ucell, - orb_, - this->gd, - &(this->pv), - *(this->psi), - dynamic_cast*>(this->pelec)->get_DM(), - p_ham_deepks, - iter, - conv_esolver, - GlobalV::MY_RANK, - GlobalV::ofs_running); - } - } + this->ld, #endif +#ifdef __EXX + *this->exd, + *this->exc, +#endif + iter, istep, conv_esolver, this->scf_ene_thr); + } template class ESolver_KS_LCAO; diff --git a/source/source_esolver/esolver_ks_pw.cpp b/source/source_esolver/esolver_ks_pw.cpp index 6bb78b5c25..97d4f61533 100644 --- a/source/source_esolver/esolver_ks_pw.cpp +++ b/source/source_esolver/esolver_ks_pw.cpp @@ -22,11 +22,12 @@ #include "source_base/kernels/dsp/dsp_connector.h" #endif -#include +//#include #include "source_pw/module_pwdft/setup_pot.h" // mohan add 20250929 #include "source_estate/setup_estate_pw.h" // mohan add 20251005 #include "source_io/ctrl_output_pw.h" // mohan add 20250927 +#include "source_estate/module_charge/chgmixing.h" // use charge mixing, mohan add 20251006 namespace ModuleESolver { @@ -185,50 +186,8 @@ void ESolver_KS_PW::iter_init(UnitCell& ucell, const int istep, const // Call iter_init() of ESolver_KS ESolver_KS::iter_init(ucell, istep, iter); - if (iter == 1) - { - this->p_chgmix->init_mixing(); - this->p_chgmix->mixing_restart_step = PARAM.inp.scf_nmax + 1; - } - - // For mixing restart - if (iter == this->p_chgmix->mixing_restart_step && PARAM.inp.mixing_restart > 0.0) - { - this->p_chgmix->init_mixing(); - this->p_chgmix->mixing_restart_count++; - - if (PARAM.inp.dft_plus_u) - { - auto* dftu = ModuleDFTU::DFTU::get_instance(); - if (dftu->uramping > 0.01 && !dftu->u_converged()) - { - this->p_chgmix->mixing_restart_step = PARAM.inp.scf_nmax + 1; - } - if (dftu->uramping > 0.01) - { - bool do_uramping = true; - if (PARAM.inp.sc_mag_switch) - { - spinconstrain::SpinConstrain>& sc - = spinconstrain::SpinConstrain>::getScInstance(); - if (!sc.mag_converged()) // skip uramping if mag not converged - { - do_uramping = false; - } - } - if (do_uramping) - { - dftu->uramping_update(); // update U by uramping if uramping > 0.01 - std::cout << " U-Ramping! Current U = "; - for (int i = 0; i < dftu->U0.size(); i++) - { - std::cout << dftu->U[i] * ModuleBase::Ry_to_eV << " "; - } - std::cout << " eV " << std::endl; - } - } - } - } + // perform charge mixing + module_charge::chgmixing_ks_pw(iter, this->p_chgmix, PARAM.inp); // mohan move harris functional to here, 2012-06-05 // use 'rho(in)' and 'v_h and v_xc'(in) diff --git a/source/source_esolver/lcao_after_scf.cpp b/source/source_esolver/lcao_after_scf.cpp index 30a4cd6997..bf819f8a75 100644 --- a/source/source_esolver/lcao_after_scf.cpp +++ b/source/source_esolver/lcao_after_scf.cpp @@ -1,5 +1,5 @@ #include "esolver_ks_lcao.h" -#include "source_io/ctrl_output_lcao.h" +#include "source_io/ctrl_scf_lcao.h" namespace ModuleESolver { @@ -23,7 +23,12 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep, const if(!estate) { - ModuleBase::WARNING_QUIT("ModuleIO::ctrl_output_lcao","pelec does not exist"); + ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::after_scf","pelec does not exist"); + } + + if(!hamilt_lcao) + { + ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::after_scf","p_hamilt does not exist"); } @@ -49,23 +54,13 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep, const if (out_flag) { - ModuleIO::ctrl_output_lcao(ucell, - PARAM.inp, - this->kv, - estate, - this->pv, - this->gd, - this->psi, - hamilt_lcao, - this->two_center_bundle_, - this->GK, - this->orb_, - this->pw_wfc, - this->pw_rho, - this->GridT, - this->pw_big, - this->sf, - this->rdmft_solver, + ModuleIO::ctrl_scf_lcao(ucell, + PARAM.inp, this->kv, estate, this->pv, + this->gd, this->psi, hamilt_lcao, + this->two_center_bundle_, this->GK, + this->orb_, this->pw_wfc, this->pw_rho, + this->GridT, this->pw_big, this->sf, + this->rdmft_solver, #ifdef __MLALGO this->ld, #endif diff --git a/source/source_estate/CMakeLists.txt b/source/source_estate/CMakeLists.txt index e22dbfc3b0..f1bd096009 100644 --- a/source/source_estate/CMakeLists.txt +++ b/source/source_estate/CMakeLists.txt @@ -17,6 +17,7 @@ list(APPEND objects module_pot/potential_new.cpp module_pot/potential_types.cpp module_pot/pot_sep.cpp + module_charge/chgmixing.cpp module_charge/charge.cpp module_charge/charge_init.cpp module_charge/charge_mpi.cpp diff --git a/source/source_estate/module_charge/chgmixing.cpp b/source/source_estate/module_charge/chgmixing.cpp new file mode 100644 index 0000000000..87496ced87 --- /dev/null +++ b/source/source_estate/module_charge/chgmixing.cpp @@ -0,0 +1,111 @@ +#include "source_estate/module_charge/chgmixing.h" +#include "source_lcao/module_dftu/dftu.h" +#include "source_lcao/module_deltaspin/spin_constrain.h" + +void module_charge::chgmixing_ks_pw(const int iter, // scf iteration number + Charge_Mixing* p_chgmix, // charge mixing class + const Input_para& inp) // input parameters +{ + ModuleBase::TITLE("module_charge", "chgmixing_ks_pw"); + + if (iter == 1) + { + p_chgmix->init_mixing(); + p_chgmix->mixing_restart_step = inp.scf_nmax + 1; + } + + // For mixing restart + if (iter == p_chgmix->mixing_restart_step && inp.mixing_restart > 0.0) + { + p_chgmix->init_mixing(); + p_chgmix->mixing_restart_count++; + + if (inp.dft_plus_u) + { + auto* dftu = ModuleDFTU::DFTU::get_instance(); + if (dftu->uramping > 0.01 && !dftu->u_converged()) + { + p_chgmix->mixing_restart_step = inp.scf_nmax + 1; + } + if (dftu->uramping > 0.01) + { + bool do_uramping = true; + if (inp.sc_mag_switch) + { + spinconstrain::SpinConstrain>& sc + = spinconstrain::SpinConstrain>::getScInstance(); + if (!sc.mag_converged()) // skip uramping if mag not converged + { + do_uramping = false; + } + } + if (do_uramping) + { + dftu->uramping_update(); // update U by uramping if uramping > 0.01 + std::cout << " U-Ramping! Current U = "; + for (int i = 0; i < dftu->U0.size(); i++) + { + std::cout << dftu->U[i] * ModuleBase::Ry_to_eV << " "; + } + std::cout << " eV " << std::endl; + } + } + } + } + + return; +} + +void module_charge::chgmixing_ks_lcao(const int iter, // scf iteration number + Charge_Mixing* p_chgmix, // charge mixing class + const int nnr, // dimension of density matrix + const Input_para& inp) // input parameters +{ + ModuleBase::TITLE("module_charge", "chgmixing_ks_lcao"); + + if (iter == 1) + { + p_chgmix->mix_reset(); // init mixing + p_chgmix->mixing_restart_step = inp.scf_nmax + 1; + p_chgmix->mixing_restart_count = 0; + // this output will be removed once the feeature is stable + if (GlobalC::dftu.uramping > 0.01) + { + std::cout << " U-Ramping! Current U = "; + for (int i = 0; i < GlobalC::dftu.U0.size(); i++) + { + std::cout << GlobalC::dftu.U[i] * ModuleBase::Ry_to_eV << " "; + } + std::cout << " eV " << std::endl; + } + } + + // for mixing restart + if (iter == p_chgmix->mixing_restart_step && inp.mixing_restart > 0.0) + { + p_chgmix->init_mixing(); + p_chgmix->mixing_restart_count++; + if (inp.dft_plus_u) + { + GlobalC::dftu.uramping_update(); // update U by uramping if uramping > 0.01 + if (GlobalC::dftu.uramping > 0.01) + { + std::cout << " U-Ramping! Current U = "; + for (int i = 0; i < GlobalC::dftu.U0.size(); i++) + { + std::cout << GlobalC::dftu.U[i] * ModuleBase::Ry_to_eV << " "; + } + std::cout << " eV " << std::endl; + } + if (GlobalC::dftu.uramping > 0.01 && !GlobalC::dftu.u_converged()) + { + p_chgmix->mixing_restart_step = inp.scf_nmax + 1; + } + } + if (inp.mixing_dmr) // for mixing_dmr + { + // allocate memory for dmr_mdata + p_chgmix->allocate_mixing_dmr(nnr); + } + } +} diff --git a/source/source_estate/module_charge/chgmixing.h b/source/source_estate/module_charge/chgmixing.h new file mode 100644 index 0000000000..2e962ed8e5 --- /dev/null +++ b/source/source_estate/module_charge/chgmixing.h @@ -0,0 +1,22 @@ +#ifndef CHGMIXING_H +#define CHGMIXING_H + +#include "source_estate/module_charge/charge_mixing.h" +#include "source_io/module_parameter/input_parameter.h" // use Input_para + +namespace module_charge +{ + +void chgmixing_ks_pw(const int iter, + Charge_Mixing* p_chgmix, + const Input_para& inp); // input parameters + +void chgmixing_ks_lcao(const int iter, // scf iteration number + Charge_Mixing* p_chgmix, // charge mixing class + const int nnr, // dimension of density matrix + const Input_para& inp); // input parameters + +} + + +#endif diff --git a/source/source_io/CMakeLists.txt b/source/source_io/CMakeLists.txt index e511ffb35c..850a124402 100644 --- a/source/source_io/CMakeLists.txt +++ b/source/source_io/CMakeLists.txt @@ -83,7 +83,9 @@ if(ENABLE_LCAO) single_R_io.cpp cal_r_overlap_R.cpp output_mat_sparse.cpp - ctrl_output_lcao.cpp + ctrl_scf_lcao.cpp + ctrl_runner_lcao.cpp + ctrl_iter_lcao.cpp ) endif() diff --git a/source/source_io/ctrl_iter_lcao.cpp b/source/source_io/ctrl_iter_lcao.cpp new file mode 100644 index 0000000000..d14ba00f63 --- /dev/null +++ b/source/source_io/ctrl_iter_lcao.cpp @@ -0,0 +1,156 @@ +#include "source_io/ctrl_iter_lcao.h" // use ctrl_iter_lcao() +#include "source_pw/module_pwdft/global.h" // use GlobalC::restart + +#ifdef __MLALGO +#include "source_lcao/module_deepks/LCAO_deepks.h" +#include "source_lcao/module_deepks/LCAO_deepks_interface.h" +#endif + +namespace ModuleIO +{ + +template +void ctrl_iter_lcao(UnitCell& ucell, // unit cell * + const Input_para& inp, // input parameters * + K_Vectors& kv, // k points * + elecstate::ElecStateLCAO* pelec, // electronic info * + Parallel_Orbitals& pv, // parallel orbital info * + Grid_Driver& gd, // adjacent atom info * + psi::Psi* psi, // wave functions * + Charge &chr, // charge density * + Charge_Mixing* p_chgmix, // charge mixing * + hamilt::HamiltLCAO* p_hamilt, // hamiltonian * + LCAO_Orbitals &orb, // orbital info * +#ifdef __MLALGO + LCAO_Deepks& ld, +#endif +#ifdef __EXX + Exx_LRI_Interface& exd, + Exx_LRI_Interface>& exc, +#endif + int &iter, + const int istep, + bool &conv_esolver, + const double &scf_ene_thr) +{ + ModuleBase::TITLE("ModuleIO", "ctrl_iter_lcao"); + ModuleBase::timer::tick("ModuleIO", "ctrl_iter_lcao"); + + // save charge density + // Peize Lin add 2020.04.04 + if (GlobalC::restart.info_save.save_charge) + { + for (int is = 0; is < inp.nspin; ++is) + { + GlobalC::restart.save_disk("charge", is, chr.nrxx, chr.rho[is]); + } + } + +#ifdef __EXX + // save exx matrix + if (inp.calculation != "nscf") + { + if (GlobalC::exx_info.info_global.cal_exx) + { + GlobalC::exx_info.info_ri.real_number ? + exd.exx_iter_finish(kv, ucell, *p_hamilt, *pelec, + *p_chgmix, scf_ene_thr, iter, istep, conv_esolver) : + exc.exx_iter_finish(kv, ucell, *p_hamilt, *pelec, + *p_chgmix, scf_ene_thr, iter, istep, conv_esolver); + } + } +#endif + + // for deepks, output labels during electronic steps (after conv_esolver is renewed) +#ifdef __MLALGO + if (inp.deepks_out_labels >0 && inp.deepks_out_freq_elec) + { + if (iter % inp.deepks_out_freq_elec == 0 ) + { + std::shared_ptr> ld_shared_ptr(&ld, [](LCAO_Deepks*) {}); + LCAO_Deepks_Interface deepks_interface(ld_shared_ptr); + + deepks_interface.out_deepks_labels(pelec->f_en.etot, kv.get_nks(), + ucell.nat, PARAM.globalv.nlocal, pelec->ekb, kv.kvec_d, + ucell, orb, gd, &pv, *psi, pelec->get_DM(), + p_hamilt, iter, conv_esolver, GlobalV::MY_RANK, GlobalV::ofs_running); + } + } +#endif + + ModuleBase::timer::tick("ModuleIO", "ctrl_iter_lcao"); +} + +// TK: double TR: double +template void ctrl_iter_lcao(UnitCell& ucell, // unit cell * + const Input_para& inp, // input parameters * + K_Vectors& kv, // k points * + elecstate::ElecStateLCAO* pelec, // electronic info * + Parallel_Orbitals& pv, // parallel orbital info * + Grid_Driver& gd, // adjacent atom info * + psi::Psi* psi, // wave functions * + Charge &chr, // charge density * + Charge_Mixing* p_chgmix, // charge mixing * + hamilt::HamiltLCAO* p_hamilt, // hamiltonian * + LCAO_Orbitals &orb, // orbital info * +#ifdef __MLALGO + LCAO_Deepks& ld, +#endif +#ifdef __EXX + Exx_LRI_Interface& exd, + Exx_LRI_Interface>& exc, +#endif + int &iter, + const int istep, + bool &conv_esolver, + const double &scf_ene_thr); + +// TK: complex TR: double +template void ctrl_iter_lcao, double>(UnitCell& ucell, // unit cell * + const Input_para& inp, // input parameters * + K_Vectors& kv, // k points * + elecstate::ElecStateLCAO>* pelec, // electronic info * + Parallel_Orbitals& pv, // parallel orbital info * + Grid_Driver& gd, // adjacent atom info * + psi::Psi>* psi, // wave functions * + Charge &chr, // charge density * + Charge_Mixing* p_chgmix, // charge mixing * + hamilt::HamiltLCAO, double>* p_hamilt, // hamiltonian * + LCAO_Orbitals &orb, // orbital info * +#ifdef __MLALGO + LCAO_Deepks>& ld, +#endif +#ifdef __EXX + Exx_LRI_Interface, double>& exd, + Exx_LRI_Interface, std::complex>& exc, +#endif + int &iter, + const int istep, + bool &conv_esolver, + const double &scf_ene_thr); + +// TK: complex TR: complex +template void ctrl_iter_lcao, std::complex>(UnitCell& ucell, // unit cell * + const Input_para& inp, // input parameters * + K_Vectors& kv, // k points * + elecstate::ElecStateLCAO>* pelec, // electronic info * + Parallel_Orbitals& pv, // parallel orbital info * + Grid_Driver& gd, // adjacent atom info * + psi::Psi>* psi, // wave functions * + Charge &chr, // charge density * + Charge_Mixing* p_chgmix, // charge mixing * + hamilt::HamiltLCAO, std::complex>* p_hamilt, // hamiltonian * + LCAO_Orbitals &orb, // orbital info * +#ifdef __MLALGO + LCAO_Deepks>& ld, +#endif +#ifdef __EXX + Exx_LRI_Interface, double>& exd, + Exx_LRI_Interface, std::complex>& exc, +#endif + int &iter, + const int istep, + bool &conv_esolver, + const double &scf_ene_thr); + +} // end ModuleIO diff --git a/source/source_io/ctrl_iter_lcao.h b/source/source_io/ctrl_iter_lcao.h new file mode 100644 index 0000000000..7b6dc32b47 --- /dev/null +++ b/source/source_io/ctrl_iter_lcao.h @@ -0,0 +1,43 @@ +#ifndef CTRL_ITER_LCAO_H +#define CTRL_ITER_LCAO_H + +#include "source_cell/unitcell.h" // use UnitCell +#include "source_cell/klist.h" // use K_Vectors +#include "source_estate/elecstate_lcao.h" // use elecstate::ElecStateLCAO +#include "source_psi/psi.h" // use Psi +#include "source_estate/module_charge/charge.h" // use charge +#include "source_estate/module_charge/charge_mixing.h" // use charge mixing +#include "source_lcao/hamilt_lcao.h" // use hamilt::HamiltLCAO +#ifdef __EXX +#include "source_lcao/module_ri/Exx_LRI_interface.h" // use EXX codes +#endif + +namespace ModuleIO +{ + +template +void ctrl_iter_lcao(UnitCell& ucell, // unit cell * + const Input_para& inp, // input parameters * + K_Vectors& kv, // k points * + elecstate::ElecStateLCAO* pelec, // electronic info * + Parallel_Orbitals& pv, // parallel orbital info * + Grid_Driver& gd, // adjacent atom info * + psi::Psi* psi, // wave functions * + Charge &chr, // charge density * + Charge_Mixing* p_chgmix, // charge mixing * + hamilt::HamiltLCAO* p_hamilt, // hamiltonian * + LCAO_Orbitals &orb, // orbital info * +#ifdef __MLALGO + LCAO_Deepks& ld, +#endif +#ifdef __EXX + Exx_LRI_Interface& exd, + Exx_LRI_Interface>& exc, +#endif + int &iter, + const int istep, + bool &conv_esolver, + const double &scf_ene_thr); + +} +#endif diff --git a/source/source_io/ctrl_runner_lcao.cpp b/source/source_io/ctrl_runner_lcao.cpp new file mode 100644 index 0000000000..d92f08ff3d --- /dev/null +++ b/source/source_io/ctrl_runner_lcao.cpp @@ -0,0 +1,219 @@ +#include "source_io/ctrl_runner_lcao.h" // use ctrl_runner_lcao() + +#include "source_estate/elecstate_lcao.h" // use elecstate::ElecState +#include "source_lcao/hamilt_lcao.h" // use hamilt::HamiltLCAO + +#include "source_io/write_proj_band_lcao.h" // projcted band structure +#include "source_io/cal_ldos.h" // cal LDOS +#include "source_io/write_eband_terms.hpp" +#include "source_io/write_vxc.hpp" +#include "source_io/write_vxc_r.hpp" + +namespace ModuleIO +{ + +template +void ctrl_runner_lcao(UnitCell& ucell, // unitcell + const Input_para &inp, // input + K_Vectors &kv, // k-point + elecstate::ElecStateLCAO* pelec,// electronic info + Parallel_Orbitals &pv, // orbital info + Parallel_Grid &pgrid, // grid info + Grid_Driver &gd, // search for adjacent atoms + psi::Psi* psi, // wave function + Charge &chr, // charge density + hamilt::HamiltLCAO* p_hamilt, // hamiltonian + TwoCenterBundle &two_center_bundle, // use two-center integration + Gint_Gamma &gg, // gint for Gamma-only + Gint_k &gk, // gint for multi k-points + LCAO_Orbitals &orb, // LCAO orbitals + ModulePW::PW_Basis* pw_rho, // charge density + ModulePW::PW_Basis* pw_rhod, // dense charge density + Structure_Factor &sf, // structure factor + ModuleBase::matrix &vloc, // local pseudopotential +#ifdef __EXX + std::shared_ptr> exd, + std::shared_ptr>> exc, +#endif + surchem &solvent) // solvent model +{ + ModuleBase::TITLE("ModuleIO", "ctrl_runner_lcao"); + ModuleBase::timer::tick("ModuleIO", "ctrl_runner_lcao"); + + // 1) write projected band structure + if (inp.out_proj_band) + { + ModuleIO::write_proj_band_lcao(psi, pv, pelec, kv, ucell, p_hamilt); + } + + // 2) out ldos + if (inp.out_ldos[0]) + { + ModuleIO::Cal_ldos::cal_ldos_lcao(pelec, psi[0], pgrid, ucell); + } + + // 3) print out exchange-correlation potential + if (inp.out_mat_xc) + { + ModuleIO::write_Vxc(inp.nspin, + PARAM.globalv.nlocal, + GlobalV::DRANK, + &pv, + *psi, + ucell, + sf, + solvent, + *pw_rho, + *pw_rhod, + vloc, + chr, + gg, + gk, + kv, + orb.cutoffs(), + pelec->wg, + gd +#ifdef __EXX + , + exd ? &exd->get_Hexxs() : nullptr, + exc ? &exc->get_Hexxs() : nullptr +#endif + ); + } + + if (inp.out_mat_xc2) + { + ModuleIO::write_Vxc_R(inp.nspin, + &pv, + ucell, + sf, + solvent, + *pw_rho, + *pw_rhod, + vloc, + chr, + gg, + gk, + kv, + orb.cutoffs(), + gd +#ifdef __EXX + , + exd ? &exd->get_Hexxs() : nullptr, + exc ? &exc->get_Hexxs() : nullptr +#endif + ); + } + + + // write eband terms + if (inp.out_eband_terms) + { + ModuleIO::write_eband_terms(inp.nspin, + PARAM.globalv.nlocal, + GlobalV::DRANK, + &pv, + *psi, + ucell, + sf, + solvent, + *pw_rho, + *pw_rhod, + vloc, + chr, + gg, + gk, + kv, + pelec->wg, + gd, + orb.cutoffs(), + two_center_bundle +#ifdef __EXX + , + exd ? &exd->get_Hexxs() : nullptr, + exc ? &exc->get_Hexxs() : nullptr +#endif + ); + } + +} + + + + +// TK: double TR: double +template void ModuleIO::ctrl_runner_lcao(UnitCell& ucell, // unitcell + const Input_para &inp, // input + K_Vectors &kv, // k-point + elecstate::ElecStateLCAO* pelec,// electronic info + Parallel_Orbitals &pv, // orbital info + Parallel_Grid &pgrid, // grid info + Grid_Driver &gd, // search for adjacent atoms + psi::Psi* psi, // wave function + Charge &chr, // charge density + hamilt::HamiltLCAO* p_hamilt, // hamiltonian + TwoCenterBundle &two_center_bundle, // use two-center integration + Gint_Gamma &gg, // gint for Gamma-only + Gint_k &gk, // gint for multi k-points + LCAO_Orbitals &orb, // LCAO orbitals + ModulePW::PW_Basis* pw_rho, // charge density + ModulePW::PW_Basis* pw_rhod, // dense charge density + Structure_Factor &sf, // structure factor + ModuleBase::matrix &vloc, // local pseudopotential +#ifdef __EXX + std::shared_ptr> exd, + std::shared_ptr>> exc, +#endif + surchem &solvent); // solvent model + +// TK: complex TR: double +template void ctrl_runner_lcao, double>(UnitCell& ucell, // unitcell + const Input_para &inp, // input + K_Vectors &kv, // k-point + elecstate::ElecStateLCAO>* pelec,// electronic info + Parallel_Orbitals &pv, // orbital info + Parallel_Grid &pgrid, // grid info + Grid_Driver &gd, // search for adjacent atoms + psi::Psi>* psi, // wave function + Charge &chr, // charge density + hamilt::HamiltLCAO, double>* p_hamilt, // hamiltonian + TwoCenterBundle &two_center_bundle, // use two-center integration + Gint_Gamma &gg, // gint for Gamma-only + Gint_k &gk, // gint for multi k-points + LCAO_Orbitals &orb, // LCAO orbitals + ModulePW::PW_Basis* pw_rho, // charge density + ModulePW::PW_Basis* pw_rhod, // dense charge density + Structure_Factor &sf, // structure factor + ModuleBase::matrix &vloc, // local pseudopotential +#ifdef __EXX + std::shared_ptr, double>> exd, + std::shared_ptr, std::complex>> exc, +#endif + surchem &solvent); // solvent model + +// TK: complex TR: complex +template void ctrl_runner_lcao, std::complex>(UnitCell& ucell, // unitcell + const Input_para &inp, // input + K_Vectors &kv, // k-point + elecstate::ElecStateLCAO>* pelec,// electronic info + Parallel_Orbitals &pv, // orbital info + Parallel_Grid &pgrid, // grid info + Grid_Driver &gd, // search for adjacent atoms + psi::Psi>* psi, // wave function + Charge &chr, // charge density + hamilt::HamiltLCAO, std::complex>* p_hamilt, // hamiltonian + TwoCenterBundle &two_center_bundle, // use two-center integration + Gint_Gamma &gg, // gint for Gamma-only + Gint_k &gk, // gint for multi k-points + LCAO_Orbitals &orb, // LCAO orbitals + ModulePW::PW_Basis* pw_rho, // charge density + ModulePW::PW_Basis* pw_rhod, // dense charge density + Structure_Factor &sf, // structure factor + ModuleBase::matrix &vloc, // local pseudopotential +#ifdef __EXX + std::shared_ptr, double>> exd, + std::shared_ptr, std::complex>> exc, +#endif + surchem &solvent); // solvent model + +} // end namespace diff --git a/source/source_io/ctrl_runner_lcao.h b/source/source_io/ctrl_runner_lcao.h new file mode 100644 index 0000000000..20318bba94 --- /dev/null +++ b/source/source_io/ctrl_runner_lcao.h @@ -0,0 +1,45 @@ +#ifndef CTRL_RUNNER_LCAO_H +#define CTRL_RUNNER_LCAO_H + +#include "source_cell/unitcell.h" // use UnitCell +#include "source_cell/klist.h" // use K_Vectors +#include "source_estate/elecstate_lcao.h" // use elecstate::ElecStateLCAO +#include "source_psi/psi.h" // use Psi +#include "source_lcao/hamilt_lcao.h" // use hamilt::HamiltLCAO +#include "source_basis/module_nao/two_center_bundle.h" // use TwoCenterBundle +#include "source_lcao/module_gint/gint_k.h" // use Gint_k +#ifdef __EXX +#include "source_lcao/module_ri/Exx_LRI_interface.h" // use EXX codes +#endif + +namespace ModuleIO +{ + +template +void ctrl_runner_lcao(UnitCell& ucell, // unitcell + const Input_para &inp, // input + K_Vectors &kv, // k-point + elecstate::ElecStateLCAO* pelec,// electronic info + Parallel_Orbitals &pv, // orbital info + Parallel_Grid &pgrid, // grid info + Grid_Driver &gd, // search for adjacent atoms + psi::Psi* psi, // wave function + Charge &chr, // charge density + hamilt::HamiltLCAO* p_hamilt, // hamiltonian + TwoCenterBundle &two_center_bundle, // use two-center integration + Gint_Gamma &gg, // gint for Gamma-only + Gint_k &gk, // gint for multi k-points + LCAO_Orbitals &orb, // LCAO orbitals + ModulePW::PW_Basis* pw_rho, // charge density + ModulePW::PW_Basis* pw_rhod, // dense charge density + Structure_Factor &sf, // structure factor + ModuleBase::matrix &vloc, // local pseudopotential +#ifdef __EXX + std::shared_ptr> exd, + std::shared_ptr>> exc, +#endif + surchem &solvent); // solvent model + +} + +#endif diff --git a/source/source_io/ctrl_output_lcao.cpp b/source/source_io/ctrl_scf_lcao.cpp similarity index 96% rename from source/source_io/ctrl_output_lcao.cpp rename to source/source_io/ctrl_scf_lcao.cpp index 5722265307..c981e631e3 100644 --- a/source/source_io/ctrl_output_lcao.cpp +++ b/source/source_io/ctrl_scf_lcao.cpp @@ -1,7 +1,7 @@ #include #include "source_estate/elecstate_lcao.h" // use elecstate::ElecState -#include "source_io/ctrl_output_lcao.h" // use ctrl_output_lcao() +#include "source_io/ctrl_scf_lcao.h" // use ctrl_scf_lcao() #include "source_lcao/hamilt_lcao.h" // use hamilt::HamiltLCAO #include "source_hamilt/hamilt.h" // use Hamilt @@ -34,7 +34,7 @@ namespace ModuleIO { template -void ctrl_output_lcao(UnitCell& ucell, +void ctrl_scf_lcao(UnitCell& ucell, const Input_para& inp, K_Vectors& kv, elecstate::ElecStateLCAO* pelec, @@ -60,8 +60,8 @@ void ctrl_output_lcao(UnitCell& ucell, #endif const int istep) { - ModuleBase::TITLE("ModuleIO", "ctrl_output_lcao"); - ModuleBase::timer::tick("ModuleIO", "ctrl_output_lcao"); + ModuleBase::TITLE("ModuleIO", "ctrl_scf_lcao"); + ModuleBase::timer::tick("ModuleIO", "ctrl_scf_lcao"); const bool out_app_flag = inp.out_app_flag; const bool gamma_only = PARAM.globalv.gamma_only_local; @@ -434,14 +434,14 @@ void ctrl_output_lcao(UnitCell& ucell, } - ModuleBase::timer::tick("ModuleIO", "ctrl_output_lcao"); + ModuleBase::timer::tick("ModuleIO", "ctrl_scf_lcao"); } } // End ModuleIO // For gamma only -template void ModuleIO::ctrl_output_lcao(UnitCell& ucell, +template void ModuleIO::ctrl_scf_lcao(UnitCell& ucell, const Input_para& inp, K_Vectors& kv, elecstate::ElecStateLCAO* pelec, @@ -468,7 +468,7 @@ template void ModuleIO::ctrl_output_lcao(UnitCell& ucell, const int istep); // For multiple k-points -template void ModuleIO::ctrl_output_lcao, double>(UnitCell& ucell, +template void ModuleIO::ctrl_scf_lcao, double>(UnitCell& ucell, const Input_para& inp, K_Vectors& kv, elecstate::ElecStateLCAO>* pelec, @@ -494,7 +494,7 @@ template void ModuleIO::ctrl_output_lcao, double>(UnitCell& #endif const int istep); -template void ModuleIO::ctrl_output_lcao, std::complex>(UnitCell& ucell, +template void ModuleIO::ctrl_scf_lcao, std::complex>(UnitCell& ucell, const Input_para& inp, K_Vectors& kv, elecstate::ElecStateLCAO>* pelec, diff --git a/source/source_io/ctrl_output_lcao.h b/source/source_io/ctrl_scf_lcao.h similarity index 94% rename from source/source_io/ctrl_output_lcao.h rename to source/source_io/ctrl_scf_lcao.h index a223e5381d..8e3dcf4d29 100644 --- a/source/source_io/ctrl_output_lcao.h +++ b/source/source_io/ctrl_scf_lcao.h @@ -1,5 +1,5 @@ -#ifndef CTRL_OUTPUT_LCAO_H -#define CTRL_OUTPUT_LCAO_H +#ifndef CTRL_SCF_LCAO_H +#define CTRL_SCF_LCAO_H #include @@ -21,7 +21,7 @@ namespace ModuleIO { // in principle, we need to add const for all of the variables, mohan note 2025-06-05 template - void ctrl_output_lcao(UnitCell& ucell, + void ctrl_scf_lcao(UnitCell& ucell, const Input_para& inp, K_Vectors& kv, elecstate::ElecStateLCAO* pelec, diff --git a/source/source_pw/module_pwdft/setup_pwrho.h b/source/source_pw/module_pwdft/setup_pwrho.h index 820edb9d2d..4e11848d62 100644 --- a/source/source_pw/module_pwdft/setup_pwrho.h +++ b/source/source_pw/module_pwdft/setup_pwrho.h @@ -1,7 +1,6 @@ #ifndef SETUP_PWRHO_H #define SETUP_PWRHO_H -#include "source_base/module_device/device.h" // use Device #include "source_cell/unitcell.h" // use UnitCell #include "source_pw/module_pwdft/structure_factor.h" // use Structure_Factor #include "source_basis/module_pw/pw_basis.h" // use PW_Basis