esolver_ks_lcao.cpp

#include "esolver_ks_lcao.h"
#include "source_base/module_external/blacs_connector.h"
#include "source_cell/module_neighbor/sltk_atom_arrange.h"
#include "source_estate/elecstate_tools.h"
#include "source_lcao/module_deltaspin/spin_constrain.h"
#include "source_lcao/module_deltaspin/deltaspin_lcao.h"
#include "source_lcao/setup_dftu_lcao.h"
#include "source_hamilt/hs_matrix_k.h"
#include "source_estate/module_charge/symm_rho.h"
#include "source_lcao/lcao_domain.h" // need DeePKS_init
#include "source_lcao/force_stress_lcao.h"
#include "source_hamilt/module_gint/gint.h"
#include "source_estate/elecstate_lcao.h"
#include "source_lcao/hamilt_lcao.h"
#include "source_hsolver/hsolver_lcao.h"
#ifdef __EXX
#include "../source_lcao/module_ri/exx_opt_orb.h"
#endif
#include "source_lcao/module_rdmft/rdmft.h"
#include "source_estate/module_charge/chgmixing.h" // use charge mixing, mohan add 20251006
#include "source_estate/module_dm/init_dm.h" // init dm from electronic wave functions
#include "source_io/module_ctrl/ctrl_runner_lcao.h" // use ctrl_runner_lcao() 
#include "source_io/module_ctrl/ctrl_iter_lcao.h" // use ctrl_iter_lcao() 
#include "source_io/module_ctrl/ctrl_scf_lcao.h" // use ctrl_scf_lcao()
#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
{

template <typename TK, typename TR>
ESolver_KS_LCAO<TK, TR>::ESolver_KS_LCAO()
{
    this->classname = "ESolver_KS_LCAO";
    this->basisname = "LCAO";
}

template <typename TK, typename TR>
ESolver_KS_LCAO<TK, TR>::~ESolver_KS_LCAO()
{
    //****************************************************
    // do not add any codes in this deconstructor funcion
    //****************************************************
    Setup_Psi<TK>::deallocate_psi(this->psi);
}

template <typename TK, typename TR>
void ESolver_KS_LCAO<TK, TR>::before_all_runners(BaseCell& basecell, const Input_para& inp)
{
    basecell.require_kind(BaseCell::Kind::unit_cell, __FUNCTION__);
    UnitCell& ucell = static_cast<UnitCell&>(basecell);

    ModuleBase::TITLE("ESolver_KS_LCAO", "before_all_runners");
    ModuleBase::timer::start("ESolver_KS_LCAO", "before_all_runners");

    // 1) before_all_runners in ESolver_KS (includes init_general_exx_info)
    ESolver_KS::before_all_runners(ucell, inp);

    // 2) init full Exx_Info for LCAO (includes info_ri, info_opt_abfs, info_lip)
    init_exx_info(this->exx_info_, inp);

    // 3) init EXX NAO - must be after init_exx_info
    this->exx_nao.init(ucell, this->exx_info_);

    // 3) autoset nbands in ElecState before init_basis (for Psi 2d division)
    if (this->pelec == nullptr)
    {
        // TK stands for double and std::complex<double>?
        this->pelec = new elecstate::ElecStateLCAO<TK>(&(this->chr), &(this->kv),
          this->kv.get_nks(), this->pw_big);
    }

    // 4) read 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_);

    // 5) setup EXX calculations
    if (inp.calculation == "gen_opt_abfs")
    {
#ifdef __EXX
        Exx_Opt_Orb exx_opt_orb;
        exx_opt_orb.generate_matrix(exx_info_.info_opt_abfs, this->kv, ucell, this->orb_);
#else
        ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::before_all_runners", "calculation=gen_opt_abfs must compile __EXX");
#endif
        return;
    }

    LCAO_domain::set_psi_occ_dm_chg<TK>(this->kv, this->psi, this->pv, this->pelec,
      this->dmat, this->chr, inp);

    LCAO_domain::set_pot<TK>(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, this->exx_info_);

    //! if kpar is not divisible by nks, print a warning
    ModuleIO::print_kpar(this->kv.get_nks(), PARAM.globalv.kpar_lcao);

    //! init rdmft, added by jghan
    if (inp.rdmft == true)
    {
        rdmft_solver.init(this->pv, ucell,
          this->gd, this->kv, *(this->pelec), this->orb_,
          two_center_bundle_, inp.dft_functional, inp.rdmft_power_alpha, this->exx_info_);
    }

    ModuleBase::timer::end("ESolver_KS_LCAO", "before_all_runners");
    return;
}


template <typename TK, typename TR>
void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
{
    ModuleBase::TITLE("ESolver_KS_LCAO", "before_scf");
    ModuleBase::timer::start("ESolver_KS_LCAO", "before_scf");

    //! 1) call before_scf() of ESolver_KS.
    ESolver_KS::before_scf(ucell, istep);

    //! 2) find search radius
    double search_radius = atom_arrange::set_sr_NL(GlobalV::ofs_running,
      this->inp_->out_level, orb_.get_rcutmax_Phi(), ucell.infoNL->get_rcutmax_Beta(),
      PARAM.globalv.gamma_only_local);

    //! 3) use search_radius to search adj atoms
    atom_arrange::search(PARAM.globalv.search_pbc, GlobalV::ofs_running,
      this->gd, ucell, search_radius, this->inp_->test_atom_input);

    //! 4) initialize NAO basis set
    // here new is a unique pointer, which will be deleted automatically
    gint_info_.reset(
        new ModuleGint::GintInfo(
        this->pw_big->nbx, this->pw_big->nby, this->pw_big->nbz,
        this->pw_rho->nx, this->pw_rho->ny, this->pw_rho->nz,
        0, 0, this->pw_big->nbzp_start,
        this->pw_big->nbx, this->pw_big->nby, this->pw_big->nbzp,
        orb_.Phi, ucell, this->gd));
    ModuleGint::Gint::set_gint_info(gint_info_.get());

    // 7) For each atom, calculate the adjacent atoms in different cells
    // and allocate the space for H(R) and S(R).
    // If k point is used here, allocate HlocR after atom_arrange.
    this->RA.for_2d(ucell, this->gd, this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs());

    // 8) initialize the Hamiltonian operators
    // if atom moves, then delete old pointer and add a new one
    if (this->p_hamilt != nullptr)
    {
        delete this->p_hamilt;
        this->p_hamilt = nullptr;
    }
    if (this->p_hamilt == nullptr)
    {
        this->p_hamilt = new hamilt::HamiltLCAO<TK, TR>(
            ucell, this->gd, &this->pv, this->pelec->pot, this->kv,
            two_center_bundle_, orb_, this->dmat.dm, &this->dftu, this->deepks, istep, exx_nao, this->exx_info_);
    }

    // 9) for each ionic step, the overlap <phi|alpha> must be rebuilt
    // since it depends on ionic positions.
    // overlap_orb_alpha is only built when DeePKS is enabled (descriptor
    // orbitals); guard the dereference so non-DeePKS runs don't form a
    // reference from a null unique_ptr (undefined behaviour).
    if (two_center_bundle_.overlap_orb_alpha)
    {
        this->deepks.build_overlap(ucell, orb_, pv, gd, *(two_center_bundle_.overlap_orb_alpha), *this->inp_);
    }

    // 10) prepare sc calculation
    init_deltaspin_lcao<TK>(ucell, *this->inp_, &(this->pv), this->kv, this->p_hamilt, this->psi, this->dmat.dm, this->pelec);

    // 11) set xc type before the first cal of xc in pelec->init_scf, Peize Lin add 2016-12-03
    this->exx_nao.before_scf(ucell, this->kv, orb_, this->p_chgmix, istep, *this->inp_, this->exx_info_);

    // 12) initalize DM(R), which has the same size with Hamiltonian(R)
    auto* hamilt_lcao = dynamic_cast<hamilt::HamiltLCAO<TK, TR>*>(this->p_hamilt);

    if(!hamilt_lcao)
    {
        ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::before_scf","p_hamilt does not exist");
    }
    this->dmat.dm->init_DMR(*hamilt_lcao->getHR());

    // 13.1) decide the strategy for initializing DMR and HR
    if(istep == 0)//if the first scf step, readin DMR from file,
    {
        //calculate or readin the density matrix DMR
        if(this->inp_->init_chg == "dm" || this->inp_->init_chg == "dm_no_renormalize")
        {
            //! 13.1.1) init charge density from density matrix file
            LCAO_domain::init_chg_dm<TK>(PARAM.globalv.global_readin_dir, this->inp_->nspin,
                this->dmat, ucell, &(this->pv), this->pelec->charge);
        }
        if(this->inp_->init_chg == "hr")
        {
            //! 13.1.2) init charge density from Hamiltonian matrix file
            LCAO_domain::init_chg_hr<TK, TR>(PARAM.globalv.global_readin_dir, this->inp_->nspin,
                static_cast<hamilt::Hamilt<TK>*>(this->p_hamilt), ucell, &(this->pv), this->psi[0], this->pelec, *this->dmat.dm,
                this->chr, this->inp_->ks_solver);
        }
    }
    else if(this->inp_->esolver_type!="tddft")//if not, use the DMR calculated from last step
    {
        // 13.1.2) two cases are considered:
        // 1. DMK in DensityMatrix is not empty (istep > 0), then DMR is initialized by DMK
        // 2. DMK in DensityMatrix is empty (istep == 0), then DMR is initialized by zeros
        this->dmat.dm->cal_DMR();
    }
    // 13.2) init_scf, should be before_scf? mohan add 2025-03-10
    elecstate::init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric,
                          istep, PARAM.globalv.global_out_dir, *this->inp_, this->pelec);

#ifdef __MLALGO
    // 14) initialize DM2(R) of DeePKS, the DM2(R) is different from DM(R)
    this->deepks.ld.init_DMR(ucell, orb_, this->pv, this->gd);
#endif

    // 16) the electron charge density should be symmetrized,
    Symmetry_rho::symmetrize_rho(this->inp_->nspin, this->chr, this->pw_rho, ucell.symm);

    // 17) update of RDMFT, added by jghan
    if (this->inp_->rdmft == true)
    {
        rdmft_solver.update_ion(ucell, *(this->pw_rho), this->locpp.vloc, this->sf.strucFac);
    }

    ModuleBase::timer::end("ESolver_KS_LCAO", "before_scf");
    return;
}


