deepmodeling · ESROAMER · Jun 27, 2025 · Jun 27, 2025 · Jun 27, 2025 · Jun 27, 2025
diff --git a/source/Makefile.Objects b/source/Makefile.Objects
@@ -354,6 +354,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\
@@ -623,8 +624,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_tddft/td_current.h"
+#include "source_lcao/module_tddft/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)
@@ -183,4 +183,4 @@ namespace LR
     }
 }
 template class LR::LR_Spectrum<double>;
-template class LR::LR_Spectrum<std::complex<double>>;
+template class LR::LR_Spectrum<std::complex<double>>;
diff --git a/source/module_parameter/input_parameter.h b/source/module_parameter/input_parameter.h
@@ -283,7 +283,8 @@ struct Input_para
     double bessel_descriptor_sigma = 0.1;           ///< spherical bessel smearing_sigma
 
     // ==============   #Parameters (9.rt-tddft) ===========================
-    double td_force_dt = 0.02; ///<"fs"
+    double td_dt = -1.0;            ///< time step for propagation
+    int estep_per_md = 1;           ///< number of electronic steps per MD step
     bool td_vext = false;      ///< add extern potential or not
     // std::string td_vext_dire = "1";   ///< vext direction
     std::vector<int> td_vext_dire = {1}; ///< vector of vext direction

diff --git a/source/module_parameter/system_parameter.h b/source/module_parameter/system_parameter.h
@@ -43,6 +43,7 @@ struct System_para
     std::string global_readin_dir = "";     ///< global readin directory
     std::string global_stru_dir = "";       ///< global structure directory
     std::string global_matrix_dir = "";     ///< global matrix directory
+    std::string global_wfc_dir = "";       ///< global wavefunction directory
     std::string global_mlkedf_descriptor_dir = ""; ///< global ML KEDF descriptor directory
     std::string global_deepks_label_elec_dir = ""; ///< global directory for DeePKS labels during electronic steps 
     std::string log_file = "log";           ///< log file name
@@ -65,4 +66,4 @@ struct System_para
 
     bool search_pbc = true;                 ///< whether to search for periodic boundary conditions, force set to true
 };
-#endif
+#endif
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,46 @@ 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/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_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_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,