diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index d0430de33c..743698459d 100644 --- a/docs/advanced/input_files/input-main.md +++ b/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 diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp index 4276cb5fb3..7e4f9ad8d0 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.cpp +++ b/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>(GlobalC::kv.nks, GlobalV::NBANDS, GlobalV::NLOCAL, nullptr); #endif - if (this->Hk_laststep == nullptr){ - this->Hk_laststep = new std::complex*[GlobalC::kv.nks]; - for(int ik = 0; ik < GlobalC::kv.nks; ++ik){ - this->Hk_laststep[ik] = new std::complex[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*[GlobalC::kv.nks]; + for (int ik = 0; ik < GlobalC::kv.nks; ++ik) + { + this->Hk_laststep[ik] = new std::complex[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> 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> 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* tmp2 = new complex[this->LOC.ParaV->nloc]; complex* tmp3 = new complex[this->LOC.ParaV->nloc]; complex* tmp4 = new complex[this->LOC.ParaV->nloc]; - complex* tmp5 = new complex[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 one_float = {1.0, 0.0}, zero_float = {0.0, 0.0}; const complex 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,19 +742,7 @@ 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; @@ -743,7 +750,6 @@ void ESolver_KS_LCAO_TDDFT::cal_edm_tddft() 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); diff --git a/source/module_hamilt_lcao/module_tddft/ELEC_evolve.cpp b/source/module_hamilt_lcao/module_tddft/ELEC_evolve.cpp index 3cfdd30b1d..4b37b54197 100644 --- a/source/module_hamilt_lcao/module_tddft/ELEC_evolve.cpp +++ b/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>* psi, psi::Psi>* psi_laststep, std::complex** 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 diff --git a/source/module_hamilt_lcao/module_tddft/ELEC_evolve.h b/source/module_hamilt_lcao/module_tddft/ELEC_evolve.h index eb647b83c3..dfa2b904eb 100644 --- a/source/module_hamilt_lcao/module_tddft/ELEC_evolve.h +++ b/source/module_hamilt_lcao/module_tddft/ELEC_evolve.h @@ -45,7 +45,9 @@ class ELEC_evolve psi::Psi>* psi, psi::Psi>* psi_laststep, std::complex** Hk_laststep, - ModuleBase::matrix& ekb); + ModuleBase::matrix& ekb, + int htype, + int propagator); }; #endif diff --git a/source/module_hamilt_lcao/module_tddft/LCAO_evolve.cpp b/source/module_hamilt_lcao/module_tddft/LCAO_evolve.cpp index 49b190f6af..cd7649f536 100644 --- a/source/module_hamilt_lcao/module_tddft/LCAO_evolve.cpp +++ b/source/module_hamilt_lcao/module_tddft/LCAO_evolve.cpp @@ -1,12 +1,12 @@ #include "LCAO_evolve.h" -#include "module_io/input.h" -#include "module_base/scalapack_connector.h" -#include "module_base/lapack_connector.h" -#include "module_io/write_HS.h" -#include "module_hamilt_pw/hamilt_pwdft/global.h" #include "ELEC_evolve.h" +#include "module_base/lapack_connector.h" +#include "module_base/scalapack_connector.h" #include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h" +#include "module_hamilt_pw/hamilt_pwdft/global.h" +#include "module_io/input.h" +#include "module_io/write_HS.h" #include // fuxiang add 2016-10-28 @@ -27,8 +27,10 @@ void Evolve_LCAO_Matrix::evolve_complex_matrix(const int& ik, hamilt::Hamilt* p_hamilt, psi::Psi>* psi_k, psi::Psi>* psi_k_laststep, - std::complex* Hk_laststep, - double* ekb) const + std::complex* H_laststep, + double* ekb, + int htype, + int propagator) const { ModuleBase::TITLE("Evolve_LCAO_Matrix", "evolve_complex_matrix"); time_t time_start = time(NULL); @@ -37,7 +39,15 @@ void Evolve_LCAO_Matrix::evolve_complex_matrix(const int& ik, if (GlobalV::ESOLVER_TYPE == "tddft") { #ifdef __MPI - this->using_ScaLAPACK_complex(GlobalV::NBANDS, GlobalV::NLOCAL, psi_k_laststep[0].get_pointer(), Hk_laststep, p_hamilt, psi_k[0].get_pointer(), ekb); + this->using_ScaLAPACK_complex(GlobalV::NBANDS, + GlobalV::NLOCAL, + psi_k_laststep[0].get_pointer(), + H_laststep, + p_hamilt, + psi_k[0].get_pointer(), + ekb, + htype, + propagator); #else this->using_LAPACK_complex(ik, p_hamilt, psi_k[0].get_pointer(), psi_k_laststep[0].get_pointer(), ekb); #endif @@ -341,14 +351,15 @@ void Evolve_LCAO_Matrix::using_LAPACK_complex(const int& ik, #ifdef __MPI -void Evolve_LCAO_Matrix::using_ScaLAPACK_complex( - const int nband, - const int nlocal, - const std::complex* psi_k_laststep, - const std::complex* Hk_laststep, - hamilt::Hamilt* p_hamilt, - std::complex* psi_k, - double* ekb) const +void Evolve_LCAO_Matrix::using_ScaLAPACK_complex(const int nband, + const int nlocal, + const std::complex* psi_k_laststep, + const std::complex* H_laststep, + hamilt::Hamilt* p_hamilt, + std::complex* psi_k, + double* ekb, + int htype, + int propagator) const { ModuleBase::TITLE("Evolve_LCAO_Matrix", "using_ScaLAPACK_complex"); @@ -364,83 +375,68 @@ void Evolve_LCAO_Matrix::using_ScaLAPACK_complex( ModuleBase::GlobalFunc::ZEROS(Htmp, this->ParaV->nloc); BlasConnector::copy(this->ParaV->nloc, h_mat.p, 1, Htmp, 1); - //print test - /*GlobalV::ofs_running << " H(t) matrix :" << endl; - for (int i = 0; i < this->ParaV->nrow; i++) - { - for (int j = 0; j < this->ParaV->ncol; j++) - { - GlobalV::ofs_running << Hk_laststep[i * this->ParaV->ncol + j].real() << "+" << Hk_laststep[i * this->ParaV->ncol + j].imag() << "i "; - } - GlobalV::ofs_running << endl;} - GlobalV::ofs_running << " H(t+dt) matrix :" << endl; - for (int i = 0; i < this->ParaV->nrow; i++) - { - for (int j = 0; j < this->ParaV->ncol; j++) - { - GlobalV::ofs_running << Htmp[i * this->ParaV->ncol + j].real() << "+" << Htmp[i * this->ParaV->ncol + j].imag() << "i "; - } - GlobalV::ofs_running << endl;} - // -*/ - - complex mixing = {0.5, 0.0}; - ScalapackConnector::geadd( - 'N', - nlocal, - nlocal, - mixing, - Hk_laststep, - 1, - 1, - this->ParaV->desc, - mixing, - Htmp, - 1, - 1, - this->ParaV->desc - ); - /*GlobalV::ofs_running << " H(t+dt/2) matrix :" << endl; - for (int i = 0; i < this->ParaV->nrow; i++) - { - for (int j = 0; j < this->ParaV->ncol; j++) - { - GlobalV::ofs_running << Htmp[i * this->ParaV->ncol + j].real() << "+" << Htmp[i * this->ParaV->ncol + j].imag() << "i "; - } - GlobalV::ofs_running << endl;}*/ + complex* Hold = new complex[this->ParaV->nloc]; + ModuleBase::GlobalFunc::ZEROS(Hold, this->ParaV->nloc); + BlasConnector::copy(this->ParaV->nloc, h_mat.p, 1, Hold, 1); + complex* U_operator = new complex[this->ParaV->nloc]; ModuleBase::GlobalFunc::ZEROS(U_operator, this->ParaV->nloc); -// (1)->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // (1)->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + + /// @brief compute H(t+dt/2) + /// @input H_laststep, Htmp, print_matrix + /// @output Htmp + if (htype == 1 && propagator!=2) + { + half_Hmatrix(nband, nlocal, Htmp, H_laststep, print_matrix); + } + + // (2)->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> /// @brief compute U_operator /// @input Stmp, Htmp, print_matrix /// @output U_operator - compute_U_operator(nband, nlocal, Stmp, Htmp, U_operator, print_matrix); + switch (propagator) + { + case 0: + compute_U_operator_CN2(nband, nlocal, Stmp, Htmp, U_operator, print_matrix); + break; + + case 1: + compute_U_operator_taylor(nband, nlocal, Stmp, Htmp, U_operator, print_matrix); + break; + + case 2: + compute_U_operator_etrs(nband, nlocal, Stmp, Htmp, H_laststep, U_operator, print_matrix); + + break; -// (2)->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + default: + std::cout << "method of propagator is wrong" << endl; + break; + } + + // (3)->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> /// @brief apply U_operator to the wave function of the previous step for new wave function - /// @input U_operator, psi_k_laststep + /// @input U_operator, psi_k_laststep, print_matrix /// @output psi_k - U_to_wfc(nband, nlocal, U_operator, psi_k_laststep, psi_k); - - + U_to_wfc(nband, nlocal, U_operator, psi_k_laststep, psi_k, print_matrix); -// (3)->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // (4)->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> /// @brief normalize psi_k /// @input Stmp, psi_not_norm, psi_k, print_matrix /// @output psi_k norm_wfc(nband, nlocal, Stmp, psi_k, print_matrix); + // (5)->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> -// (4)->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> - /// @brief compute ekb /// @input Htmp, psi_k /// @output ekb - compute_ekb(nband, nlocal, Htmp, psi_k, ekb); + compute_ekb(nband, nlocal, Hold, psi_k, ekb); delete[] Stmp; delete[] Htmp; @@ -448,14 +444,80 @@ void Evolve_LCAO_Matrix::using_ScaLAPACK_complex( return; } +void Evolve_LCAO_Matrix::half_Hmatrix(const int nband, + const int nlocal, + std::complex* Htmp, + const std::complex* H_laststep, + const int print_matrix) const +{ + if (print_matrix) + { + GlobalV::ofs_running << std::setprecision(10); + GlobalV::ofs_running << endl; + GlobalV::ofs_running << " H(t+dt) :" << endl; + for (int i = 0; i < this->ParaV->nrow; i++) + { + for (int j = 0; j < this->ParaV->ncol; j++) + { + GlobalV::ofs_running << Htmp[i * this->ParaV->ncol + j].real() << "+" + << Htmp[i * this->ParaV->ncol + j].imag() << "i "; + } + GlobalV::ofs_running << endl; + } + GlobalV::ofs_running << endl; + GlobalV::ofs_running << endl; + GlobalV::ofs_running << " H(t):" << endl; + for (int i = 0; i < this->ParaV->nrow; i++) + { + for (int j = 0; j < this->ParaV->ncol; j++) + { + GlobalV::ofs_running << H_laststep[i * this->ParaV->ncol + j].real() << "+" + << H_laststep[i * this->ParaV->ncol + j].imag() << "i "; + } + GlobalV::ofs_running << endl; + } + GlobalV::ofs_running << endl; + } + + complex alpha = {0.5, 0.0}; + complex beta = {0.5, 0.0}; + ScalapackConnector::geadd('N', + nlocal, + nlocal, + alpha, + H_laststep, + 1, + 1, + this->ParaV->desc, + beta, + Htmp, + 1, + 1, + this->ParaV->desc); -void Evolve_LCAO_Matrix::compute_U_operator( - const int nband, - const int nlocal, - const std::complex* Stmp, - const std::complex* Htmp, - std::complex* U_operator, - const int print_matrix) const + if (print_matrix) + { + GlobalV::ofs_running << endl; + GlobalV::ofs_running << " H (t+dt/2) :" << endl; + for (int i = 0; i < this->ParaV->nrow; i++) + { + for (int j = 0; j < this->ParaV->ncol; j++) + { + GlobalV::ofs_running << Htmp[i * this->ParaV->ncol + j].real() << "+" + << Htmp[i * this->ParaV->ncol + j].imag() << "i "; + } + GlobalV::ofs_running << endl; + } + GlobalV::ofs_running << endl; + } +} + +void Evolve_LCAO_Matrix::compute_U_operator_CN2(const int nband, + const int nlocal, + const std::complex* Stmp, + const std::complex* Htmp, + std::complex* U_operator, + const int print_matrix) const { // (1) copy Htmp to Numerator & Denominator complex* Numerator = new complex[this->ParaV->nloc]; @@ -474,7 +536,8 @@ void Evolve_LCAO_Matrix::compute_U_operator( { for (int j = 0; j < this->ParaV->ncol; j++) { - GlobalV::ofs_running << Stmp[i * this->ParaV->ncol + j].real() << "+" << Stmp[i * this->ParaV->ncol + j].imag() << "i "; + GlobalV::ofs_running << Stmp[i * this->ParaV->ncol + j].real() << "+" + << Stmp[i * this->ParaV->ncol + j].imag() << "i "; } GlobalV::ofs_running << endl; } @@ -485,14 +548,15 @@ void Evolve_LCAO_Matrix::compute_U_operator( { for (int j = 0; j < this->ParaV->ncol; j++) { - GlobalV::ofs_running << Numerator[i * this->ParaV->ncol + j].real() << "+" << Numerator[i * this->ParaV->ncol + j].imag() << "i "; + GlobalV::ofs_running << Numerator[i * this->ParaV->ncol + j].real() << "+" + << Numerator[i * this->ParaV->ncol + j].imag() << "i "; } GlobalV::ofs_running << endl; } GlobalV::ofs_running << endl; } -// ->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // ->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> // (2) compute Numerator & Denominator by GEADD // Numerator = Stmp - i*para * Htmp; beta1 = - para = -0.25 * INPUT.mdp.md_dt // Denominator = Stmp + i*para * Htmp; beta2 = para = 0.25 * INPUT.mdp.md_dt @@ -500,128 +564,111 @@ void Evolve_LCAO_Matrix::compute_U_operator( complex beta1 = {0.0, -0.25 * INPUT.mdp.md_dt}; complex beta2 = {0.0, 0.25 * INPUT.mdp.md_dt}; - ScalapackConnector::geadd( - 'N', - nlocal, - nlocal, - alpha, - Stmp, - 1, - 1, - this->ParaV->desc, - beta1, - Numerator, - 1, - 1, - this->ParaV->desc - ); - ScalapackConnector::geadd( - 'N', - nlocal, - nlocal, - alpha, - Stmp, - 1, - 1, - this->ParaV->desc, - beta2, - Denominator, - 1, - 1, - this->ParaV->desc - ); + ScalapackConnector::geadd('N', + nlocal, + nlocal, + alpha, + Stmp, + 1, + 1, + this->ParaV->desc, + beta1, + Numerator, + 1, + 1, + this->ParaV->desc); + ScalapackConnector::geadd('N', + nlocal, + nlocal, + alpha, + Stmp, + 1, + 1, + this->ParaV->desc, + beta2, + Denominator, + 1, + 1, + this->ParaV->desc); if (print_matrix) { + GlobalV::ofs_running << " beta=" << beta1 << endl; GlobalV::ofs_running << " fenmu:" << endl; for (int i = 0; i < this->ParaV->nrow; i++) { for (int j = 0; j < this->ParaV->ncol; j++) { - GlobalV::ofs_running << Denominator[i * this->ParaV->ncol + j].real() << "+" << Denominator[i * this->ParaV->ncol + j].imag() << "i "; + GlobalV::ofs_running << Denominator[i * this->ParaV->ncol + j].real() << "+" + << Denominator[i * this->ParaV->ncol + j].imag() << "i "; } GlobalV::ofs_running << endl; } GlobalV::ofs_running << endl; } -//->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + //->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> // (3) Next, invert Denominator int* ipiv = new int[this->ParaV->nloc]; int info = 0; // (3.1) compute ipiv - ScalapackConnector::getrf( - nlocal, - nlocal, - Denominator, - 1, - 1, - this->ParaV->desc, - ipiv, - &info - ); + ScalapackConnector::getrf(nlocal, nlocal, Denominator, 1, 1, this->ParaV->desc, ipiv, &info); int lwork = -1; int liwotk = -1; std::vector> work(1, 0); std::vector iwork(1, 0); // (3.2) compute work - ScalapackConnector::getri( - nlocal, - Denominator, - 1, - 1, - this->ParaV->desc, - ipiv, - work.data(), - &lwork, - iwork.data(), - &liwotk, - &info - ); + ScalapackConnector::getri(nlocal, + Denominator, + 1, + 1, + this->ParaV->desc, + ipiv, + work.data(), + &lwork, + iwork.data(), + &liwotk, + &info); lwork = work[0].real(); work.resize(lwork, 0); liwotk = iwork[0]; iwork.resize(liwotk, 0); - // (3.3) compute inverse matrix of matrix_A - ScalapackConnector::getri( - nlocal, - Denominator, - 1, - 1, - this->ParaV->desc, - ipiv, - work.data(), - &lwork, - iwork.data(), - &liwotk, - &info - ); + // (3.3) compute inverse matrix of Denominator + ScalapackConnector::getri(nlocal, + Denominator, + 1, + 1, + this->ParaV->desc, + ipiv, + work.data(), + &lwork, + iwork.data(), + &liwotk, + &info); assert(0 == info); -//->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + //->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> // (4) U_operator = Denominator * Numerator; - ScalapackConnector::gemm( - 'N', - 'N', - nlocal, - nlocal, - nlocal, - 1.0, - Denominator, - 1, - 1, - this->ParaV->desc, - Numerator, - 1, - 1, - this->ParaV->desc, - 0.0, - U_operator, - 1, - 1, - this->ParaV->desc - ); + ScalapackConnector::gemm('N', + 'N', + nlocal, + nlocal, + nlocal, + 1.0, + Denominator, + 1, + 1, + this->ParaV->desc, + Numerator, + 1, + 1, + this->ParaV->desc, + 0.0, + U_operator, + 1, + 1, + this->ParaV->desc); if (print_matrix) { @@ -630,7 +677,8 @@ void Evolve_LCAO_Matrix::compute_U_operator( { for (int j = 0; j < this->ParaV->ncol; j++) { - GlobalV::ofs_running << Denominator[i * this->ParaV->ncol + j].real() << "+" << Denominator[i * this->ParaV->ncol + j].imag() << "i "; + GlobalV::ofs_running << Denominator[i * this->ParaV->ncol + j].real() << "+" + << Denominator[i * this->ParaV->ncol + j].imag() << "i "; } GlobalV::ofs_running << endl; } @@ -640,7 +688,8 @@ void Evolve_LCAO_Matrix::compute_U_operator( { for (int j = 0; j < this->ParaV->ncol; j++) { - GlobalV::ofs_running << Numerator[i * this->ParaV->ncol + j].real() << "+" << Numerator[i * this->ParaV->ncol + j].imag() << "i "; + GlobalV::ofs_running << Numerator[i * this->ParaV->ncol + j].real() << "+" + << Numerator[i * this->ParaV->ncol + j].imag() << "i "; } GlobalV::ofs_running << endl; } @@ -663,103 +712,513 @@ void Evolve_LCAO_Matrix::compute_U_operator( } } -// cout << "U_operator Success!!!" <* U_operator, - const std::complex* psi_k_laststep, - std::complex* psi_k) const +void Evolve_LCAO_Matrix::compute_U_operator_taylor(const int nband, + const int nlocal, + const std::complex* Stmp, + const std::complex* Htmp, + std::complex* U_operator, + const int print_matrix) const { + ModuleBase::GlobalFunc::ZEROS(U_operator, this->ParaV->nloc); + complex* A_matrix = new complex[this->ParaV->nloc]; + ModuleBase::GlobalFunc::ZEROS(A_matrix, this->ParaV->nloc); + complex* rank0 = new complex[this->ParaV->nloc]; + ModuleBase::GlobalFunc::ZEROS(rank0, this->ParaV->nloc); + complex* rank2 = new complex[this->ParaV->nloc]; + ModuleBase::GlobalFunc::ZEROS(rank2, this->ParaV->nloc); + complex* rank3 = new complex[this->ParaV->nloc]; + ModuleBase::GlobalFunc::ZEROS(rank3, this->ParaV->nloc); + complex* rank4 = new complex[this->ParaV->nloc]; + ModuleBase::GlobalFunc::ZEROS(rank4, this->ParaV->nloc); + complex* tmp1 = new complex[this->ParaV->nloc]; + ModuleBase::GlobalFunc::ZEROS(tmp1, this->ParaV->nloc); + complex* tmp2 = new complex[this->ParaV->nloc]; + ModuleBase::GlobalFunc::ZEROS(tmp2, this->ParaV->nloc); + complex* Sinv = new complex[this->ParaV->nloc]; + ModuleBase::GlobalFunc::ZEROS(Sinv, this->ParaV->nloc); + BlasConnector::copy(this->ParaV->nloc, Stmp, 1, Sinv, 1); + + if (print_matrix) + { + GlobalV::ofs_running << endl; + GlobalV::ofs_running << " S matrix :" << endl; + for (int i = 0; i < this->ParaV->nrow; i++) + { + for (int j = 0; j < this->ParaV->ncol; j++) + { + GlobalV::ofs_running << Stmp[i * this->ParaV->ncol + j].real() << "+" + << Stmp[i * this->ParaV->ncol + j].imag() << "i "; + } + GlobalV::ofs_running << endl; + } + GlobalV::ofs_running << endl; + GlobalV::ofs_running << endl; + GlobalV::ofs_running << " H matrix :" << endl; + for (int i = 0; i < this->ParaV->nrow; i++) + { + for (int j = 0; j < this->ParaV->ncol; j++) + { + GlobalV::ofs_running << Htmp[i * this->ParaV->ncol + j].real() << "+" + << Htmp[i * this->ParaV->ncol + j].imag() << "i "; + } + GlobalV::ofs_running << endl; + } + GlobalV::ofs_running << endl; + } + + // set rank0 + int info; + int myid; + MPI_Comm_rank(this->ParaV->comm_2D, &myid); + int naroc[2]; // maximum number of row or column + + for (int iprow = 0; iprow < this->ParaV->dim0; ++iprow) + { + for (int ipcol = 0; ipcol < this->ParaV->dim1; ++ipcol) + { + const int coord[2] = {iprow, ipcol}; + int src_rank; + info = MPI_Cart_rank(this->ParaV->comm_2D, coord, &src_rank); + if (myid == src_rank) + { + naroc[0] = this->ParaV->nrow; + naroc[1] = this->ParaV->ncol; + for (int j = 0; j < naroc[1]; ++j) + { + int igcol = globalIndex(j, this->ParaV->nb, this->ParaV->dim1, ipcol); + if (igcol >= nlocal) + continue; + for (int i = 0; i < naroc[0]; ++i) + { + int igrow = globalIndex(i, this->ParaV->nb, this->ParaV->dim0, iprow); + if (igrow >= nlocal) + continue; + if (igcol == igrow) + { + rank0[j * naroc[0] + i] = {1.0, 0.0}; + } + else + { + rank0[j * naroc[0] + i] = {0.0, 0.0}; + } + } + } + } + } // loop ipcol + } // loop iprow + + // complex beta = {0.0, -0.5 * INPUT.mdp.md_dt}; + complex beta = {0.0, -0.25 * INPUT.mdp.md_dt}; // for ETRS + + //->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // invert Stmp + int* ipiv = new int[this->ParaV->nloc]; + // (3.1) compute ipiv + ScalapackConnector::getrf(nlocal, nlocal, Sinv, 1, 1, this->ParaV->desc, ipiv, &info); + int lwork = -1; + int liwotk = -1; + std::vector> work(1, 0); + std::vector iwork(1, 0); + // (3.2) compute work + ScalapackConnector::getri(nlocal, + Sinv, + 1, + 1, + this->ParaV->desc, + ipiv, + work.data(), + &lwork, + iwork.data(), + &liwotk, + &info); + lwork = work[0].real(); + work.resize(lwork, 0); + liwotk = iwork[0]; + iwork.resize(liwotk, 0); + ScalapackConnector::getri(nlocal, + Sinv, + 1, + 1, + this->ParaV->desc, + ipiv, + work.data(), + &lwork, + iwork.data(), + &liwotk, + &info); + assert(0 == info); + + //->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + + // A_matrix = - idt S^-1 H ; + ScalapackConnector::gemm('N', + 'N', + nlocal, + nlocal, + nlocal, + beta, + Sinv, + 1, + 1, + this->ParaV->desc, + Htmp, + 1, + 1, + this->ParaV->desc, + 0.0, + U_operator, + 1, + 1, + this->ParaV->desc); + + // rank2 = A^2 ; + ScalapackConnector::gemm('N', + 'N', + nlocal, + nlocal, + nlocal, + 1.0, + U_operator, + 1, + 1, + this->ParaV->desc, + U_operator, + 1, + 1, + this->ParaV->desc, + 0.0, + rank2, + 1, + 1, + this->ParaV->desc); + + // rank3 = A^3 ; + ScalapackConnector::gemm('N', + 'N', + nlocal, + nlocal, + nlocal, + 1.0, + U_operator, + 1, + 1, + this->ParaV->desc, + rank2, + 1, + 1, + this->ParaV->desc, + 0.0, + rank3, + 1, + 1, + this->ParaV->desc); + + // rank4 = A^4 ; + ScalapackConnector::gemm('N', + 'N', + nlocal, + nlocal, + nlocal, + 1.0, + U_operator, + 1, + 1, + this->ParaV->desc, + rank3, + 1, + 1, + this->ParaV->desc, + 0.0, + rank4, + 1, + 1, + this->ParaV->desc); + + complex p1 = {1.0, 0.0}; + complex p2 = {1.0 / 2.0, 0.0}; + complex p3 = {1.0 / 6.0, 0.0}; + complex p4 = {1.0 / 24.0, 0.0}; + + ScalapackConnector::geadd('N', + nlocal, + nlocal, + p1, + rank0, + 1, + 1, + this->ParaV->desc, + p1, + U_operator, + 1, + 1, + this->ParaV->desc); + + ScalapackConnector::geadd('N', + nlocal, + nlocal, + p2, + rank2, + 1, + 1, + this->ParaV->desc, + p1, + U_operator, + 1, + 1, + this->ParaV->desc); + + ScalapackConnector::geadd('N', + nlocal, + nlocal, + p3, + rank3, + 1, + 1, + this->ParaV->desc, + p1, + U_operator, + 1, + 1, + this->ParaV->desc); + + ScalapackConnector::geadd('N', + nlocal, + nlocal, + p4, + rank4, + 1, + 1, + this->ParaV->desc, + p1, + U_operator, + 1, + 1, + this->ParaV->desc); + + if (print_matrix) + { + GlobalV::ofs_running << " A_matrix:" << endl; + for (int i = 0; i < this->ParaV->nrow; i++) + { + for (int j = 0; j < this->ParaV->ncol; j++) + { + GlobalV::ofs_running << A_matrix[i * this->ParaV->ncol + j].real() << "+" + << A_matrix[i * this->ParaV->ncol + j].imag() << "i "; + } + GlobalV::ofs_running << endl; + } + GlobalV::ofs_running << endl; + GlobalV::ofs_running << " U operator:" << endl; + for (int i = 0; i < this->ParaV->nrow; i++) + { + for (int j = 0; j < this->ParaV->ncol; j++) + { + double aa, bb; + aa = U_operator[i * this->ParaV->ncol + j].real(); + bb = U_operator[i * this->ParaV->ncol + j].imag(); + if (abs(aa) < 1e-8) + aa = 0.0; + if (abs(bb) < 1e-8) + bb = 0.0; + GlobalV::ofs_running << aa << "+" << bb << "i "; + } + GlobalV::ofs_running << endl; + } + } - ScalapackConnector::gemm( - 'N', - 'N', - nlocal, - nband, - nlocal, - 1.0, - U_operator, - 1, - 1, - this->ParaV->desc, - psi_k_laststep, - 1, - 1, - this->ParaV->desc_wfc, - 0.0, - psi_k, - 1, - 1, - this->ParaV->desc_wfc - ); + // cout << "U_operator Success!!!" <* Stmp, + const std::complex* Htmp, + const std::complex* H_laststep, + std::complex* U_operator, + const int print_matrix) const +{ + complex* U1 = new complex[this->ParaV->nloc]; + complex* U2 = new complex[this->ParaV->nloc]; + ModuleBase::GlobalFunc::ZEROS(U1, this->ParaV->nloc); + ModuleBase::GlobalFunc::ZEROS(U2, this->ParaV->nloc); + compute_U_operator_taylor(nband, nlocal, Stmp, Htmp, U1, print_matrix); + compute_U_operator_taylor(nband, nlocal, Stmp, H_laststep, U2, print_matrix); + ScalapackConnector::gemm('N', + 'N', + nlocal, + nlocal, + nlocal, + 1.0, + U1, + 1, + 1, + this->ParaV->desc, + U2, + 1, + 1, + this->ParaV->desc, + 0.0, + U_operator, + 1, + 1, + this->ParaV->desc); +} -void Evolve_LCAO_Matrix::norm_wfc( - const int nband, - const int nlocal, - const std::complex* Stmp, - std::complex* psi_k, - const int print_matrix) const +void Evolve_LCAO_Matrix::U_to_wfc(const int nband, + const int nlocal, + const std::complex* U_operator, + const std::complex* psi_k_laststep, + std::complex* psi_k, + const int print_matrix) const +{ + + ScalapackConnector::gemm('N', + 'N', + nlocal, + nband, + nlocal, + 1.0, + U_operator, + 1, + 1, + this->ParaV->desc, + psi_k_laststep, + 1, + 1, + this->ParaV->desc_wfc, + 0.0, + psi_k, + 1, + 1, + this->ParaV->desc_wfc); + + if (print_matrix) + { + GlobalV::ofs_running << endl; + GlobalV::ofs_running << " psi_k:" << endl; + for (int i = 0; i < this->ParaV->ncol_bands; i++) + { + for (int j = 0; j < this->ParaV->ncol; j++) + { + double aa, bb; + aa = psi_k[i * this->ParaV->ncol + j].real(); + bb = psi_k[i * this->ParaV->ncol + j].imag(); + if (abs(aa) < 1e-8) + aa = 0.0; + if (abs(bb) < 1e-8) + bb = 0.0; + GlobalV::ofs_running << aa << "+" << bb << "i "; + } + GlobalV::ofs_running << endl; + } + GlobalV::ofs_running << endl; + GlobalV::ofs_running << " psi_k_laststep:" << endl; + for (int i = 0; i < this->ParaV->ncol_bands; i++) + { + for (int j = 0; j < this->ParaV->ncol; j++) + { + double aa, bb; + aa = psi_k_laststep[i * this->ParaV->ncol + j].real(); + bb = psi_k_laststep[i * this->ParaV->ncol + j].imag(); + if (abs(aa) < 1e-8) + aa = 0.0; + if (abs(bb) < 1e-8) + bb = 0.0; + GlobalV::ofs_running << aa << "+" << bb << "i "; + } + GlobalV::ofs_running << endl; + } + GlobalV::ofs_running << endl; + } +} + +void Evolve_LCAO_Matrix::norm_wfc(const int nband, + const int nlocal, + const std::complex* Stmp, + std::complex* psi_k, + const int print_matrix) const { complex* tmp1 = new complex[this->ParaV->nloc_wfc]; ModuleBase::GlobalFunc::ZEROS(tmp1, this->ParaV->nloc_wfc); - + complex* Cij = new complex[this->ParaV->nloc]; ModuleBase::GlobalFunc::ZEROS(Cij, this->ParaV->nloc); - - ScalapackConnector::gemm( - 'N', - 'N', - nlocal, - nband, - nlocal, - 1.0, - Stmp, - 1, - 1, - this->ParaV->desc, - psi_k, - 1, - 1, - this->ParaV->desc_wfc, - 0.0, - tmp1, - 1, - 1, - this->ParaV->desc_wfc - ); - - ScalapackConnector::gemm( - 'C', - 'N', - nband, - nband, - nlocal, - 1.0, - psi_k, - 1, - 1, - this->ParaV->desc_wfc, - tmp1, - 1, - 1, - this->ParaV->desc_wfc, - 0.0, - Cij, - 1, - 1, - this->ParaV->desc - ); + ScalapackConnector::gemm('N', + 'N', + nlocal, + nband, + nlocal, + 1.0, + Stmp, + 1, + 1, + this->ParaV->desc, + psi_k, + 1, + 1, + this->ParaV->desc_wfc, + 0.0, + tmp1, + 1, + 1, + this->ParaV->desc_wfc); + + ScalapackConnector::gemm('C', + 'N', + nband, + nband, + nlocal, + 1.0, + psi_k, + 1, + 1, + this->ParaV->desc_wfc, + tmp1, + 1, + 1, + this->ParaV->desc_wfc, + 0.0, + Cij, + 1, + 1, + this->ParaV->desc_Eij); + + if (print_matrix) + { + GlobalV::ofs_running << "original Cij :" << endl; + for (int i = 0; i < this->ParaV->ncol; i++) + { + for (int j = 0; j < this->ParaV->nrow; j++) + { + double aa, bb; + aa = Cij[i * this->ParaV->ncol + j].real(); + bb = Cij[i * this->ParaV->ncol + j].imag(); + if (abs(aa) < 1e-8) + aa = 0.0; + if (abs(bb) < 1e-8) + bb = 0.0; + GlobalV::ofs_running << aa << "+" << bb << "i "; + } + GlobalV::ofs_running << endl; + } + GlobalV::ofs_running << endl; + } int info; int myid; @@ -801,32 +1260,27 @@ void Evolve_LCAO_Matrix::norm_wfc( } // loop ipcol } // loop iprow - // std::cout << "nlocal" << nlocal << std::endl; - // std::cout << "GlobalV::NLOCAL" << GlobalV::NLOCAL << std::endl; BlasConnector::copy(this->ParaV->nloc_wfc, psi_k, 1, tmp1, 1); - ScalapackConnector::gemm( - 'N', - 'N', - nlocal, - nband, - nband, - 1.0, - tmp1, - 1, - 1, - this->ParaV->desc_wfc, - Cij, - 1, - 1, - this->ParaV->desc, - 0.0, - psi_k, - 1, - 1, - this->ParaV->desc_wfc - ); - + ScalapackConnector::gemm('N', + 'N', + nlocal, + nband, + nband, + 1.0, + tmp1, + 1, + 1, + this->ParaV->desc_wfc, + Cij, + 1, + 1, + this->ParaV->desc_Eij, + 0.0, + psi_k, + 1, + 1, + this->ParaV->desc_wfc); if (print_matrix) { @@ -835,7 +1289,8 @@ void Evolve_LCAO_Matrix::norm_wfc( { for (int j = 0; j < this->ParaV->nrow; j++) { - GlobalV::ofs_running << Cij[i * this->ParaV->ncol + j].real() << "+" << Cij[i * this->ParaV->ncol_bands + j].imag() << "i "; + GlobalV::ofs_running << Cij[i * this->ParaV->ncol + j].real() << "+" + << Cij[i * this->ParaV->ncol + j].imag() << "i "; } GlobalV::ofs_running << endl; } @@ -858,7 +1313,7 @@ void Evolve_LCAO_Matrix::norm_wfc( GlobalV::ofs_running << endl; } GlobalV::ofs_running << endl; - GlobalV::ofs_running << " psi_k nlocal*nlocal:" << endl; + GlobalV::ofs_running << " psi_k before normalization:" << endl; for (int i = 0; i < this->ParaV->ncol; i++) { for (int j = 0; j < this->ParaV->ncol; j++) @@ -878,19 +1333,15 @@ void Evolve_LCAO_Matrix::norm_wfc( GlobalV::ofs_running << endl; } - delete[] tmp1; delete[] Cij; - } - -void Evolve_LCAO_Matrix::compute_ekb( - const int nband, - const int nlocal, - const std::complex* Htmp, - const std::complex* psi_k, - double* ekb) const +void Evolve_LCAO_Matrix::compute_ekb(const int nband, + const int nlocal, + const std::complex* Htmp, + const std::complex* psi_k, + double* ekb) const { complex* tmp1 = new complex[this->ParaV->nloc_wfc]; @@ -899,50 +1350,45 @@ void Evolve_LCAO_Matrix::compute_ekb( complex* Eij = new complex[this->ParaV->nloc]; ModuleBase::GlobalFunc::ZEROS(Eij, this->ParaV->nloc); - ScalapackConnector::gemm( - 'N', - 'N', - nlocal, - nband, - nlocal, - 1.0, - Htmp, - 1, - 1, - this->ParaV->desc, - psi_k, - 1, - 1, - this->ParaV->desc_wfc, - 0.0, - tmp1, - 1, - 1, - this->ParaV->desc_wfc - ); - - ScalapackConnector::gemm( - 'C', - 'N', - nband, - nband, - nlocal, - 1.0, - psi_k, - 1, - 1, - this->ParaV->desc_wfc, - tmp1, - 1, - 1, - this->ParaV->desc_wfc, - 0.0, - Eij, - 1, - 1, - this->ParaV->desc - ); - + ScalapackConnector::gemm('N', + 'N', + nlocal, + nband, + nlocal, + 1.0, + Htmp, + 1, + 1, + this->ParaV->desc, + psi_k, + 1, + 1, + this->ParaV->desc_wfc, + 0.0, + tmp1, + 1, + 1, + this->ParaV->desc_wfc); + + ScalapackConnector::gemm('C', + 'N', + nband, + nband, + nlocal, + 1.0, + psi_k, + 1, + 1, + this->ParaV->desc_wfc, + tmp1, + 1, + 1, + this->ParaV->desc_wfc, + 0.0, + Eij, + 1, + 1, + this->ParaV->desc_Eij); if (ELEC_evolve::td_print_eij > 0.0) { @@ -950,9 +1396,9 @@ void Evolve_LCAO_Matrix::compute_ekb( << "------------------------------------------------------------------------------------------------" << endl; GlobalV::ofs_running << " Eij:" << endl; - for (int i = 0; i < this->ParaV->ncol; i++) + for (int i = 0; i < this->ParaV->nrow_bands; i++) { - for (int j = 0; j < this->ParaV->nrow; j++) + for (int j = 0; j < this->ParaV->ncol_bands; j++) { double aa, bb; aa = Eij[i * this->ParaV->ncol + j].real(); @@ -973,12 +1419,11 @@ void Evolve_LCAO_Matrix::compute_ekb( << endl; } - int info; int myid; int naroc[2]; MPI_Comm_rank(this->ParaV->comm_2D, &myid); - + double* Eii = new double[nband]; ModuleBase::GlobalFunc::ZEROS(Eii, nband); for (int iprow = 0; iprow < this->ParaV->dim0; ++iprow) @@ -1012,11 +1457,10 @@ void Evolve_LCAO_Matrix::compute_ekb( } // loop ipcol } // loop iprow info = MPI_Allreduce(Eii, ekb, nband, MPI_DOUBLE, MPI_SUM, this->ParaV->comm_2D); - + delete[] tmp1; delete[] Eij; delete[] Eii; } - #endif diff --git a/source/module_hamilt_lcao/module_tddft/LCAO_evolve.h b/source/module_hamilt_lcao/module_tddft/LCAO_evolve.h index 1f1f4e1fc0..bcf470a771 100644 --- a/source/module_hamilt_lcao/module_tddft/LCAO_evolve.h +++ b/source/module_hamilt_lcao/module_tddft/LCAO_evolve.h @@ -22,8 +22,10 @@ class Evolve_LCAO_Matrix hamilt::Hamilt* p_hamilt, psi::Psi>* psi_k, psi::Psi>* psi_k_laststep, - std::complex* Hk_laststep, - double* ekb) const; + std::complex* H_laststep, + double* ekb, + int htype, + int propagator) const; private: // LCAO_Matrix* LM; @@ -39,16 +41,42 @@ class Evolve_LCAO_Matrix const int nband, const int nlocal, const std::complex* psi_k_laststep, - const std::complex* Hk_laststep, + const std::complex* H_laststep, hamilt::Hamilt* p_hamilt, std::complex* psi_k, - double* ekb) const; + double* ekb, + int htype, + int propagator) const; + + void half_Hmatrix( + const int nband, + const int nlocal, + std::complex* Htmp, + const std::complex* H_laststep, + const int print_matrix) const; + + void compute_U_operator_CN2( + const int nband, + const int nlocal, + const std::complex* Stmp, + const std::complex* Htmp, + std::complex* U_operator, + const int print_matrix) const; - void compute_U_operator( + void compute_U_operator_taylor( + const int nband, + const int nlocal, + const std::complex* Stmp, + const std::complex* Htmp, + std::complex* U_operator, + const int print_matrix) const; + + void compute_U_operator_etrs( const int nband, const int nlocal, const std::complex* Stmp, const std::complex* Htmp, + const std::complex* H_laststep, std::complex* U_operator, const int print_matrix) const; @@ -57,7 +85,8 @@ class Evolve_LCAO_Matrix const int nlocal, const std::complex* U_operator, const std::complex* psi_k_laststep, - std::complex* psi_k) const; + std::complex* psi_k, + const int print_matrix) const; void norm_wfc( const int nband, diff --git a/source/module_io/input.cpp b/source/module_io/input.cpp index 1f1eaac91c..200e7a3ccf 100644 --- a/source/module_io/input.cpp +++ b/source/module_io/input.cpp @@ -408,6 +408,9 @@ void Input::Default(void) td_vext = false; td_vext_dire = "1"; + td_htype = 0; + propagator = 0; + out_dipole = false; out_efield = false; @@ -1545,6 +1548,14 @@ bool Input::Read(const std::string &fn) { read_value(ifs, td_edm); } + else if (strcmp("td_htype", word) == 0) + { + read_value(ifs, td_htype); + } + else if (strcmp("propagator", word) == 0) + { + read_value(ifs, propagator); + } else if (strcmp("td_stype", word) == 0) { read_value(ifs, td_stype); @@ -2957,6 +2968,8 @@ void Input::Bcast() Parallel_Common::bcast_double(td_force_dt); Parallel_Common::bcast_bool(td_vext); Parallel_Common::bcast_string(td_vext_dire); + Parallel_Common::bcast_int(td_htype); + Parallel_Common::bcast_int(propagator); Parallel_Common::bcast_int(td_stype); Parallel_Common::bcast_string(td_ttype); Parallel_Common::bcast_int(td_tstart); diff --git a/source/module_io/input.h b/source/module_io/input.h index f9002da4f6..fc2960c1a9 100644 --- a/source/module_io/input.h +++ b/source/module_io/input.h @@ -380,6 +380,9 @@ class Input double td_print_eij; // threshold to output Eij elements int td_edm; //0: new edm method 1: old edm method + int td_htype; // type of Hamiltonian to calculate propagator + int propagator; // method of propagator + int td_stype ; //type of space domain 0 : length gauge 1: velocity gauge std::string td_ttype ; //type of time domain diff --git a/source/module_io/test/input_test.cpp b/source/module_io/test/input_test.cpp index b594b63e8d..6d61502863 100644 --- a/source/module_io/test/input_test.cpp +++ b/source/module_io/test/input_test.cpp @@ -252,6 +252,8 @@ TEST_F(InputTest, Default) EXPECT_DOUBLE_EQ(INPUT.td_force_dt,0.02); EXPECT_FALSE(INPUT.td_vext); EXPECT_EQ(INPUT.td_vext_dire,"1"); + EXPECT_EQ(INPUT.td_htype,0); + EXPECT_EQ(INPUT.propagator,0); EXPECT_EQ(INPUT.td_stype,0); EXPECT_EQ(INPUT.td_ttype,"0"); EXPECT_EQ(INPUT.td_tstart,1); @@ -584,6 +586,8 @@ TEST_F(InputTest, Read) EXPECT_DOUBLE_EQ(INPUT.td_force_dt,0.02); EXPECT_EQ(INPUT.td_vext,0); // EXPECT_EQ(INPUT.td_vext_dire,"1"); + EXPECT_EQ(INPUT.td_htype,0); + EXPECT_EQ(INPUT.propagator,0); EXPECT_EQ(INPUT.td_stype,0); // EXPECT_EQ(INPUT.td_ttype,"0"); EXPECT_EQ(INPUT.td_tstart,1); diff --git a/source/module_io/test/input_test_para.cpp b/source/module_io/test/input_test_para.cpp index b4b3afbcba..af1f917b0d 100644 --- a/source/module_io/test/input_test_para.cpp +++ b/source/module_io/test/input_test_para.cpp @@ -258,6 +258,8 @@ TEST_F(InputParaTest,Bcast) EXPECT_DOUBLE_EQ(INPUT.td_force_dt,0.02); EXPECT_FALSE(INPUT.td_vext); EXPECT_EQ(INPUT.td_vext_dire,"1"); + EXPECT_EQ(INPUT.td_htype,0); + EXPECT_EQ(INPUT.propagator,0); EXPECT_EQ(INPUT.td_stype,0); EXPECT_EQ(INPUT.td_ttype,"0"); EXPECT_EQ(INPUT.td_tstart,1); diff --git a/source/module_io/test/support/INPUT b/source/module_io/test/support/INPUT index 91b4fc9de4..92e7e6e1cc 100644 --- a/source/module_io/test/support/INPUT +++ b/source/module_io/test/support/INPUT @@ -284,6 +284,8 @@ exx_opt_orb_tolerence 0 # td_force_dt 0.02 #time of force change td_vext 0 #add extern potential or not td_vext_dire 1 #extern potential direction +td_htype 0 # type of Hamiltonian to calculate propagator +propagator 0 # method of propagator td_stype 0 #space domain type td_ttype 0 #time domain type td_tstart 1 #the start step of electric field diff --git a/source/module_io/test/support/witestfile b/source/module_io/test/support/witestfile index 67ce58da05..ad0a1e5a37 100644 --- a/source/module_io/test/support/witestfile +++ b/source/module_io/test/support/witestfile @@ -278,6 +278,8 @@ exx_opt_orb_tolerence 0 # td_force_dt 0.02 #time of force change td_vext 0 #add extern potential or not td_vext_dire 1 #extern potential direction +td_htype 0 # type of Hamiltonian to calculate propagator +propagator 0 # method of propagator td_stype 0 #space domain type td_ttype 0 #time domain type td_tstart 1 #the start step of electric field diff --git a/tests/integrate/601_NO_TDDFT_H2_etrs/INPUT b/tests/integrate/601_NO_TDDFT_H2_etrs/INPUT new file mode 100755 index 0000000000..d0ea674f1b --- /dev/null +++ b/tests/integrate/601_NO_TDDFT_H2_etrs/INPUT @@ -0,0 +1,35 @@ +INPUT_PARAMETERS +#Parameters (General) +suffix autotest +pseudo_dir ../../PP_ORB +orbital_dir ../../PP_ORB + +nbands 5 +calculation md +esolver_type tddft + +#Parameter (Accuracy) +ecutwfc 100 +scf_nmax 50 + +ks_solver scalapack_gvx +basis_type lcao +gamma_only 0 +md_nstep 2 + +mixing_type pulay +mixing_beta 0.7 +scf_thr 1.0e-6 + +cal_stress 1 +stress_thr 1e-6 +cal_force 1 +force_thr_ev 1.0e-3 + +propagator 2 + +md_type nve +md_dt 0.05 +init_vel 1 +ocp 1 +ocp_set 1*1 1*1 3*0 diff --git a/tests/integrate/601_NO_TDDFT_H2_etrs/KPT b/tests/integrate/601_NO_TDDFT_H2_etrs/KPT new file mode 100755 index 0000000000..c289c0158a --- /dev/null +++ b/tests/integrate/601_NO_TDDFT_H2_etrs/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/tests/integrate/601_NO_TDDFT_H2_etrs/STRU b/tests/integrate/601_NO_TDDFT_H2_etrs/STRU new file mode 100755 index 0000000000..af32239121 --- /dev/null +++ b/tests/integrate/601_NO_TDDFT_H2_etrs/STRU @@ -0,0 +1,22 @@ +ATOMIC_SPECIES +H 1.00794 H.LDA.UPF + +NUMERICAL_ORBITAL +H_lda_8.0au_100Ry_2s1p.orb + +LATTICE_CONSTANT +15 + +LATTICE_VECTORS +1 0 0 #latvec1 +0 1 0 #latvec2 +0 0 1 #latvec3 + +ATOMIC_POSITIONS +Cartesian + +H #label +0 #magnetism +2 #number of atoms +0.502095122795 0.474071446432 0.509122850526 m 1 1 1 v -1.58317459417e-05 0.000196314331857 0.000681467807833 +0.497904877205 0.525928553568 0.592277149473 m 1 1 1 v 1.58317459417e-05 -0.000196314331857 -0.000681467807833 diff --git a/tests/integrate/601_NO_TDDFT_H2_etrs/result.ref b/tests/integrate/601_NO_TDDFT_H2_etrs/result.ref new file mode 100644 index 0000000000..799a64080a --- /dev/null +++ b/tests/integrate/601_NO_TDDFT_H2_etrs/result.ref @@ -0,0 +1,5 @@ +etotref +etotperatomref +totalforceref 0.714836 +totalstressref 1.295806 +totaltimeref diff --git a/tests/integrate/601_NO_TDDFT_H2_halfH/INPUT b/tests/integrate/601_NO_TDDFT_H2_halfH/INPUT new file mode 100755 index 0000000000..901858802f --- /dev/null +++ b/tests/integrate/601_NO_TDDFT_H2_halfH/INPUT @@ -0,0 +1,36 @@ +INPUT_PARAMETERS +#Parameters (General) +suffix autotest +pseudo_dir ../../PP_ORB +orbital_dir ../../PP_ORB + +nbands 5 +calculation md +esolver_type tddft + +#Parameter (Accuracy) +ecutwfc 100 +scf_nmax 50 + +ks_solver scalapack_gvx +basis_type lcao +gamma_only 0 +md_nstep 2 + +mixing_type pulay +mixing_beta 0.7 +scf_thr 1.0e-6 + +cal_stress 1 +stress_thr 1e-6 +cal_force 1 +force_thr_ev 1.0e-3 + +td_htype 1 +propagator 0 + +md_type nve +md_dt 0.05 +init_vel 1 +ocp 1 +ocp_set 1*1 1*1 3*0 diff --git a/tests/integrate/601_NO_TDDFT_H2_halfH/KPT b/tests/integrate/601_NO_TDDFT_H2_halfH/KPT new file mode 100755 index 0000000000..c289c0158a --- /dev/null +++ b/tests/integrate/601_NO_TDDFT_H2_halfH/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/tests/integrate/601_NO_TDDFT_H2_halfH/STRU b/tests/integrate/601_NO_TDDFT_H2_halfH/STRU new file mode 100755 index 0000000000..af32239121 --- /dev/null +++ b/tests/integrate/601_NO_TDDFT_H2_halfH/STRU @@ -0,0 +1,22 @@ +ATOMIC_SPECIES +H 1.00794 H.LDA.UPF + +NUMERICAL_ORBITAL +H_lda_8.0au_100Ry_2s1p.orb + +LATTICE_CONSTANT +15 + +LATTICE_VECTORS +1 0 0 #latvec1 +0 1 0 #latvec2 +0 0 1 #latvec3 + +ATOMIC_POSITIONS +Cartesian + +H #label +0 #magnetism +2 #number of atoms +0.502095122795 0.474071446432 0.509122850526 m 1 1 1 v -1.58317459417e-05 0.000196314331857 0.000681467807833 +0.497904877205 0.525928553568 0.592277149473 m 1 1 1 v 1.58317459417e-05 -0.000196314331857 -0.000681467807833 diff --git a/tests/integrate/601_NO_TDDFT_H2_halfH/result.ref b/tests/integrate/601_NO_TDDFT_H2_halfH/result.ref new file mode 100644 index 0000000000..d1095ae3ec --- /dev/null +++ b/tests/integrate/601_NO_TDDFT_H2_halfH/result.ref @@ -0,0 +1,5 @@ +etotref -18.11864981293731 +etotperatomref -9.0593249065 +totalforceref 44.341410 +totalstressref 78.497930 +totaltimeref +3.9966 diff --git a/tests/integrate/601_NO_TDDFT_H2_taylor/INPUT b/tests/integrate/601_NO_TDDFT_H2_taylor/INPUT new file mode 100755 index 0000000000..9246c5724e --- /dev/null +++ b/tests/integrate/601_NO_TDDFT_H2_taylor/INPUT @@ -0,0 +1,36 @@ +INPUT_PARAMETERS +#Parameters (General) +suffix autotest +pseudo_dir ../../PP_ORB +orbital_dir ../../PP_ORB + +nbands 5 +calculation md +esolver_type tddft + +#Parameter (Accuracy) +ecutwfc 100 +scf_nmax 50 + +ks_solver scalapack_gvx +basis_type lcao +gamma_only 0 +md_nstep 2 + +mixing_type pulay +mixing_beta 0.7 +scf_thr 1.0e-6 + +cal_stress 1 +stress_thr 1e-6 +cal_force 1 +force_thr_ev 1.0e-3 + +td_htype 1 +propagator 1 + +md_type nve +md_dt 0.05 +init_vel 1 +ocp 1 +ocp_set 1*1 1*1 3*0 diff --git a/tests/integrate/601_NO_TDDFT_H2_taylor/KPT b/tests/integrate/601_NO_TDDFT_H2_taylor/KPT new file mode 100755 index 0000000000..c289c0158a --- /dev/null +++ b/tests/integrate/601_NO_TDDFT_H2_taylor/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/tests/integrate/601_NO_TDDFT_H2_taylor/STRU b/tests/integrate/601_NO_TDDFT_H2_taylor/STRU new file mode 100755 index 0000000000..af32239121 --- /dev/null +++ b/tests/integrate/601_NO_TDDFT_H2_taylor/STRU @@ -0,0 +1,22 @@ +ATOMIC_SPECIES +H 1.00794 H.LDA.UPF + +NUMERICAL_ORBITAL +H_lda_8.0au_100Ry_2s1p.orb + +LATTICE_CONSTANT +15 + +LATTICE_VECTORS +1 0 0 #latvec1 +0 1 0 #latvec2 +0 0 1 #latvec3 + +ATOMIC_POSITIONS +Cartesian + +H #label +0 #magnetism +2 #number of atoms +0.502095122795 0.474071446432 0.509122850526 m 1 1 1 v -1.58317459417e-05 0.000196314331857 0.000681467807833 +0.497904877205 0.525928553568 0.592277149473 m 1 1 1 v 1.58317459417e-05 -0.000196314331857 -0.000681467807833 diff --git a/tests/integrate/601_NO_TDDFT_H2_taylor/result.ref b/tests/integrate/601_NO_TDDFT_H2_taylor/result.ref new file mode 100644 index 0000000000..ab82b76038 --- /dev/null +++ b/tests/integrate/601_NO_TDDFT_H2_taylor/result.ref @@ -0,0 +1,5 @@ +etotref -18.06759546636300 +etotperatomref -9.0337977332 +totalforceref 46.058216 +totalstressref 81.553158 +totaltimeref +3.8555 diff --git a/tests/integrate/CASES b/tests/integrate/CASES index ccf4623199..78face74b9 100644 --- a/tests/integrate/CASES +++ b/tests/integrate/CASES @@ -184,7 +184,7 @@ 315_NO_sol_H2O 316_NO_scan_Si2 320_NO_GO_MD_MSST -320_NO_GO_MD_NHC +320_NO_GO_MD_NVT 345_NO_GO_BS 360_NO_15_GO_PU_AF #381_NO_GO_HSE @@ -192,6 +192,9 @@ #401_NP_KP_sp #401_NP_KP_spd 601_NO_TDDFT_H2 +601_NO_TDDFT_H2_halfH +601_NO_TDDFT_H2_taylor +601_NO_TDDFT_H2_etrs 601_NO_TDDFT_H2_kpoint 601_NO_TDDFT_H2_oldedm 601_NO_TDDFT_CO