Simplify matrix multiplication and add new propagator method

docs/advanced/input_files/input-main.md

@@ -282,6 +282,8 @@
 		- [td\_edm](#td_edm)
 		- [td\_print\_eij](#td_print_eij)
 		- [td\_force\_dt](#td_force_dt)
+		- [td\_htype](#td_htype)
+		- [propagator](#propagator)
 		- [td\_vext](#td_vext)
 		- [td\_vext\_dire](#td_vext_dire)
 		- [td\_stype](#td_stype)
@@ -2321,6 +2323,25 @@ These variables are used to control berry phase and wannier90 interface paramete
 - **Description**: Time-dependent evolution force changes time step. (fs)
 - **Default**: 0.02
 
+### td_htype
+
+- **Type**: Integer
+- **Description**:
+  type of Hamiltonian to calculate propagator (for propagator = 0 or 1).
+  - 0: H(t+dt).
+  - 1: H(t+dt/2).
+- **Default**: 0
+
+### propagator
+
+- **Type**: Integer
+- **Description**:
+  method of propagator.
+  - 0: Crank-Nicolson.
+  - 1: 4th Taylor expansions of exponential.
+  - 2: enforced time-reversal symmetry (ETRS).
+- **Default**: 0
+
 ### td_vext
 
 - **Type**: Boolean

source/module_esolver/esolver_ks_lcao_tddft.cpp

@@ -1,11 +1,11 @@
 #include "esolver_ks_lcao_tddft.h"
 
 #include "module_io/cal_r_overlap_R.h"
+#include "module_io/dipole_io.h"
 #include "module_io/dm_io.h"
 #include "module_io/rho_io.h"
-#include "module_io/dipole_io.h"
-#include "module_io/write_HS_R.h"
 #include "module_io/write_HS.h"
+#include "module_io/write_HS_R.h"
 
 //--------------temporary----------------------------
 #include "module_base/blas_connector.h"
@@ -41,9 +41,14 @@ ESolver_KS_LCAO_TDDFT::~ESolver_KS_LCAO_TDDFT()
 {
     // this->orb_con.clear_after_ions(GlobalC::UOT, GlobalC::ORB, GlobalV::deepks_setorb, GlobalC::ucell.infoNL.nproj);
     delete psi_laststep;
-    for(int ik = 0; ik < GlobalC::kv.nks; ++ik){
-        delete Hk_laststep[ik];}
+    if (Hk_laststep != nullptr)
+    {
+        for (int ik = 0; ik < GlobalC::kv.nks; ++ik)
+        {
+            delete Hk_laststep[ik];
+        }
         delete Hk_laststep;
+    }
 }
 
 void ESolver_KS_LCAO_TDDFT::Init(Input& inp, UnitCell& ucell)
@@ -191,7 +196,15 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density(int istep, int iter, double ethr)
 
     if (GlobalV::ESOLVER_TYPE == "tddft" && istep >= 2 && !GlobalV::GAMMA_ONLY_LOCAL)
     {
-        ELEC_evolve::evolve_psi(istep, this->p_hamilt, this->LOWF, this->psi, this->psi_laststep, this->Hk_laststep, this->pelec_td->ekb);
+        ELEC_evolve::evolve_psi(istep,
+                                this->p_hamilt,
+                                this->LOWF,
+                                this->psi,
+                                this->psi_laststep,
+                                this->Hk_laststep,
+                                this->pelec_td->ekb,
+                                INPUT.td_htype,
+                                INPUT.propagator);
         this->pelec_td->psiToRho_td(this->psi[0]);
         // this->pelec_td->psiToRho(this->psi[0]);
     }
@@ -277,7 +290,7 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density(int istep, int iter, double ethr)
 
 void ESolver_KS_LCAO_TDDFT::updatepot(const int istep, const int iter)
 {
-    //print Hamiltonian and Overlap matrix
+    // print Hamiltonian and Overlap matrix
     if (this->conv_elec)
     {
         if (!GlobalV::GAMMA_ONLY_LOCAL)
@@ -286,7 +299,7 @@ void ESolver_KS_LCAO_TDDFT::updatepot(const int istep, const int iter)
         }
         for (int ik = 0; ik < GlobalC::kv.nks; ++ik)
         {
-            if(hsolver::HSolverLCAO::out_mat_hs) 
+            if (hsolver::HSolverLCAO::out_mat_hs)
             {
                 this->p_hamilt->updateHk(ik);
             }
@@ -345,7 +358,8 @@ void ESolver_KS_LCAO_TDDFT::updatepot(const int istep, const int iter)
 
     if (!this->conv_elec)
     {
-        if(GlobalV::NSPIN==4) GlobalC::ucell.cal_ux();
+        if (GlobalV::NSPIN == 4)
+            GlobalC::ucell.cal_ux();
         this->pelec->pot->update_from_charge(this->pelec->charge, &GlobalC::ucell);
         GlobalC::en.delta_escf(this->pelec);
     }
@@ -367,11 +381,17 @@ void ESolver_KS_LCAO_TDDFT::updatepot(const int istep, const int iter)
             this->psi_laststep
                 = new psi::Psi<std::complex<double>>(GlobalC::kv.nks, GlobalV::NBANDS, GlobalV::NLOCAL, nullptr);
 #endif
-        if (this->Hk_laststep == nullptr){
-            this->Hk_laststep = new std::complex<double>*[GlobalC::kv.nks];
-            for(int ik = 0; ik < GlobalC::kv.nks; ++ik){
-                this->Hk_laststep[ik] = new std::complex<double>[this->LOC.ParaV->nloc];
-                ModuleBase::GlobalFunc::ZEROS(Hk_laststep[ik], this->LOC.ParaV->nloc);
+
+        if (INPUT.td_htype == 1)
+        {
+            if (this->Hk_laststep == nullptr)
+            {
+                this->Hk_laststep = new std::complex<double>*[GlobalC::kv.nks];
+                for (int ik = 0; ik < GlobalC::kv.nks; ++ik)
+                {
+                    this->Hk_laststep[ik] = new std::complex<double>[this->LOC.ParaV->nloc];
+                    ModuleBase::GlobalFunc::ZEROS(Hk_laststep[ik], this->LOC.ParaV->nloc);
+                }
             }
         }
 
@@ -382,14 +402,16 @@ void ESolver_KS_LCAO_TDDFT::updatepot(const int istep, const int iter)
             int size0 = psi->get_nbands() * psi->get_nbasis();
             for (int index = 0; index < size0; ++index)
                 psi_laststep[0].get_pointer()[index] = psi[0].get_pointer()[index];
-            //use for new propagator
-            this->p_hamilt->updateHk(ik);
-            //store Hk laststep
-            hamilt::MatrixBlock<complex<double>> h_mat, s_mat;
-            this->p_hamilt->matrix(h_mat, s_mat);
-            BlasConnector::copy(this->LOC.ParaV->nloc, h_mat.p, 1, Hk_laststep[ik], 1);
+
+            // store Hamiltonian
+            if (INPUT.td_htype == 1)
+            {
+                this->p_hamilt->updateHk(ik);
+                hamilt::MatrixBlock<complex<double>> h_mat, s_mat;
+                this->p_hamilt->matrix(h_mat, s_mat);
+                BlasConnector::copy(this->LOC.ParaV->nloc, h_mat.p, 1, Hk_laststep[ik], 1);
+            }
         }
-
 
         if (istep > 1 && ELEC_evolve::td_edm == 0)
             this->cal_edm_tddft();
@@ -433,7 +455,7 @@ void ESolver_KS_LCAO_TDDFT::afterscf(const int istep)
         const int precision = 3;
         std::stringstream ssc;
         ssc << GlobalV::global_out_dir << "SPIN" << is + 1 << "_CHG.cube";
-        double& ef_tmp = GlobalC::en.get_ef(is,GlobalV::TWO_EFERMI);
+        double& ef_tmp = GlobalC::en.get_ef(is, GlobalV::TWO_EFERMI);
         ModuleIO::write_rho(
 #ifdef __MPI
             GlobalC::bigpw->bz,
@@ -522,7 +544,7 @@ void ESolver_KS_LCAO_TDDFT::afterscf(const int istep)
             std::cout << " !! CONVERGENCE HAS NOT BEEN ACHIEVED !!" << std::endl;
     }
 
-    if( GlobalV::CALCULATION != "md" || (istep % hsolver::HSolverLCAO::out_hsR_interval == 0))
+    if (GlobalV::CALCULATION != "md" || (istep % hsolver::HSolverLCAO::out_hsR_interval == 0))
     {
         if (hsolver::HSolverLCAO::out_mat_hsR)
         {
@@ -552,7 +574,6 @@ void ESolver_KS_LCAO_TDDFT::afterscf(const int istep)
             else
             {
                 r_matrix.out_rR(istep);
-
             }
         }
     }
@@ -572,14 +593,12 @@ void ESolver_KS_LCAO_TDDFT::cal_edm_tddft()
         complex<double>* tmp2 = new complex<double>[this->LOC.ParaV->nloc];
         complex<double>* tmp3 = new complex<double>[this->LOC.ParaV->nloc];
         complex<double>* tmp4 = new complex<double>[this->LOC.ParaV->nloc];
-        complex<double>* tmp5 = new complex<double>[this->LOC.ParaV->nloc];
         ModuleBase::GlobalFunc::ZEROS(Htmp, this->LOC.ParaV->nloc);
         ModuleBase::GlobalFunc::ZEROS(Sinv, this->LOC.ParaV->nloc);
         ModuleBase::GlobalFunc::ZEROS(tmp1, this->LOC.ParaV->nloc);
         ModuleBase::GlobalFunc::ZEROS(tmp2, this->LOC.ParaV->nloc);
         ModuleBase::GlobalFunc::ZEROS(tmp3, this->LOC.ParaV->nloc);
         ModuleBase::GlobalFunc::ZEROS(tmp4, this->LOC.ParaV->nloc);
-        ModuleBase::GlobalFunc::ZEROS(tmp5, this->LOC.ParaV->nloc);
         const int inc = 1;
         int nrow = this->LOC.ParaV->nrow;
         int ncol = this->LOC.ParaV->ncol;
@@ -629,8 +648,8 @@ void ESolver_KS_LCAO_TDDFT::cal_edm_tddft()
         const char N_char = 'N', T_char = 'T';
         const complex<double> one_float = {1.0, 0.0}, zero_float = {0.0, 0.0};
         const complex<double> half_float = {0.5, 0.0};
-        pzgemm_(&T_char,
-                &T_char,
+        pzgemm_(&N_char,
+                &N_char,
                 &GlobalV::NLOCAL,
                 &GlobalV::NLOCAL,
                 &GlobalV::NLOCAL,
@@ -670,7 +689,7 @@ void ESolver_KS_LCAO_TDDFT::cal_edm_tddft()
                 this->LOC.ParaV->desc);
 
         pzgemm_(&N_char,
-                &T_char,
+                &N_char,
                 &GlobalV::NLOCAL,
                 &GlobalV::NLOCAL,
                 &GlobalV::NLOCAL,
@@ -690,7 +709,7 @@ void ESolver_KS_LCAO_TDDFT::cal_edm_tddft()
                 this->LOC.ParaV->desc);
 
         pzgemm_(&N_char,
-                &T_char,
+                &N_char,
                 &GlobalV::NLOCAL,
                 &GlobalV::NLOCAL,
                 &GlobalV::NLOCAL,
@@ -723,27 +742,14 @@ void ESolver_KS_LCAO_TDDFT::cal_edm_tddft()
                  &one_int,
                  this->LOC.ParaV->desc);
 
-        pztranu_(&GlobalV::NLOCAL,
-                 &GlobalV::NLOCAL,
-                 &one_float,
-                 tmp4,
-                 &one_int,
-                 &one_int,
-                 this->LOC.ParaV->desc,
-                 &zero_float,
-                 tmp5,
-                 &one_int,
-                 &one_int,
-                 this->LOC.ParaV->desc);
-        zcopy_(&this->LOC.ParaV->nloc, tmp5, &inc, this->LOC.edm_k_tddft[ik].c, &inc);
+        zcopy_(&this->LOC.ParaV->nloc, tmp4, &inc, this->LOC.edm_k_tddft[ik].c, &inc);
 
         delete[] Htmp;
         delete[] Sinv;
         delete[] tmp1;
         delete[] tmp2;
         delete[] tmp3;
         delete[] tmp4;
-        delete[] tmp5;
         delete[] ipiv;
 #else
         this->LOC.edm_k_tddft[ik].create(this->LOC.ParaV->ncol, this->LOC.ParaV->nrow);

source/module_hamilt_lcao/module_tddft/ELEC_evolve.cpp

@@ -1,12 +1,12 @@
 #include "ELEC_evolve.h"
 
-#include "module_base/timer.h"
+#include "LCAO_evolve.h"
 #include "module_base/parallel_reduce.h"
-#include "module_hamilt_pw/hamilt_pwdft/global.h"
+#include "module_base/timer.h"
 #include "module_elecstate/module_charge/symmetry_rho.h"
-#include "LCAO_evolve.h"
-#include "module_hamilt_lcao/module_dftu/dftu.h"
 #include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h"
+#include "module_hamilt_lcao/module_dftu/dftu.h"
+#include "module_hamilt_pw/hamilt_pwdft/global.h"
 
 ELEC_evolve::ELEC_evolve(){};
 ELEC_evolve::~ELEC_evolve(){};
@@ -26,7 +26,9 @@ void ELEC_evolve::evolve_psi(const int& istep,
                              psi::Psi<std::complex<double>>* psi,
                              psi::Psi<std::complex<double>>* psi_laststep,
                              std::complex<double>** Hk_laststep,
-                             ModuleBase::matrix& ekb)
+                             ModuleBase::matrix& ekb,
+                             int htype,
+                             int propagator)
 {
     ModuleBase::TITLE("ELEC_evolve", "eveolve_psi");
     ModuleBase::timer::tick("ELEC_evolve", "evolve_psi");
@@ -40,7 +42,14 @@ void ELEC_evolve::evolve_psi(const int& istep,
         Evolve_LCAO_Matrix ELM(lowf.ParaV);
         psi->fix_k(ik);
         psi_laststep->fix_k(ik);
-        ELM.evolve_complex_matrix(ik, phm, psi, psi_laststep, Hk_laststep[ik], &(ekb(ik, 0)));
+        if (Hk_laststep == nullptr)
+        {
+            ELM.evolve_complex_matrix(ik, phm, psi, psi_laststep, nullptr, &(ekb(ik, 0)), htype, propagator);
+        }
+        else
+        {
+            ELM.evolve_complex_matrix(ik, phm, psi, psi_laststep, Hk_laststep[ik], &(ekb(ik, 0)), htype, propagator);
+        }
         ModuleBase::timer::tick("Efficience", "evolve_k");
     } // end k
 

source/module_hamilt_lcao/module_tddft/ELEC_evolve.h

@@ -45,7 +45,9 @@ class ELEC_evolve
                            psi::Psi<std::complex<double>>* psi,
                            psi::Psi<std::complex<double>>* psi_laststep,
                            std::complex<double>** Hk_laststep,
-                           ModuleBase::matrix& ekb);
+                           ModuleBase::matrix& ekb,
+                           int htype,
+                           int propagator);
 };
 
 #endif