Upload

source/Makefile.Objects

@@ -329,6 +329,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\
@@ -598,8 +599,8 @@ OBJS_LCAO=evolve_elec.o\
       propagator_cn2.o\
       propagator_taylor.o\
       propagator_etrs.o\
-      td_velocity.o\
-      td_current.o\
+      td_info.o\
+      velocity_op.o\
       snap_psibeta_half_tddft.o\
       solve_propagation.o\
       upsi.o\

source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt

@@ -13,6 +13,7 @@ if(ENABLE_LCAO)
         operator_lcao/nonlocal_new.cpp
         operator_lcao/td_ekinetic_lcao.cpp
         operator_lcao/td_nonlocal_lcao.cpp
+        operator_lcao/td_pot_hybrid.cpp
         operator_lcao/dspin_lcao.cpp
         operator_lcao/dftu_lcao.cpp
         pulay_force_stress_center2.cpp

source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp

@@ -38,6 +38,7 @@
 #include "operator_lcao/overlap_new.h"
 #include "operator_lcao/td_ekinetic_lcao.h"
 #include "operator_lcao/td_nonlocal_lcao.h"
+#include "operator_lcao/td_pot_hybrid.h"
 #include "operator_lcao/veff_lcao.h"
 
 namespace hamilt
@@ -344,12 +345,8 @@ HamiltLCAO<TK, TR>::HamiltLCAO(Gint_Gamma* GG_in,
         }
 #endif
         // TDDFT_velocity_gauge
