diff --git a/source/Makefile.Objects b/source/Makefile.Objects
index b76a83cc82..0a730089a2 100644
--- a/source/Makefile.Objects
+++ b/source/Makefile.Objects
@@ -496,6 +496,7 @@ OBJS_LCAO=DM_gamma.o\
       middle_hamilt.o\
       norm_psi.o\
       propagator.o\
+      td_velocity.o\
       upsi.o\
       FORCE_STRESS.o\
       FORCE_gamma.o\
diff --git a/source/module_elecstate/potentials/H_TDDFT_pw.h b/source/module_elecstate/potentials/H_TDDFT_pw.h
index 7e244ef0bc..ea6b4f0ef9 100644
--- a/source/module_elecstate/potentials/H_TDDFT_pw.h
+++ b/source/module_elecstate/potentials/H_TDDFT_pw.h
@@ -133,4 +133,4 @@ class H_TDDFT_pw : public PotBase
 
 } // namespace elecstate
 
-#endif
\ No newline at end of file
+#endif
diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp
index 3626a5ecf1..995e5b5e08 100644
--- a/source/module_esolver/esolver_ks_lcao_tddft.cpp
+++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp
@@ -15,6 +15,7 @@
 #include "module_elecstate/module_charge/symmetry_rho.h"
 #include "module_elecstate/occupy.h"
 #include "module_hamilt_lcao/module_tddft/evolve_elec.h"
+#include "module_hamilt_lcao/module_tddft/td_velocity.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
 #include "module_io/print_info.h"
 
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.cpp
index acc6328acc..4102aae997 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.cpp
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.cpp
@@ -5,6 +5,7 @@
 #include "module_hamilt_lcao/hamilt_lcaodft/center2_orb-orb11.h"
 #include "module_elecstate/potentials/H_TDDFT_pw.h"
 #include "module_hamilt_lcao/module_hcontainer/hcontainer_funcs.h"
+#include "module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h"
 
 #include "module_cell/module_neighbor/sltk_grid_driver.h"
 #include "module_base/libm/libm.h"
@@ -29,6 +30,10 @@ TDEkinetic<OperatorLCAO<TK, TR>>::TDEkinetic(LCAO_Matrix* LM_in,
         this->init_td();
         // initialize HR to get adjs info.
         this->initialize_HR(Grid,this->LM->ParaV);
+        if(TD_Velocity::out_mat_R == true)
+        {
+            out_mat_R = true;
+        }
     }
 template <typename TK, typename TR>
 TDEkinetic<OperatorLCAO<TK, TR>>::~TDEkinetic()
@@ -37,6 +42,7 @@ TDEkinetic<OperatorLCAO<TK, TR>>::~TDEkinetic()
     {
         delete this->hR_tmp;
     }
+    TD_Velocity::td_vel_op = nullptr;
 }
 //term A^2*S
 template <typename TK, typename TR>