template <typename TK, typename TR>
double ESolver_KS_LCAO<TK, TR>::cal_energy()
{
    return this->pelec->f_en.etot;
}

template <typename TK, typename TR>
void ESolver_KS_LCAO<TK, TR>::cal_force(BaseCell& basecell, ModuleBase::matrix& force)
{
    basecell.require_kind(BaseCell::Kind::unit_cell, __FUNCTION__);
    UnitCell& ucell = static_cast<UnitCell&>(basecell);

    ModuleBase::TITLE("ESolver_KS_LCAO", "cal_force");
    ModuleBase::timer::start("ESolver_KS_LCAO", "cal_force");

    Force_Stress_LCAO<TK> fsl(this->RA, ucell.nat);

    deepks.dpks_out_type = "tot";  // for deepks method

    fsl.getForceStress(ucell, this->get_vdw_result(), this->inp_->cal_force, this->inp_->cal_stress,
                       this->inp_->test_force, this->inp_->test_stress,
                       this->gd, this->pv, this->pelec, this->dmat, this->psi,
                       two_center_bundle_, orb_, force, this->scs,
                       this->locpp, this->sf, this->kv,
                       this->pw_rho, this->solvent, this->dftu, this->deepks,
                       this->exx_nao, &ucell.symm, this->exx_info_, this->inp_->td_stype,
                       static_cast<hamilt::Hamilt<TK>*>(this->p_hamilt));

    // delete RA after cal_force
    this->RA.delete_grid();

    this->have_force = true;

    ModuleBase::timer::end("ESolver_KS_LCAO", "cal_force");
}

