fix(rt): correct hybrid gauge input density for RI exchange

source/source_esolver/esolver_ks_lcao.cpp

@@ -431,13 +431,17 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
     // calculate exact-exchange
     if (this->inp_->calculation != "nscf")
     {
+        const auto* td_info = this->get_td_info();
+        const auto* field = td_info != nullptr ? &td_info->field() : nullptr;
         if (exx_info_.info_ri.real_number)
         {
-            this->exx_nao.exd->exx_eachiterinit(istep, ucell, *this->dmat.dm, this->kv, this->pv, iter);
+            this->exx_nao.exd->exx_eachiterinit(istep, ucell, *this->dmat.dm, this->kv, this->pv, iter,
+                field);
         }
         else
         {
-            this->exx_nao.exc->exx_eachiterinit(istep, ucell, *this->dmat.dm, this->kv, this->pv, iter);
+            this->exx_nao.exc->exx_eachiterinit(istep, ucell, *this->dmat.dm, this->kv, this->pv, iter,
+                field);
         }
     }
 #endif

source/source_estate/module_pot/potential_new.h

@@ -80,6 +80,9 @@ class Potential : public PotBase
 
     PotBase* get_pot_type(const std::string& pot_type);
 
+    /// @brief Read the shared RT-TDDFT field state.
+    const TDFieldManager* td_field_manager() const { return td_field_manager_.get(); }
+
     /**
      * @brief Inject the shared RT-TDDFT field state before potential setup.
      *

source/source_hamilt/module_xc/exx_info.cpp

@@ -218,6 +218,11 @@ bool init_exx_info(Exx_Info& exx_info, const Input_para& inp)
             { exx_info.info_lip.lambda = std::stod(inp.exx_fock_lambda[0]); }
 
         exx_info.info_ri.real_number = std::stoi(inp.exx_real_number);
+        if (exx_info.info_global.cal_exx && inp.basis_type == "lcao"
+            && inp.esolver_type == "tddft" && inp.td_stype == 2 && exx_info.info_ri.real_number)
+        {
+            throw std::invalid_argument("td_stype=2 EXX requires exx_real_number=0");
+        }
         exx_info.info_ri.pca_threshold = inp.exx_pca_threshold;
         exx_info.info_ri.C_threshold = inp.exx_c_threshold;
         exx_info.info_ri.V_threshold = inp.exx_v_threshold;

source/source_io/module_parameter/read_inp_exx_dftu.cpp

@@ -444,7 +444,8 @@ void ReadInput::item_exx()
         item.category = "Exact Exchange (LCAO)";
         item.type = "String";
         item.description = R"(* True: Enforce LibRI to use double data type.
-* False: Enforce LibRI to use complex data type. Setting it to True can effectively improve the speed of self-consistent calculations with hybrid functionals.)";
+* False: Enforce LibRI to use complex data type. Setting it to True can effectively improve the speed of self-consistent calculations with hybrid functionals.
+* LCAO RT-TDDFT exact exchange with td_stype=2 requires exx_real_number=0 to represent the input density phase.)";
         item.default_value = "depends on the gamma_only option";
         item.unit = "";
         read_sync_string(input.exx_real_number);

source/source_lcao/module_ri/exx_lri_interface.h

@@ -4,6 +4,7 @@
 #include "exx_lri.h"
 #include "source_estate/module_charge/chg_mix.h"
 #include "source_lcao/module_ri/mix_dmk_2d.h"
+#include "source_lcao/module_ri/hybrid_gauge_phase.h"
 #include "source_lcao/module_ri/module_exx_symmetry/symm_rotation.h"
 #include "source_estate/module_dm/density_matrix.h" // mohan add 2025-11-04
 #include "source_hamilt/hamilt.h"
@@ -15,6 +16,7 @@ class Charge_Mixing;
 namespace elecstate
 {
     class ElecState;
+    class TDFieldManager;
 }
 namespace module_dm
 {
@@ -105,7 +107,8 @@ class Exx_LRI_Interface
                           const module_dm::DensityMatrix<T, double>& dm/**< double should be Tdata if complex-PBE-DM is supported*/,
                           const K_Vectors& kv,
                           const Parallel_Orbitals& pv,
