Unify RT-TDDFT EXX with the general LCAO H(R) path (#7864)

* Refactor LCAO EXX accumulation into H(R)

* Fix value in current_tot.txt.ref

Updated the third row's last value in current_tot.txt.ref.

* Update current_tot.txt.ref

* Update current_tot.txt.ref

* Fix reference values in current_tot.txt.ref

* Update current_tot.txt.ref

* Update current_tot.txt.ref

* Update current_tot.txt.ref

* Fix RT-TDDFT EXX cell mapping

Author: ESROAMER

source/source_esolver/esolver_factory.cpp

@@ -4,6 +4,7 @@
 #include "esolver_ks_pw.h"
 #include "esolver_sdft_pw.h"
 #include "source_base/module_device/device.h"
+#include "source_hamilt/module_xc/general_exx_info.h"
 #include "source_io/module_parameter/parameter.h"
 #ifdef __LCAO
 #include "esolver_dm2rho.h"
@@ -254,7 +255,13 @@ ESolver* init_esolver(const Input_para& inp)
     }
     else if (esolver_type == "ksdft_lcao_tddft")
     {
-        if (inp.nspin < 4)
+        // Hybrid RT-TDDFT stores the complete Hamiltonian in complex H(R),
+        // even for collinear spin. The final operator-chain fold then applies
+        // the same TD gauge phase to local, non-local, and EXX terms.
+        General_Exx_Info exx_info;
+        init_general_exx_info(exx_info, inp);
+        const bool use_complex_hr = inp.nspin >= 4 || exx_info.cal_exx;
+        if (!use_complex_hr)
         {
 #if ((defined __CUDA) /* || (defined __ROCM) */)
             if (inp.device == "gpu")

source/source_io/module_current/td_current_io_comm.cpp

@@ -198,37 +198,8 @@ void ModuleIO::sum_HR(const UnitCell& ucell,
 
     // init complex full_hR
     init_from_hR(hR, full_hR);
-#ifdef __EXX
-    const bool use_cell_nearest = (ModuleBase::Vector3<double>(std::fmod(kv.get_koffset(0), 1.0),
-                                                               std::fmod(kv.get_koffset(1), 1.0),
-                                                               std::fmod(kv.get_koffset(2), 1.0))
-                                       .norm()
-                                   < 1e-10);
-    RI::Cell_Nearest<int, int, 3, double, 3> cell_nearest;
-    // reallocate full_hR for BvK used in EXX
-    if (exx_info.info_global.cal_exx)
-    {
-        const std::array<int, 3> Rs_period = {kv.nmp[0], kv.nmp[1], kv.nmp[2]};
-        if (use_cell_nearest)
-        {
-            // set cell_nearest
-            std::map<int, std::array<double, 3>> atoms_pos;
-            for (int iat = 0; iat < ucell.nat; ++iat)
-            {
-                atoms_pos[iat] = RI_Util::Vector3_to_array3(ucell.atoms[ucell.iat2it[iat]].tau[ucell.iat2ia[iat]]);
-            }
-            const std::array<std::array<double, 3>, 3> latvec
-                = {RI_Util::Vector3_to_array3(ucell.a1), RI_Util::Vector3_to_array3(ucell.a2), RI_Util::Vector3_to_array3(ucell.a3)};
-            cell_nearest.init(atoms_pos, latvec, Rs_period);
-            hamilt::reallocate_hcontainer(ucell.nat, full_hR, Rs_period, &cell_nearest);
-        }
-        else
-        {
-            hamilt::reallocate_hcontainer(ucell.nat, full_hR, Rs_period);
-        }
-    }
-#endif
-    // add other hR
+    // The complete H(R) already contains exact exchange. Copy it once into
+    // full_hR; rebuilding BvK cells and adding HexxR here would double count.
     add_HR(hR, full_hR);
     // add velocity complex hR
     if (PARAM.inp.td_stype == 1)
@@ -240,36 +211,6 @@ void ModuleIO::sum_HR(const UnitCell& ucell,
         const hamilt::HContainer<std::complex<double>>* velocity_hR = TD_info::td_vel_op->get_velocity_HR_pointer();
         add_HR(velocity_hR, full_hR);
     }
-#ifdef __EXX
-    // add HexxR to complex full_hR
-    if (exx_info.info_global.cal_exx)
-    {
-        for (size_t is = 0; is != PARAM.inp.nspin; ++is)
-        {
-            if (use_cell_nearest)
-            {
-                RI_2D_Comm::add_HexxR(is,
-                                      exx_info.info_global.hybrid_alpha,
-                                      exx_nao.exc->get_Hexxs(),
-                                      pv,
-                                      PARAM.globalv.npol,
-                                      *full_hR,
-                                      &cell_nearest);
-            }
-            else
-            {
-                RI_2D_Comm::add_HexxR(is,
-                                      exx_info.info_global.hybrid_alpha,
-                                      exx_nao.exc->get_Hexxs(),
-                                      pv,
-                                      PARAM.globalv.npol,
-                                      *full_hR,
-                                      nullptr);
-            }
-        }
-    }
-#endif
-
     ModuleBase::timer::end("ModuleIO", "sum_HR");
 }
 

source/source_lcao/hamilt_lcao.cpp

@@ -422,35 +422,20 @@ HamiltLCAO<TK, TR>::HamiltLCAO(const UnitCell& ucell,
         // Peize Lin add 2016-12-03
         // set xc type before the first cal of xc in pelec->init_scf
         // and calculate Cs, Vs
-        Operator<TK>* exx;
-        if (PARAM.inp.esolver_type == "tddft")
-        {
-            exx = new OperatorEXX<OperatorLCAO<TK, TR>>(this->hsk,
-                                                        this->hR,
-                                                        ucell,
-                                                        *this->kv,
-                                                        exx_nao.exd.get(),
-                                                        exx_nao.exc.get(),
-                                                        exx_info,
-                                                        Add_Hexx_Type::k,
-                                                        istep,
-                                                        !GlobalC::restart.info_load.restart_exx
-                                                            && GlobalC::restart.info_load.load_H);
-        }
-        else
-        {
-            exx = new OperatorEXX<OperatorLCAO<TK, TR>>(this->hsk,
-                                                        this->hR,
-                                                        ucell,
-                                                        *kv,
-                                                        exx_nao.exd.get(),
-                                                        exx_nao.exc.get(),
-                                                        exx_info,
-                                                        Add_Hexx_Type::R,
-                                                        istep,
-                                                        !GlobalC::restart.info_load.restart_exx
-                                                            && GlobalC::restart.info_load.load_H);
-        }
+        // Keep exact exchange in H(R) for every workflow. For RT-TDDFT the
+        // factory selects complex H(R) when EXX is active, so the operator
+        // chain folds the complete Hamiltonian with one common TD phase.
+        Operator<TK>* exx = new OperatorEXX<OperatorLCAO<TK, TR>>(this->hsk,
+                                                                  this->hR,
+                                                                  ucell,
+                                                                  *this->kv,
+                                                                  exx_nao.exd.get(),
+                                                                  exx_nao.exc.get(),
+                                                                  exx_info,
+                                                                  Add_Hexx_Type::R,
+                                                                  istep,
+                                                                  !GlobalC::restart.info_load.restart_exx
+                                                                      && GlobalC::restart.info_load.load_H);
         this->getOperator()->add(exx);
     }
 #endif

source/source_lcao/module_operator_lcao/op_exx_lcao.cpp

@@ -9,7 +9,6 @@
 #include "source_hamilt/module_hcontainer/read_hcontainer.h"
 #include "source_lcao/module_ri/exx_lri_interface.h"
 #include "source_lcao/module_ri/ri_2d_comm.h"
-#include "source_lcao/module_rt/td_info.h"
 
 namespace hamilt
 {
@@ -312,11 +311,15 @@ OperatorEXX<OperatorLCAO<TK, TR>>::OperatorEXX(HS_Matrix_K<TK>* hsk_in,
         if (this->add_hexx_type == Add_Hexx_Type::R)
         {
             // if k points has no shift, use cell_nearest to reduce the memory cost
-            this->use_cell_nearest = (ModuleBase::Vector3<double>(std::fmod(this->kv.get_koffset(0), 1.0),
-                                                                  std::fmod(this->kv.get_koffset(1), 1.0),
-                                                                  std::fmod(this->kv.get_koffset(2), 1.0))
-                                          .norm()
-                                      < 1e-10);
+            // In the hybrid gauge, BvK-equivalent cells can carry different
+            // finite-field phases, so preserve the original cell indices.
+            const bool hybrid_gauge_rt = PARAM.inp.esolver_type == "tddft" && PARAM.inp.td_stype == 2;
+            this->use_cell_nearest = !hybrid_gauge_rt
+                                     && (ModuleBase::Vector3<double>(std::fmod(this->kv.get_koffset(0), 1.0),
+                                                                     std::fmod(this->kv.get_koffset(1), 1.0),
+                                                                     std::fmod(this->kv.get_koffset(2), 1.0))
+                                             .norm()
+                                         < 1e-10);
 
             const std::array<int, 3> Rs_period = {this->kv.nmp[0], this->kv.nmp[1], this->kv.nmp[2]};
             if (this->use_cell_nearest)
@@ -460,14 +463,23 @@ void OperatorEXX<OperatorLCAO<TK, TR>>::contributeHR()
         const int two_level_step
             = exx_info_ptr->info_ri.real_number ? this->exd->get_two_level_step() : this->exc->get_two_level_step();
 
+        // Remember that the initial GGA-only stage has completed. In
+        // RT-TDDFT (and in subsequent ionic steps) LibRI may reset
+        // two_level_step to zero; without this persistent state the HR path
+        // would incorrectly skip EXX again.
+        if (two_level_step > 0)
+        {
+            this->initial_gga_done = true;
+        }
+
         // Check if we are in the pre-convergence stage of the two-level SCF (i.e., the pure GGA loop)
         bool in_gga_pre_loop = (two_level_step == 0);
 
         // Check if a high-quality initial guess is missing (neither reading wavefunctions from a file nor restarting)
         bool lacks_good_guess = (PARAM.inp.init_wfc != "file" && !this->restart);
 
         // If in the pre-convergence loop and lacking a good initial guess, skip adding the EXX contribution
-        if (in_gga_pre_loop && lacks_good_guess)
+        if (in_gga_pre_loop && lacks_good_guess && !this->initial_gga_done)
         {
             return; // In the non-EXX loop, skip adding EXX contribution
         }
@@ -513,52 +525,25 @@ template <typename TK, typename TR>
 void OperatorEXX<OperatorLCAO<TK, TR>>::contributeHk(int ik)
 {
     ModuleBase::TITLE("OperatorEXX", "constributeHk");
-    const bool has_workflow = exx_info_ptr->info_ri.real_number ? (this->exd != nullptr) : (this->exc != nullptr);
-    int two_level_step = 0;
-    if (has_workflow)
-    {
-        two_level_step
-            = exx_info_ptr->info_ri.real_number ? this->exd->get_two_level_step() : this->exc->get_two_level_step();
-    }
-
-    // Peize Lin add 2016-12-03
 
-    // Taoni Bao add 2026-05-15
-    // In RT-TDDFT, contributeHk is used, but two_level_step is reset to 0 at each ionic step.
-    // In order to add EXX correctly in for istep > 0, this->istep == 0 is needed to avoid skipping EXX calculation.
-    // 1. For NSCF
-    if (PARAM.inp.calculation == "nscf" || !has_workflow)
+    // The main LCAO path stores EXX in H(R) and lets the final base operator
+    // fold H(R) into H(k). Keep this override only for one-shot k-space EXX
+    // operators used by write_Vxc and RDMFT.
+    if (this->add_hexx_type == Add_Hexx_Type::R)
     {
-        // Do nothing here, allow the code to proceed and calculate EXX.
+        return;
     }
-    // 2. For the first ionic step:
-    else if (this->istep == 0)
-    {
-        // If EXX is once turned on (two_level_step > 0), let OperatorEXX remember this
-        if (two_level_step > 0)
-        {
-            this->initial_gga_done = true;
-        }
 
-        // Check if we are in the pre-convergence stage of the two-level SCF (i.e., the pure GGA loop)
-        bool in_gga_pre_loop = (two_level_step == 0);
-
-        // Check if a high-quality initial guess is missing
-        bool lacks_good_guess = (!this->restart);
-
-        // If in the pre-convergence loop and lacking a good initial guess, skip adding the EXX contribution
-        // Taoni Bao add 2026-05-18, only skip EXX if initial GGA loop is not done
-        // Fix RT-TDDFT EXX missing problem in the evolution
-        if (in_gga_pre_loop && lacks_good_guess && !this->initial_gga_done)
-        {
-            return; // In the non-EXX loop, skip adding EXX contribution
-        }
+    // The restart path may still use a full EXX workflow. For one-shot
+    // operators exd/exc are null and the value remains zero.
+    int two_level_step = 0;
+    if (exx_info_ptr->info_ri.real_number && this->exd != nullptr)
+    {
+        two_level_step = this->exd->get_two_level_step();
     }
-    // 3. For subsequent ionic steps (istep > 0), add EXX normally
-
-    if (this->add_hexx_type == Add_Hexx_Type::R)
+    else if (!exx_info_ptr->info_ri.real_number && this->exc != nullptr)
     {
-        OperatorLCAO<TK, TR>::contributeHk(ik);
+        two_level_step = this->exc->get_two_level_step();
     }
 
     if (XC_Functional::get_func_type() == 4 || XC_Functional::get_func_type() == 5)
@@ -584,41 +569,25 @@ void OperatorEXX<OperatorLCAO<TK, TR>>::contributeHk(int ik)
                 }
             }
         }
-        // cal H(k) from H(R) normally
-        if (PARAM.inp.esolver_type == "tddft" && PARAM.inp.td_stype == 2)
+        if (exx_info_ptr->info_ri.real_number)
         {
-            RI_2D_Comm::add_Hexx_td(ucell,
-                                    this->kv,
-                                    ik,
-                                    exx_info_ptr->info_global.hybrid_alpha,
-                                    *this->Hexxc,
-                                    *this->hR->get_paraV(),
-                                    TD_info::td_vel_op->cart_At,
-                                    TD_info::td_vel_op->get_phase_hybrid(),
-                                    this->hsk->get_hk());
+            RI_2D_Comm::add_Hexx(ucell,
+                                 this->kv,
+                                 ik,
+                                 exx_info_ptr->info_global.hybrid_alpha,
+                                 *this->Hexxd,
+                                 *this->hR->get_paraV(),
+                                 this->hsk->get_hk());
         }
         else
         {
-            if (exx_info_ptr->info_ri.real_number)
-            {
-                RI_2D_Comm::add_Hexx(ucell,
-                                     this->kv,
-                                     ik,
-                                     exx_info_ptr->info_global.hybrid_alpha,
-                                     *this->Hexxd,
-                                     *this->hR->get_paraV(),
-                                     this->hsk->get_hk());
-            }
-            else
-            {
-                RI_2D_Comm::add_Hexx(ucell,
-                                     this->kv,
-                                     ik,
-                                     exx_info_ptr->info_global.hybrid_alpha,
-                                     *this->Hexxc,
-                                     *this->hR->get_paraV(),
-                                     this->hsk->get_hk());
-            }
+            RI_2D_Comm::add_Hexx(ucell,
+                                 this->kv,
+                                 ik,
+                                 exx_info_ptr->info_global.hybrid_alpha,
+                                 *this->Hexxc,
+                                 *this->hR->get_paraV(),
+                                 this->hsk->get_hk());
         }
     }
 }
@@ -748,4 +717,4 @@ template void reallocate_hcontainer<std::complex<double>>(
 } // namespace hamilt
 
 // End content migrated from op_exx_lcao.hpp
-#endif
+#endif

source/source_lcao/module_operator_lcao/op_exx_lcao.h

@@ -62,6 +62,8 @@ class OperatorEXX<OperatorLCAO<TK, TR>> : public OperatorLCAO<TK, TR>
         const Exx_Info* exx_info,
         Add_Hexx_Type add_hexx_type_in);
 
+    // Retained for the one-shot Add_Hexx_Type::k path used by write_Vxc/RDMFT;
+    // the main LCAO path uses contributeHR().
     virtual void contributeHk(int ik) override;
     virtual void contributeHR() override;
 

source/source_lcao/module_operator_lcao/operator_lcao.cpp

@@ -178,19 +178,12 @@ void OperatorLCAO<TK, TR>::init(const int ik_in) {
         }
         case calculation_type::lcao_exx:
         {
-            //update HR first
-            if (!this->hr_done && PARAM.inp.esolver_type != "tddft")
+            // EXX is accumulated in H(R); the last operator-chain node folds
+            // the complete H(R) into H(k), including the TD gauge phase.
+            if (!this->hr_done)
             {
                 this->contributeHR();
             }
-            else if(PARAM.inp.esolver_type == "tddft")
-            {
-                this->contributeHk(ik_in);
-            }
-
-            //update HK next
-            //in cal_type=lcao_exx, HK only need to update from one node
-            // this->contributeHk(ik_in);
 
             break;
         }