template <typename TK, typename TR>
void ESolver_KS_LCAO<TK, TR>::cal_stress(BaseCell& basecell, ModuleBase::matrix& stress)
{
    basecell.require_kind(BaseCell::Kind::unit_cell, __FUNCTION__);
    UnitCell& ucell = static_cast<UnitCell&>(basecell);

    ModuleBase::TITLE("ESolver_KS_LCAO", "cal_stress");
    ModuleBase::timer::start("ESolver_KS_LCAO", "cal_stress");

    if (!this->have_force)
    {
        ModuleBase::matrix fcs;
        this->cal_force(ucell, fcs);
    }

    // the stress has been calculated in 'cal_force'
    stress = this->scs;
    this->have_force = false;

    ModuleBase::timer::end("ESolver_KS_LCAO", "cal_stress");
}

template <typename TK, typename TR>
void ESolver_KS_LCAO<TK, TR>::after_all_runners(BaseCell& basecell)
{
    basecell.require_kind(BaseCell::Kind::unit_cell, __FUNCTION__);
    UnitCell& ucell = static_cast<UnitCell&>(basecell);

    ModuleBase::TITLE("ESolver_KS_LCAO", "after_all_runners");
    ModuleBase::timer::start("ESolver_KS_LCAO", "after_all_runners");

    ESolver_KS::after_all_runners(ucell);

    auto* hamilt_lcao = dynamic_cast<hamilt::HamiltLCAO<TK, TR>*>(this->p_hamilt);
    if(!hamilt_lcao)
    {
	    ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::after_all_runners","p_hamilt does not exist");
    }

    ModuleIO::ctrl_runner_lcao<TK, TR>(ucell,
		    *this->inp_, this->kv, this->pelec, this->dmat, this->pv, this->Pgrid, 
		    this->gd, this->psi, this->chr, hamilt_lcao,
		    this->two_center_bundle_,
		    this->orb_, this->pw_rho, this->pw_rhod,
		    this->sf, this->locpp.vloc, this->exx_nao, this->exx_info_, this->solvent);


#ifdef __MPI
#ifdef __LCAO
    // Exit BLACS environment for LCAO calculations
    Cblacs_exit(1);
#endif
#endif

    ModuleBase::timer::end("ESolver_KS_LCAO", "after_all_runners");
}

