diff --git a/source/source_esolver/esolver_ks.cpp b/source/source_esolver/esolver_ks.cpp index d9544209d3..1718f4b342 100644 --- a/source/source_esolver/esolver_ks.cpp +++ b/source/source_esolver/esolver_ks.cpp @@ -24,14 +24,13 @@ #include "source_io/json_output/output_info.h" #include "source_estate/update_pot.h" // mohan add 20251016 +#include "source_estate/module_charge/chgmixing.h" // mohan add 20251018 namespace ModuleESolver { template -ESolver_KS::ESolver_KS() -{ -} +ESolver_KS::ESolver_KS(){} template @@ -274,31 +273,15 @@ void ESolver_KS::iter_init(UnitCell& ucell, const int istep, const in if (PARAM.inp.esolver_type == "ksdft") { - diag_ethr = hsolver::set_diagethr_ks(PARAM.inp.basis_type, - PARAM.inp.esolver_type, - PARAM.inp.calculation, - PARAM.inp.init_chg, - PARAM.inp.precision, - istep, - iter, - drho, - PARAM.inp.pw_diag_thr, - diag_ethr, - PARAM.inp.nelec); + diag_ethr = hsolver::set_diagethr_ks(PARAM.inp.basis_type, PARAM.inp.esolver_type, + PARAM.inp.calculation, PARAM.inp.init_chg, PARAM.inp.precision, istep, iter, + drho, PARAM.inp.pw_diag_thr, diag_ethr, PARAM.inp.nelec); } else if (PARAM.inp.esolver_type == "sdft") { - diag_ethr = hsolver::set_diagethr_sdft(PARAM.inp.basis_type, - PARAM.inp.esolver_type, - PARAM.inp.calculation, - PARAM.inp.init_chg, - istep, - iter, - drho, - PARAM.inp.pw_diag_thr, - diag_ethr, - PARAM.inp.nbands, - esolver_KS_ne); + diag_ethr = hsolver::set_diagethr_sdft(PARAM.inp.basis_type, PARAM.inp.esolver_type, + PARAM.inp.calculation, PARAM.inp.init_chg, istep, iter, drho, + PARAM.inp.pw_diag_thr, diag_ethr, PARAM.inp.nbands, esolver_KS_ne); } // save input charge density (rho) @@ -309,9 +292,7 @@ template void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& iter, bool &conv_esolver) { - //---------------------------------------------------------------- - // 1) print out band gap - //---------------------------------------------------------------- + // 1.1) print out band gap if (!PARAM.globalv.two_fermi) { this->pelec->cal_bandgap(); @@ -321,125 +302,30 @@ void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& i this->pelec->cal_bandgap_updw(); } + // 1.2) print out eigenvalues and occupations if(iter % PARAM.inp.out_freq_elec == 0) { - //---------------------------------------------------------------- - // 2) print out eigenvalues and occupations - //---------------------------------------------------------------- if (PARAM.inp.out_band[0] || iter == PARAM.inp.scf_nmax || conv_esolver) { ModuleIO::write_eig_iter(this->pelec->ekb,this->pelec->wg,*this->pelec->klist); } } - //---------------------------------------------------------------- - // 2) compute magnetization, only for LSDA(spin==2) - //---------------------------------------------------------------- + // 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()); - //---------------------------------------------------------------- - // 3) charge mixing - //---------------------------------------------------------------- - if (PARAM.globalv.ks_run) - { - // mixing will restart at this->p_chgmix->mixing_restart steps - if (drho <= PARAM.inp.mixing_restart && PARAM.inp.mixing_restart > 0.0 - && this->p_chgmix->mixing_restart_step > iter) - { - this->p_chgmix->mixing_restart_step = iter + 1; - } - - if (PARAM.inp.scf_os_stop) // if oscillation is detected, SCF will stop - { - this->oscillate_esolver - = this->p_chgmix->if_scf_oscillate(iter, drho, PARAM.inp.scf_os_ndim, PARAM.inp.scf_os_thr); - } - - // drho will be 0 at this->p_chgmix->mixing_restart step, which is - // not ground state - bool not_restart_step = !(iter == this->p_chgmix->mixing_restart_step && PARAM.inp.mixing_restart > 0.0); - // SCF will continue if U is not converged for uramping calculation - bool is_U_converged = true; - // to avoid unnecessary dependence on dft+u, refactor is needed -#ifdef __LCAO - if (PARAM.inp.dft_plus_u) - { - is_U_converged = GlobalC::dftu.u_converged(); - } -#endif - - conv_esolver = (drho < this->scf_thr && not_restart_step && is_U_converged); - - // add energy threshold for SCF convergence - if (this->scf_ene_thr > 0.0) - { - // calculate energy of output charge density - elecstate::update_pot(ucell, this->pelec, this->chr, conv_esolver); - this->pelec->cal_energies(2); // 2 means Kohn-Sham functional - // now, etot_old is the energy of input density, while etot is the energy of output density - this->pelec->f_en.etot_delta = this->pelec->f_en.etot - this->pelec->f_en.etot_old; - // output etot_delta - GlobalV::ofs_running << " DeltaE_womix = " << this->pelec->f_en.etot_delta * ModuleBase::Ry_to_eV << " eV" - << std::endl; - if (iter > 1 && conv_esolver == 1) // only check when density is converged - { - // update the convergence flag - conv_esolver - = (std::abs(this->pelec->f_en.etot_delta * ModuleBase::Ry_to_eV) < this->scf_ene_thr); - } - } + // 2.2) charge mixing + module_charge::chgmixing_ks(iter, ucell, this->pelec, this->chr, this->p_chgmix, + this->pw_rhod->nrxx, this->drho, this->oscillate_esolver, conv_esolver, hsolver_error, + this->scf_thr, this->scf_ene_thr, PARAM.inp); - // If drho < hsolver_error in the first iter or drho < scf_thr, we - // do not change rho. - if (drho < hsolver_error || conv_esolver || PARAM.inp.calculation == "nscf") - { - if (drho < hsolver_error) - { - GlobalV::ofs_warning << " drho < hsolver_error, keep " - "charge density unchanged." - << std::endl; - } - } - else - { - //----------charge mixing--------------- - // mixing will restart after this->p_chgmix->mixing_restart - // steps - if (PARAM.inp.mixing_restart > 0 && iter == this->p_chgmix->mixing_restart_step - 1 - && drho <= PARAM.inp.mixing_restart) - { - // do not mix charge density - } - else - { - p_chgmix->mix_rho(&this->chr); // update chr->rho by mixing - } - if (PARAM.inp.scf_thr_type == 2) - { - this->chr.renormalize_rho(); // renormalize rho in R-space would - // induce a error in K-space - } - //----------charge mixing done----------- - } - } - -#ifdef __MPI - MPI_Bcast(&drho, 1, MPI_DOUBLE, 0, BP_WORLD); - - // change MPI_DOUBLE to MPI_C_BOOL, mohan 2025-04-13 - MPI_Bcast(&conv_esolver, 1, MPI_C_BOOL, 0, BP_WORLD); - MPI_Bcast(this->chr.rho[0], this->pw_rhod->nrxx, MPI_DOUBLE, 0, BP_WORLD); -#endif - - // 4) Update potentials (should be done every SF iter) + // 2.3) Update potentials (should be done every SF iter) elecstate::update_pot(ucell, this->pelec, this->chr, conv_esolver); - // 5) calculate energies - // 1 means Harris-Foulkes functional - // 2 means Kohn-Sham functional - this->pelec->cal_energies(1); - this->pelec->cal_energies(2); + // 3.1) calculate energies + this->pelec->cal_energies(1); // Harris-Foulkes functional + this->pelec->cal_energies(2); // Kohn-Sham functional if (iter == 1) { @@ -448,10 +334,18 @@ void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& i this->pelec->f_en.etot_delta = this->pelec->f_en.etot - this->pelec->f_en.etot_old; this->pelec->f_en.etot_old = this->pelec->f_en.etot; + // 4) get meta-GGA related parameters + double dkin = 0.0; // for meta-GGA + if (XC_Functional::get_ked_flag()) + { + dkin = p_chgmix->get_dkin(&this->chr, PARAM.inp.nelec); + } - //---------------------------------------------------------------- - // 6) time and meta-GGA - //---------------------------------------------------------------- + // Iter finish + ESolver_FP::iter_finish(ucell, istep, iter, conv_esolver); + + + // the end, print time #ifdef __MPI double duration = (double)(MPI_Wtime() - iter_time); #else @@ -460,38 +354,19 @@ void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& i / static_cast(1e6); #endif - // get mtaGGA related parameters - double dkin = 0.0; // for meta-GGA - if (XC_Functional::get_ked_flag()) - { - dkin = p_chgmix->get_dkin(&this->chr, PARAM.inp.nelec); - } - - // pint energy - elecstate::print_etot(ucell.magnet, *pelec,conv_esolver, iter, drho, + // print energies + elecstate::print_etot(ucell.magnet, *pelec, conv_esolver, iter, drho, dkin, duration, diag_ethr); #ifdef __RAPIDJSON - // 7) add Json of scf mag + // add Json of scf mag Json::add_output_scf_mag(ucell.magnet.tot_mag, ucell.magnet.abs_mag, this->pelec->f_en.etot * ModuleBase::Ry_to_eV, this->pelec->f_en.etot_delta * ModuleBase::Ry_to_eV, drho, duration); #endif //__RAPIDJSON - - // 7) SCF restart information - if (PARAM.inp.mixing_restart > 0 - && iter == this->p_chgmix->mixing_restart_step - 1 - && iter != PARAM.inp.scf_nmax) - { - this->p_chgmix->mixing_restart_last = iter; - std::cout << " SCF restart after this step!" << std::endl; - } - - // 8) Iter finish - ESolver_FP::iter_finish(ucell, istep, iter, conv_esolver); } //! Something to do after SCF iterations when SCF is converged or comes to the max iter step. @@ -534,13 +409,9 @@ void ESolver_KS::after_scf(UnitCell& ucell, const int istep, const bo ss << ".txt"; const double eshift = 0.0; - ModuleIO::nscf_band(is, - ss.str(), - PARAM.inp.nbands, - eshift, - PARAM.inp.out_band[1], // precision - this->pelec->ekb, - this->kv); + ModuleIO::nscf_band(is, ss.str(), PARAM.inp.nbands, + eshift, PARAM.inp.out_band[1], // precision + this->pelec->ekb, this->kv); } } } diff --git a/source/source_estate/module_charge/chgmixing.cpp b/source/source_estate/module_charge/chgmixing.cpp index 87496ced87..59467fc010 100644 --- a/source/source_estate/module_charge/chgmixing.cpp +++ b/source/source_estate/module_charge/chgmixing.cpp @@ -1,7 +1,133 @@ #include "source_estate/module_charge/chgmixing.h" +#include "source_estate/update_pot.h" #include "source_lcao/module_dftu/dftu.h" #include "source_lcao/module_deltaspin/spin_constrain.h" +void module_charge::chgmixing_ks(const int iter, // scf iteration number + UnitCell& ucell, + elecstate::ElecState* pelec, + Charge &chr, // charge density + Charge_Mixing* p_chgmix, // charge mixing class + const int nrxx, // charge density + double &drho, // charge density deviation + bool &oscillate_esolver, // whether the esolver has oscillation of charge density + bool &conv_esolver, + const double &hsolver_error, + const double &scf_thr, + const double &scf_ene_thr, + const Input_para& inp) // input parameters +{ + + if (PARAM.globalv.ks_run) + { + // mixing will restart at p_chgmix->mixing_restart steps + if (drho <= inp.mixing_restart && inp.mixing_restart > 0.0 + && p_chgmix->mixing_restart_step > iter) + { + p_chgmix->mixing_restart_step = iter + 1; + } + + if (inp.scf_os_stop) // if oscillation is detected, SCF will stop + { + oscillate_esolver = p_chgmix->if_scf_oscillate(iter, drho, + inp.scf_os_ndim, inp.scf_os_thr); + } + + // drho will be 0 at p_chgmix->mixing_restart step, which is + // not ground state + bool not_restart_step = !(iter == p_chgmix->mixing_restart_step && inp.mixing_restart > 0.0); + // SCF will continue if U is not converged for uramping calculation + bool is_U_converged = true; + // to avoid unnecessary dependence on dft+u, refactor is needed +#ifdef __LCAO + if (inp.dft_plus_u) + { + is_U_converged = GlobalC::dftu.u_converged(); + } +#endif + + conv_esolver = (drho < scf_thr && not_restart_step && is_U_converged); + + // add energy threshold for SCF convergence + if (scf_ene_thr > 0.0) + { + // calculate energy of output charge density + elecstate::update_pot(ucell, pelec, chr, conv_esolver); + pelec->cal_energies(2); // 2 means Kohn-Sham functional + // now, etot_old is the energy of input density, while etot is the energy of output density + pelec->f_en.etot_delta = pelec->f_en.etot - pelec->f_en.etot_old; + // output etot_delta + GlobalV::ofs_running << " DeltaE_womix = " << pelec->f_en.etot_delta * ModuleBase::Ry_to_eV << " eV" + << std::endl; + if (iter > 1 && conv_esolver == 1) // only check when density is converged + { + // update the convergence flag + conv_esolver + = (std::abs(pelec->f_en.etot_delta * ModuleBase::Ry_to_eV) < scf_ene_thr); + } + } + + + + // If drho < hsolver_error in the first iter or drho < scf_thr, we + // do not change rho. + if (drho < hsolver_error || conv_esolver || inp.calculation == "nscf") + { + if (drho < hsolver_error) + { + GlobalV::ofs_warning << " drho < hsolver_error, keep " + "charge density unchanged." + << std::endl; + } + } + else + { + //----------charge mixing--------------- + // mixing will restart after p_chgmix->mixing_restart + // steps + if (inp.mixing_restart > 0 && iter == p_chgmix->mixing_restart_step - 1 + && drho <= inp.mixing_restart) + { + // do not mix charge density + } + else + { + p_chgmix->mix_rho(&chr); // update chr->rho by mixing + } + if (inp.scf_thr_type == 2) + { + chr.renormalize_rho(); // renormalize rho in R-space would + // induce a error in K-space + } + //----------charge mixing done----------- + } + } + +#ifdef __MPI + MPI_Bcast(&drho, 1, MPI_DOUBLE, 0, BP_WORLD); + + // change MPI_DOUBLE to MPI_C_BOOL, mohan 2025-04-13 + MPI_Bcast(&conv_esolver, 1, MPI_C_BOOL, 0, BP_WORLD); + + assert(nrxx>=0); // mohan add 2025-10-18 + MPI_Bcast(chr.rho[0], nrxx, MPI_DOUBLE, 0, BP_WORLD); +#endif + + // mohan move the following code here, 2025-10-18 + // SCF restart information + if (PARAM.inp.mixing_restart > 0 + && iter == p_chgmix->mixing_restart_step - 1 + && iter != PARAM.inp.scf_nmax) + { + p_chgmix->mixing_restart_last = iter; + std::cout << " SCF restart after this step!" << std::endl; + } + + return; +} + + + 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 diff --git a/source/source_estate/module_charge/chgmixing.h b/source/source_estate/module_charge/chgmixing.h index 2e962ed8e5..3be8660013 100644 --- a/source/source_estate/module_charge/chgmixing.h +++ b/source/source_estate/module_charge/chgmixing.h @@ -1,12 +1,29 @@ #ifndef CHGMIXING_H #define CHGMIXING_H -#include "source_estate/module_charge/charge_mixing.h" +#include "source_estate/elecstate.h" // use pelec +#include "source_estate/module_charge/charge.h" // use chr +#include "source_estate/module_charge/charge_mixing.h" // use p_chgmix #include "source_io/module_parameter/input_parameter.h" // use Input_para +#include "source_cell/unitcell.h" namespace module_charge { +void chgmixing_ks(const int iter, // scf iteration number + UnitCell& ucell, + elecstate::ElecState* pelec, + Charge &chr, // charge density + Charge_Mixing* p_chgmix, // charge mixing class + const int nrxx, // charge density + double &drho, // charge density deviation + bool &oscillate_esolver, // whether the esolver has oscillation of charge density + bool &conv_esolver, + const double &hsolver_error, + const double &scf_thr, + const double &scf_ene_thr, + const Input_para& inp); // input parameters + void chgmixing_ks_pw(const int iter, Charge_Mixing* p_chgmix, const Input_para& inp); // input parameters