Fix HR output in tddft

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\

source/module_elecstate/potentials/H_TDDFT_pw.h

@@ -133,4 +133,4 @@ class H_TDDFT_pw : public PotBase
 
 } // namespace elecstate
 
-#endif
+#endif

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"
 

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>>>;
 
-}
+}

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

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;
     }

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
 {

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)
         {

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, 

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(

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);
+}