template <typename TK, typename TR>
void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const int iter)
{
    ModuleBase::TITLE("ESolver_KS_LCAO", "iter_init");

    // call iter_init() of ESolver_KS
    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(), *this->inp_); 

    if (iter == 1)
    {
        this->gint_precision_controller_.set_mode(this->inp_->gint_precision);
        this->gint_precision_controller_.reset_for_new_scf();
        this->gint_info_->set_exec_precision(this->gint_precision_controller_.current_precision());
        if (this->inp_->gint_precision == "mix")
        {
            GlobalV::ofs_running << "\n >> Gint mixed-precision mode: starting SCF with fp32"
                                 << " (will switch to fp64 when drho is small enough)" << std::endl;
            std::cout << " >> NOTICE: Gint grid-integration starts with fp32 (mixed-precision mode)" << std::endl;
        }
        else if (this->inp_->gint_precision == "single")
        {
            GlobalV::ofs_running << "\n >> Gint single-precision mode: using fp32 throughout SCF" << std::endl;
            std::cout << " >> NOTICE: Gint grid-integration uses fp32 throughout SCF (single-precision mode)" << std::endl;
        }
    }

    // mohan update 2012-06-05
    this->pelec->f_en.deband_harris = this->pelec->cal_delta_eband(ucell);

    if (istep == 0 && this->inp_->init_wfc == "file")
	{
		int exx_two_level_step = 0;
#ifdef __EXX
		if (exx_info_.info_global.cal_exx)
		{
			// the following steps are only needed in the first outer exx loop
			exx_two_level_step
				= exx_info_.info_ri.real_number ? 
                  this->exx_nao.exd->two_level_step : this->exx_nao.exc->two_level_step;
		}
#endif
		elecstate::init_dm<TK>(ucell, this->pelec, this->dmat, this->psi, this->chr, iter, exx_two_level_step);
	}

#ifdef __EXX
    // calculate exact-exchange
    if (this->inp_->calculation != "nscf")
    {
        if (exx_info_.info_ri.real_number)
        {
            this->exx_nao.exd->exx_eachiterinit(istep, ucell, *this->dmat.dm, this->kv, iter);
        }
        else
        {
            this->exx_nao.exc->exx_eachiterinit(istep, ucell, *this->dmat.dm, this->kv, iter);
        }
    }
#endif

    init_dftu_lcao<TK>(istep, iter, this->inp_->dft_plus_u, &(this->dftu), this->dmat.dm, ucell, this->chr.rho, this->pw_rho->nrxx);