-                          const int& iter);
+                          const int& iter,
+                          const elecstate::TDFieldManager* field);
 
     /// @brief in hamilt2rho: calculate Hexx and Eexx
     void exx_hamilt2rho(elecstate::ElecState& elec, const Parallel_Orbitals& pv, const int iter);
@@ -132,7 +135,8 @@ class Exx_LRI_Interface
                             int& iter,
                             const int& istep,
                             const double& etot,
-                            const double& scf_ene_thr);
+                            const double& scf_ene_thr,
+                            const elecstate::TDFieldManager* field);
 
     /// @brief  the step of the outer loop.
     /// nullptr: no dependence on the number of two_level_step, contributeHk will do enerything normally.
@@ -150,6 +154,17 @@ class Exx_LRI_Interface
 private:
 
     Mix_DMk_2D<T> mix_DMk_2D;
+    Mix_DMk_2D<std::complex<double>> mix_common_density_;
+    RI_2D_Comm::HybridDensityLayout common_density_layout_;
+
+    std::vector<std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>> prepare_exx_density(
+        const std::vector<std::vector<T>>& dmk, const UnitCell& ucell,
+        const K_Vectors& kv, const Parallel_Orbitals& pv, const int nspin,
+        const bool restart, const ModuleBase::Vector3<double>* source_A);
+    void mix_common_density(RI_2D_Comm::HybridGaugeComplexDensity& density,
+        const UnitCell& ucell, const ModuleBase::Vector3<double>& source_A, const bool restart);
+    void mix_common_density(std::vector<std::map<TA, std::map<TAC, RI::Tensor<double>>>>& density,
+        const UnitCell& ucell, const ModuleBase::Vector3<double>& source_A, const bool restart);
 
     Exx_Info_Global info_global;
     size_t hybrid_step_ = 1;

source/source_lcao/module_ri/exx_lri_interface.hpp

@@ -1,6 +1,7 @@
 #ifndef EXX_LRI_INTERFACE_HPP
 #define EXX_LRI_INTERFACE_HPP
 #include "exx_lri_interface.h"
+#include "source_estate/module_pot/td_field_manager.h"
 #include "source_cell/unitcell.h"
 #include "source_base/formatter.h"
 #include "source_base/parallel_common.h"
@@ -79,6 +80,64 @@ void Exx_LRI_Interface<T, Tdata>::cal_exx_stress(const double& omega, const doub
     this->flag_finish.stress = true;
 }
 