-        if (TD_Velocity::tddft_velocity)
+        if (PARAM.inp.esolver_type == "tddft" && PARAM.inp.td_stype == 1)
         {
-            if (!TD_Velocity::init_vecpot_file)
-            {
-                elecstate::H_TDDFT_pw::update_At();
-            }
             Operator<TK>* td_ekinetic = new TDEkinetic<OperatorLCAO<TK, TR>>(this->hsk,
                                                                              this->hR,
                                                                              this->kv,
@@ -359,10 +356,27 @@ HamiltLCAO<TK, TR>::HamiltLCAO(Gint_Gamma* GG_in,
                                                                              two_center_bundle.overlap_orb.get());
             this->getOperator()->add(td_ekinetic);
 
-            Operator<TK>* td_nonlocal
-                = new TDNonlocal<OperatorLCAO<TK, TR>>(this->hsk, this->kv->kvec_d, this->hR, &ucell, orb, &grid_d);
+            Operator<TK>* td_nonlocal = new TDNonlocal<OperatorLCAO<TK, TR>>(this->hsk,
+                                                                             this->kv->kvec_d,
+                                                                             this->hR,
+                                                                             &ucell,
+                                                                             orb,
+                                                                             &grid_d);
             this->getOperator()->add(td_nonlocal);
         }
+        if (PARAM.inp.esolver_type == "tddft" && PARAM.inp.td_stype == 2)
+        {
+            Operator<TK>* td_pot_hybrid = new TD_pot_hybrid<OperatorLCAO<TK, TR>>(this->hsk,
+                                                                           this->kv,
+                                                                           this->hR,
+                                                                           this->sR,
+                                                                           orb,
+                                                                           &ucell,
+                                                                           orb.cutoffs(),
+                                                                           &grid_d,
+                                                                           two_center_bundle.kinetic_orb.get());
+            this->getOperator()->add(td_pot_hybrid);
+        }
         if (PARAM.inp.dft_plus_u)
         {
             Operator<TK>* dftu = nullptr;

source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/CMakeLists.txt

@@ -11,6 +11,7 @@ add_library(
     nonlocal_new.cpp
     td_ekinetic_lcao.cpp
     td_nonlocal_lcao.cpp
+    td_pot_hybrid.cpp
     dspin_lcao.cpp
     dftu_lcao.cpp
 )

source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.cpp

@@ -12,6 +12,8 @@
 #include "source_hsolver/diago_elpa_native.h"
 #endif
 
+#include "module_hamilt_lcao/module_tddft/td_info.h"
+
 namespace hamilt {
 
 template <>
@@ -200,7 +202,7 @@ void OperatorLCAO<TK, TR>::init(const int ik_in) {
 
         break;
     }
-    case calculation_type::lcao_tddft_velocity: {
+    case calculation_type::lcao_tddft_periodic: {
         if (!this->hr_done) {
             // in cal_type=lcao_fixed, HR should be updated by each sub-chain
             // nodes
@@ -240,8 +242,8 @@ void OperatorLCAO<TK, TR>::init(const int ik_in) {
 }
 
 // contributeHk()
-template <typename TK, typename TR>
-void OperatorLCAO<TK, TR>::contributeHk(int ik) {
+template <>
+void OperatorLCAO<double, double>::contributeHk(int ik) {
     ModuleBase::TITLE("OperatorLCAO", "contributeHk");
     ModuleBase::timer::tick("OperatorLCAO", "contributeHk");
     if(ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER(PARAM.inp.ks_solver))
@@ -256,6 +258,37 @@ void OperatorLCAO<TK, TR>::contributeHk(int ik) {
     }
     ModuleBase::timer::tick("OperatorLCAO", "contributeHk");
 }
+// contributeHk()
+template <typename TK, typename TR>
+void OperatorLCAO<TK, TR>::contributeHk(int ik) {
+    ModuleBase::TITLE("OperatorLCAO", "contributeHk");
+    ModuleBase::timer::tick("OperatorLCAO", "contributeHk");
+    if(ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER(PARAM.inp.ks_solver))
+    {
+        const int nrow = this->hsk->get_pv()->get_row_size();
+        if(PARAM.inp.td_stype == 2)
+        {
+            TD_info::td_vel_op->folding_HR_td(*this->hR, this->hsk->get_hk(), this->kvec_d[ik], nrow, 1);
+        }
+        else
+        {
+            hamilt::folding_HR(*this->hR, this->hsk->get_hk(), this->kvec_d[ik], nrow, 1);
+        }
+    }
+    else
+    {
+        const int ncol = this->hsk->get_pv()->get_col_size();
+        if(PARAM.inp.td_stype == 2)
+        {
+            TD_info::td_vel_op->folding_HR_td(*this->hR, this->hsk->get_hk(), this->kvec_d[ik], ncol, 0);
+        }
+        else
+        {
+            hamilt::folding_HR(*this->hR, this->hsk->get_hk(), this->kvec_d[ik], ncol, 0);
+        }
+    }
+    ModuleBase::timer::tick("OperatorLCAO", "contributeHk");
+}
 
 template class OperatorLCAO<double, double>;
 template class OperatorLCAO<std::complex<double>, double>;

source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/overlap_new.cpp

@@ -7,6 +7,7 @@
 #include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h"
 #include "module_hamilt_lcao/module_hcontainer/hcontainer_funcs.h"
 #include <vector>
+#include "module_hamilt_lcao/module_tddft/td_info.h"
 
 template <typename TK, typename TR>
 hamilt::OverlapNew<hamilt::OperatorLCAO<TK, TR>>::OverlapNew(HS_Matrix_K<TK>* hsk_in,
@@ -188,11 +189,11 @@ void hamilt::OverlapNew<hamilt::OperatorLCAO<TK, TR>>::contributeHR()
 }
 
 // contributeHk()
-template <typename TK, typename TR>
-void hamilt::OverlapNew<hamilt::OperatorLCAO<TK, TR>>::contributeHk(int ik)
+template <>
+void hamilt::OverlapNew<hamilt::OperatorLCAO<double, double>>::contributeHk(int ik)
 {
     //! if k vector is not changed, then do nothing and return, only for gamma_only case
-    if (this->kvec_d[ik] == this->kvec_d_old && std::is_same<TK, double>::value)
+    if (this->kvec_d[ik] == this->kvec_d_old)
     {
         return;
     }
@@ -217,7 +218,44 @@ void hamilt::OverlapNew<hamilt::OperatorLCAO<TK, TR>>::contributeHk(int ik)
 
     ModuleBase::timer::tick("OverlapNew", "contributeHk");
 }
+template <typename TK, typename TR>
+void hamilt::OverlapNew<hamilt::OperatorLCAO<TK, TR>>::contributeHk(int ik)
+{
+    ModuleBase::TITLE("OverlapNew", "contributeHk");
+    ModuleBase::timer::tick("OverlapNew", "contributeHk");
+
+    //! set SK to zero and then calculate SK for each k vector
+    this->hsk->set_zero_sk();
+    if (ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER(PARAM.inp.ks_solver))
+    {
+        const int nrow = this->SR->get_atom_pair(0).get_paraV()->get_row_size();
+        if(PARAM.inp.td_stype == 2)
+        {
+            TD_info::td_vel_op->folding_HR_td(*this->SR, this->hsk->get_sk(), this->kvec_d[ik], nrow, 1);
+        }
+        else
+        {
+            hamilt::folding_HR(*this->SR, this->hsk->get_sk(), this->kvec_d[ik], nrow, 1);
+        }
+    }
+    else
+    {
+        const int ncol = this->SR->get_atom_pair(0).get_paraV()->get_col_size();
+        if(PARAM.inp.td_stype == 2)
+        {
+            TD_info::td_vel_op->folding_HR_td(*this->SR, this->hsk->get_sk(), this->kvec_d[ik], ncol, 0);
+        }
+        else
+        {
+            hamilt::folding_HR(*this->SR, this->hsk->get_sk(), this->kvec_d[ik], ncol, 0);
+        }
+    }
+
+    // update kvec_d_old
+    this->kvec_d_old = this->kvec_d[ik];
 
+    ModuleBase::timer::tick("OverlapNew", "contributeHk");
+}
 template <typename TK, typename TR>
 TK* hamilt::OverlapNew<hamilt::OperatorLCAO<TK, TR>>::getSk()
 {