@@ -185,22 +191,24 @@ void TDEkinetic<OperatorLCAO<TK, TR>>::cal_HR_IJR(const int& iat1,
 template <typename TK, typename TR>
 void TDEkinetic<OperatorLCAO<TK, TR>>::init_td(void)
 {
+    TD_Velocity::td_vel_op = &td_velocity;
     //calculate At in cartesian coorinates.
-	double l_norm[3]={GlobalC::ucell.a1.norm() ,GlobalC::ucell.a2.norm() ,GlobalC::ucell.a3.norm()};
+	double l_norm[3]={this->ucell->a1.norm() ,this->ucell->a2.norm() ,this->ucell->a3.norm()};
     double (&A)[3] = elecstate::H_TDDFT_pw::At;
-	cart_At = GlobalC::ucell.a1*A[0]/l_norm[0] + GlobalC::ucell.a2*A[1]/l_norm[1] + GlobalC::ucell.a3*A[2]/l_norm[2];
+	cart_At = this->ucell->a1*A[0]/l_norm[0] + this->ucell->a2*A[1]/l_norm[1] + this->ucell->a3*A[2]/l_norm[2];
     std::cout << "cart_At: " << cart_At[0] << " " <<cart_At[1]<< " " << cart_At[2] << std::endl;
-	
+
     //init MOT,MGT
+    const LCAO_Orbitals& orb = LCAO_Orbitals::get_const_instance();
     this->MOT.allocate(
-		GlobalC::ORB.get_ntype(),	// number of atom types
-		GlobalC::ORB.get_lmax(),	// max L used to calculate overlap
-		static_cast<int>(GlobalC::ORB.get_kmesh()) | 1,				// kpoints, for integration in k space
-		GlobalC::ORB.get_Rmax(),				// max value of radial table
-		GlobalC::ORB.get_dR(),								// delta R, for making radial table
-		GlobalC::ORB.get_dk());											// Peize Lin change 2017-04-16
+		orb.get_ntype(),	                                // number of atom types
+		orb.get_lmax(),	                                    // max L used to calculate overlap
+		static_cast<int>(orb.get_kmesh()) | 1,				// kpoints, for integration in k space
+		orb.get_Rmax(),				                        // max value of radial table
+		orb.get_dR(),								        // delta R, for making radial table
+		orb.get_dk());                                      // Peize Lin change 2017-04-16
 	int Lmax_used, Lmax;
-	this->MOT.init_Table_Spherical_Bessel (2, 1, Lmax_used, Lmax, 1, GlobalC::ORB, GlobalC::ucell.infoNL.Beta);
+	this->MOT.init_Table_Spherical_Bessel (2, 1, Lmax_used, Lmax, 1, orb, this->ucell->infoNL.Beta);
 
 	//=========================================
 	// (2) init Ylm Coef
@@ -214,17 +222,17 @@ void TDEkinetic<OperatorLCAO<TK, TR>>::init_td(void)
 	this->MGT.init_Gaunt( Lmax );
 
     //init_radial table
-    for( size_t TA=0; TA!=GlobalC::ORB.get_ntype(); ++TA )
-		for( size_t TB=0; TB!=GlobalC::ORB.get_ntype(); ++TB )
-			for( int LA=0; LA<=GlobalC::ORB.Phi[TA].getLmax(); ++LA )
-				for( size_t NA=0; NA!=GlobalC::ORB.Phi[TA].getNchi(LA); ++NA )
-					for( int LB=0; LB<=GlobalC::ORB.Phi[TB].getLmax(); ++LB )
-						for( size_t NB=0; NB!=GlobalC::ORB.Phi[TB].getNchi(LB); ++NB )
+    for( size_t TA=0; TA!=orb.get_ntype(); ++TA )
+		for( size_t TB=0; TB!=orb.get_ntype(); ++TB )
+			for( int LA=0; LA<=orb.Phi[TA].getLmax(); ++LA )
+				for( size_t NA=0; NA!=orb.Phi[TA].getNchi(LA); ++NA )
+					for( int LB=0; LB<=orb.Phi[TB].getLmax(); ++LB )
+						for( size_t NB=0; NB!=orb.Phi[TB].getNchi(LB); ++NB )
 							center2_orb11_s[TA][TB][LA][NA][LB].insert(
 								std::make_pair(NB, 
                                 Center2_Orb::Orb11(
-									GlobalC::ORB.Phi[TA].PhiLN(LA,NA),
-									GlobalC::ORB.Phi[TB].PhiLN(LB,NB),
+									orb.Phi[TA].PhiLN(LA,NA),
+									orb.Phi[TB].PhiLN(LB,NB),
 									this->MOT, this->MGT)));
 	for( auto &coA : center2_orb11_s )
 		for( auto &coB : coA.second )
@@ -361,22 +369,39 @@ void TDEkinetic<OperatorLCAO<TK, TR>>::contributeHk(int ik)
 template<>
 void TDEkinetic<OperatorLCAO<std::complex<double>, double>>::contributeHk(int ik)
 {
-    if (GlobalV::ESOLVER_TYPE != "tddft" || elecstate::H_TDDFT_pw::stype != 1)
+    if (TD_Velocity::tddft_velocity == false)
     {
         return;
     }
     else{        
         ModuleBase::TITLE("TDEkinetic", "contributeHk");
         ModuleBase::timer::tick("TDEkinetic", "contributeHk");
+        const Parallel_Orbitals* paraV = this->hR_tmp->get_atom_pair(0).get_paraV();
+        //save HR data for output
+        int spin_tot = paraV->nspin;
+        if(spin_tot==4);
+        else if(!output_hR_done && out_mat_R)
+        {
+            for(int spin_now = 0;spin_now < spin_tot;spin_now++)
+            {
+                sparse_format::cal_HContainer_cd(
+                    *(paraV),
+                    spin_now, 
+                    1e-10, 
+                    *hR_tmp, 
+                    td_velocity.HR_sparse_td_vel[spin_now]);
+            }
+            output_hR_done = true;
+        }
         //folding inside HR to HK
         if(ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER())
         {
-            const int nrow = this->LM->ParaV->get_row_size();
+            const int nrow = paraV->get_row_size();
             hamilt::folding_HR(*this->hR_tmp, this->hK->data(), this->kvec_d[ik], nrow, 1);
         }
         else
         {
-            const int ncol = this->LM->ParaV->get_col_size();
+            const int ncol = paraV->get_col_size();
             hamilt::folding_HR(*this->hR_tmp, this->hK->data(), this->kvec_d[ik], ncol, 0);
         }
         
@@ -388,4 +413,4 @@ template class TDEkinetic<hamilt::OperatorLCAO<double, double>>;
 template class TDEkinetic<hamilt::OperatorLCAO<std::complex<double>, double>>;
 template class TDEkinetic<hamilt::OperatorLCAO<std::complex<double>, std::complex<double>>>;
 
-}
\ No newline at end of file
+}
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.h
index 4f83268bff..70fa1f08c5 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.h
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.h
@@ -8,6 +8,7 @@
 #include "module_hamilt_lcao/module_hcontainer/hcontainer.h"
 #include "module_basis/module_ao/ORB_table_phi.h"
 #include "module_basis/module_ao/ORB_gaunt_table.h"
+#include "module_hamilt_lcao/module_tddft/td_velocity.h"
 
 
 namespace hamilt
@@ -75,6 +76,8 @@ class TDEkinetic<OperatorLCAO<TK,TR>> : public OperatorLCAO<TK, TR>
     void calculate_HR(void);
     virtual void set_HR_fixed(void*)override;
 
+    TD_Velocity td_velocity;
+
   private:
     const UnitCell* ucell = nullptr;
     
@@ -109,6 +112,8 @@ class TDEkinetic<OperatorLCAO<TK,TR>> : public OperatorLCAO<TK, TR>
 
     bool hR_tmp_done = false;
     bool allocated = false;
+    bool output_hR_done = false;
+    bool out_mat_R = false;
 };
 
 } // namespace hamilt
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_nonlocal_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_nonlocal_lcao.cpp
index 1e29913f21..e1e0bcede9 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_nonlocal_lcao.cpp
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_nonlocal_lcao.cpp
@@ -47,9 +47,9 @@ template <typename TK, typename TR>
 void hamilt::TDNonlocal<hamilt::OperatorLCAO<TK, TR>>::init_td(void)
 {
     //calculate At in cartesian coorinates.
-	double l_norm[3]={GlobalC::ucell.a1.norm() ,GlobalC::ucell.a2.norm() ,GlobalC::ucell.a3.norm()};
+	double l_norm[3]={this->ucell->a1.norm() ,this->ucell->a2.norm() ,this->ucell->a3.norm()};
     double (&A)[3] = elecstate::H_TDDFT_pw::At;
-	cart_At = -(GlobalC::ucell.a1*A[0]/l_norm[0] + GlobalC::ucell.a2*A[1]/l_norm[1] + GlobalC::ucell.a3*A[2]/l_norm[2]);
+	cart_At = -(this->ucell->a1*A[0]/l_norm[0] + this->ucell->a2*A[1]/l_norm[1] + this->ucell->a3*A[2]/l_norm[2]);
     std::cout << "cart_At: " << cart_At[0] << " " <<cart_At[1]<< " " << cart_At[2] << std::endl;
 }
 // initialize_HR()
