deepmodeling · mohanchen · Jul 10, 2025 · Jul 9, 2025 · Jul 9, 2025 · Jul 9, 2025
diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md
@@ -375,6 +375,8 @@
     - [out\_wannier\_unk](#out_wannier_unk)
     - [out\_wannier\_wvfn\_formatted](#out_wannier_wvfn_formatted)
   - [rt-TDDFT: Real-time time dependent density functional theory](#tddft-time-dependent-density-functional-theory)
+    - [estep\_per\_md](#estep_per_md)
+    - [td\_dt](#td_dt)
     - [td\_edm](#td_edm)
     - [td\_print\_eij](#td_print_eij)
     - [td\_propagator](#td_propagator)
@@ -1728,7 +1730,7 @@ These variables are used to control the output of properties.
 
 - **Type**: Integer
 - **Availability**: Numerical atomic orbital basis
-- **Description**: Whether to output the electronic wavefunction coefficients into files and store them in the folder `OUT.${suffix}`. The files are named as `wf{s}{spin index}{k(optional)}{k-point index}{g(optional)}{geometry index1}{_nao} + {".txt"/".dat"}`. Here, 's' refers to spin, where s1 means spin up channel while s2 means spin down channel, and 's12' refer to spinor wave functions that contains both spin channels with spin-orbital coupling or noncollinear calculations enabled. In addition, if 'gamma_only' is set to 0, then the optinoal k-point sampling index appears with the k-point index attached to the electronic wave function file names. Finally, if [out_app_flag](#out_app_flag) is set to false, the file name contains the optional 'g' index for each ionic step that may have different geometries, and if [out_app_flag](#out_app_flag) is set to true, the wave functions accumulate during ionic steps.
+- **Description**: Whether to output the electronic wavefunction coefficients into files and store them in the folder `OUT.${suffix}`. The files are named as `wf{s}{spin index}{k(optional)}{k-point index}{g(optional)}{geometry index1}{_nao} + {".txt"/".dat"}`. Here, 's' refers to spin, where s1 means spin up channel while s2 means spin down channel, and 's12' refer to spinor wave functions that contains both spin channels with spin-orbital coupling or noncollinear calculations enabled. In addition, if 'gamma_only' is set to 0, then the optinoal k-point sampling index appears with the k-point index attached to the electronic wave function file names. Finally, if [out_app_flag](#out_app_flag) is set to false, the file name contains the optional 'g' index for each ionic step that may have different geometries, and if [out_app_flag](#out_app_flag) is set to true, the wave functions accumulate during ionic steps. If the out_app_flag is set to false, a new folder named WFC will be created, and the wave function files will be saved into it.
   - 0: no output
   - 1: (txt format) 
     - gamma-only: `wfs1_nao.txt` or `wfs2_nao.txt`, ...;
@@ -3687,6 +3689,18 @@ These variables are used to control berry phase and wannier90 interface paramete
 
 ## TDDFT: time dependent density functional theory
 
+### estep_per_md
+
+- **Type**: Integer
+- **Description**: The number of electron propagation steps between two ionic steps.
+- **Default**: 1
+
+### td_dt
+
+- **Type**: Real
+- **Description**: The time step used in electron propagation. Setting td_dt will reset the md_dt value to td_dt * estep_per_md.
+- **Default**: md_dt/estep_per_md
+
 ### td_edm
 
 - **Type**: Integer
@@ -3739,6 +3753,7 @@ These variables are used to control berry phase and wannier90 interface paramete
   Type of electric field in space domain
   - 0: length gauge.
   - 1: velocity gauge.
+  - 2: hybrid gauge.
 - **Default**: 0
 
 ### td_ttype

diff --git a/source/Makefile.Objects b/source/Makefile.Objects
@@ -356,6 +356,7 @@ OBJS_HAMILT_LCAO=hamilt_lcao.o\
     overlap_new.o\
     td_ekinetic_lcao.o\
     td_nonlocal_lcao.o\
+    td_pot_hybrid.o\
     veff_lcao.o\
     meta_lcao.o\
     op_dftu_lcao.o\
@@ -624,8 +625,9 @@ OBJS_LCAO=evolve_elec.o\
       propagator_cn2.o\
       propagator_taylor.o\
       propagator_etrs.o\
-      td_velocity.o\
-      td_current.o\
+      td_folding.o\
+      td_info.o\
+      velocity_op.o\
       snap_psibeta_half_tddft.o\
       solve_propagation.o\
       upsi.o\

diff --git a/source/module_lr/lr_spectrum.h b/source/module_lr/lr_spectrum.h
@@ -5,7 +5,7 @@
 #include "source_estate/module_dm/density_matrix.h"
 #include "module_lr/utils/lr_util.h"
 #include "source_basis/module_nao/two_center_bundle.h"
-#include "source_lcao/module_rt/td_current.h"
+#include "source_lcao/module_rt/velocity_op.h"
 namespace LR
 {
     template<typename T>
@@ -52,8 +52,8 @@ namespace LR
         /// calculate the transition dipole of all states in length gauge
         void cal_transition_dipoles_length();
         /// calculate the transition dipole of state S in velocity gauge: $i(\sum_{iak}X^S_{iak}<ik|v|ak>)/\Omega_S$
-        ModuleBase::Vector3<T> cal_transition_dipole_istate_velocity_R(const int istate, const TD_current& vR);
-        ModuleBase::Vector3<T> cal_transition_dipole_istate_velocity_k(const int istate, const TD_current& vR);
+        ModuleBase::Vector3<T> cal_transition_dipole_istate_velocity_R(const int istate, const Velocity_op<std::complex<double>>& vR);
+        ModuleBase::Vector3<T> cal_transition_dipole_istate_velocity_k(const int istate, const Velocity_op<std::complex<double>>& vR);
         /// calculate the transition dipole of all states in velocity gauge
         void cal_transition_dipoles_velocity();
         double cal_mean_squared_dipole(ModuleBase::Vector3<T> dipole);

diff --git a/source/module_lr/lr_spectrum_velocity.cpp b/source/module_lr/lr_spectrum_velocity.cpp
@@ -7,17 +7,17 @@
 namespace LR
 {
     /// get the velocity matrix v(R)
-    inline TD_current get_velocity_matrix_R(const UnitCell& ucell,
+    inline Velocity_op<std::complex<double>> get_velocity_matrix_R(const UnitCell& ucell,
         const Grid_Driver& gd,
         const Parallel_Orbitals& pmat,
         const TwoCenterBundle& two_center_bundle)
     {
-        // convert the orbital object to the old class for TD_current
+        // convert the orbital object to the old class for Velocity_op
         LCAO_Orbitals orb;
         const auto& inp = PARAM.inp;
         two_center_bundle.to_LCAO_Orbitals(orb, inp.lcao_ecut, inp.lcao_dk, inp.lcao_dr, inp.lcao_rmax);
         // actually this class calculates the velocity matrix v(R) at A=0
-        TD_current vR(&ucell, &gd, &pmat, orb, two_center_bundle.overlap_orb.get());
+        Velocity_op<std::complex<double>> vR(&ucell, &gd, &pmat, orb, two_center_bundle.overlap_orb.get());
         vR.calculate_vcomm_r(); // $<\mu, 0|[Vnl, r]|\nu, R>$
         vR.calculate_grad_term();   // $<\mu, 0|\nabla|\nu, R>$
         return vR;
@@ -47,7 +47,7 @@ namespace LR
 
     /// this algorithm has bug in multi-k cases, just for test
     template<typename T>
-    ModuleBase::Vector3<T> LR::LR_Spectrum<T>::cal_transition_dipole_istate_velocity_R(const int istate, const TD_current& vR)
+    ModuleBase::Vector3<T> LR::LR_Spectrum<T>::cal_transition_dipole_istate_velocity_R(const int istate, const Velocity_op<std::complex<double>>& vR)
     {
         // transition density matrix D(R)
         const elecstate::DensityMatrix<T, T>& DM_trans = this->cal_transition_density_matrix(istate);
@@ -69,7 +69,7 @@ namespace LR
 
     // this algorithm is actually in use
     template<typename T>
-    ModuleBase::Vector3<T> LR::LR_Spectrum<T>::cal_transition_dipole_istate_velocity_k(const int istate, const TD_current& vR)
+    ModuleBase::Vector3<T> LR::LR_Spectrum<T>::cal_transition_dipole_istate_velocity_k(const int istate, const Velocity_op<std::complex<double>>& vR)
     {
         // transition density matrix D(R)
         const elecstate::DensityMatrix<T, T>& DM_trans = this->cal_transition_density_matrix(istate, this->X, false);
@@ -97,7 +97,7 @@ namespace LR
     template<typename T>
     void LR::LR_Spectrum<T>::cal_transition_dipoles_velocity()
     {
-        const TD_current& vR = get_velocity_matrix_R(ucell, gd_, pmat, two_center_bundle_);     // velocity matrix v(R)
+        const Velocity_op<std::complex<double>>& vR = get_velocity_matrix_R(ucell, gd_, pmat, two_center_bundle_);     // velocity matrix v(R)
         transition_dipole_.resize(nstate);
         this->mean_squared_transition_dipole_.resize(nstate);
         for (int istate = 0;istate < nstate;++istate)
@@ -148,7 +148,7 @@ namespace LR
     void LR::LR_Spectrum<T>::test_transition_dipoles_velocity_ks(const double* const ks_eig)
     {
         // velocity matrix v(R)
-        const TD_current& vR = get_velocity_matrix_R(ucell, gd_, pmat, two_center_bundle_);
+        const Velocity_op<std::complex<double>>& vR = get_velocity_matrix_R(ucell, gd_, pmat, two_center_bundle_);
         //  (e_c-e_v) of KS eigenvalues
         std::vector<double> eig_ks_diff(this->ldim);
         for (int is = 0;is < this->nspin_x;++is)

diff --git a/source/source_base/global_file.cpp b/source/source_base/global_file.cpp
@@ -25,6 +25,7 @@ void ModuleBase::Global_File::make_dir_out(
     const std::string &suffix,
 	const std::string &calculation,
     const bool &out_dir,
+    const bool &out_wfc_dir,
     const int rank,
     const bool &restart,
     const bool out_alllog)
@@ -153,6 +154,45 @@ void ModuleBase::Global_File::make_dir_out(
 #endif
     }
 
+    if(out_wfc_dir)
+    {
+        int make_dir_wfc = 0;
+        std::string command1 =  "test -d " + PARAM.globalv.global_wfc_dir + " || mkdir " + PARAM.globalv.global_wfc_dir;
+
+        times = 0;
+        while(times<GlobalV::NPROC)
+        {
+            if(rank==times)
+            {
+                if ( system( command1.c_str() ) == 0 )
+                {
+                    std::cout << " MAKE THE WFC DIR    : " << PARAM.globalv.global_wfc_dir << std::endl;
+                    make_dir_wfc = 1;
+                }
+                else
+                {
+                    std::cout << " PROC " << rank << " CAN NOT MAKE THE WFC DIR !!! " << std::endl;
+                    make_dir_wfc = 0;
+                }
+            }
+#ifdef __MPI
+            Parallel_Reduce::reduce_all(make_dir_wfc);
+#endif
+            if(make_dir_wfc>0) { break;
+}
+            ++times;
+        }
+
+#ifdef __MPI
+        if(make_dir_wfc==0)
+        {
+            std::cout << " CAN NOT MAKE THE WFC DIR......." << std::endl;
+            ModuleBase::QUIT();
+        }
+        MPI_Barrier(MPI_COMM_WORLD);
+#endif
+    }
+
     if(PARAM.inp.of_ml_gene_data == 1)
     {
         int make_dir_descrip = 0;

diff --git a/source/source_base/global_file.h b/source/source_base/global_file.h
@@ -23,6 +23,7 @@ namespace Global_File
 	void make_dir_out(const std::string &suffix,
 		const std::string &calculation,
         const bool &out_dir,
+		const bool &out_wfc_dir,
 		const int rank,
         const bool &restart,
 		const bool out_alllog = false); 

diff --git a/source/source_base/module_mixing/test/mixing_test.cpp b/source/source_base/module_mixing/test/mixing_test.cpp
diff --git a/source/source_base/test/global_file_test.cpp b/source/source_base/test/global_file_test.cpp
@@ -33,7 +33,7 @@ TEST_F(GlobalFile,mkdirout)
 {
 	    std::string output;
 	    testing::internal::CaptureStdout();
-	    ModuleBase::Global_File::make_dir_out("Si","m",false,0,true,true);
+	    ModuleBase::Global_File::make_dir_out("Si","m",false,false,0,true,true);
 	    output = testing::internal::GetCapturedStdout();
 	    EXPECT_THAT(output,testing::HasSubstr("MAKE THE DIR"));
 	    GlobalV::ofs_warning.close();
@@ -43,7 +43,7 @@ TEST_F(GlobalFile,mkdirout)
         remove(dd.c_str());
 
 	    testing::internal::CaptureStdout();
-	    ModuleBase::Global_File::make_dir_out("Si","md",false,0,true,false);
+	    ModuleBase::Global_File::make_dir_out("Si","md",false,false,0,true,false);
 	    output = testing::internal::GetCapturedStdout();
 	    EXPECT_THAT(output,testing::HasSubstr("MAKE THE STRU DIR"));
 	    EXPECT_TRUE(GlobalV::ofs_running.is_open());
@@ -53,7 +53,7 @@ TEST_F(GlobalFile,mkdirout)
         remove(bb.c_str());
 
 	    testing::internal::CaptureStdout();
-	    ModuleBase::Global_File::make_dir_out("Si","md",true,0,true,true);
+	    ModuleBase::Global_File::make_dir_out("Si","md",true,false,0,true,true);
 	    output = testing::internal::GetCapturedStdout();
 	    EXPECT_THAT(output,testing::HasSubstr("MAKE THE MATRIX DIR"));
 	    EXPECT_TRUE(GlobalV::ofs_running.is_open());

diff --git a/source/source_base/test_parallel/parallel_2d_test.sh b/source/source_base/test_parallel/parallel_2d_test.sh
diff --git a/source/source_base/test_parallel/parallel_common_test.sh b/source/source_base/test_parallel/parallel_common_test.sh
diff --git a/source/source_base/test_parallel/parallel_global_test.sh b/source/source_base/test_parallel/parallel_global_test.sh
diff --git a/source/source_base/test_parallel/parallel_reduce_test.sh b/source/source_base/test_parallel/parallel_reduce_test.sh
diff --git a/source/source_basis/module_ao/ORB_nonlocal_lm.h b/source/source_basis/module_ao/ORB_nonlocal_lm.h
@@ -45,6 +45,9 @@ class Numerical_Nonlocal_Lm
 	const double& getKpoint(const int &ik) const { return this->k_radial[ik]; }
 	const double* getBeta_k() const { return this->beta_k; }
 	const double& getBeta_k(const int &ik) const { return this->beta_k[ik]; }
+
+	const int& getNk() const { return nk; }
+	const double& getDruniform() const { return dr_uniform; }
 
     // enables deep copy
 	Numerical_Nonlocal_Lm& operator= (const Numerical_Nonlocal_Lm& nol );

diff --git a/source/source_basis/module_ao/test/parallel_orbitals_test.sh b/source/source_basis/module_ao/test/parallel_orbitals_test.sh
diff --git a/source/source_basis/module_nao/test/I_sphbes_coeff.txt b/source/source_basis/module_nao/test/I_sphbes_coeff.txt
diff --git a/source/source_basis/module_pw/test/test_gnu.sh b/source/source_basis/module_pw/test/test_gnu.sh
diff --git a/source/source_basis/module_pw/test/test_intel.sh b/source/source_basis/module_pw/test/test_intel.sh
diff --git a/source/source_cell/test/bcast_atom_pseudo_test.sh b/source/source_cell/test/bcast_atom_pseudo_test.sh
diff --git a/source/source_cell/test/bcast_atom_spec_test.sh b/source/source_cell/test/bcast_atom_spec_test.sh
diff --git a/source/source_cell/test/klist_test_para.sh b/source/source_cell/test/klist_test_para.sh
diff --git a/source/source_cell/test/parallel_kpoints_test.sh b/source/source_cell/test/parallel_kpoints_test.sh
diff --git a/source/source_cell/test/unitcell_test_parallel.sh b/source/source_cell/test/unitcell_test_parallel.sh
diff --git a/source/source_cell/test_pw/unitcell_test_pw_para.sh b/source/source_cell/test_pw/unitcell_test_pw_para.sh
diff --git a/source/source_esolver/esolver.cpp b/source/source_esolver/esolver.cpp
@@ -230,13 +230,26 @@ ESolver* init_esolver(const Input_para& inp, UnitCell& ucell)
     }
     else if (esolver_type == "ksdft_lcao_tddft")
     {
-#if ((defined __CUDA) /* || (defined __ROCM) */)
-        if (PARAM.inp.device == "gpu")
+        if (PARAM.inp.nspin < 4)
         {
-            return new ESolver_KS_LCAO_TDDFT<base_device::DEVICE_GPU>();
+#if ((defined __CUDA) /* || (defined __ROCM) */)
+            if (PARAM.inp.device == "gpu")
+            {
+                return new ESolver_KS_LCAO_TDDFT<double, base_device::DEVICE_GPU>();
+            }
+#endif
+            return new ESolver_KS_LCAO_TDDFT<double, base_device::DEVICE_CPU>();
         }
+        else
+        {
+#if ((defined __CUDA) /* || (defined __ROCM) */)
+            if (PARAM.inp.device == "gpu")
+            {
+                return new ESolver_KS_LCAO_TDDFT<std::complex<double>, base_device::DEVICE_GPU>();
+            }
 #endif
-        return new ESolver_KS_LCAO_TDDFT<base_device::DEVICE_CPU>();
+            return new ESolver_KS_LCAO_TDDFT<std::complex<double>, base_device::DEVICE_CPU>();
+        }
     }
     else if (esolver_type == "lr_lcao")
     {

diff --git a/source/source_esolver/esolver_ks.cpp b/source/source_esolver/esolver_ks.cpp
@@ -287,7 +287,10 @@ void ESolver_KS<T, Device>::before_scf(UnitCell& ucell, const int istep)
 template <typename T, typename Device>
 void ESolver_KS<T, Device>::iter_init(UnitCell& ucell, const int istep, const int iter)
 {
-    ModuleIO::write_head(GlobalV::ofs_running, istep, iter, this->basisname);
+    if(PARAM.inp.esolver_type != "tddft")
+    {
+        ModuleIO::write_head(GlobalV::ofs_running, istep, iter, this->basisname);
+    }
 
 #ifdef __MPI
     iter_time = MPI_Wtime();

diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp
@@ -161,7 +161,7 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
     }
 
     // 5) read psi from file
-    if (PARAM.inp.init_wfc == "file")
+    if (PARAM.inp.init_wfc == "file"&& PARAM.inp.esolver_type != "tddft")
 	{
 		if (!ModuleIO::read_wfc_nao(PARAM.globalv.global_readin_dir, 
 					this->pv,