Fix the issue of inconsistent behavior between estep and md_nstep. (#7198)

* update tddft subloop

* Add files via upload

* Update result.ref

* Add files via upload

* Add files via upload

* Add files via upload

* Delete tests/05_rtTDDFT/19_NO_hyb_exx_TDDFT directory

* Update CASES_CPU.txt

* Update reference values in result.ref file

Author: ESROAMER

source/source_esolver/esolver_ks_lcao.cpp

@@ -185,7 +185,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
                 this->chr, PARAM.inp.ks_solver);
         }
     }
-    else //if not, use the DMR calculated from last step
+    else if(PARAM.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

source/source_esolver/esolver_ks_lcao_tddft.cpp

@@ -94,6 +94,14 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::runner(UnitCell& ucell, const int istep)
     // 1) before_scf (electronic iteration loops)
     //----------------------------------------------------------------
     this->before_scf(ucell, istep); // From ESolver_KS_LCAO
+    if (PARAM.inp.td_stype == 2)
+    {
+        this->dmat.dm->cal_DMR_td(ucell, TD_info::cart_At);
+    }
+    else
+    {
+        this->dmat.dm->cal_DMR();
+    }
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT SCF");
 
     // Initialize velocity operator for current calculation
@@ -151,7 +159,7 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::runner(UnitCell& ucell, const int istep)
                                         &this->sf,
                                         GlobalV::ofs_running,
                                         GlobalV::ofs_warning);
-            // need to test if correct when estep>0
+            this->exx_nao.before_scf(ucell, this->kv, this->orb_, this->p_chgmix, totstep, PARAM.inp);
             this->pelec->init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm);
 
             if (totstep <= PARAM.inp.td_tend + 1)

source/source_estate/module_dm/density_matrix.cpp

@@ -352,7 +352,11 @@ void DensityMatrix_Tools::cal_DMR_td(
     }
     ModuleBase::timer::end("DensityMatrix", "cal_DMR_td");
 }
-
+template <>
+void DensityMatrix<double, double>::cal_DMR_td(const UnitCell& ucell, const ModuleBase::Vector3<double> At, const int ik_in)
+{
+    return;
+}
 template <>
 void DensityMatrix<std::complex<double>, double>::cal_DMR_td(const UnitCell& ucell, const ModuleBase::Vector3<double> At, const int ik_in)
 {

source/source_estate/module_dm/init_dm.cpp

@@ -3,6 +3,7 @@
 #include "source_estate/elecstate_tools.h"
 #include "source_estate/cal_ux.h"
 #include "source_lcao/rho_tau_lcao.h" // mohan add 2025-11-12
+#include "source_lcao/module_rt/td_info.h"
 
 template <typename TK>
 void elecstate::init_dm(UnitCell& ucell,
@@ -22,7 +23,14 @@ void elecstate::init_dm(UnitCell& ucell,
 		elecstate::calEBand(pelec->ekb, pelec->wg, pelec->f_en);
 
 		elecstate::cal_dm_psi(dmat.dm->get_paraV_pointer(), pelec->wg, *psi, *dmat.dm);
-		dmat.dm->cal_DMR();
+		if (PARAM.inp.esolver_type!="tddft" && PARAM.inp.td_stype == 2)
+		{
+			dmat.dm->cal_DMR_td(ucell, TD_info::cart_At);
+		}
+		else
+		{
+			dmat.dm->cal_DMR();
+		}
 
         // mohan add 2025-11-12, use density matrix to calculate the charge density
         LCAO_domain::dm2rho(dmat.dm->get_DMR_vector(), PARAM.inp.nspin, &chr);

tests/05_rtTDDFT/CASES_CPU.txt

@@ -15,4 +15,4 @@
 15_NO_TRI_TDDFT
 16_NO_vel_TDDFT
 17_NO_vel_TDDFT
-18_NO_hyb_TDDFT
+18_NO_hyb_TDDFT

tests/08_EXX/14_NO_TDDFT_PBE0/INPUT

@@ -0,0 +1,50 @@
+INPUT_PARAMETERS
+ 
+
+                
+# general information
+calculation       md
+esolver_type      tddft
+md_type           nve
+md_nstep          1
+estep_per_md      5
+td_dt             0.005
+md_tfirst         0
+
+                
+# rt-TDDFT parameters
+td_vext           1 # add time-dependent external potential
+td_vext_dire      3 # direction along z
+td_stype          2 # 2: hybrid gauge
+td_ttype          3 # Heaviside type potential
+td_tstart         1 # the step electric field starts
+td_tend           10 # the step electric field ends
+td_heavi_t0       3
+td_heavi_amp      0.05
+
+                
+# print out current information
+out_current       1
+
+                
+suffix            autotest
+pseudo_dir        ../../PP_ORB
+orbital_dir       ../../PP_ORB
+basis_type        lcao
+gamma_only        0
+
+                
+# electronic structure calculations
+ecutwfc           20
+scf_nmax          50
+scf_thr           1e-6
+ks_solver         scalapack_gvx
+
+                
+# charge mixing
+mixing_type       broyden
+mixing_beta       0.7
+mixing_gg0        0.0
+
+# functional
+dft_functional        pbe0

tests/08_EXX/14_NO_TDDFT_PBE0/KPT

@@ -0,0 +1,4 @@
+K_POINTS
+0
+Gamma
+2 1 1 0 0 0

tests/08_EXX/14_NO_TDDFT_PBE0/STRU

@@ -0,0 +1,22 @@
+ATOMIC_SPECIES
+Si  28.085  Si_ONCV_PBE-1.0.upf
+
+NUMERICAL_ORBITAL
+Si_gga_6au_100Ry_2s2p1d.orb
+
+LATTICE_CONSTANT
+10.2  // add lattice constant
+
+LATTICE_VECTORS
+0.0 0.5 0.5
+0.5 0.0 0.5
+0.5 0.5 0.0
+
+ATOMIC_POSITIONS
+Cartesian //Cartesian or Direct coordinate.
+
+Si	// Element type	
+0.0	// magnetism
+2	// number of atoms
+0.00 0.00 0.00 0 0 0
+0.25 0.25 0.25 1 1 1

tests/08_EXX/14_NO_TDDFT_PBE0/current_tot.txt.ref

@@ -0,0 +1,6 @@
+1 -2.4064298950993433e-19 2.4067857486819294e-19 2.4048227372973307e-19
+2 1.3583754132510189e-08 -1.3583754049585592e-08 2.4306592013239020e-06
+3 3.8454029971142809e-08 -3.8454029713369349e-08 7.2734897464115248e-06
+4 4.3394109439952714e-08 -4.3394109112344272e-08 9.6350822407378473e-06
+5 2.6144935011100098e-08 -2.6144935077756982e-08 9.4987415747287647e-06
+6 6.1114518547755749e-10 -6.1114527641090651e-10 9.2969143096824206e-06

tests/08_EXX/14_NO_TDDFT_PBE0/result.ref

@@ -0,0 +1,4 @@
+etotref -170.67031862
+etotperatomref -85.33515931
+CompareCurrent_pass 0
+totaltimeref 16.38