#ifdef __MLALGO
    // the density matrixes of DeePKS have been updated in each iter
    this->deepks.ld.set_hr_cal(true);

    // HR in HamiltLCAO should be recalculate
    if (this->inp_->deepks_scf)
    {
        this->p_hamilt->refresh();
    }
#endif

    if (this->inp_->vl_in_h)
    {
        // update real space Hamiltonian
        this->p_hamilt->refresh();
    }

    // save density matrix DMR for mixing
    if (this->inp_->mixing_restart > 0 && this->inp_->mixing_dmr && this->p_chgmix->mixing_restart_count > 0)
    {
        this->dmat.dm->save_DMR();
    }
}

template <typename TK, typename TR>
void ESolver_KS_LCAO<TK, TR>::hamilt2rho_single(UnitCell& ucell, int istep, int iter, double ethr)
{
    ModuleBase::TITLE("ESolver_KS_LCAO", "hamilt2rho_single");

    // 1) reset energy
    this->pelec->f_en.eband = 0.0;
    this->pelec->f_en.demet = 0.0;
    bool skip_charge = this->inp_->calculation == "nscf" ? true : false;

    // 2) run the inner lambda loop to contrain atomic moments with the DeltaSpin method
    bool skip_solve = false;
    if (this->inp_->sc_mag_switch)
    {
        spinconstrain::SpinConstrain<TK>& sc = spinconstrain::SpinConstrain<TK>::getScInstance();
        if (this->inp_->sc_lambda_strategy == "linear_scan")
        {
            sc.run_lambda_linear_scan(iter - 1);
            skip_solve = true;
        }
        else if (!sc.mag_converged() && this->drho > 0 && this->drho < this->inp_->sc_scf_thr)
        {
            sc.run_lambda_loop(iter - 1);
            sc.set_mag_converged(true);
            skip_solve = true;
        }
        else if (sc.mag_converged())
        {
            sc.run_lambda_loop(iter - 1);
            skip_solve = true;
        }
    }

    // 3) run Hsolver
    if (!skip_solve)
    {
        hsolver::HSolverLCAO<TK> hsolver_lcao_obj(&(this->pv),
                                                  this->inp_->ks_solver,
                                                  PARAM.globalv.kpar_lcao,
                                                  PARAM.globalv.nlocal,
                                                  this->inp_->nbands,
                                                  this->inp_->nelec,
                                                  this->inp_->device == "gpu");
        hsolver_lcao_obj.solve(static_cast<hamilt::Hamilt<TK>*>(this->p_hamilt), this->psi[0], this->pelec, *this->dmat.dm, 
          this->chr, this->inp_->nspin, skip_charge);
    }
    else
    {
        // Lambda loop updated the density matrix (DM) but not the real-space charge density.
        // HSolver was skipped, so we need to sync rho from DM manually.
        LCAO_domain::dm2rho(this->dmat.dm->get_DMR_vector(), this->inp_->nspin, &this->chr);
    }

    // 4) EXX
#ifdef __EXX
    if (this->inp_->calculation != "nscf")
    {
        if (exx_info_.info_ri.real_number)
        {
            this->exx_nao.exd->exx_hamilt2rho(*this->pelec, this->pv, iter);
        }
        else
        {
            this->exx_nao.exc->exx_hamilt2rho(*this->pelec, this->pv, iter);
        }
    }
#endif

    // 5) symmetrize the charge density
    Symmetry_rho::symmetrize_rho(this->inp_->nspin, this->chr, this->pw_rho, ucell.symm);

    // 6) calculate delta energy
    this->pelec->f_en.deband = this->pelec->cal_delta_eband(ucell);
}