+template<typename T, typename Tdata>
+void Exx_LRI_Interface<T, Tdata>::mix_common_density(
+    RI_2D_Comm::HybridGaugeComplexDensity& density,
+    const UnitCell& ucell,
+    const ModuleBase::Vector3<double>& source_A,
+    const bool restart)
+{
+    RI_2D_Comm::mix_hybrid_density(ucell.latvec, ucell.lat0, source_A, restart,
+                                 this->mix_common_density_, this->common_density_layout_, density);
+}
+
+template<typename T, typename Tdata>
+void Exx_LRI_Interface<T, Tdata>::mix_common_density(
+    std::vector<std::map<TA, std::map<TAC, RI::Tensor<double>>>>& density,
+    const UnitCell& ucell,
+    const ModuleBase::Vector3<double>& source_A,
+    const bool restart)
+{
+    throw std::invalid_argument("sHG EXX density requires complex RI tensors (exx_real_number 0)");
+}
+
+template<typename T, typename Tdata>
+std::vector<std::map<typename Exx_LRI_Interface<T, Tdata>::TA,
+    std::map<typename Exx_LRI_Interface<T, Tdata>::TAC, RI::Tensor<Tdata>>>>
+Exx_LRI_Interface<T, Tdata>::prepare_exx_density(
+    const std::vector<std::vector<T>>& dmk,
+    const UnitCell& ucell,
+    const K_Vectors& kv,
+    const Parallel_Orbitals& pv,
+    const int nspin,
+    const bool restart,
+    const ModuleBase::Vector3<double>* source_A)
+{
+    std::vector<std::vector<T>> restored_storage;
+    const std::vector<std::vector<T>>* restored = &dmk;
+    if (this->exx_spacegroup_symmetry)
+    {
+        restored_storage = this->symrot_.restore_dm(kv, dmk, pv);
+        restored = &restored_storage;
+    }
+    if (source_A == nullptr)
+    {
+        this->mix_DMk_2D.mix(*restored, restart);
+        return RI_2D_Comm::split_m2D_ktoR<Tdata>(
+            ucell, *this->exx_ptr->p_kv, this->mix_DMk_2D.get_DMk_out(), pv, nspin,
+            this->exx_spacegroup_symmetry);
+    }
+    std::vector<const std::vector<T>*> raw;
+    for (const auto& block : *restored)
+    {
+        raw.push_back(&block);
+    }
+    auto density = RI_2D_Comm::split_m2D_ktoR<Tdata>(
+        ucell, *this->exx_ptr->p_kv, raw, pv, nspin, this->exx_spacegroup_symmetry);
+    this->mix_common_density(density, ucell, *source_A, restart);
+    return density;
+}
+
 template<typename T, typename Tdata>
 void Exx_LRI_Interface<T, Tdata>::exx_before_all_runners(
     const K_Vectors& kv,
@@ -163,9 +222,11 @@ void Exx_LRI_Interface<T, Tdata>::exx_beforescf(const int istep,
             { this->mix_DMk_2D.set_nks(kv.get_nks()); }
 
         if (this->info_global.separate_loop)
-            { this->mix_DMk_2D.set_mixing_plain(this->info_global.mixing_beta_for_loop1); }
+            { this->mix_DMk_2D.set_mixing_plain(this->info_global.mixing_beta_for_loop1);
+              this->mix_common_density_.set_mixing_plain(this->info_global.mixing_beta_for_loop1); }
         else
-            { this->mix_DMk_2D.set_mixing(chgmix.get_mixing()); }
+            { this->mix_DMk_2D.set_mixing(chgmix.get_mixing());
+              this->mix_common_density_.set_mixing(chgmix.get_mixing()); }
 
         // remember chgmix so exx_eachiterinit can re-borrow its mixing pointer after
         // mixing_restart's init_mixing() has reallocated it (else use-after-free -> SIGSEGV)
@@ -183,7 +244,8 @@ void Exx_LRI_Interface<T, Tdata>::exx_eachiterinit(const int istep,
                                                    const module_dm::DensityMatrix<T, double>& dm,
                                                    const K_Vectors& kv,
                                                    const Parallel_Orbitals& pv,
-                                                   const int& iter)
+                                                   const int& iter,
+    const elecstate::TDFieldManager* field)
 {
     ModuleBase::TITLE("Exx_LRI_Interface","exx_eachiterinit");
     if (this->info_global.cal_exx)
@@ -212,32 +274,23 @@ void Exx_LRI_Interface<T, Tdata>::exx_eachiterinit(const int istep,
                     { flag_restart = true; }   
                 this->last_borrowed_mixing_ = cur_mixing;
                 this->mix_DMk_2D.set_mixing(this->p_chgmix_->get_mixing());
+                this->mix_common_density_.set_mixing(this->p_chgmix_->get_mixing());
             }
 
-            auto cal = [this, &ucell,&kv, &pv, &flag_restart](const module_dm::DensityMatrix<T, double>& dm_in)
+            const auto density_time = RI_2D_Comm::hybrid_density_time(istep, iter);
+            const bool use_left_density = density_time == RI_2D_Comm::HybridDensityTime::Left;
+            const auto& source_density = use_left_density ? *this->dm_last_step : dm;
+            const ModuleBase::Vector3<double>* source_A = nullptr;
+            if (field != nullptr && field->gauge() == 2)
             {
-                if (this->exx_spacegroup_symmetry)
-                    { this->mix_DMk_2D.mix(symrot_.restore_dm(kv, dm_in.get_dmk_vec(), pv), flag_restart); }
-                else
-                    { this->mix_DMk_2D.mix(dm_in.get_dmk_vec(), flag_restart); }
-                const std::vector<std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>> Ds =
-                    RI_2D_Comm::split_m2D_ktoR<Tdata>(
-                        ucell,
-                        *this->exx_ptr->p_kv,
-                        this->mix_DMk_2D.get_DMk_out(),
-                        pv,
-                        PARAM.inp.nspin,
-                        this->exx_spacegroup_symmetry);
-                if(this->exx_spacegroup_symmetry && this->exx_ptr->info.exx_symmetry_realspace)
-                    { this->cal_exx_elec(Ds, ucell, pv, &this->symrot_); }
-                else
-                    { this->cal_exx_elec(Ds, ucell, pv); }
-            };
-
-            if(istep > 0 && flag_restart)
-                { cal(*this->dm_last_step); }
+                source_A = use_left_density ? &field->A_left_ha() : &field->A_right_ha();
+            }
+            const auto Ds = this->prepare_exx_density(source_density.get_dmk_vec(), ucell, kv, pv,
+                                                      PARAM.inp.nspin, flag_restart, source_A);
+            if (this->exx_spacegroup_symmetry && this->exx_ptr->info.exx_symmetry_realspace)
+                { this->cal_exx_elec(Ds, ucell, pv, &this->symrot_); }
             else
-                { cal(dm); }
+                { this->cal_exx_elec(Ds, ucell, pv); }
         }
     }
 }
