diff --git a/source/source_esolver/esolver_ks_lcao_tddft.cpp b/source/source_esolver/esolver_ks_lcao_tddft.cpp index 66d7717f14f..fd47b15fc53 100644 --- a/source/source_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/source_esolver/esolver_ks_lcao_tddft.cpp @@ -101,7 +101,7 @@ void ESolver_KS_LCAO_TDDFT::runner(UnitCell& ucell, const int istep) // 1) before_scf (electronic iteration loops) //---------------------------------------------------------------- this->before_scf(ucell, istep); // From ESolver_KS_LCAO - + td_p->initialize_phase_hybrid(ucell, dynamic_cast, TR>*>(this->p_hamilt)->getHR()); // Initialize the moving spatial gauge if (use_td_moving_gauge && this->td_mg_ == nullptr) { @@ -113,7 +113,7 @@ void ESolver_KS_LCAO_TDDFT::runner(UnitCell& ucell, const int istep) if (PARAM.inp.td_stype == 2) { - this->dmat.dm->cal_DMR_td(ucell, TD_info::cart_At); + this->dmat.dm->cal_DMR_td(td_p->get_phase_hybrid(),TD_info::cart_At); } else { @@ -594,7 +594,7 @@ void ESolver_KS_LCAO_TDDFT::weight_dm_rho(const UnitCell& ucell) elecstate::cal_dm_psi(this->dmat.dm->get_paraV_pointer(), this->pelec->wg, this->psi[0], *this->dmat.dm); if (PARAM.inp.td_stype == 2) { - this->dmat.dm->cal_DMR_td(ucell, TD_info::cart_At); + this->dmat.dm->cal_DMR_td(td_p->get_phase_hybrid(), TD_info::cart_At); } else { diff --git a/source/source_estate/module_dm/cal_dm_psi.cpp b/source/source_estate/module_dm/cal_dm_psi.cpp index 0445a10a16b..5ad2144c5bc 100644 --- a/source/source_estate/module_dm/cal_dm_psi.cpp +++ b/source/source_estate/module_dm/cal_dm_psi.cpp @@ -70,11 +70,11 @@ void cal_dm_psi(const Parallel_Orbitals* ParaV, return; } - +template void cal_dm_psi(const Parallel_Orbitals* ParaV, const ModuleBase::matrix& wg, const psi::Psi>& wfc, - elecstate::DensityMatrix, double>& DM) + elecstate::DensityMatrix, TR>& DM) { ModuleBase::TITLE("elecstate", "cal_dm_psi"); ModuleBase::timer::start("elecstate", "cal_dm_psi"); @@ -268,5 +268,14 @@ void psiMulPsi(const psi::Psi>& psi1, dm_out, nlocal); } - +template +void cal_dm_psi(const Parallel_Orbitals* ParaV, + const ModuleBase::matrix& wg, + const psi::Psi>& wfc, + elecstate::DensityMatrix, std::complex>& DM); +template +void cal_dm_psi(const Parallel_Orbitals* ParaV, + const ModuleBase::matrix& wg, + const psi::Psi>& wfc, + elecstate::DensityMatrix, double>& DM); } // namespace elecstate diff --git a/source/source_estate/module_dm/cal_dm_psi.h b/source/source_estate/module_dm/cal_dm_psi.h index 09ab3d3974a..a32bc8c22c3 100644 --- a/source/source_estate/module_dm/cal_dm_psi.h +++ b/source/source_estate/module_dm/cal_dm_psi.h @@ -11,7 +11,8 @@ namespace elecstate void cal_dm_psi(const Parallel_Orbitals* ParaV, const ModuleBase::matrix& wg, const psi::Psi& wfc, elecstate::DensityMatrix& DM); // for Multi-k case where DMK is std::complex - void cal_dm_psi(const Parallel_Orbitals* ParaV, const ModuleBase::matrix& wg, const psi::Psi>& wfc, elecstate::DensityMatrix, double>& DM); + template + void cal_dm_psi(const Parallel_Orbitals* ParaV, const ModuleBase::matrix& wg, const psi::Psi>& wfc, elecstate::DensityMatrix, TR>& DM); // for Gamma-Only case with MPI void psiMulPsiMpi(const psi::Psi& psi1, diff --git a/source/source_estate/module_dm/density_matrix.cpp b/source/source_estate/module_dm/density_matrix.cpp index b0734e6d065..7414e00fd3b 100644 --- a/source/source_estate/module_dm/density_matrix.cpp +++ b/source/source_estate/module_dm/density_matrix.cpp @@ -214,7 +214,7 @@ template void DensityMatrix_Tools::cal_DMR_td( const DensityMatrix &dm, std::vector*> &dmR_out, - const UnitCell& ucell, + const std::map, std::complex>& phase_hybrid, const ModuleBase::Vector3 At, const int ik_in) { @@ -262,19 +262,21 @@ void DensityMatrix_Tools::cal_DMR_td( } #endif target_DMR_mat_vec[iR] = target_mat->get_pointer(); - //cal tddft phase for hybrid gauge - const ModuleBase::Vector3 dtau = ucell.cal_dtau(iat1, iat2, R_index); - const double arg_td = At * dtau * ucell.lat0; for(int ik = 0; ik < dm._nk; ++ik) { if(ik_in >= 0 && ik_in != ik) { continue; } // cal k_phase // if TK==std::complex, kphase is e^{ikR} const ModuleBase::Vector3 dR(R_index[0], R_index[1], R_index[2]); - const double arg = (dm._kvec_d[ik] * dR) * ModuleBase::TWO_PI + arg_td; + const double arg = (dm._kvec_d[ik] * dR) * ModuleBase::TWO_PI; double sinp, cosp; ModuleBase::libm::sincos(arg, &sinp, &cosp); kphase_vec[ik][iR] = TK(cosp, sinp); + if(PARAM.inp.td_stype==2) + { + //phase for hybrid gauge tddft + kphase_vec[ik][iR] *= phase_hybrid.at(R_index); + } } } @@ -353,20 +355,20 @@ void DensityMatrix_Tools::cal_DMR_td( ModuleBase::timer::end("DensityMatrix", "cal_DMR_td"); } template <> -void DensityMatrix::cal_DMR_td(const UnitCell& ucell, const ModuleBase::Vector3 At, const int ik_in) +void DensityMatrix::cal_DMR_td(const std::map, std::complex>& phase_hybrid, const ModuleBase::Vector3 At, const int ik_in) { return; } template <> -void DensityMatrix, double>::cal_DMR_td(const UnitCell& ucell, const ModuleBase::Vector3 At, const int ik_in) +void DensityMatrix, double>::cal_DMR_td(const std::map, std::complex>& phase_hybrid, const ModuleBase::Vector3 At, const int ik_in) { - DensityMatrix_Tools::cal_DMR_td(*this, this->_DMR, ucell, At, ik_in); + DensityMatrix_Tools::cal_DMR_td(*this, this->_DMR, phase_hybrid, At, ik_in); } template <> -void DensityMatrix, std::complex>::cal_DMR_td(const UnitCell& ucell, const ModuleBase::Vector3 At, const int ik_in) +void DensityMatrix, std::complex>::cal_DMR_td(const std::map, std::complex>& phase_hybrid, const ModuleBase::Vector3 At, const int ik_in) { - DensityMatrix_Tools::cal_DMR_td(*this, this->_DMR, ucell, At, ik_in); + DensityMatrix_Tools::cal_DMR_td(*this, this->_DMR, phase_hybrid, At, ik_in); } diff --git a/source/source_estate/module_dm/density_matrix.h b/source/source_estate/module_dm/density_matrix.h index a8b0e1c4ecb..90c51e9c061 100644 --- a/source/source_estate/module_dm/density_matrix.h +++ b/source/source_estate/module_dm/density_matrix.h @@ -47,7 +47,7 @@ namespace DensityMatrix_Tools extern void cal_DMR_td( const DensityMatrix &dm, std::vector*> &dmR_out, - const UnitCell& ucell, + const std::map, std::complex>& phase_hybrid, const ModuleBase::Vector3 At, const int ik_in); @@ -225,7 +225,7 @@ class DensityMatrix * if ik_in < 0, calculate all k-points * if ik_in >= 0, calculate only one k-point without summing over k-points */ - void cal_DMR_td(const UnitCell& ucell, const ModuleBase::Vector3 At, const int ik_in = -1); + void cal_DMR_td(const std::map, std::complex>& phase_hybrid, const ModuleBase::Vector3 At, const int ik_in = -1); /** * @brief calculate complex density matrix DMR with both real and imaginary part for noncollinear-spin calculation @@ -327,7 +327,7 @@ class DensityMatrix TR* dmr_tmp_ = nullptr; friend void DensityMatrix_Tools::cal_DMR(const DensityMatrix &dm, std::vector*> &dmR_out, const int ik_in); - friend void DensityMatrix_Tools::cal_DMR_td(const DensityMatrix &dm, std::vector*> &dmR_out, const UnitCell& ucell, const ModuleBase::Vector3 At, const int ik_in); + friend void DensityMatrix_Tools::cal_DMR_td(const DensityMatrix &dm, std::vector*> &dmR_out, const std::map, std::complex>& phase_hybrid, const ModuleBase::Vector3 At, const int ik_in); friend void DensityMatrix_Tools::cal_DMR_full(const DensityMatrix &dm, hamilt::HContainer>* dmR_out, const int ik_in); }; diff --git a/source/source_estate/module_dm/init_dm.cpp b/source/source_estate/module_dm/init_dm.cpp index ac7a9b19dd2..20b7fa1217d 100644 --- a/source/source_estate/module_dm/init_dm.cpp +++ b/source/source_estate/module_dm/init_dm.cpp @@ -25,7 +25,7 @@ void elecstate::init_dm(UnitCell& ucell, elecstate::cal_dm_psi(dmat.dm->get_paraV_pointer(), pelec->wg, *psi, *dmat.dm); if (PARAM.inp.esolver_type!="tddft" && PARAM.inp.td_stype == 2) { - dmat.dm->cal_DMR_td(ucell, TD_info::cart_At); + dmat.dm->cal_DMR_td(TD_info::td_vel_op->get_phase_hybrid(), TD_info::cart_At); } else { diff --git a/source/source_io/module_ctrl/ctrl_output_td.cpp b/source/source_io/module_ctrl/ctrl_output_td.cpp index d1d2a9aeb5c..1a5c37dc89c 100644 --- a/source/source_io/module_ctrl/ctrl_output_td.cpp +++ b/source/source_io/module_ctrl/ctrl_output_td.cpp @@ -42,16 +42,16 @@ void ctrl_output_td(const UnitCell& ucell, { if (TD_info::out_current_k) { - ModuleIO::write_current_eachk(ucell, istep, psi, pelec, kv, intor, pv, orb, velocity_mat, RA); + ModuleIO::write_current_eachk(ucell, istep, psi, pelec, kv, intor, pv, orb, velocity_mat, td_p, RA); } else { - ModuleIO::write_current(ucell, istep, psi, pelec, kv, intor, pv, orb, velocity_mat, RA); + ModuleIO::write_current(ucell, istep, psi, pelec, kv, intor, pv, orb, velocity_mat, td_p, RA); } } else if(TD_info::out_current==2) { - ModuleIO::write_current(ucell, grid, istep, psi, pelec, kv, pv, orb, td_p->r_calculator, p_hamilt->getSR(), p_hamilt->getHR(), exx_nao); + ModuleIO::write_current(ucell, grid, istep, psi, pelec, kv, pv, orb, td_p, p_hamilt->getSR(), p_hamilt->getHR(), exx_nao); } // (3) Output file for restart if (PARAM.inp.out_freq_td > 0) // default value of out_freq_td is 0 diff --git a/source/source_io/module_current/td_current_io.cpp b/source/source_io/module_current/td_current_io.cpp index b5a59cbe270..b8e160bb3b4 100644 --- a/source/source_io/module_current/td_current_io.cpp +++ b/source/source_io/module_current/td_current_io.cpp @@ -10,26 +10,9 @@ #include "source_estate/module_dm/cal_dm_psi.h" #include "source_estate/module_pot/H_TDDFT_pw.h" #include "source_lcao/LCAO_domain.h" -#include "source_lcao/module_rt/td_info.h" #include "source_io/module_parameter/parameter.h" #ifdef __LCAO -void ModuleIO::cal_tmp_DM(const UnitCell& ucell, - elecstate::DensityMatrix, double>& DM_real, - elecstate::DensityMatrix, double>& DM_imag, - int nspin_dm) -{ - ModuleBase::TITLE("ModuleIO", "cal_tmp_DM"); - ModuleBase::timer::start("ModuleIO", "cal_tmp_DM"); - for (int is = 1; is <= nspin_dm; ++is) - { - for (int ik = 0; ik < DM_real.get_DMK_nks() / nspin_dm; ++ik) - { - cal_tmp_DM_k(ucell, DM_real, DM_imag, ik, nspin_dm, is, false); - } - } - ModuleBase::timer::end("ModuleIO", "cal_tmp_DM"); -} template void ModuleIO::write_current(const UnitCell& ucell, const int istep, @@ -40,6 +23,7 @@ void ModuleIO::write_current(const UnitCell& ucell, const Parallel_Orbitals* pv, const LCAO_Orbitals& orb, const Velocity_op* cal_current, + TD_info* td_p, Record_adj& ra) { @@ -70,15 +54,21 @@ void ModuleIO::write_current(const UnitCell& ucell, // be refactored in the future. const int nspin0 = PARAM.inp.nspin; const int nspin_dm = std::map({ {1,1},{2,2},{4,1} })[nspin0]; - elecstate::DensityMatrix, double> DM_real(pv, nspin_dm, kv.kvec_d, kv.get_nks() / nspin_dm); - elecstate::DensityMatrix, double> DM_imag(pv, nspin_dm, kv.kvec_d, kv.get_nks() / nspin_dm); + elecstate::DensityMatrix, std::complex> tmp_dm(pv, nspin_dm, kv.kvec_d, kv.get_nks() / nspin_dm); // calculate DMK - elecstate::cal_dm_psi(DM_real.get_paraV_pointer(), pelec->wg, psi[0], DM_real); + elecstate::cal_dm_psi(pv, pelec->wg, psi[0], tmp_dm); // init DMR - DM_real.init_DMR(ra, &ucell); - DM_imag.init_DMR(ra, &ucell); - cal_tmp_DM(ucell, DM_real, DM_imag, nspin_dm); + tmp_dm.init_DMR(ra, &ucell); + + if(PARAM.inp.td_stype!=2) + { + tmp_dm.cal_DMR(); + } + else + { + tmp_dm.cal_DMR_td(td_p->get_phase_hybrid(),TD_info::cart_At); + } //DM_real.sum_DMR_spin(); //DM_imag.sum_DMR_spin(); @@ -120,10 +110,8 @@ void ModuleIO::write_current(const UnitCell& ucell, double Rz = ra.info[iat][cb][2]; //std::cout<< "iat1: " << iat1 << " iat2: " << iat2 << " Rx: " << Rx << " Ry: " << Ry << " Rz:" << Rz << std::endl; // get BaseMatrix - hamilt::BaseMatrix* tmp_matrix_real - = DM_real.get_DMR_pointer(1)->find_matrix(iat1, iat2, Rx, Ry, Rz); - hamilt::BaseMatrix* tmp_matrix_imag - = DM_imag.get_DMR_pointer(1)->find_matrix(iat1, iat2, Rx, Ry, Rz); + hamilt::BaseMatrix>* tmp_matrix + = tmp_dm.get_DMR_pointer(1)->find_matrix(iat1, iat2, Rx, Ry, Rz); // refactor hamilt::BaseMatrix>* tmp_m_rvx = current_term[0]->find_matrix(iat1, iat2, Rx, Ry, Rz); @@ -131,7 +119,7 @@ void ModuleIO::write_current(const UnitCell& ucell, = current_term[1]->find_matrix(iat1, iat2, Rx, Ry, Rz); hamilt::BaseMatrix>* tmp_m_rvz = current_term[2]->find_matrix(iat1, iat2, Rx, Ry, Rz); - if (tmp_matrix_real == nullptr) + if (tmp_matrix == nullptr) { continue; } @@ -142,8 +130,7 @@ void ModuleIO::write_current(const UnitCell& ucell, { for (int nu = 0; nu < pv->get_col_size(iat2); ++nu) { - double dm2d1_real = tmp_matrix_real->get_value(mu, nu); - double dm2d1_imag = tmp_matrix_imag->get_value(mu, nu); + std::complex dm2d1 = tmp_matrix->get_value(mu, nu); std::complex rvx = {0, 0}; std::complex rvy = {0, 0}; @@ -158,9 +145,9 @@ void ModuleIO::write_current(const UnitCell& ucell, //std::cout<<"mu: "<< mu <<" nu: "<< nu << std::endl; // std::cout<<"dm2d1_real: "<< dm2d1_real << " dm2d1_imag: "<< dm2d1_imag << std::endl; //std::cout<<"rvz: "<< rvz.real() << " " << rvz.imag() << std::endl; - local_current[0] -= dm2d1_real * rvx.real() - dm2d1_imag * rvx.imag(); - local_current[1] -= dm2d1_real * rvy.real() - dm2d1_imag * rvy.imag(); - local_current[2] -= dm2d1_real * rvz.real() - dm2d1_imag * rvz.imag(); + local_current[0] -= dm2d1.real() * rvx.real() - dm2d1.imag() * rvx.imag(); + local_current[1] -= dm2d1.real() * rvy.real() - dm2d1.imag() * rvy.imag(); + local_current[2] -= dm2d1.real() * rvz.real() - dm2d1.imag() * rvz.imag(); } // end kk } // end jj } // end cb @@ -193,211 +180,6 @@ void ModuleIO::write_current(const UnitCell& ucell, ModuleBase::timer::end("ModuleIO", "write_current"); return; } -void ModuleIO::cal_tmp_DM_k(const UnitCell& ucell, - elecstate::DensityMatrix, double>& DM_real, - elecstate::DensityMatrix, double>& DM_imag, - const int ik, - const int nspin, - const int is, - const bool reset) -{ - ModuleBase::TITLE("ModuleIO", "cal_tmp_DM_k"); - ModuleBase::timer::start("ModuleIO", "cal_tmp_DM_k"); - int ld_hk = DM_real.get_paraV_pointer()->nrow; - int ld_hk2 = 2 * ld_hk; - // tmp for is - int ik_begin = DM_real.get_DMK_nks() / nspin * (is - 1); // jump nk for spin_down if nspin==2 - //sum spin up and down into up - hamilt::HContainer* tmp_DMR_real = DM_real.get_DMR_vector()[0]; - hamilt::HContainer* tmp_DMR_imag = DM_imag.get_DMR_vector()[0]; - if(reset) - { - tmp_DMR_real->set_zero(); - tmp_DMR_imag->set_zero(); - } -#ifdef _OPENMP -#pragma omp parallel for -#endif - for (int i = 0; i < tmp_DMR_real->size_atom_pairs(); ++i) - { - hamilt::AtomPair& tmp_ap_real = tmp_DMR_real->get_atom_pair(i); - hamilt::AtomPair& tmp_ap_imag = tmp_DMR_imag->get_atom_pair(i); - int iat1 = tmp_ap_real.get_atom_i(); - int iat2 = tmp_ap_real.get_atom_j(); - // get global indexes of whole matrix for each atom in this process - int row_ap = DM_real.get_paraV_pointer()->atom_begin_row[iat1]; - int col_ap = DM_real.get_paraV_pointer()->atom_begin_col[iat2]; - // SOC - std::vector> tmp_DMR; - if (PARAM.inp.nspin == 4) - { - tmp_DMR.resize(tmp_ap_real.get_size()); - } - for (int ir = 0; ir < tmp_ap_real.get_R_size(); ++ir) - { - const ModuleBase::Vector3 r_index = tmp_ap_real.get_R_index(ir); - hamilt::BaseMatrix* tmp_matrix_real = tmp_ap_real.find_matrix(r_index); - hamilt::BaseMatrix* tmp_matrix_imag = tmp_ap_imag.find_matrix(r_index); -#ifdef __DEBUG - if (tmp_matrix_real == nullptr) - { - std::cout << "tmp_matrix is nullptr" << std::endl; - continue; - } -#endif - // only ik - if (PARAM.inp.nspin != 4) - { - double arg_td = 0.0; - if(elecstate::H_TDDFT_pw::stype == 2) - { - //cal tddft phase for hybrid gauge - const int iat1 = tmp_ap_real.get_atom_i(); - const int iat2 = tmp_ap_real.get_atom_j(); - ModuleBase::Vector3 dtau = ucell.cal_dtau(iat1, iat2, r_index); - double& tmp_lat0 = ucell.lat0; - arg_td = TD_info::cart_At * dtau * tmp_lat0; - } - // cal k_phase - // if TK==std::complex, kphase is e^{ikR} - const ModuleBase::Vector3 dR(r_index.x, r_index.y, r_index.z); - const double arg = (DM_real.get_kvec_d()[ik] * dR) * ModuleBase::TWO_PI + arg_td; - double sinp, cosp; - ModuleBase::libm::sincos(arg, &sinp, &cosp); - std::complex kphase = std::complex(cosp, sinp); - // set DMR element - double* tmp_DMR_real_pointer = tmp_matrix_real->get_pointer(); - double* tmp_DMR_imag_pointer = tmp_matrix_imag->get_pointer(); - std::complex* tmp_DMK_pointer = DM_real.get_DMK_pointer(ik + ik_begin); - double* DMK_real_pointer = nullptr; - double* DMK_imag_pointer = nullptr; - // jump DMK to fill DMR - // DMR is row-major, DMK is column-major - tmp_DMK_pointer += col_ap * DM_real.get_paraV_pointer()->nrow + row_ap; - for (int mu = 0; mu < DM_real.get_paraV_pointer()->get_row_size(iat1); ++mu) - { - DMK_real_pointer = (double*)tmp_DMK_pointer; - DMK_imag_pointer = DMK_real_pointer + 1; - // calculate real part - BlasConnector::axpy(DM_real.get_paraV_pointer()->get_col_size(iat2), - -kphase.imag(), - DMK_imag_pointer, - ld_hk2, - tmp_DMR_real_pointer, - 1); - BlasConnector::axpy(DM_real.get_paraV_pointer()->get_col_size(iat2), - kphase.real(), - DMK_real_pointer, - ld_hk2, - tmp_DMR_real_pointer, - 1); - // calculate imag part - BlasConnector::axpy(DM_imag.get_paraV_pointer()->get_col_size(iat2), - kphase.imag(), - DMK_real_pointer, - ld_hk2, - tmp_DMR_imag_pointer, - 1); - BlasConnector::axpy(DM_imag.get_paraV_pointer()->get_col_size(iat2), - kphase.real(), - DMK_imag_pointer, - ld_hk2, - tmp_DMR_imag_pointer, - 1); - tmp_DMK_pointer += 1; - tmp_DMR_real_pointer += DM_real.get_paraV_pointer()->get_col_size(iat2); - tmp_DMR_imag_pointer += DM_imag.get_paraV_pointer()->get_col_size(iat2); - } - } - // treat DMR as pauli matrix when NSPIN=4 - if (PARAM.inp.nspin == 4) - { - tmp_DMR.assign(tmp_ap_real.get_size(), std::complex(0.0, 0.0)); - { - // cal k_phase - // if TK==std::complex, kphase is e^{ikR} - const ModuleBase::Vector3 dR(r_index.x, r_index.y, r_index.z); - double arg_td = 0.0; - if(elecstate::H_TDDFT_pw::stype == 2) - { - //new - //cal tddft phase for mixing gauge - const int iat1 = tmp_ap_real.get_atom_i(); - const int iat2 = tmp_ap_real.get_atom_j(); - ModuleBase::Vector3 dtau = ucell.cal_dtau(iat1, iat2, r_index); - double& tmp_lat0 = ucell.lat0; - arg_td = TD_info::cart_At * dtau * tmp_lat0; - } - const double arg = (DM_real.get_kvec_d()[ik] * dR) * ModuleBase::TWO_PI + arg_td; - double sinp, cosp; - ModuleBase::libm::sincos(arg, &sinp, &cosp); - std::complex kphase = std::complex(cosp, sinp); - // set DMR element - std::complex* tmp_DMR_pointer = tmp_DMR.data(); - std::complex* tmp_DMK_pointer = DM_real.get_DMK_pointer(ik + ik_begin);; - double* DMK_real_pointer = nullptr; - double* DMK_imag_pointer = nullptr; - // jump DMK to fill DMR - // DMR is row-major, DMK is column-major - tmp_DMK_pointer += col_ap * DM_real.get_paraV_pointer()->nrow + row_ap; - for (int mu = 0; mu < tmp_ap_real.get_row_size(); ++mu) - { - BlasConnector::axpy(tmp_ap_real.get_col_size(), - kphase, - tmp_DMK_pointer, - ld_hk, - tmp_DMR_pointer, - 1); - tmp_DMK_pointer += 1; - tmp_DMR_pointer += tmp_ap_real.get_col_size(); - } - } - int npol = 2; - // step_trace = 0 for NSPIN=1,2; ={0, 1, local_col, local_col+1} for NSPIN=4 - int step_trace[4]; - for (int is = 0; is < npol; is++) - { - for (int is2 = 0; is2 < npol; is2++) - { - step_trace[is * npol + is2] = tmp_ap_real.get_col_size() * is + is2; - } - } - std::complex tmp[4]; - double* target_DMR_real = tmp_matrix_real->get_pointer(); - double* target_DMR_imag = tmp_matrix_imag->get_pointer(); - std::complex* tmp_DMR_pointer = tmp_DMR.data(); - for (int irow = 0; irow < tmp_ap_real.get_row_size(); irow += 2) - { - for (int icol = 0; icol < tmp_ap_real.get_col_size(); icol += 2) - { - // catch the 4 spin component value of one orbital pair - tmp[0] = tmp_DMR_pointer[icol + step_trace[0]]; - tmp[1] = tmp_DMR_pointer[icol + step_trace[1]]; - tmp[2] = tmp_DMR_pointer[icol + step_trace[2]]; - tmp[3] = tmp_DMR_pointer[icol + step_trace[3]]; - // transfer to Pauli matrix and save the real part - // save them back to the tmp_matrix - target_DMR_real[icol + step_trace[0]] += tmp[0].real() + tmp[3].real(); - target_DMR_real[icol + step_trace[1]] += tmp[1].real() + tmp[2].real(); - target_DMR_real[icol + step_trace[2]] - += -tmp[1].imag() + tmp[2].imag(); // (i * (rho_updown - rho_downup)).real() - target_DMR_real[icol + step_trace[3]] += tmp[0].real() - tmp[3].real(); - //imag part - target_DMR_imag[icol + step_trace[0]] += tmp[0].imag() + tmp[3].imag(); - target_DMR_imag[icol + step_trace[1]] += tmp[1].imag() + tmp[2].imag(); - target_DMR_imag[icol + step_trace[2]] - += tmp[1].real() - tmp[2].real(); // (i * (rho_updown - rho_downup)).real() - target_DMR_imag[icol + step_trace[3]] += tmp[0].imag() - tmp[3].imag(); - } - tmp_DMR_pointer += tmp_ap_real.get_col_size() * 2; - target_DMR_real += tmp_ap_real.get_col_size() * 2; - target_DMR_imag += tmp_ap_real.get_col_size() * 2; - } - } - } - } - ModuleBase::timer::end("ModuleIO", "cal_tmp_DM_k"); -} template void ModuleIO::write_current_eachk(const UnitCell& ucell, const int istep, @@ -408,6 +190,7 @@ void ModuleIO::write_current_eachk(const UnitCell& ucell, const Parallel_Orbitals* pv, const LCAO_Orbitals& orb, const Velocity_op* cal_current, + TD_info* td_p, Record_adj& ra) { @@ -440,22 +223,30 @@ void ModuleIO::write_current_eachk(const UnitCell& ucell, const int nspin0 = PARAM.inp.nspin; const int nspin_dm = std::map({ {1,1},{2,2},{4,1} })[nspin0]; - elecstate::DensityMatrix, double> DM_real(pv, nspin_dm, kv.kvec_d, kv.get_nks() / nspin_dm); - elecstate::DensityMatrix, double> DM_imag(pv, nspin_dm, kv.kvec_d, kv.get_nks() / nspin_dm); + elecstate::DensityMatrix, std::complex> tmp_dm(pv, nspin_dm, kv.kvec_d, kv.get_nks() / nspin_dm); + //elecstate::DensityMatrix, double> DM_real(pv, nspin_dm, kv.kvec_d, kv.get_nks() / nspin_dm); + //elecstate::DensityMatrix, double> DM_imag(pv, nspin_dm, kv.kvec_d, kv.get_nks() / nspin_dm); // calculate DMK - elecstate::cal_dm_psi(DM_real.get_paraV_pointer(), pelec->wg, psi[0], DM_real); + elecstate::cal_dm_psi(pv, pelec->wg, psi[0], tmp_dm); // init DMR - DM_real.init_DMR(ra, &ucell); - DM_imag.init_DMR(ra, &ucell); + tmp_dm.init_DMR(ra, &ucell); - int nks = DM_real.get_DMK_nks() / nspin_dm; + int nks = tmp_dm.get_DMK_nks() / nspin_dm; double current_total[3] = {0.0, 0.0, 0.0}; for (int is = 1; is <= nspin_dm; ++is) { for (int ik = 0; ik < nks; ++ik) { - cal_tmp_DM_k(ucell, DM_real, DM_imag, ik, nspin_dm, is); + if(PARAM.inp.td_stype!=2) + { + tmp_dm.cal_DMR(ik); + } + else + { + tmp_dm.cal_DMR_td(td_p->get_phase_hybrid(),TD_info::cart_At,ik); + } + // check later double current_ik[3] = {0.0, 0.0, 0.0}; #ifdef _OPENMP @@ -497,10 +288,8 @@ void ModuleIO::write_current_eachk(const UnitCell& ucell, double Rz = ra.info[iat][cb][2]; //std::cout<< "iat1: " << iat1 << " iat2: " << iat2 << " Rx: " << Rx << " Ry: " << Ry << " Rz:" << Rz << std::endl; // get BaseMatrix - hamilt::BaseMatrix* tmp_matrix_real - = DM_real.get_DMR_pointer(is)->find_matrix(iat1, iat2, Rx, Ry, Rz); - hamilt::BaseMatrix* tmp_matrix_imag - = DM_imag.get_DMR_pointer(is)->find_matrix(iat1, iat2, Rx, Ry, Rz); + hamilt::BaseMatrix>* tmp_matrix + = tmp_dm.get_DMR_pointer(is)->find_matrix(iat1, iat2, Rx, Ry, Rz); // refactor hamilt::BaseMatrix>* tmp_m_rvx = current_term[0]->find_matrix(iat1, iat2, Rx, Ry, Rz); @@ -508,7 +297,7 @@ void ModuleIO::write_current_eachk(const UnitCell& ucell, = current_term[1]->find_matrix(iat1, iat2, Rx, Ry, Rz); hamilt::BaseMatrix>* tmp_m_rvz = current_term[2]->find_matrix(iat1, iat2, Rx, Ry, Rz); - if (tmp_matrix_real == nullptr) + if (tmp_matrix == nullptr) { continue; } @@ -519,8 +308,7 @@ void ModuleIO::write_current_eachk(const UnitCell& ucell, { for (int nu = 0; nu < pv->get_col_size(iat2); ++nu) { - double dm2d1_real = tmp_matrix_real->get_value(mu, nu); - double dm2d1_imag = tmp_matrix_imag->get_value(mu, nu); + std::complex dm2d1 = tmp_matrix->get_value(mu, nu); std::complex rvx = {0, 0}; std::complex rvy = {0, 0}; @@ -535,9 +323,9 @@ void ModuleIO::write_current_eachk(const UnitCell& ucell, // std::cout<<"mu: "<< mu <<" nu: "<< nu << std::endl; // std::cout<<"dm2d1_real: "<< dm2d1_real << " dm2d1_imag: "<< dm2d1_imag << std::endl; // std::cout<<"rvz: "<< rvz.real() << " " << rvz.imag() << std::endl; - local_current_ik[0] -= dm2d1_real * rvx.real() - dm2d1_imag * rvx.imag(); - local_current_ik[1] -= dm2d1_real * rvy.real() - dm2d1_imag * rvy.imag(); - local_current_ik[2] -= dm2d1_real * rvz.real() - dm2d1_imag * rvz.imag(); + local_current_ik[0] -= dm2d1.real() * rvx.real() - dm2d1.imag() * rvx.imag(); + local_current_ik[1] -= dm2d1.real() * rvy.real() - dm2d1.imag() * rvy.imag(); + local_current_ik[2] -= dm2d1.real() * rvz.real() - dm2d1.imag() * rvz.imag(); } // end kk } // end jj } // end cb @@ -601,6 +389,7 @@ void ModuleIO::write_current_eachk( const Parallel_Orbitals* pv, const LCAO_Orbitals& orb, const Velocity_op* cal_current, + TD_info* td_p, Record_adj& ra); template void ModuleIO::write_current_eachk>(const UnitCell& ucell, @@ -612,6 +401,7 @@ void ModuleIO::write_current_eachk>(const UnitCell& ucell, const Parallel_Orbitals* pv, const LCAO_Orbitals& orb, const Velocity_op>* cal_current, + TD_info* td_p, Record_adj& ra); template void ModuleIO::write_current(const UnitCell& ucell, @@ -623,6 +413,7 @@ void ModuleIO::write_current(const UnitCell& ucell, const Parallel_Orbitals* pv, const LCAO_Orbitals& orb, const Velocity_op* cal_current, + TD_info* td_p, Record_adj& ra); template void ModuleIO::write_current>(const UnitCell& ucell, @@ -634,6 +425,7 @@ void ModuleIO::write_current>(const UnitCell& ucell, const Parallel_Orbitals* pv, const LCAO_Orbitals& orb, const Velocity_op>* cal_current, + TD_info* td_p, Record_adj& ra); #endif //__LCAO diff --git a/source/source_io/module_current/td_current_io.h b/source/source_io/module_current/td_current_io.h index 2a64362b032..c0ac2d6f409 100644 --- a/source/source_io/module_current/td_current_io.h +++ b/source/source_io/module_current/td_current_io.h @@ -6,6 +6,7 @@ #include "source_psi/psi.h" #include "source_lcao/module_rt/velocity_op.h" #include "source_lcao/setup_exx.h" +#include "source_lcao/module_rt/td_info.h" #ifdef __EXX #include #endif @@ -24,6 +25,7 @@ void write_current_eachk(const UnitCell& ucell, const Parallel_Orbitals* pv, const LCAO_Orbitals& orb, const Velocity_op* cal_current, + TD_info* td_p, Record_adj& ra); template void write_current(const UnitCell& ucell, @@ -35,6 +37,7 @@ void write_current(const UnitCell& ucell, const Parallel_Orbitals* pv, const LCAO_Orbitals& orb, const Velocity_op* cal_current, + TD_info* td_p, Record_adj& ra); /// @brief func to output current calculated using i[r,H] directly template @@ -47,24 +50,11 @@ void write_current( const K_Vectors& kv, const Parallel_Orbitals* pv, const LCAO_Orbitals& orb, - cal_r_overlap_R& r_calculator, + TD_info* td_p, const hamilt::HContainer* sR, const hamilt::HContainer* hR, const Exx_NAO>& exx_nao ); -/// @brief calculate sum_n[𝜌_(𝑛𝑘,𝜇𝜈)] for current calculation -void cal_tmp_DM_k(const UnitCell& ucell, - elecstate::DensityMatrix, double>& DM_real, - elecstate::DensityMatrix, double>& DM_imag, - const int ik, - const int nspin, - const int is, - const bool reset = true); - -void cal_tmp_DM(const UnitCell& ucell, - elecstate::DensityMatrix, double>& DM_real, - elecstate::DensityMatrix, double>& DM_imag, - const int nspin); void set_rR_from_hR(const UnitCell& ucell, const Grid_Driver& GridD, const LCAO_Orbitals& orb, @@ -95,6 +85,7 @@ template void init_from_hR(const hamilt::HContainer* hR, hamilt::HContainer* aimR); template void cal_velocity_basis_k(const UnitCell& ucell, + TD_info* td_p, const LCAO_Orbitals& orb, const Parallel_Orbitals* pv, const K_Vectors& kv, @@ -114,7 +105,7 @@ void cal_current_comm_k(const UnitCell& ucell, const LCAO_Orbitals& orb, const Parallel_Orbitals* pv, const K_Vectors& kv, - cal_r_overlap_R& r_calculator, + TD_info* td_p, const hamilt::HContainer& sR, const hamilt::HContainer>& hR, const psi::Psi>* psi, diff --git a/source/source_io/module_current/td_current_io_comm.cpp b/source/source_io/module_current/td_current_io_comm.cpp index 66392a80394..2311c08079f 100644 --- a/source/source_io/module_current/td_current_io_comm.cpp +++ b/source/source_io/module_current/td_current_io_comm.cpp @@ -321,6 +321,7 @@ void ModuleIO::add_HR(const hamilt::HContainer* hR, hamilt::HContainer void ModuleIO::cal_velocity_basis_k(const UnitCell& ucell, + TD_info* td_p, const LCAO_Orbitals& orb, const Parallel_Orbitals* pv, const K_Vectors& kv, @@ -373,7 +374,7 @@ void ModuleIO::cal_velocity_basis_k(const UnitCell& ucell, const int nrow = pv->get_row_size(); if (elecstate::H_TDDFT_pw::stype == 2) { - module_rt::folding_HR_td(ucell, hR, hk, kv.kvec_d[ik], TD_info::cart_At, nrow, 1); + module_rt::folding_HR_td(hR, hk, kv.kvec_d[ik], TD_info::cart_At, td_p->get_phase_hybrid(), nrow, 1); } else { @@ -383,7 +384,7 @@ void ModuleIO::cal_velocity_basis_k(const UnitCell& ucell, ModuleBase::GlobalFunc::ZEROS(sk, pv->nloc); if (elecstate::H_TDDFT_pw::stype == 2) { - module_rt::folding_HR_td(ucell, sR, sk, kv.kvec_d[ik], TD_info::cart_At, nrow, 1); + module_rt::folding_HR_td(sR, sk, kv.kvec_d[ik], TD_info::cart_At, td_p->get_phase_hybrid(), nrow, 1); } else { @@ -443,6 +444,7 @@ void ModuleIO::cal_velocity_basis_k(const UnitCell& ucell, partial_hk, kv.kvec_d[ik], TD_info::cart_At, + td_p->get_phase_hybrid(), i_alpha, nrow, 1); @@ -460,6 +462,7 @@ void ModuleIO::cal_velocity_basis_k(const UnitCell& ucell, partial_sk, kv.kvec_d[ik], TD_info::cart_At, + td_p->get_phase_hybrid(), i_alpha, nrow, 1); @@ -490,7 +493,7 @@ void ModuleIO::cal_velocity_basis_k(const UnitCell& ucell, // folding_rR(rR[i_alpha], partial_sk, rk, pv, kv.kvec_d[ik], nrow, 1); if (elecstate::H_TDDFT_pw::stype == 2) { - module_rt::folding_HR_td(ucell, *rR[i_alpha], rk, kv.kvec_d[ik], TD_info::cart_At, nrow, 1); + module_rt::folding_HR_td(*rR[i_alpha], rk, kv.kvec_d[ik], TD_info::cart_At, td_p->get_phase_hybrid(), nrow, 1); } else { @@ -801,7 +804,7 @@ void ModuleIO::cal_current_comm_k(const UnitCell& ucell, const LCAO_Orbitals& orb, const Parallel_Orbitals* pv, const K_Vectors& kv, - cal_r_overlap_R& r_calculator, + TD_info* td_p, const hamilt::HContainer& sR, const hamilt::HContainer>& hR, const psi::Psi>* psi, @@ -830,9 +833,9 @@ void ModuleIO::cal_current_comm_k(const UnitCell& ucell, } } // set rR - set_rR_from_hR(ucell, GridD, orb, pv, r_calculator, &hR, rR); + set_rR_from_hR(ucell, GridD, orb, pv, td_p->r_calculator, &hR, rR); // set velocity_basis_k - cal_velocity_basis_k(ucell, orb, pv, kv, rR, sR, hR, velocity_basis_k); + cal_velocity_basis_k(ucell, td_p, orb, pv, kv, rR, sR, hR, velocity_basis_k); // set velocity_k cal_velocity_matrix(psi, pv, kv, velocity_basis_k, velocity_k); @@ -866,7 +869,7 @@ void ModuleIO::write_current(const UnitCell& ucell, const K_Vectors& kv, const Parallel_Orbitals* pv, const LCAO_Orbitals& orb, - cal_r_overlap_R& r_calculator, + TD_info* td_p, const hamilt::HContainer* sR, const hamilt::HContainer* hR, const Exx_NAO>& exx_nao) @@ -880,7 +883,7 @@ void ModuleIO::write_current(const UnitCell& ucell, full_hR = new hamilt::HContainer>(pv); current_k.resize(kv.get_nks()); sum_HR(ucell, *pv, kv, hR, full_hR, exx_nao); - cal_current_comm_k(ucell, GridD, orb, pv, kv, r_calculator, *sR, *full_hR, psi, pelec, current_k); + cal_current_comm_k(ucell, GridD, orb, pv, kv, td_p, *sR, *full_hR, psi, pelec, current_k); delete full_hR; int nspin0 = 1; @@ -940,7 +943,7 @@ template void ModuleIO::write_current(const UnitCell& ucell, const K_Vectors& kv, const Parallel_Orbitals* pv, const LCAO_Orbitals& orb, - cal_r_overlap_R& r_calculator, + TD_info* td_p, const hamilt::HContainer* sR, const hamilt::HContainer* hR, const Exx_NAO>& exx_nao); @@ -953,7 +956,7 @@ template void ModuleIO::write_current>(const UnitCell& ucel const K_Vectors& kv, const Parallel_Orbitals* pv, const LCAO_Orbitals& orb, - cal_r_overlap_R& r_calculator, + TD_info* td_p, const hamilt::HContainer>* sR, const hamilt::HContainer>* hR, const Exx_NAO>& exx_nao); diff --git a/source/source_lcao/module_operator_lcao/op_exx_lcao.hpp b/source/source_lcao/module_operator_lcao/op_exx_lcao.hpp index 7c1ad176f75..8f6f0f832cd 100644 --- a/source/source_lcao/module_operator_lcao/op_exx_lcao.hpp +++ b/source/source_lcao/module_operator_lcao/op_exx_lcao.hpp @@ -477,6 +477,7 @@ void OperatorEXX>::contributeHk(int ik) *this->Hexxc, *this->hR->get_paraV(), TD_info::td_vel_op->cart_At, + TD_info::td_vel_op->get_phase_hybrid(), this->hsk->get_hk()); } else diff --git a/source/source_lcao/module_operator_lcao/operator_lcao.cpp b/source/source_lcao/module_operator_lcao/operator_lcao.cpp index 1e4a5f728bd..b44d8ae67e7 100644 --- a/source/source_lcao/module_operator_lcao/operator_lcao.cpp +++ b/source/source_lcao/module_operator_lcao/operator_lcao.cpp @@ -276,7 +276,7 @@ void OperatorLCAO::contributeHk(int ik) { const int nrow = this->hsk->get_pv()->get_row_size(); if(PARAM.inp.td_stype == 2) { - module_rt::folding_HR_td(*(TD_info::td_vel_op->get_ucell()), *this->hR, this->hsk->get_hk(), this->kvec_d[ik], TD_info::cart_At, nrow, 1); + module_rt::folding_HR_td(*this->hR, this->hsk->get_hk(), this->kvec_d[ik], TD_info::cart_At, TD_info::td_vel_op->get_phase_hybrid(), nrow, 1); } else { @@ -288,7 +288,7 @@ void OperatorLCAO::contributeHk(int ik) { const int ncol = this->hsk->get_pv()->get_col_size(); if(PARAM.inp.td_stype == 2) { - module_rt::folding_HR_td(*(TD_info::td_vel_op->get_ucell()), *this->hR, this->hsk->get_hk(), this->kvec_d[ik], TD_info::cart_At, ncol, 0); + module_rt::folding_HR_td(*this->hR, this->hsk->get_hk(), this->kvec_d[ik], TD_info::cart_At, TD_info::td_vel_op->get_phase_hybrid(), ncol, 0); } else { diff --git a/source/source_lcao/module_operator_lcao/overlap.cpp b/source/source_lcao/module_operator_lcao/overlap.cpp index 4e031c277aa..5355d0f14dd 100644 --- a/source/source_lcao/module_operator_lcao/overlap.cpp +++ b/source/source_lcao/module_operator_lcao/overlap.cpp @@ -288,7 +288,7 @@ void hamilt::Overlap>::contributeHk(int ik) const int nrow = this->SR->get_atom_pair(0).get_paraV()->get_row_size(); if(PARAM.inp.td_stype == 2) { - module_rt::folding_HR_td(*ucell, *this->SR, this->hsk->get_sk(), this->kvec_d[ik], TD_info::cart_At, nrow, 1); + module_rt::folding_HR_td(*this->SR, this->hsk->get_sk(), this->kvec_d[ik], TD_info::cart_At, TD_info::td_vel_op->get_phase_hybrid(), nrow, 1); } else { @@ -300,7 +300,7 @@ void hamilt::Overlap>::contributeHk(int ik) const int ncol = this->SR->get_atom_pair(0).get_paraV()->get_col_size(); if(PARAM.inp.td_stype == 2) { - module_rt::folding_HR_td(*ucell, *this->SR, this->hsk->get_sk(), this->kvec_d[ik], TD_info::cart_At, ncol, 0); + module_rt::folding_HR_td(*this->SR, this->hsk->get_sk(), this->kvec_d[ik], TD_info::cart_At, TD_info::td_vel_op->get_phase_hybrid(), ncol, 0); } else { diff --git a/source/source_lcao/module_operator_lcao/td_pot_hybrid.cpp b/source/source_lcao/module_operator_lcao/td_pot_hybrid.cpp index ba16cfa0325..015af41fec1 100644 --- a/source/source_lcao/module_operator_lcao/td_pot_hybrid.cpp +++ b/source/source_lcao/module_operator_lcao/td_pot_hybrid.cpp @@ -129,12 +129,13 @@ void hamilt::TD_pot_hybrid>::calculate_HR() const int iat2 = ucell->itia2iat(T2, I2); const ModuleBase::Vector3& R_index2 = adjs.box[ad]; ModuleBase::Vector3 dtau = this->ucell->cal_dtau(iat1, iat2, R_index2); + ModuleBase::Vector3 dR = this->ucell->cal_dtau(0, 0, R_index2); hamilt::BaseMatrix* tmp = this->HR_fixed->find_matrix(iat1, iat2, R_index2); hamilt::BaseMatrix* tmp_overlap = this->SR->find_matrix(iat1, iat2, R_index2); if (tmp != nullptr) { - this->cal_HR_IJR(iat1, iat2, paraV, dtau, tmp->get_pointer(), tmp_overlap->get_pointer()); + this->cal_HR_IJR(iat1, iat2, paraV, dtau, dR, tmp->get_pointer(), tmp_overlap->get_pointer()); } else { @@ -152,6 +153,7 @@ void hamilt::TD_pot_hybrid>::cal_HR_IJR(const int& const int& iat2, const Parallel_Orbitals* paraV, const ModuleBase::Vector3& dtau, + const ModuleBase::Vector3& dR, TR* hr_mat_p, TR* sr_p) { @@ -211,7 +213,7 @@ void hamilt::TD_pot_hybrid>::cal_HR_IJR(const int& for (int ipol = 0; ipol < npol; ipol++) { hr_mat_p[ipol * step_trace] += tmp_r * Et; - hr_mat_p[ipol * step_trace] -= ((dtau + tau1) * Et) * sr_p[ipol * step_trace] * this->ucell->lat0; + hr_mat_p[ipol * step_trace] -= (dR * Et) * sr_p[ipol * step_trace] * this->ucell->lat0; } hr_mat_p += npol; sr_p += npol; diff --git a/source/source_lcao/module_operator_lcao/td_pot_hybrid.h b/source/source_lcao/module_operator_lcao/td_pot_hybrid.h index 751e94f3fe0..7b86218e099 100644 --- a/source/source_lcao/module_operator_lcao/td_pot_hybrid.h +++ b/source/source_lcao/module_operator_lcao/td_pot_hybrid.h @@ -114,6 +114,7 @@ class TD_pot_hybrid> : public OperatorLCAO const int& iat2, const Parallel_Orbitals* paraV, const ModuleBase::Vector3& dtau, + const ModuleBase::Vector3& dR, TR* hr_mat_p, TR* sr_p); diff --git a/source/source_lcao/module_operator_lcao/test/tmp_mocks.cpp b/source/source_lcao/module_operator_lcao/test/tmp_mocks.cpp index 2c27539199a..1a362575aed 100644 --- a/source/source_lcao/module_operator_lcao/test/tmp_mocks.cpp +++ b/source/source_lcao/module_operator_lcao/test/tmp_mocks.cpp @@ -213,8 +213,17 @@ void Numerical_Orbital::set_orbital_info(const int&, // mock of TD_info class TD_info { public: - TD_info() {} + TD_info(): phase_hybrid{} {td_vel_op = this;} ~TD_info() {} static ModuleBase::Vector3 cart_At; + static TD_info* td_vel_op; + std::map, std::complex> phase_hybrid; + const std::map, std::complex>& get_phase_hybrid() const + { + return this->phase_hybrid; + } }; -ModuleBase::Vector3 TD_info::cart_At(0.0, 0.0, 0.0); \ No newline at end of file +ModuleBase::Vector3 TD_info::cart_At(0.0, 0.0, 0.0); +TD_info* TD_info::td_vel_op = nullptr; + +static TD_info mock_td_info; diff --git a/source/source_lcao/module_ri/RI_2D_Comm.h b/source/source_lcao/module_ri/RI_2D_Comm.h index 88b0e3421cc..001f8969b57 100644 --- a/source/source_lcao/module_ri/RI_2D_Comm.h +++ b/source/source_lcao/module_ri/RI_2D_Comm.h @@ -78,6 +78,7 @@ extern std::vector>>> split_m2D_kto const std::vector>>>& Hs, const Parallel_Orbitals& pv, const ModuleBase::Vector3& At, + const std::map, std::complex>& phase_hybrid, TK* hk); template diff --git a/source/source_lcao/module_ri/RI_2D_Comm.hpp b/source/source_lcao/module_ri/RI_2D_Comm.hpp index 7df8d8bb5ac..9d322672375 100644 --- a/source/source_lcao/module_ri/RI_2D_Comm.hpp +++ b/source/source_lcao/module_ri/RI_2D_Comm.hpp @@ -11,6 +11,7 @@ #include "source_base/timer.h" #include "source_lcao/LCAO_domain.h" #include "source_io/module_parameter/parameter.h" +#include "source_lcao/module_rt/td_info.h" #include #include #include @@ -321,6 +322,7 @@ void RI_2D_Comm::add_Hexx_td( const std::vector>>>& Hs, const Parallel_Orbitals& pv, const ModuleBase::Vector3& At, + const std::map, std::complex>& phase_hybrid, TK* hk) { ModuleBase::TITLE("RI_2D_Comm", "add_Hexx_td"); @@ -340,13 +342,10 @@ void RI_2D_Comm::add_Hexx_td( const TA& iat1 = Hs_tmpB.first.first; const TC& cell1 = Hs_tmpB.first.second; const ModuleBase::Vector3 r_index = RI_Util::array3_to_Vector3(cell1); - const ModuleBase::Vector3 dtau = ucell.cal_dtau(iat0, iat1, r_index); - const double arg_td = At * dtau * ucell.lat0; - const std::complex frac = alpha * std::exp(ModuleBase::IMAG_UNIT - * ((ModuleBase::TWO_PI * kv.kvec_c[ik] * (r_index * ucell.latvec)) + arg_td)); + * (ModuleBase::TWO_PI * kv.kvec_c[ik] * (r_index * ucell.latvec))) * phase_hybrid.at(r_index); const RI::Tensor& H = Hs_tmpB.second; for (size_t iw0_b = 0; iw0_b < H.shape[0]; ++iw0_b) diff --git a/source/source_lcao/module_rt/td_folding.cpp b/source/source_lcao/module_rt/td_folding.cpp index ed42b7524e1..bd3b4cd24d0 100644 --- a/source/source_lcao/module_rt/td_folding.cpp +++ b/source/source_lcao/module_rt/td_folding.cpp @@ -2,11 +2,11 @@ #include "source_base/libm/libm.h" namespace module_rt{ template -void folding_HR_td(const UnitCell& ucell, - const hamilt::HContainer& hR, +void folding_HR_td(const hamilt::HContainer& hR, std::complex* hk, const ModuleBase::Vector3& kvec_d_in, const ModuleBase::Vector3& cart_At, + const std::map, std::complex>& phase_hybrid, const int ncol, const int hk_type) { @@ -19,22 +19,14 @@ void folding_HR_td(const UnitCell& ucell, for(int ir = 0;ir < tmp.get_R_size(); ++ir ) { const ModuleBase::Vector3 r_index = tmp.get_R_index(ir); - - //new - //cal tddft phase for hybrid gauge - const int iat1 = tmp.get_atom_i(); - const int iat2 = tmp.get_atom_j(); - ModuleBase::Vector3 dtau = ucell.cal_dtau(iat1, iat2, r_index); - const double arg_td = cart_At * dtau * ucell.lat0; - //new - // cal k_phase // if TK==std::complex, kphase is e^{ikR} const ModuleBase::Vector3 dR(r_index.x, r_index.y, r_index.z); - const double arg = (kvec_d_in * dR) * ModuleBase::TWO_PI + arg_td; + const double arg = (kvec_d_in * dR) * ModuleBase::TWO_PI; double sinp = 0.0, cosp = 0.0; ModuleBase::libm::sincos(arg, &sinp, &cosp); std::complex kphase = std::complex(cosp, sinp); + kphase *= phase_hybrid.at(r_index); tmp.find_R(r_index); tmp.add_to_matrix(hk, ncol, kphase, hk_type); @@ -79,6 +71,7 @@ void folding_partial_HR_td(const UnitCell& ucell, std::complex* hk, const ModuleBase::Vector3& kvec_d_in, const ModuleBase::Vector3& cart_At, + const std::map, std::complex>& phase_hybrid, const int ix, const int ncol, const int hk_type) @@ -93,21 +86,15 @@ void folding_partial_HR_td(const UnitCell& ucell, { const ModuleBase::Vector3 r_index = tmp.get_R_index(ir); - //new - //cal tddft phase for mixing gague - const int iat1 = tmp.get_atom_i(); - const int iat2 = tmp.get_atom_j(); - ModuleBase::Vector3 dtau = ucell.cal_dtau(iat1, iat2, r_index); - const double arg_td = cart_At * dtau * ucell.lat0; - //new // cal k_phase // if TK==std::complex, kphase is e^{ikR} const ModuleBase::Vector3 dR(r_index.x, r_index.y, r_index.z); - const double arg = (kvec_d_in * dR) * ModuleBase::TWO_PI + arg_td; + const double arg = (kvec_d_in * dR) * ModuleBase::TWO_PI; double sinp = 0.0, cosp = 0.0; ModuleBase::libm::sincos(arg, &sinp, &cosp); std::complex kphase = std::complex(cosp, sinp); + kphase *= phase_hybrid.at(r_index); const ModuleBase::Vector3 dR_car = dR * ucell.latvec * ucell.lat0; tmp.find_R(r_index); @@ -116,21 +103,21 @@ void folding_partial_HR_td(const UnitCell& ucell, } } template -void folding_HR_td(const UnitCell& ucell, - const hamilt::HContainer& hR, - std::complex* hk, - const ModuleBase::Vector3& kvec_d_in, - const ModuleBase::Vector3& At, - const int ncol, - const int hk_type); +void folding_HR_td(const hamilt::HContainer& hR, + std::complex* hk, + const ModuleBase::Vector3& kvec_d_in, + const ModuleBase::Vector3& At, + const std::map, std::complex>& phase_hybrid, + const int ncol, + const int hk_type); template -void folding_HR_td>(const UnitCell& ucell, - const hamilt::HContainer>& hR, - std::complex* hk, - const ModuleBase::Vector3& kvec_d_in, - const ModuleBase::Vector3& At, - const int ncol, - const int hk_type); +void folding_HR_td>(const hamilt::HContainer>& hR, + std::complex* hk, + const ModuleBase::Vector3& kvec_d_in, + const ModuleBase::Vector3& At, + const std::map, std::complex>& phase_hybrid, + const int ncol, + const int hk_type); template void folding_partial_HR>(const UnitCell& ucell, const hamilt::HContainer>& hR, @@ -153,6 +140,7 @@ void folding_partial_HR_td>(const UnitCell& ucell, std::complex* hk, const ModuleBase::Vector3& kvec_d_in, const ModuleBase::Vector3& cart_At, + const std::map, std::complex>& phase_hybrid, const int ix, const int ncol, const int hk_type); @@ -162,6 +150,7 @@ void folding_partial_HR_td(const UnitCell& ucell, std::complex* hk, const ModuleBase::Vector3& kvec_d_in, const ModuleBase::Vector3& cart_At, + const std::map, std::complex>& phase_hybrid, const int ix, const int ncol, const int hk_type); diff --git a/source/source_lcao/module_rt/td_folding.h b/source/source_lcao/module_rt/td_folding.h index 486b65af5c8..d61d72488e2 100644 --- a/source/source_lcao/module_rt/td_folding.h +++ b/source/source_lcao/module_rt/td_folding.h @@ -5,13 +5,13 @@ namespace module_rt{ // folding HR to hk, for hybrid gauge template -void folding_HR_td(const UnitCell& ucell, - const hamilt::HContainer& hR, - std::complex* hk, - const ModuleBase::Vector3& kvec_d_in, - const ModuleBase::Vector3& At, - const int ncol, - const int hk_type); +void folding_HR_td(const hamilt::HContainer& hR, + std::complex* hk, + const ModuleBase::Vector3& kvec_d_in, + const ModuleBase::Vector3& At, + const std::map, std::complex>& phase_hybrid, + const int ncol, + const int hk_type); template void folding_partial_HR(const UnitCell& ucell, const hamilt::HContainer& hR, @@ -26,6 +26,7 @@ void folding_partial_HR_td(const UnitCell& ucell, std::complex* hk, const ModuleBase::Vector3& kvec_d_in, const ModuleBase::Vector3& cart_At, + const std::map, std::complex>& phase_hybrid, const int ix, const int ncol, const int hk_type); diff --git a/source/source_lcao/module_rt/td_info.cpp b/source/source_lcao/module_rt/td_info.cpp index a51840cb66f..f4251cc3731 100644 --- a/source/source_lcao/module_rt/td_info.cpp +++ b/source/source_lcao/module_rt/td_info.cpp @@ -1,5 +1,6 @@ #include "td_info.h" +#include "source_base/libm/libm.h" #include "source_estate/module_pot/H_TDDFT_pw.h" #include "source_io/module_parameter/parameter.h" @@ -20,7 +21,6 @@ std::vector> TD_info::At_from_file; TD_info::TD_info(const UnitCell* ucell_in,const Parallel_Orbitals& pv, const LCAO_Orbitals& orb) { - this->ucell = ucell_in; if (init_vecpot_file && istep == -1) { this->read_cart_At(); @@ -41,7 +41,7 @@ TD_info::TD_info(const UnitCell* ucell_in,const Parallel_Orbitals& pv, const LCA this->istep += estep_shift; if(out_current==2||elecstate::H_TDDFT_pw::stype == 2) { - r_calculator.init(*ucell, pv, orb); + r_calculator.init(*ucell_in, pv, orb); } return; } @@ -109,6 +109,19 @@ void TD_info::cal_cart_At(const ModuleBase::Vector3& At) { this->output_cart_At(PARAM.globalv.global_out_dir); } + // update hybrid gauge phase + if(elecstate::H_TDDFT_pw::stype == 2) + { + for(const auto& phase_pair : phase_hybrid) + { + const ModuleBase::Vector3& r_index = phase_pair.first; + ModuleBase::Vector3 dR = double(r_index.x) * a1 + double(r_index.y) * a2 + double(r_index.z) * a3; + const double arg_td = cart_At * dR * lat0; + double sinp, cosp; + ModuleBase::libm::sincos(arg_td, &sinp, &cosp); + phase_hybrid[r_index] = std::complex(cosp, sinp); + } + } } void TD_info::read_cart_At(void) @@ -184,7 +197,30 @@ void TD_info::out_restart_info(const int nstep, return; } +template +void TD_info::initialize_phase_hybrid(const UnitCell& ucell, const hamilt::HContainer* hR) +{ + this->a1 = ucell.a1; + this->a2 = ucell.a2; + this->a3 = ucell.a3; + this->lat0 = ucell.lat0; + for (int i = 0; i < hR->size_atom_pairs(); ++i) + { + hamilt::AtomPair& tmp = hR->get_atom_pair(i); + for(int ir = 0;ir < tmp.get_R_size(); ++ir ) + { + const ModuleBase::Vector3 r_index = tmp.get_R_index(ir); + if(phase_hybrid.count(r_index))continue; + + ModuleBase::Vector3 dR = double(r_index.x) * a1 + double(r_index.y) * a2 + double(r_index.z) * a3; + const double arg_td = cart_At * dR * lat0; + double sinp, cosp; + ModuleBase::libm::sincos(arg_td, &sinp, &cosp); + phase_hybrid[r_index] = std::complex(cosp, sinp); + } + } +} void TD_info::initialize_current_term(const hamilt::HContainer>* HR, const Parallel_Orbitals* paraV) { @@ -230,3 +266,9 @@ void TD_info::destroy_HS_R_td_sparse(void) HR_sparse_td_vel[0].swap(empty_HR_sparse_td_vel_up); HR_sparse_td_vel[1].swap(empty_HR_sparse_td_vel_down); } + + +template +void TD_info::initialize_phase_hybrid>(const UnitCell& ucell, const hamilt::HContainer>* hR); +template +void TD_info::initialize_phase_hybrid(const UnitCell& ucell, const hamilt::HContainer* hR); \ No newline at end of file diff --git a/source/source_lcao/module_rt/td_info.h b/source/source_lcao/module_rt/td_info.h index 92f91ae638b..bccb1f27c49 100644 --- a/source/source_lcao/module_rt/td_info.h +++ b/source/source_lcao/module_rt/td_info.h @@ -55,6 +55,15 @@ class TD_info { return this->current_term[i]; } + // allocate memory for phase_hybrid. + template + void initialize_phase_hybrid(const UnitCell& ucell, const hamilt::HContainer* hR); + + const std::map, std::complex>& get_phase_hybrid() const + { + return this->phase_hybrid; + } + // set velocity HR. void set_velocity_HR(hamilt::HContainer>* HR) { @@ -69,12 +78,6 @@ class TD_info { return istep; } - - const UnitCell* get_ucell() - { - return this->ucell; - } - // For TDDFT velocity gauge, to fix the output of HR std::map, std::map>>> HR_sparse_td_vel[2]; @@ -82,8 +85,12 @@ class TD_info cal_r_overlap_R r_calculator; private: - /// @brief pointer to the unit cell - const UnitCell* ucell = nullptr; + /// @brief lattice vectors, used to calculate the extra phase for hybrid gauge + ModuleBase::Vector3a1, a2, a3; + double lat0; + + /// @brief store time-dependent phase for hybrid gauge + std::map, std::complex> phase_hybrid; /// @brief read At from output file void read_cart_At(); diff --git a/tests/05_rtTDDFT/18_NO_hyb_TDDFT/current_tot.txt.ref b/tests/05_rtTDDFT/18_NO_hyb_TDDFT/current_tot.txt.ref index e19d79a92be..8b2e64bf219 100644 --- a/tests/05_rtTDDFT/18_NO_hyb_TDDFT/current_tot.txt.ref +++ b/tests/05_rtTDDFT/18_NO_hyb_TDDFT/current_tot.txt.ref @@ -1,11 +1,11 @@ -1 -1.2080528880846261e-19 1.2080528880846066e-19 1.2080528880845773e-19 -2 9.6988649203065955e-09 -9.6988647723290744e-09 2.4283263685858898e-06 -3 2.6772041590862983e-08 -2.6772038966231164e-08 7.2673271979666816e-06 -4 2.7319060127512691e-08 -2.7319057936239790e-08 9.6256951598233973e-06 -5 9.5656716023055945e-09 -9.5656714606956783e-09 9.4882668471607465e-06 -6 -1.5805949584040751e-08 1.5805949767647706e-08 9.2880215292310497e-06 -7 -4.7582174004461010e-08 4.7582174205200823e-08 9.0311994708317000e-06 -8 -8.4545887860098823e-08 8.4545888087037147e-08 8.7249297042081756e-06 -9 -1.2566554788492410e-07 1.2566554814454782e-07 8.3765586452439330e-06 -10 -1.7019855732619671e-07 1.7019855762999082e-07 7.9931205849008677e-06 -11 -2.1769144339022118e-07 2.1769144375450180e-07 7.5810111633604005e-06 +1 -1.2080528880845217e-19 1.2080528880844668e-19 1.2080528880847506e-19 +2 9.7416687175202817e-09 -9.7416676366009574e-09 2.4283018795250870e-06 +3 2.6892878119640875e-08 -2.6892868192289296e-08 7.2672237757325802e-06 +4 2.7470324656261670e-08 -2.7470316658463090e-08 9.6255586558291574e-06 +5 9.6933238461489889e-09 -9.6933251237513404e-09 9.4881134038928501e-06 +6 -1.5708609086986498e-08 1.5708607992448572e-08 9.2878533914560557e-06 +7 -4.7517340777625543e-08 4.7517340077239463e-08 9.0310281059803680e-06 +8 -8.4512110564646110e-08 8.4512110408745541e-08 8.7247731748824823e-06 +9 -1.2565959390615339e-07 1.2565959434062147e-07 8.3764366928988363e-06 +10 -1.7021765236158024e-07 1.7021765331323308e-07 7.9930490640377748e-06 +11 -2.1773512779931505e-07 2.1773512905882795e-07 7.5809976831549012e-06 diff --git a/tests/05_rtTDDFT/README b/tests/05_rtTDDFT/README index 1b1ea8caa4f..f898d311361 100644 --- a/tests/05_rtTDDFT/README +++ b/tests/05_rtTDDFT/README @@ -1,4 +1,4 @@ In order to run these tests, you need to read ../integrate/Autotest.sh and set abacus directory in the file. To start running these rt-TDDFT (in LCAO basis) you just need to type: -../integrate/Autotst.sh +../integrate/Autotest.sh diff --git a/tests/08_EXX/14_NO_TDDFT_PBE0/current_tot.txt.ref b/tests/08_EXX/14_NO_TDDFT_PBE0/current_tot.txt.ref index 1248e06af2b..7a529b6b59f 100644 --- a/tests/08_EXX/14_NO_TDDFT_PBE0/current_tot.txt.ref +++ b/tests/08_EXX/14_NO_TDDFT_PBE0/current_tot.txt.ref @@ -1,6 +1,6 @@ -1 -2.4064812749552146e-19 2.4058918971736318e-19 2.4056235004085640e-19 -2 1.3401891378677289e-08 -1.3401891246959429e-08 2.4369369671292067e-06 -3 3.9314244407525817e-08 -3.9314246669204259e-08 7.3211566892522424e-06 -4 4.7934170184878929e-08 -4.7934175958037426e-08 9.7956029317437335e-06 -5 3.4376283041277918e-08 -3.4376284895741181e-08 9.8560550504479175e-06 -6 5.9843794461236954e-09 -5.9843813295237610e-09 9.9145073758045465e-06 +1 -2.4061105780900972e-19 2.4057746355074182e-19 2.4064137152788421e-19 +2 1.3424215192380810e-08 -1.3424215730033235e-08 2.4369257186711371e-06 +3 3.9266797908203167e-08 -3.9266803842385155e-08 7.3211720761214893e-06 +4 4.7558646233794483e-08 -4.7558654016584407e-08 9.7958584306165467e-06 +5 3.3345693326762254e-08 -3.3345695331858395e-08 9.8568251466587576e-06 +6 4.0507435957670880e-09 -4.0507458176112040e-09 9.9160372634246860e-06 diff --git a/tests/08_EXX/14_NO_TDDFT_PBE0/result.ref b/tests/08_EXX/14_NO_TDDFT_PBE0/result.ref index 1b0730251e2..fac54b62e9b 100644 --- a/tests/08_EXX/14_NO_TDDFT_PBE0/result.ref +++ b/tests/08_EXX/14_NO_TDDFT_PBE0/result.ref @@ -1,4 +1,4 @@ -etotref -203.7710904354790 -etotperatomref -101.8855452177 +etotref -203.77109235 +etotperatomref -101.88554617 CompareCurrent_pass 0 totaltimeref 16.38 diff --git a/tests/15_rtTDDFT_GPU/18_NO_hyb_TDDFT_GPU/current_tot.txt.ref b/tests/15_rtTDDFT_GPU/18_NO_hyb_TDDFT_GPU/current_tot.txt.ref index e19d79a92be..8b2e64bf219 100644 --- a/tests/15_rtTDDFT_GPU/18_NO_hyb_TDDFT_GPU/current_tot.txt.ref +++ b/tests/15_rtTDDFT_GPU/18_NO_hyb_TDDFT_GPU/current_tot.txt.ref @@ -1,11 +1,11 @@ -1 -1.2080528880846261e-19 1.2080528880846066e-19 1.2080528880845773e-19 -2 9.6988649203065955e-09 -9.6988647723290744e-09 2.4283263685858898e-06 -3 2.6772041590862983e-08 -2.6772038966231164e-08 7.2673271979666816e-06 -4 2.7319060127512691e-08 -2.7319057936239790e-08 9.6256951598233973e-06 -5 9.5656716023055945e-09 -9.5656714606956783e-09 9.4882668471607465e-06 -6 -1.5805949584040751e-08 1.5805949767647706e-08 9.2880215292310497e-06 -7 -4.7582174004461010e-08 4.7582174205200823e-08 9.0311994708317000e-06 -8 -8.4545887860098823e-08 8.4545888087037147e-08 8.7249297042081756e-06 -9 -1.2566554788492410e-07 1.2566554814454782e-07 8.3765586452439330e-06 -10 -1.7019855732619671e-07 1.7019855762999082e-07 7.9931205849008677e-06 -11 -2.1769144339022118e-07 2.1769144375450180e-07 7.5810111633604005e-06 +1 -1.2080528880845217e-19 1.2080528880844668e-19 1.2080528880847506e-19 +2 9.7416687175202817e-09 -9.7416676366009574e-09 2.4283018795250870e-06 +3 2.6892878119640875e-08 -2.6892868192289296e-08 7.2672237757325802e-06 +4 2.7470324656261670e-08 -2.7470316658463090e-08 9.6255586558291574e-06 +5 9.6933238461489889e-09 -9.6933251237513404e-09 9.4881134038928501e-06 +6 -1.5708609086986498e-08 1.5708607992448572e-08 9.2878533914560557e-06 +7 -4.7517340777625543e-08 4.7517340077239463e-08 9.0310281059803680e-06 +8 -8.4512110564646110e-08 8.4512110408745541e-08 8.7247731748824823e-06 +9 -1.2565959390615339e-07 1.2565959434062147e-07 8.3764366928988363e-06 +10 -1.7021765236158024e-07 1.7021765331323308e-07 7.9930490640377748e-06 +11 -2.1773512779931505e-07 2.1773512905882795e-07 7.5809976831549012e-06