diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp
index 76e9d49ba3..5e9da976d0 100644
--- a/source/module_esolver/esolver_ks.cpp
+++ b/source/module_esolver/esolver_ks.cpp
@@ -609,7 +609,7 @@ void ESolver_KS<T, Device>::runner(const int istep, UnitCell& ucell)
         this->update_pot(istep, iter);
 
         // 10) finish scf iterations
-        this->iter_finish(iter);
+        this->iter_finish(istep, iter);
 #ifdef __MPI
         double duration = (double)(MPI_Wtime() - iterstart);
 #else
@@ -663,7 +663,7 @@ void ESolver_KS<T, Device>::runner(const int istep, UnitCell& ucell)
 };
 
 template <typename T, typename Device>
-void ESolver_KS<T, Device>::iter_finish(int& iter)
+void ESolver_KS<T, Device>::iter_finish(const int istep, int& iter)
 {
     // 1 means Harris-Foulkes functional
     // 2 means Kohn-Sham functional
diff --git a/source/module_esolver/esolver_ks.h b/source/module_esolver/esolver_ks.h
index 26a2535885..5cc6ccbd55 100644
--- a/source/module_esolver/esolver_ks.h
+++ b/source/module_esolver/esolver_ks.h
@@ -58,7 +58,7 @@ class ESolver_KS : public ESolver_FP
 		virtual void iter_init(const int istep, const int iter) {};
 
 		//! Something to do after hamilt2density function in each iter loop.
-        virtual void iter_finish(int& iter);
+        virtual void iter_finish(const int istep, int& iter);
 
         //! Something to do after SCF iterations when SCF is converged or comes to the max iter step.
         virtual void after_scf(const int istep) override;
diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp
index f028f9f601..b12fda5916 100644
--- a/source/module_esolver/esolver_ks_lcao.cpp
+++ b/source/module_esolver/esolver_ks_lcao.cpp
@@ -625,13 +625,15 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(const int istep, const int iter)
     // calculate exact-exchange
     if (GlobalC::exx_info.info_ri.real_number)
     {
-        this->exd->exx_eachiterinit(*dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM(),
+        this->exd->exx_eachiterinit(istep,
+                                    *dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM(),
                                     this->kv,
                                     iter);
     }
     else
     {
-        this->exc->exx_eachiterinit(*dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM(),
+        this->exc->exx_eachiterinit(istep,
+                                    *dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM(),
                                     this->kv,
                                     iter);
     }
@@ -914,12 +916,12 @@ void ESolver_KS_LCAO<TK, TR>::update_pot(const int istep, const int iter)
 //! 5) cal_MW? (why put it here?)
 //------------------------------------------------------------------------------
 template <typename TK, typename TR>
-void ESolver_KS_LCAO<TK, TR>::iter_finish(int& iter)
+void ESolver_KS_LCAO<TK, TR>::iter_finish(const int istep, int& iter)
 {
     ModuleBase::TITLE("ESolver_KS_LCAO", "iter_finish");
 
     // call iter_finish() of ESolver_KS
-    ESolver_KS<TK>::iter_finish(iter);
+    ESolver_KS<TK>::iter_finish(istep, iter);
 
     // 1) mix density matrix if mixing_restart + mixing_dmr + not first
     // mixing_restart at every iter
@@ -944,8 +946,8 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(int& iter)
     if (GlobalC::exx_info.info_global.cal_exx)
     {
         GlobalC::exx_info.info_ri.real_number ?
-            this->exd->exx_iter_finish(this->kv, GlobalC::ucell, *this->p_hamilt, *this->pelec, *this->p_chgmix, this->scf_ene_thr, iter, this->conv_esolver) :
-            this->exc->exx_iter_finish(this->kv, GlobalC::ucell, *this->p_hamilt, *this->pelec, *this->p_chgmix, this->scf_ene_thr, iter, this->conv_esolver);
+            this->exd->exx_iter_finish(this->kv, GlobalC::ucell, *this->p_hamilt, *this->pelec, *this->p_chgmix, this->scf_ene_thr, iter, istep, this->conv_esolver) :
+            this->exc->exx_iter_finish(this->kv, GlobalC::ucell, *this->p_hamilt, *this->pelec, *this->p_chgmix, this->scf_ene_thr, iter, istep, this->conv_esolver);
     }
 #endif
 
diff --git a/source/module_esolver/esolver_ks_lcao.h b/source/module_esolver/esolver_ks_lcao.h
index 4d0cc10514..5ebdd79548 100644
--- a/source/module_esolver/esolver_ks_lcao.h
+++ b/source/module_esolver/esolver_ks_lcao.h
@@ -56,7 +56,7 @@ class ESolver_KS_LCAO : public ESolver_KS<TK> {
 
     virtual void update_pot(const int istep, const int iter) override;
 
-    virtual void iter_finish(int& iter) override;
+    virtual void iter_finish(const int istep, int& iter) override;
 
     virtual void after_scf(const int istep) override;
 
diff --git a/source/module_esolver/esolver_ks_lcaopw.cpp b/source/module_esolver/esolver_ks_lcaopw.cpp
index ec95ea8761..0d868dfc34 100644
--- a/source/module_esolver/esolver_ks_lcaopw.cpp
+++ b/source/module_esolver/esolver_ks_lcaopw.cpp
@@ -191,9 +191,9 @@ namespace ModuleESolver
     }
 
     template <typename T>
-    void ESolver_KS_LIP<T>::iter_finish(int& iter)
+    void ESolver_KS_LIP<T>::iter_finish(const int istep, int& iter)
     {
-        ESolver_KS_PW<T>::iter_finish(iter);
+        ESolver_KS_PW<T>::iter_finish(istep, iter);
 
 #ifdef __EXX
         if (GlobalC::exx_info.info_global.cal_exx && this->conv_esolver)
diff --git a/source/module_esolver/esolver_ks_lcaopw.h b/source/module_esolver/esolver_ks_lcaopw.h
index e5202aff46..c133791dcb 100644
--- a/source/module_esolver/esolver_ks_lcaopw.h
+++ b/source/module_esolver/esolver_ks_lcaopw.h
@@ -28,7 +28,7 @@ namespace ModuleESolver
 
     protected:
       virtual void iter_init(const int istep, const int iter) override;
-      virtual void iter_finish(int& iter) override;
+      virtual void iter_finish(const int istep, int& iter) override;
 
       virtual void allocate_hamilt() override;
       virtual void deallocate_hamilt() override;
diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp
index a245a44358..17ec60b6bb 100644
--- a/source/module_esolver/esolver_ks_pw.cpp
+++ b/source/module_esolver/esolver_ks_pw.cpp
@@ -434,10 +434,10 @@ void ESolver_KS_PW<T, Device>::update_pot(const int istep, const int iter)
 }
 
 template <typename T, typename Device>
-void ESolver_KS_PW<T, Device>::iter_finish(int& iter)
+void ESolver_KS_PW<T, Device>::iter_finish(const int istep, int& iter)
 {
     // call iter_finish() of ESolver_KS
-    ESolver_KS<T, Device>::iter_finish(iter);
+    ESolver_KS<T, Device>::iter_finish(istep, iter);
 
     // liuyu 2023-10-24
     // D in uspp need vloc, thus needs update when veff updated
diff --git a/source/module_esolver/esolver_ks_pw.h b/source/module_esolver/esolver_ks_pw.h
index f2e65a4ef3..dc586a67a5 100644
--- a/source/module_esolver/esolver_ks_pw.h
+++ b/source/module_esolver/esolver_ks_pw.h
@@ -46,7 +46,7 @@ class ESolver_KS_PW : public ESolver_KS<T, Device>
 
     virtual void update_pot(const int istep, const int iter) override;
 
-    virtual void iter_finish(int& iter) override;
+    virtual void iter_finish(const int istep, int& iter) override;
 
     virtual void after_scf(const int istep) override;
 
diff --git a/source/module_esolver/esolver_sdft_pw.cpp b/source/module_esolver/esolver_sdft_pw.cpp
index d7e5a5b317..2e21355208 100644
--- a/source/module_esolver/esolver_sdft_pw.cpp
+++ b/source/module_esolver/esolver_sdft_pw.cpp
@@ -153,10 +153,10 @@ void ESolver_SDFT_PW<T, Device>::before_scf(const int istep)
 }
 
 template <typename T, typename Device>
-void ESolver_SDFT_PW<T, Device>::iter_finish(int& iter)
+void ESolver_SDFT_PW<T, Device>::iter_finish(const int istep, int& iter)
 {
     // call iter_finish() of ESolver_KS
-    ESolver_KS<T, Device>::iter_finish(iter);
+    ESolver_KS<T, Device>::iter_finish(istep, iter);
 }
 
 template <typename T, typename Device>
diff --git a/source/module_esolver/esolver_sdft_pw.h b/source/module_esolver/esolver_sdft_pw.h
index 8f544febe6..71e5670c44 100644
--- a/source/module_esolver/esolver_sdft_pw.h
+++ b/source/module_esolver/esolver_sdft_pw.h
@@ -39,7 +39,7 @@ class ESolver_SDFT_PW : public ESolver_KS_PW<T, Device>
 
     virtual void others(const int istep) override;
 
-    virtual void iter_finish(int& iter) override;
+    virtual void iter_finish(const int istep, int& iter) override;
 
     virtual void after_scf(const int istep) override;
 
diff --git a/source/module_esolver/lcao_before_scf.cpp b/source/module_esolver/lcao_before_scf.cpp
index b50722c883..2062f5afda 100644
--- a/source/module_esolver/lcao_before_scf.cpp
+++ b/source/module_esolver/lcao_before_scf.cpp
@@ -115,6 +115,7 @@ void ESolver_KS_LCAO<TK, TR>::beforesolver(const int istep)
             orb_,
             DM
 #ifdef __EXX
+            , istep
             , GlobalC::exx_info.info_ri.real_number ? &this->exd->two_level_step : &this->exc->two_level_step
             , GlobalC::exx_info.info_ri.real_number ? &exx_lri_double->Hexxs : nullptr
             , GlobalC::exx_info.info_ri.real_number ? nullptr : &exx_lri_complex->Hexxs
@@ -210,11 +211,11 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(const int istep)
 #ifdef __EXX // set xc type before the first cal of xc in pelec->init_scf
     if (GlobalC::exx_info.info_ri.real_number)
     {
-        this->exd->exx_beforescf(this->kv, *this->p_chgmix, GlobalC::ucell, orb_);
+        this->exd->exx_beforescf(istep, this->kv, *this->p_chgmix, GlobalC::ucell, orb_);
     }
     else
     {
-        this->exc->exx_beforescf(this->kv, *this->p_chgmix, GlobalC::ucell, orb_);
+        this->exc->exx_beforescf(istep, this->kv, *this->p_chgmix, GlobalC::ucell, orb_);
     }
 #endif // __EXX
 
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp
index bff12cb2ca..d03fe2ec86 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp
@@ -73,6 +73,7 @@ HamiltLCAO<TK, TR>::HamiltLCAO(Gint_Gamma* GG_in,
     const LCAO_Orbitals& orb,
     elecstate::DensityMatrix<TK, double>* DM_in
 #ifdef __EXX
+    , const int istep
     , int* exx_two_level_step
     , std::vector<std::map<int, std::map<TAC, RI::Tensor<double>>>>* Hexxd
     , std::vector<std::map<int, std::map<TAC, RI::Tensor<std::complex<double>>>>>* Hexxc
@@ -397,6 +398,7 @@ HamiltLCAO<TK, TR>::HamiltLCAO(Gint_Gamma* GG_in,
             Hexxd,
             Hexxc,
             Add_Hexx_Type::R,
+            istep,
             exx_two_level_step,
             !GlobalC::restart.info_load.restart_exx
             && GlobalC::restart.info_load.load_H);
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h
index c11d419c3d..0cd597b5b0 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h
@@ -37,6 +37,7 @@ class HamiltLCAO : public Hamilt<TK>
           const LCAO_Orbitals& orb,
           elecstate::DensityMatrix<TK, double>* DM_in
 #ifdef __EXX
+          , const int istep
           , int* exx_two_level_step = nullptr
           , std::vector<std::map<int, std::map<TAC, RI::Tensor<double>>>>* Hexxd = nullptr
           , std::vector<std::map<int, std::map<TAC, RI::Tensor<std::complex<double>>>>>* Hexxc = nullptr
@@ -105,6 +106,8 @@ class HamiltLCAO : public Hamilt<TK>
     /// current_spin for NSPIN=2, 0: hamiltonian for spin up, 1: hamiltonian for spin down
     int current_spin = 0;
 
+    const int istep = 0;
+
     // sk and hk will be refactored to HamiltLCAO later
     // std::vector<TK> sk;
     // std::vector<TK> hk;
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.h
index afe27a7733..9e1653be1a 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.h
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.h
@@ -32,6 +32,7 @@ class OperatorEXX<OperatorLCAO<TK, TR>> : public OperatorLCAO<TK, TR>
         std::vector<std::map<int, std::map<TAC, RI::Tensor<double>>>>* Hexxd_in = nullptr,
         std::vector<std::map<int, std::map<TAC, RI::Tensor<std::complex<double>>>>>* Hexxc_in = nullptr,
         Add_Hexx_Type add_hexx_type_in = Add_Hexx_Type::R,
+        const int istep_in = 0,
         int* two_level_step_in = nullptr,
         const bool restart_in = false);
 
@@ -54,6 +55,8 @@ class OperatorEXX<OperatorLCAO<TK, TR>> : public OperatorLCAO<TK, TR>
       /// @brief if restart, read and save Hexx, and directly use it during the first outer loop.
       bool restart = false;
 
+      const int istep = 0; // the ion step
+
       void add_loaded_Hexx(const int ik);
 
       const K_Vectors& kv;
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.hpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.hpp
index da72d37c03..c4c01f9cde 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.hpp
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.hpp
@@ -85,13 +85,15 @@ OperatorEXX<OperatorLCAO<TK, TR>>::OperatorEXX(HS_Matrix_K<TK>* hsk_in,
 	std::vector<std::map<int, std::map<TAC, RI::Tensor<double>>>>* Hexxd_in,
 	std::vector<std::map<int, std::map<TAC, RI::Tensor<std::complex<double>>>>>* Hexxc_in,
     Add_Hexx_Type add_hexx_type_in,
-	int* two_level_step_in,
+    const int istep,
+    int* two_level_step_in,
 	const bool restart_in)
     : OperatorLCAO<TK, TR>(hsk_in, kv_in.kvec_d, hR_in),
     kv(kv_in),
     Hexxd(Hexxd_in),
     Hexxc(Hexxc_in),
     add_hexx_type(add_hexx_type_in),
+    istep(istep),
     two_level_step(two_level_step_in),
     restart(restart_in)
 {
@@ -201,7 +203,7 @@ void OperatorEXX<OperatorLCAO<TK, TR>>::contributeHR()
 {
     ModuleBase::TITLE("OperatorEXX", "contributeHR");
     // Peize Lin add 2016-12-03
-    if (PARAM.inp.calculation != "nscf" && this->two_level_step != nullptr && *this->two_level_step == 0 && !this->restart) { return; }  //in the non-exx loop, do nothing 
+    if (this->istep == 0 && PARAM.inp.calculation != "nscf" && this->two_level_step != nullptr && *this->two_level_step == 0 && !this->restart) { return; }  //in the non-exx loop, do nothing 
     if (this->add_hexx_type == Add_Hexx_Type::k) { return; }
 
     if (XC_Functional::get_func_type() == 4 || XC_Functional::get_func_type() == 5)
diff --git a/source/module_ri/Exx_LRI_interface.h b/source/module_ri/Exx_LRI_interface.h
index 902c105b04..b3d358a925 100644
--- a/source/module_ri/Exx_LRI_interface.h
+++ b/source/module_ri/Exx_LRI_interface.h
@@ -54,10 +54,10 @@ class Exx_LRI_Interface
     void exx_before_all_runners(const K_Vectors& kv, const UnitCell& ucell, const Parallel_2D& pv);
 
     /// @brief in beforescf: set xc type, opt_orb, do DM mixing
-    void exx_beforescf(const K_Vectors& kv, const Charge_Mixing& chgmix, const UnitCell& ucell, const LCAO_Orbitals& orb);
+    void exx_beforescf(const int istep, const K_Vectors& kv, const Charge_Mixing& chgmix, const UnitCell& ucell, const LCAO_Orbitals& orb);
 
     /// @brief in eachiterinit:  do DM mixing and calculate Hexx when entering 2nd SCF
-    void exx_eachiterinit(const elecstate::DensityMatrix<T, double>& dm/**< double should be Tdata if complex-PBE-DM is supported*/,
+    void exx_eachiterinit(const int istep, const elecstate::DensityMatrix<T, double>& dm/**< double should be Tdata if complex-PBE-DM is supported*/,
         const K_Vectors& kv, const int& iter);
 
     /// @brief in hamilt2density: calculate Hexx and Eexx
@@ -66,7 +66,7 @@ class Exx_LRI_Interface
     /// @brief in iter_finish: write Hexx, do something according to whether SCF is converged
     void exx_iter_finish(const K_Vectors& kv, const UnitCell& ucell,
         hamilt::Hamilt<T>& hamilt, elecstate::ElecState& elec, Charge_Mixing& chgmix,
-        const double& scf_ene_thr, int& iter, bool& conv_esolver);
+        const double& scf_ene_thr, int& iter, const int istep, bool& conv_esolver);
     /// @brief: in do_after_converge: add exx operators; do DM mixing if seperate loop
     bool exx_after_converge(
         hamilt::Hamilt<T>& hamilt,
@@ -74,10 +74,12 @@ class Exx_LRI_Interface
         const K_Vectors& kv,
         const int& nspin,
         int& iter,
+        const int& istep,
         const double& etot,
         const double& scf_ene_thr);
     int two_level_step = 0;
     double etot_last_outer_loop = 0.0;
+    elecstate::DensityMatrix<T, double>* dm_last_step;
 private:
     std::shared_ptr<Exx_LRI<Tdata>> exx_ptr;
     Mix_DMk_2D mix_DMk_2D;
diff --git a/source/module_ri/Exx_LRI_interface.hpp b/source/module_ri/Exx_LRI_interface.hpp
index ec0cf1fa04..5f7ba2b964 100644
--- a/source/module_ri/Exx_LRI_interface.hpp
+++ b/source/module_ri/Exx_LRI_interface.hpp
@@ -53,13 +53,16 @@ void Exx_LRI_Interface<T, Tdata>::exx_before_all_runners(const K_Vectors& kv, co
 }
 
 template<typename T, typename Tdata>
-void Exx_LRI_Interface<T, Tdata>::exx_beforescf(const K_Vectors& kv, const Charge_Mixing& chgmix, const UnitCell& ucell, const LCAO_Orbitals& orb)
+void Exx_LRI_Interface<T, Tdata>::exx_beforescf(const int istep, const K_Vectors& kv, const Charge_Mixing& chgmix, const UnitCell& ucell, const LCAO_Orbitals& orb)
 {
 #ifdef __MPI
     if (GlobalC::exx_info.info_global.cal_exx)
     {
         if (GlobalC::restart.info_load.load_H_finish && !GlobalC::restart.info_load.restart_exx) { XC_Functional::set_xc_type(GlobalC::ucell.atoms[0].ncpp.xc_func);
-        } else
+        } 
+        else if (istep > 0) { XC_Functional::set_xc_type(GlobalC::ucell.atoms[0].ncpp.xc_func);
+        } 
+        else
         {
             if (ucell.atoms[0].ncpp.xc_func == "HF" || ucell.atoms[0].ncpp.xc_func == "PBE0" || ucell.atoms[0].ncpp.xc_func == "HSE")
             {
@@ -101,21 +104,28 @@ void Exx_LRI_Interface<T, Tdata>::exx_beforescf(const K_Vectors& kv, const Charg
 }
 
 template<typename T, typename Tdata>
-void Exx_LRI_Interface<T, Tdata>::exx_eachiterinit(const elecstate::DensityMatrix<T, double>& dm, const K_Vectors& kv, const int& iter)
+void Exx_LRI_Interface<T, Tdata>::exx_eachiterinit(const int istep, const elecstate::DensityMatrix<T, double>& dm, const K_Vectors& kv, const int& iter)
 {
     if (GlobalC::exx_info.info_global.cal_exx)
     {
-        if (!GlobalC::exx_info.info_global.separate_loop && this->two_level_step)
+        if (!GlobalC::exx_info.info_global.separate_loop && (this->two_level_step || istep > 0))
         {
             const bool flag_restart = (iter == 1) ? true : false;
-            if (this->exx_spacegroup_symmetry) { this->mix_DMk_2D.mix(symrot_.restore_dm(kv, dm.get_DMK_vector(), *dm.get_paraV_pointer()), flag_restart); }
-            else { this->mix_DMk_2D.mix(dm.get_DMK_vector(), flag_restart); }
-			const std::vector<std::map<int,std::map<std::pair<int, std::array<int, 3>>,RI::Tensor<Tdata>>>>
-				Ds = PARAM.globalv.gamma_only_local
-                ? RI_2D_Comm::split_m2D_ktoR<Tdata>(*this->exx_ptr->p_kv, this->mix_DMk_2D.get_DMk_gamma_out(), *dm.get_paraV_pointer(), PARAM.inp.nspin)
-                : RI_2D_Comm::split_m2D_ktoR<Tdata>(*this->exx_ptr->p_kv, this->mix_DMk_2D.get_DMk_k_out(), *dm.get_paraV_pointer(), PARAM.inp.nspin, this->exx_spacegroup_symmetry);
-            if (this->exx_spacegroup_symmetry && GlobalC::exx_info.info_global.exx_symmetry_realspace) { this->exx_ptr->cal_exx_elec(Ds, *dm.get_paraV_pointer(), &this->symrot_); }
-            else { this->exx_ptr->cal_exx_elec(Ds, *dm.get_paraV_pointer()); }
+            auto cal = [this, &kv, &flag_restart](const elecstate::DensityMatrix<T, double>& dm_in)
+            {
+                if (this->exx_spacegroup_symmetry) { this->mix_DMk_2D.mix(symrot_.restore_dm(kv,dm_in.get_DMK_vector(), *dm_in.get_paraV_pointer()), flag_restart); }
+                else { this->mix_DMk_2D.mix(dm_in.get_DMK_vector(), flag_restart); }
+			        const std::vector<std::map<int,std::map<std::pair<int, std::array<int, 3>>,RI::Tensor<Tdata>>>>
+				            Ds = PARAM.globalv.gamma_only_local
+                                ? RI_2D_Comm::split_m2D_ktoR<Tdata>(*this->exx_ptr->p_kv, this->mix_DMk_2D.get_DMk_gamma_out(), *dm_in.get_paraV_pointer(), PARAM.inp.nspin)
+                                : RI_2D_Comm::split_m2D_ktoR<Tdata>(*this->exx_ptr->p_kv, this->mix_DMk_2D.get_DMk_k_out(), *dm_in.get_paraV_pointer(), PARAM.inp.nspin, this->exx_spacegroup_symmetry);
+                if (this->exx_spacegroup_symmetry && GlobalC::exx_info.info_global.exx_symmetry_realspace) { this->exx_ptr->cal_exx_elec(Ds, *dm_in.get_paraV_pointer(), &this->symrot_); }
+                else { this->exx_ptr->cal_exx_elec(Ds, *dm_in.get_paraV_pointer()); }
+            };
+            if(istep > 0 && flag_restart)
+                cal(*dm_last_step);
+            else
+                cal(dm);
         }
     }
 }
@@ -147,9 +157,9 @@ void Exx_LRI_Interface<T, Tdata>::exx_hamilt2density(elecstate::ElecState& elec,
 template<typename T, typename Tdata>
 void Exx_LRI_Interface<T, Tdata>::exx_iter_finish(const K_Vectors& kv, const UnitCell& ucell,
     hamilt::Hamilt<T>& hamilt, elecstate::ElecState& elec, Charge_Mixing& chgmix,
-    const double& scf_ene_thr, int& iter, bool& conv_esolver)
+    const double& scf_ene_thr, int& iter, const int istep, bool& conv_esolver)
 {
-    if (GlobalC::restart.info_save.save_H && this->two_level_step > 0
+    if (GlobalC::restart.info_save.save_H && (this->two_level_step > 0 || istep > 0)
         && (!GlobalC::exx_info.info_global.separate_loop || iter == 1)) // to avoid saving the same value repeatedly
     {
         ////////// for Add_Hexx_Type::k
@@ -189,12 +199,14 @@ void Exx_LRI_Interface<T, Tdata>::exx_iter_finish(const K_Vectors& kv, const Uni
         {
             chgmix.close_kerker_gg0();
         }
+        this->dm_last_step = dynamic_cast<const elecstate::ElecStateLCAO<T>*>(&elec)->get_DM();
         conv_esolver = this->exx_after_converge(
             hamilt,
             *dynamic_cast<const elecstate::ElecStateLCAO<T>*>(&elec)->get_DM(),
             kv,
             PARAM.inp.nspin,
             iter,
+            istep,
             elec.f_en.etot,
             scf_ene_thr);
     }
@@ -208,6 +220,7 @@ bool Exx_LRI_Interface<T, Tdata>::exx_after_converge(
     const K_Vectors& kv,
     const int& nspin,
     int& iter,
+    const int& istep,
     const double& etot,
     const double& scf_ene_thr)
 {   // only called if (GlobalC::exx_info.info_global.cal_exx)
@@ -216,7 +229,7 @@ bool Exx_LRI_Interface<T, Tdata>::exx_after_converge(
             GlobalC::restart.info_load.restart_exx = true;
             this->exx_ptr->Eexx = 0;
         };
-
+        
         // no separate_loop case
         if (!GlobalC::exx_info.info_global.separate_loop)
         {
@@ -226,7 +239,7 @@ bool Exx_LRI_Interface<T, Tdata>::exx_after_converge(
             // in first scf loop, exx updated once in beginning,
             // in second scf loop, exx updated every iter
 
-            if (this->two_level_step)
+            if (this->two_level_step || istep > 0)
             {
                 restart_reset();
                 return true;
diff --git a/tests/integrate/285_NO_KP_RE_HSE/result.ref b/tests/integrate/285_NO_KP_RE_HSE/result.ref
index e09444967f..5b613afba2 100644
--- a/tests/integrate/285_NO_KP_RE_HSE/result.ref
+++ b/tests/integrate/285_NO_KP_RE_HSE/result.ref
@@ -1,3 +1,3 @@
-etotref -200.8969758707964
-etotperatomref -100.4484879354
+etotref -199.35070744
+etotperatomref -99.67535372
 totaltimeref 19.17
diff --git a/tests/integrate/286_NO_KP_CR_HSE/result.ref b/tests/integrate/286_NO_KP_CR_HSE/result.ref
index c938417ee6..de76ea26c2 100644
--- a/tests/integrate/286_NO_KP_CR_HSE/result.ref
+++ b/tests/integrate/286_NO_KP_CR_HSE/result.ref
@@ -1,4 +1,4 @@
-etotref -196.7312320266655092
-etotperatomref -98.3656160133
+etotref -192.82622412
+etotperatomref -96.41311206
 totaltimeref 56.4969
 37.12
diff --git a/tests/integrate/385_NO_GO_RE_S1_HSE/result.ref b/tests/integrate/385_NO_GO_RE_S1_HSE/result.ref
index ddb1a478be..356e231971 100644
--- a/tests/integrate/385_NO_GO_RE_S1_HSE/result.ref
+++ b/tests/integrate/385_NO_GO_RE_S1_HSE/result.ref
@@ -1,3 +1,3 @@
-etotref -465.7333468413376
-etotperatomref -155.2444489471
+etotref -465.73367480
+etotperatomref -155.24455827
 totaltimeref 16.44
diff --git a/tests/integrate/386_NO_GO_MD_S1_HSE/result.ref b/tests/integrate/386_NO_GO_MD_S1_HSE/result.ref
index 953aacd308..fb27bcbb5a 100644
--- a/tests/integrate/386_NO_GO_MD_S1_HSE/result.ref
+++ b/tests/integrate/386_NO_GO_MD_S1_HSE/result.ref
@@ -1,3 +1,3 @@
-etotref -465.8658203382321
-etotperatomref -155.2886067794
+etotref -465.86709810
+etotperatomref -155.28903270
 totaltimeref 17.49