@@ -341,7 +394,7 @@ void Exx_LRI_Interface<T, Tdata>::exx_iter_finish(const K_Vectors& kv,
             iter,
             istep,
             elec.f_en.etot,
-            scf_ene_thr);
+            scf_ene_thr, elec.pot->td_field_manager());
 
         // On the 0->1 transition, exx_after_converge switches the XC functional (PBE->hybrid), 
         // but v_eff was already set for this (converged) iteration under the OLD functional. 
@@ -371,7 +424,8 @@ bool Exx_LRI_Interface<T, Tdata>::exx_after_converge(
     int& iter,
     const int& istep,
     const double& etot,
-    const double& scf_ene_thr)
+    const double& scf_ene_thr,
+    const elecstate::TDFieldManager* field)
 {   // only called if (this->info_global.cal_exx)
     ModuleBase::TITLE("Exx_LRI_Interface","exx_after_converge");
     auto restart_reset = [this]()
@@ -432,18 +486,13 @@ bool Exx_LRI_Interface<T, Tdata>::exx_after_converge(
             // if init_wfc == "file", DM is calculated in the 1st iter of the 1st two-level step, so we mix it here
             const bool flag_restart = (this->two_level_step == 0 && PARAM.inp.init_wfc != "file") ? true : false;
 
-            if(this->exx_spacegroup_symmetry)
-                { this->mix_DMk_2D.mix(symrot_.restore_dm(kv, dm.get_dmk_vec(), pv), flag_restart); }
-            else
-                { this->mix_DMk_2D.mix(dm.get_dmk_vec(), flag_restart); }
-            const std::vector<std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>> Ds =
-                RI_2D_Comm::split_m2D_ktoR<Tdata>(
-                    ucell,
-                    *this->exx_ptr->p_kv,
-                    this->mix_DMk_2D.get_DMk_out(),
-                    pv,
-                    nspin,
-                    this->exx_spacegroup_symmetry);
+            const ModuleBase::Vector3<double>* source_A = nullptr;
+            if (field != nullptr && field->gauge() == 2)
+            {
+                source_A = &field->A_right_ha();
+            }
+            const auto Ds = this->prepare_exx_density(dm.get_dmk_vec(), ucell, kv, pv,
+                                                      nspin, flag_restart, source_A);
             if(this->exx_spacegroup_symmetry && this->exx_ptr->info.exx_symmetry_realspace)
                 { this->cal_exx_elec(Ds, ucell, pv, &this->symrot_); }
             else