@@ -372,7 +372,7 @@ void hamilt::TDNonlocal<hamilt::OperatorLCAO<TK, TR>>::contributeHk(int ik)
 template<>
 void hamilt::TDNonlocal<hamilt::OperatorLCAO<std::complex<double>, double>>::contributeHk(int ik)
 {
-    if (GlobalV::ESOLVER_TYPE != "tddft" || elecstate::H_TDDFT_pw::stype != 1)
+    if (TD_Velocity::tddft_velocity == false)
     {
         return;
     }
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_nonlocal_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_nonlocal_lcao.h
index 5d7b250c67..e485e5b66a 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_nonlocal_lcao.h
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_nonlocal_lcao.h
@@ -8,6 +8,7 @@
 #include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h"
 #include "module_hamilt_lcao/module_hcontainer/hcontainer.h"
 #include "module_elecstate/potentials/H_TDDFT_pw.h"
+#include "module_hamilt_lcao/module_tddft/td_velocity.h"
 
 namespace hamilt
 {
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp
index 5b018e9731..8829bef774 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp
@@ -3,6 +3,7 @@
 #include "spar_u.h"
 #include "spar_exx.h"
 #include "module_hamilt_lcao/module_hcontainer/hcontainer.h"
+#include "module_hamilt_lcao/module_tddft/td_velocity.h"
 
 void sparse_format::cal_HSR(
         const Parallel_Orbitals &pv,
@@ -25,12 +26,24 @@ void sparse_format::cal_HSR(
         hamilt::HamiltLCAO<std::complex<double>, double>* p_ham_lcao = 
         dynamic_cast<hamilt::HamiltLCAO<std::complex<double>, double>*>(p_ham);
 
-		sparse_format::cal_HContainer_d(
+        if(TD_Velocity::tddft_velocity)
+        {
+            sparse_format::cal_HContainer_td(
+                pv,
+                current_spin, 
+				sparse_thr, 
+				*(p_ham_lcao->getHR()), 
+				TD_Velocity::td_vel_op->HR_sparse_td_vel[current_spin]);
+        }
+        else
+        {
+            sparse_format::cal_HContainer_d(
                 pv,
                 current_spin, 
 				sparse_thr, 
 				*(p_ham_lcao->getHR()), 
 				lm.HR_sparse[current_spin]);
+        }
 
 		sparse_format::cal_HContainer_d(
 				pv,
@@ -202,6 +215,49 @@ void sparse_format::cal_HContainer_cd(
     return;
 }
 
+void sparse_format::cal_HContainer_td(
+        const Parallel_Orbitals &pv,
+		const int &current_spin, 
+		const double &sparse_thr, 
+		const hamilt::HContainer<double>& hR, 
+		std::map<Abfs::Vector3_Order<int>, 
+		std::map<size_t, std::map<size_t, std::complex<double>>>>& target)
+{
+    ModuleBase::TITLE("sparse_format","cal_HContainer_td");
+
+    auto row_indexes = pv.get_indexes_row();
+    auto col_indexes = pv.get_indexes_col();
+    for(int iap=0;iap<hR.size_atom_pairs();++iap)
+    {
+        int atom_i = hR.get_atom_pair(iap).get_atom_i();
+        int atom_j = hR.get_atom_pair(iap).get_atom_j();
+        int start_i = pv.atom_begin_row[atom_i];
+        int start_j = pv.atom_begin_col[atom_j];
+        int row_size = pv.get_row_size(atom_i);
+        int col_size = pv.get_col_size(atom_j);
+        for(int iR=0;iR<hR.get_atom_pair(iap).get_R_size();++iR)
+        {
+            auto& matrix = hR.get_atom_pair(iap).get_HR_values(iR);
+            int* r_index = hR.get_atom_pair(iap).get_R_index(iR);
+            Abfs::Vector3_Order<int> dR(r_index[0], r_index[1], r_index[2]);
+            for(int i=0;i<row_size;++i)
+            {
+                int mu = row_indexes[start_i+i];
+                for(int j=0;j<col_size;++j)
+                {
+                    int nu = col_indexes[start_j+j];
+                    const auto& value_tmp = std::complex<double>(matrix.get_value(i,j) , 0.0);
+                    if(std::abs(value_tmp)>sparse_thr)
+                    {
+                        target[dR][mu][nu] += value_tmp;
+                    }
+                }
+            }
+        }
+    }
+
+    return;
+}
 
 // in case there are elements smaller than the threshold
 void sparse_format::clear_zero_elements(
@@ -232,6 +288,28 @@ void sparse_format::clear_zero_elements(
                 }
             }
         }
+        if(TD_Velocity::tddft_velocity)
+        {
+            for (auto &R_loop : TD_Velocity::td_vel_op->HR_sparse_td_vel[current_spin])
+        {
+            for (auto &row_loop : R_loop.second)
+            {
+                auto &col_map = row_loop.second;
+                auto iter = col_map.begin();
+                while (iter != col_map.end())
+                {
+                    if (std::abs(iter->second) <= sparse_thr)
+                    {
+                        col_map.erase(iter++);
+                    }
+                    else
+                    {
+                        iter++;
+                    }
+                }
+            }
+        }
+        }
 
         for (auto &R_loop : lm.SR_sparse)
         {
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h b/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h
index 908dfb1909..d7b36f397f 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h
@@ -31,6 +31,14 @@ namespace sparse_format
 			std::map<Abfs::Vector3_Order<int>, 
 			std::map<size_t, std::map<size_t, std::complex<double>>>>& target);
 
+	void cal_HContainer_td(
+            const Parallel_Orbitals &pv,
+			const int &current_spin, 
+			const double &sparse_threshold, 
+			const hamilt::HContainer<double>& hR, 
+			std::map<Abfs::Vector3_Order<int>, 
+			std::map<size_t, std::map<size_t, std::complex<double>>>>& target);
+
 	void clear_zero_elements(
 			LCAO_Matrix &lm,
 			const int &current_spin, 
diff --git a/source/module_hamilt_lcao/module_tddft/CMakeLists.txt b/source/module_hamilt_lcao/module_tddft/CMakeLists.txt
index f41dad565d..ad921933f3 100644
--- a/source/module_hamilt_lcao/module_tddft/CMakeLists.txt
+++ b/source/module_hamilt_lcao/module_tddft/CMakeLists.txt
@@ -7,6 +7,7 @@ if(ENABLE_LCAO)
         norm_psi.cpp
         propagator.cpp
         upsi.cpp
+        td_velocity.cpp
     )
 
     add_library(
diff --git a/source/module_hamilt_lcao/module_tddft/td_velocity.cpp b/source/module_hamilt_lcao/module_tddft/td_velocity.cpp
new file mode 100644
index 0000000000..5de3d35209
--- /dev/null
+++ b/source/module_hamilt_lcao/module_tddft/td_velocity.cpp
@@ -0,0 +1,24 @@
+#include "module_base/timer.h"
+#include "td_velocity.h"
+
+
+bool TD_Velocity::tddft_velocity = false;
+bool TD_Velocity::out_mat_R = false;
+TD_Velocity* TD_Velocity::td_vel_op = nullptr;
+
+TD_Velocity::TD_Velocity()
+{
+  return;
+}
+TD_Velocity::~TD_Velocity()
+{
+  this->destroy_HS_R_td_sparse();
+}
+
+void TD_Velocity::destroy_HS_R_td_sparse(void)
+{
+  std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, std::complex<double>>>> empty_HR_sparse_td_vel_up;
+  std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, std::complex<double>>>> empty_HR_sparse_td_vel_down;
+  HR_sparse_td_vel[0].swap(empty_HR_sparse_td_vel_up);
+  HR_sparse_td_vel[1].swap(empty_HR_sparse_td_vel_down);
+}
\ No newline at end of file
diff --git a/source/module_hamilt_lcao/module_tddft/td_velocity.h b/source/module_hamilt_lcao/module_tddft/td_velocity.h
new file mode 100644
index 0000000000..0d274228dd
--- /dev/null
+++ b/source/module_hamilt_lcao/module_tddft/td_velocity.h
@@ -0,0 +1,34 @@
+#ifndef TD_VELOCITY_H
+#define TD_VELOCITY_H
+#include <map>
+
+#include "module_base/timer.h"
+#include "module_base/abfs-vector3_order.h"
+//Class to store TDDFT velocity gague infos.
+class TD_Velocity
+{
+  public:
+    TD_Velocity();
+    ~TD_Velocity();
+
+    /// @brief Judge if in tddft calculation or not 
+    static bool tddft_velocity;
+
+    /// @brief switch to control the output of HR
+    static bool out_mat_R;
+
+    /// @brief pointer to the only TD_Velocity object itself
+    static TD_Velocity* td_vel_op;
+
+    //For TDDFT velocity gague, to fix the output of HR
+    std::map<Abfs::Vector3_Order<int>, std::map<size_t, std::map<size_t, std::complex<double>>>> HR_sparse_td_vel[2];
+
+
+  private:
+
+    /// @brief destory HSR data stored
+    void destroy_HS_R_td_sparse(void);
+
+};
+
+#endif
\ No newline at end of file
diff --git a/source/module_io/input_conv.cpp b/source/module_io/input_conv.cpp
index 54de3d6e1e..0f269d336d 100644
--- a/source/module_io/input_conv.cpp
+++ b/source/module_io/input_conv.cpp
@@ -25,6 +25,7 @@
 #include "module_hamilt_lcao/hamilt_lcaodft/local_orbital_charge.h"
 #include "module_hamilt_lcao/module_dftu/dftu.h"
 #include "module_hamilt_lcao/module_tddft/evolve_elec.h"
+#include "module_hamilt_lcao/module_tddft/td_velocity.h"
 #endif
 #ifdef __PEXSI
 #include "module_hsolver/module_pexsi/pexsi_solver.h"
@@ -135,7 +136,22 @@ std::vector<double> Input_Conv::convert_units(std::string params, double c)
 void Input_Conv::read_td_efield()
 {
     elecstate::H_TDDFT_pw::stype = INPUT.td_stype;
-
+    if(INPUT.esolver_type == "tddft" && elecstate::H_TDDFT_pw::stype == 1)
+    {
+        TD_Velocity::tddft_velocity = true;
+    }
+    else
+    {
+        TD_Velocity::tddft_velocity = false;
+    }
+    if(INPUT.out_mat_hs2==1)
+    {
+        TD_Velocity::out_mat_R = true;
+    }
+    else
+    {
+        TD_Velocity::out_mat_R = false;
+    }
     parse_expression(INPUT.td_ttype, elecstate::H_TDDFT_pw::ttype);
 
     elecstate::H_TDDFT_pw::tstart = INPUT.td_tstart;
diff --git a/source/module_io/td_current_io.cpp b/source/module_io/td_current_io.cpp
index 5f69694801..4edbc2c979 100644
--- a/source/module_io/td_current_io.cpp
+++ b/source/module_io/td_current_io.cpp
@@ -8,6 +8,7 @@
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
 #include "module_elecstate/module_dm/cal_dm_psi.h"
 #include "module_base/parallel_reduce.h"
+#include "module_elecstate/potentials/H_TDDFT_pw.h"
 
 #ifdef __LCAO
 //init DSloc_R for current calculation
@@ -103,7 +104,7 @@ void ModuleIO::cal_tmp_DM(elecstate::DensityMatrix<std::complex<double>, double>
                 std::complex<double> kphase = std::complex<double>(cosp, sinp);
                 // set DMR element
                 double* tmp_DMR_pointer = tmp_matrix->get_pointer();
-                std::complex<double>* tmp_DMK_pointer = DM.get_DMK_vector()[ik + ik_begin].data();
+                std::complex<double>* tmp_DMK_pointer = DM.get_DMK_pointer(ik + ik_begin);
                 double* DMK_real_pointer = nullptr;
                 double* DMK_imag_pointer = nullptr;
                 // jump DMK to fill DMR
diff --git a/source/module_io/test/for_testing_input_conv.h b/source/module_io/test/for_testing_input_conv.h
index b012dbf183..d4449bdbf1 100644
--- a/source/module_io/test/for_testing_input_conv.h
+++ b/source/module_io/test/for_testing_input_conv.h
@@ -9,6 +9,7 @@
 #include "module_elecstate/module_charge/charge_mixing.h"
 #include "module_elecstate/occupy.h"
 #include "module_elecstate/potentials/H_TDDFT_pw.h"
+#include "module_hamilt_lcao/module_tddft/td_velocity.h"
 #include "module_elecstate/potentials/efield.h"
 #include "module_elecstate/potentials/gatefield.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h"
@@ -77,6 +78,9 @@ double elecstate::H_TDDFT_pw::lcut2;
 
 // time domain parameters
 
+bool TD_Velocity::tddft_velocity;
+bool TD_Velocity::out_mat_R;
+
 // Gauss
 int elecstate::H_TDDFT_pw::gauss_count;
 std::vector<double> elecstate::H_TDDFT_pw::gauss_omega; // time(a.u.)^-1
diff --git a/source/module_io/write_HS_sparse.cpp b/source/module_io/write_HS_sparse.cpp
index aafe937d66..c7e7d5475a 100644
--- a/source/module_io/write_HS_sparse.cpp
+++ b/source/module_io/write_HS_sparse.cpp
@@ -4,6 +4,7 @@
 #include "module_base/timer.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
 #include "single_R_io.h"
+#include "module_hamilt_lcao/module_tddft/td_velocity.h"
 
 void ModuleIO::save_HSR_sparse(
     const int &istep,
@@ -53,12 +54,26 @@ void ModuleIO::save_HSR_sparse(
         {
             for (int ispin = 0; ispin < spin_loop; ++ispin)
             {
-                auto iter = HR_sparse_ptr[ispin].find(R_coor);
-                if (iter != HR_sparse_ptr[ispin].end())
+                if(TD_Velocity::tddft_velocity)
                 {
-                    for (auto &row_loop : iter->second)
+                    auto iter = TD_Velocity::td_vel_op->HR_sparse_td_vel[ispin].find(R_coor);
+                    if (iter != TD_Velocity::td_vel_op->HR_sparse_td_vel[ispin].end())
                     {
-                        H_nonzero_num[ispin][count] += row_loop.second.size();
+                        for (auto &row_loop : iter->second)
+                        {
+                            H_nonzero_num[ispin][count] += row_loop.second.size();
+                        }
+                    }
+                }
+                else
+                {
+                    auto iter = HR_sparse_ptr[ispin].find(R_coor);
+                    if (iter != HR_sparse_ptr[ispin].end())
+                    {
+                        for (auto &row_loop : iter->second)
+                        {
+                            H_nonzero_num[ispin][count] += row_loop.second.size();
+                        }
                     }
                 }
             }
@@ -279,7 +294,14 @@ void ModuleIO::save_HSR_sparse(
             {
                 if (GlobalV::NSPIN != 4)
                 {
-                    output_single_R(g1[ispin], HR_sparse_ptr[ispin][R_coor], sparse_thr, binary, *lm.ParaV);
+                    if(TD_Velocity::tddft_velocity)
+                    {
+                        output_single_R(g1[ispin], TD_Velocity::td_vel_op->HR_sparse_td_vel[ispin][R_coor], sparse_thr, binary, *lm.ParaV);
+                    }
+                    else
+                    {
+                        output_single_R(g1[ispin], HR_sparse_ptr[ispin][R_coor], sparse_thr, binary, *lm.ParaV);
+                    }
                 }
                 else
                 {