template <typename TK, typename TR>
void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int& iter, bool& conv_esolver)
{
    ModuleBase::TITLE("ESolver_KS_LCAO", "iter_finish");

    auto* hamilt_lcao = dynamic_cast<hamilt::HamiltLCAO<TK, TR>*>(this->p_hamilt);

    if(!hamilt_lcao)
    {
        ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::iter_finish","p_hamilt does not exist");
    }

	const std::vector<std::vector<TK>>& dm_vec = this->dmat.dm->get_DMK_vector();

    // 1) calculate the local occupation number matrix and energy correction in DFT+U
    finish_dftu_lcao<TK>(iter, conv_esolver, this->inp_->dft_plus_u, this->inp_->out_chg[0], &(this->dftu), ucell, dm_vec, this->kv, this->p_chgmix->get_mixing_beta(), hamilt_lcao, PARAM.globalv.global_out_dir, this->inp_->nspin, PARAM.globalv.npol);

    // mohan add 2025-11: push DFT+U energy from Plus_U instance to ElecState.
    // Covers both dft_plus_u==1 (new method, energy accumulated by DFTU::contributeHR
    // via cal_v_of_u) and dft_plus_u==2 (old method, energy from cal_energy_correction).
    if (this->inp_->dft_plus_u)
    {
        this->pelec->set_dftu_energy(this->dftu.get_energy());
    }

    // 2) for deepks, calculate delta_e, output labels during electronic steps
    this->deepks.delta_e(ucell, this->kv, this->orb_, this->pv, this->gd, dm_vec, this->pelec->f_en, *this->inp_);

    // 3) for delta spin
    cal_mi_lcao_wrapper<TK>(iter, *this->inp_);

    // 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);
    const bool precision_switched = this->gint_precision_controller_.update_after_iteration(this->drho, this->scf_thr);
    this->gint_info_->set_exec_precision(this->gint_precision_controller_.current_precision());
    if (precision_switched)
    {
        GlobalV::ofs_running << "\n >> Gint precision switched: fp32 -> fp64 (drho = "
                             << this->drho << ")" << std::endl;
        std::cout << " >> NOTICE: Gint grid-integration precision switched from fp32 to fp64" << std::endl;
    }

    // mix density matrix if mixing_restart + mixing_dmr + not first
    // mixing_restart at every iter except the last iter
    if(iter != this->inp_->scf_nmax && !conv_esolver)
    {
        if (this->inp_->mixing_restart > 0 && this->p_chgmix->mixing_restart_count > 0 && this->inp_->mixing_dmr)
        {
            this->p_chgmix->mix_dmr(this->dmat.dm);
        }
    }

    // control the output related to the finished iteration
    ModuleIO::ctrl_iter_lcao<TK, TR>(ucell, *this->inp_, this->kv, this->pelec, *this->dmat.dm,
      this->pv, this->gd, this->psi, this->chr, this->p_chgmix, 
      hamilt_lcao, this->orb_, this->deepks, 
      this->exx_nao, this->exx_info_, iter, istep, conv_esolver, this->scf_ene_thr);
}

template <typename TK, typename TR>
void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep, const bool conv_esolver)
{
    ModuleBase::TITLE("ESolver_KS_LCAO", "after_scf");
    ModuleBase::timer::start("ESolver_KS_LCAO", "after_scf");

    auto* hamilt_lcao = dynamic_cast<hamilt::HamiltLCAO<TK, TR>*>(this->p_hamilt);

    if(!hamilt_lcao)
    {
        ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::after_scf","p_hamilt does not exist");
    }

    if (this->inp_->out_elf[0] > 0)
	{
		LCAO_domain::dm2tau(this->dmat.dm->get_DMR_vector(), this->inp_->nspin, this->pelec->charge);
	}

    //! 1) call after_scf() of ESolver_KS
    ESolver_KS::after_scf(ucell, istep, conv_esolver);

    //! 2) output of lcao every few ionic steps
    ModuleIO::ctrl_scf_lcao<TK, TR>(ucell,
            *this->inp_, this->kv, this->pelec, this->dmat.dm, this->pv,
            this->gd, this->psi, hamilt_lcao, this->dftu, this->two_center_bundle_,
            this->orb_, this->pw_wfc, this->pw_rho, this->pw_big, this->sf,
            this->pw_rhod, this->locpp.vloc, this->solvent,
            this->rdmft_solver, this->deepks, this->exx_nao, this->exx_info_,
            this->conv_esolver, this->scf_nmax_flag, istep);

    //! 3) Clean up RA, which is used to serach for adjacent atoms
    if (!this->inp_->cal_force && !this->inp_->cal_stress)
    {
        this->RA.delete_grid();
    }

    ModuleBase::timer::end("ESolver_KS_LCAO", "after_scf");
}

template class ESolver_KS_LCAO<double, double>;
template class ESolver_KS_LCAO<std::complex<double>, double>;
template class ESolver_KS_LCAO<std::complex<double>, std::complex<double>>;
} // namespace ModuleESolver