diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index e10171a1d07..c022f0520d3 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -4133,11 +4133,12 @@ These variables are used to control berry phase and wannier90 interface paramete ### out_current -- **Type**: Boolean +- **Type**: Integer - **Description**: - - True: Output current. - - False: Do not output current. -- **Default**: False + - 0: Do not output current. + - 1: Output current using the two-center integral, faster. + - 2: Output current using the matrix commutation, more precise. +- **Default**: 0 ### out_current_k diff --git a/source/Makefile.Objects b/source/Makefile.Objects index d7e5116df8a..22b8e283402 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -553,6 +553,7 @@ OBJS_IO=input_conv.o\ write_dipole.o\ write_init.o\ td_current_io.o\ + td_current_io_comm.o\ write_libxc_r.o\ output_log.o\ output_mat_sparse.o\ diff --git a/source/source_esolver/esolver_ks_lcao_tddft.cpp b/source/source_esolver/esolver_ks_lcao_tddft.cpp index 0dc9ff9e126..652918e1945 100644 --- a/source/source_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/source_esolver/esolver_ks_lcao_tddft.cpp @@ -55,7 +55,7 @@ void ESolver_KS_LCAO_TDDFT::before_all_runners(UnitCell& ucell, cons // Run before_all_runners in ESolver_KS_LCAO ESolver_KS_LCAO, TR>::before_all_runners(ucell, inp); - td_p = new TD_info(&ucell); + td_p = new TD_info(&ucell, this->pv, this->orb_); TD_info::td_vel_op = td_p; totstep += TD_info::estep_shift; @@ -90,7 +90,7 @@ void ESolver_KS_LCAO_TDDFT::runner(UnitCell& ucell, const int istep) ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT SCF"); // Initialize velocity operator for current calculation - if (PARAM.inp.td_stype != 1 && TD_info::out_current) + if (PARAM.inp.td_stype != 1 && TD_info::out_current == 1) { // initialize the velocity operator velocity_mat = new Velocity_op(&ucell, @@ -203,7 +203,7 @@ void ESolver_KS_LCAO_TDDFT::runner(UnitCell& ucell, const int istep) } } - if (PARAM.inp.td_stype != 1 && TD_info::out_current) + if(PARAM.inp.td_stype != 1 && TD_info::out_current == 1) { delete velocity_mat; } @@ -296,6 +296,12 @@ void ESolver_KS_LCAO_TDDFT::hamilt2rho_single(UnitCell& ucell, srho.begin(is, this->chr, this->pw_rho, ucell.symm); } } +#ifdef __EXX + if (GlobalC::exx_info.info_ri.real_number) + this->exx_nao.exd->exx_hamilt2rho(*this->pelec, this->pv, iter); + else + this->exx_nao.exc->exx_hamilt2rho(*this->pelec, this->pv, iter); +#endif // Calculate delta energy this->pelec->f_en.deband = this->pelec->cal_delta_eband(ucell); @@ -474,6 +480,7 @@ void ESolver_KS_LCAO_TDDFT::after_scf(UnitCell& ucell, const int ist std::cout << " Potential (Ry): " << std::setprecision(15) << this->pelec->f_en.etot << std::endl; // Output dipole, current, etc. + auto* hamilt_lcao = dynamic_cast, TR>*>(this->p_hamilt); ModuleIO::ctrl_output_td(ucell, this->chr.rho_save, this->chr.rhopw, @@ -485,8 +492,12 @@ void ESolver_KS_LCAO_TDDFT::after_scf(UnitCell& ucell, const int ist &this->pv, this->orb_, this->velocity_mat, + this->gd, + hamilt_lcao, this->RA, - this->td_p); + this->td_p, + this->exx_nao + ); ModuleBase::timer::tick(this->classname, "after_scf"); } diff --git a/source/source_io/CMakeLists.txt b/source/source_io/CMakeLists.txt index 3322add8bba..540c23df6c7 100644 --- a/source/source_io/CMakeLists.txt +++ b/source/source_io/CMakeLists.txt @@ -34,6 +34,7 @@ list(APPEND objects write_init.cpp write_mlkedf_descriptors.cpp td_current_io.cpp + td_current_io_comm.cpp write_libxc_r.cpp output_log.cpp para_json.cpp diff --git a/source/source_io/ctrl_output_td.cpp b/source/source_io/ctrl_output_td.cpp index d4b7010627d..db0fcc6ecbf 100644 --- a/source/source_io/ctrl_output_td.cpp +++ b/source/source_io/ctrl_output_td.cpp @@ -20,8 +20,12 @@ void ctrl_output_td(const UnitCell& ucell, const Parallel_Orbitals* pv, const LCAO_Orbitals& orb, const Velocity_op* velocity_mat, + const Grid_Driver& grid, + hamilt::HamiltLCAO, TR>* p_hamilt, Record_adj& RA, - TD_info* td_p) + TD_info* td_p, + const Exx_NAO>& exx_nao + ) { ModuleBase::TITLE("ModuleIO", "ctrl_output_td"); @@ -46,7 +50,7 @@ void ctrl_output_td(const UnitCell& ucell, ModuleBase::WARNING_QUIT("ModuleIO::ctrl_output_td", "Failed to cast ElecState to ElecStateLCAO"); } - if (TD_info::out_current) + if (TD_info::out_current == 1) { if (TD_info::out_current_k) { @@ -57,6 +61,10 @@ void ctrl_output_td(const UnitCell& ucell, ModuleIO::write_current(ucell, istep, psi, pelec, kv, intor, pv, orb, velocity_mat, 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); + } // (3) Output file for restart if (PARAM.inp.out_freq_td > 0) // default value of out_freq_td is 0 @@ -88,8 +96,12 @@ template void ctrl_output_td(const UnitCell&, const Parallel_Orbitals*, const LCAO_Orbitals&, const Velocity_op*, + const Grid_Driver&, + hamilt::HamiltLCAO, double>*, Record_adj&, - TD_info*); + TD_info*, + const Exx_NAO>& + ); template void ctrl_output_td>(const UnitCell&, double**, @@ -102,7 +114,11 @@ template void ctrl_output_td>(const UnitCell&, const Parallel_Orbitals*, const LCAO_Orbitals&, const Velocity_op>*, + const Grid_Driver&, + hamilt::HamiltLCAO, std::complex>*, Record_adj&, - TD_info*); + TD_info*, + const Exx_NAO>& + ); } // namespace ModuleIO \ No newline at end of file diff --git a/source/source_io/ctrl_output_td.h b/source/source_io/ctrl_output_td.h index 7f1dac4ed3f..1b57944d87e 100644 --- a/source/source_io/ctrl_output_td.h +++ b/source/source_io/ctrl_output_td.h @@ -11,6 +11,11 @@ #include "source_lcao/module_rt/velocity_op.h" #include "source_lcao/record_adj.h" #include "source_psi/psi.h" +#include "source_lcao/hamilt_lcao.h" +#include "source_lcao/setup_exx.h" +#ifdef __EXX +#include +#endif namespace ModuleIO { @@ -27,8 +32,12 @@ void ctrl_output_td(const UnitCell& ucell, const Parallel_Orbitals* pv, const LCAO_Orbitals& orb, const Velocity_op* velocity_mat, + const Grid_Driver& grid, + hamilt::HamiltLCAO, TR>* p_hamilt, Record_adj& RA, - TD_info* td_p); + TD_info* td_p, + const Exx_NAO>& exx_nao + ); } // namespace ModuleIO diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h index 3a6192bc80f..6dbae92e7f1 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -403,7 +403,7 @@ struct Input_para int out_wfc_lcao = 0; ///< output the wave functions in local basis. bool out_dipole = false; ///< output the dipole or not bool out_efield = false; ///< output the efield or not - bool out_current = false; ///< output the current or not + int out_current = 0; ///< output the current or not bool out_current_k = false; ///< output tddft current for all k points bool out_vecpot = false; ///< output the vector potential or not bool restart_save = false; ///< restart //Peize Lin add 2020-04-04 diff --git a/source/source_io/read_input_item_exx_dftu.cpp b/source/source_io/read_input_item_exx_dftu.cpp index 9889fd11560..580e4d7b4d9 100644 --- a/source/source_io/read_input_item_exx_dftu.cpp +++ b/source/source_io/read_input_item_exx_dftu.cpp @@ -164,6 +164,13 @@ void ReadInput::item_exx() Input_Item item("exx_separate_loop"); item.annotation = "if 1, a two-step method is employed, else it will " "start with a GGA-Loop, and then Hybrid-Loop"; + item.reset_value = [](const Input_Item& item, Parameter& para) { + if (para.input.esolver_type == "tddft" && para.input.exx_separate_loop) + { + GlobalV::ofs_running << "For RT-TDDFT with hybrid functionals, only exx_separate_loop = 0 is supported" << std::endl; + para.input.exx_separate_loop = false; + } + }; read_sync_bool(input.exx_separate_loop); this->add_item(item); } diff --git a/source/source_io/read_input_item_output.cpp b/source/source_io/read_input_item_output.cpp index f5e068a581f..dce309772c7 100644 --- a/source/source_io/read_input_item_output.cpp +++ b/source/source_io/read_input_item_output.cpp @@ -443,7 +443,7 @@ void ReadInput::item_output() { Input_Item item("out_current"); item.annotation = "output current or not"; - read_sync_bool(input.out_current); + read_sync_int(input.out_current); this->add_item(item); } { diff --git a/source/source_io/td_current_io.h b/source/source_io/td_current_io.h index 5158d1de26a..9ed7b6a1115 100644 --- a/source/source_io/td_current_io.h +++ b/source/source_io/td_current_io.h @@ -6,6 +6,10 @@ #include "source_estate/module_dm/density_matrix.h" #include "source_psi/psi.h" #include "source_lcao/module_rt/velocity_op.h" +#include "source_lcao/setup_exx.h" +#ifdef __EXX +#include +#endif namespace ModuleIO { @@ -33,7 +37,22 @@ void write_current(const UnitCell& ucell, const LCAO_Orbitals& orb, const Velocity_op* cal_current, Record_adj& ra); - +/// @brief func to output current calculated using i[r,H] directly +template +void write_current( + const UnitCell& ucell, + const Grid_Driver& GridD, + const int istep, + const psi::Psi>* psi, + const elecstate::ElecState* pelec, + const K_Vectors& kv, + const Parallel_Orbitals* pv, + const LCAO_Orbitals& orb, + cal_r_overlap_R& r_calculator, + 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, @@ -47,7 +66,61 @@ 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, + const Parallel_Orbitals* pv, + cal_r_overlap_R& r_calculator, + const hamilt::HContainer>* hR, + ModuleBase::Vector3*>& rR); +template +void sum_HR( + const UnitCell& ucell, + const Parallel_Orbitals& pv, + const K_Vectors& kv, + const hamilt::HContainer* hR, + hamilt::HContainer>* full_hR, + const Exx_NAO>& exx_nao +); + +template +void add_HR(const hamilt::HContainer* hR, hamilt::HContainer* full_hR); + +void init_from_adj(const UnitCell& ucell, + const Grid_Driver& GridD, + const LCAO_Orbitals& orb, + const Parallel_Orbitals* pv, + std::vector& adjs_all, + ModuleBase::Vector3*>& rR); +template +void init_from_hR(const hamilt::HContainer* hR, hamilt::HContainer* aimR); +template +void cal_velocity_basis_k(const UnitCell& ucell, + const LCAO_Orbitals& orb, + const Parallel_Orbitals* pv, + const K_Vectors& kv, + const ModuleBase::Vector3*>& rR, + const hamilt::HContainer& sR, + const hamilt::HContainer>& hR, + std::vector*>>& velocity_basis_k); +void cal_velocity_matrix(const psi::Psi>* psi, + const Parallel_Orbitals* pv, + const K_Vectors& kv, + const std::vector*>>& velocity_basis_k, + std::vector>& velocity_k); +template +void cal_current_comm_k(const UnitCell& ucell, + const Grid_Driver& GridD, + const LCAO_Orbitals& orb, + const Parallel_Orbitals* pv, + const K_Vectors& kv, + cal_r_overlap_R& r_calculator, + const hamilt::HContainer& sR, + const hamilt::HContainer>& hR, + const psi::Psi>* psi, + const elecstate::ElecState* pelec, + std::vector>& current_k); #endif // __LCAO } // namespace ModuleIO #endif diff --git a/source/source_io/td_current_io_comm.cpp b/source/source_io/td_current_io_comm.cpp new file mode 100644 index 00000000000..394a8524501 --- /dev/null +++ b/source/source_io/td_current_io_comm.cpp @@ -0,0 +1,910 @@ +#include "td_current_io.h" +#include "source_base/timer.h" +#include "source_base/global_function.h" +#include "source_base/global_variable.h" +#include "source_io/module_parameter/parameter.h" +#include "source_lcao/module_rt/td_info.h" +#include "source_base/libm/libm.h" +#include "source_base/parallel_reduce.h" +#include "source_base/tool_threading.h" +#include "source_base/vector3.h" +#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_base/module_external/lapack_connector.h" +#include "source_base/module_external/scalapack_connector.h" +#include "source_lcao/module_hcontainer/hcontainer_funcs.h" +#include "source_lcao/module_rt/td_folding.h" +#ifdef __EXX +#include "source_lcao/module_operator_lcao/op_exx_lcao.h" +#include "source_lcao/module_ri/Exx_LRI.h" +#endif +#ifdef __LCAO +template +void ModuleIO::init_from_hR(const hamilt::HContainer* hR, hamilt::HContainer* aimR) +{ + ModuleBase::TITLE("ModuleIO", "init_from_hR"); + ModuleBase::timer::tick("ModuleIO", "init_from_hR"); + for (int i = 0; i < hR->size_atom_pairs(); i++) + { + hamilt::AtomPair atom_ij = hR->get_atom_pair(i); + const int iat1 = atom_ij.get_atom_i(); + const int iat2 = atom_ij.get_atom_j(); + for (int iR = 0; iR < atom_ij.get_R_size(); iR++) + { + const ModuleBase::Vector3 r_index = atom_ij.get_R_index(iR); + hamilt::AtomPair atom_ij_ta(iat1, iat2, r_index, hR->get_paraV()); + aimR->insert_pair(atom_ij_ta); + } + } + aimR->allocate(nullptr, true); + + ModuleBase::timer::tick("ModuleIO", "init_from_hR"); +} +void ModuleIO::init_from_adj(const UnitCell& ucell, + const Grid_Driver& GridD, + const LCAO_Orbitals& orb, + const Parallel_Orbitals* pv, + std::vector& adjs_all, + ModuleBase::Vector3*>& rR) +{ + ModuleBase::TITLE("ModuleIOTD_mixing_pot", "init_from_adj"); + ModuleBase::timer::tick("ModuleIO", "init_from_adj"); + + auto orb_cutoff_ = orb.cutoffs(); + adjs_all.clear(); + adjs_all.reserve(ucell.nat); + + for (int iat1 = 0; iat1 < ucell.nat; iat1++) + { + auto tau1 = ucell.get_tau(iat1); + int T1, I1; + ucell.iat2iait(iat1, &I1, &T1); + AdjacentAtomInfo adjs; + GridD.Find_atom(ucell, tau1, T1, I1, &adjs); + std::vector is_adj(adjs.adj_num + 1, false); + for (int ad1 = 0; ad1 < adjs.adj_num + 1; ++ad1) + { + const int T2 = adjs.ntype[ad1]; + const int I2 = adjs.natom[ad1]; + const int iat2 = ucell.itia2iat(T2, I2); + if (pv->get_row_size(iat1) <= 0 || pv->get_col_size(iat2) <= 0) + { + continue; + } + const ModuleBase::Vector3& R_index2 = adjs.box[ad1]; + // choose the real adjacent atoms + // Note: the distance of atoms should less than the cutoff radius, + // When equal, the theoretical value of matrix element is zero, + // but the calculated value is not zero due to the numerical error, which would lead to result changes. + if (ucell.cal_dtau(iat1, iat2, R_index2).norm() * ucell.lat0 + < orb_cutoff_[T1] + orb_cutoff_[T2]) + { + is_adj[ad1] = true; + } + } + filter_adjs(is_adj, adjs); + adjs_all.push_back(adjs); + for (int ad = 0; ad < adjs.adj_num + 1; ++ad) + { + const int T2 = adjs.ntype[ad]; + const int I2 = adjs.natom[ad]; + int iat2 = ucell.itia2iat(T2, I2); + ModuleBase::Vector3& R_index = adjs.box[ad]; + hamilt::AtomPair tmp(iat1, iat2, R_index, pv); + for (size_t i_alpha = 0; i_alpha != 3; ++i_alpha) + { + rR[i_alpha]->insert_pair(tmp); + } + } + } + // allocate the memory of BaseMatrix in HR, and set the new values to zero + for (size_t i_alpha = 0; i_alpha != 3; ++i_alpha) + { + rR[i_alpha]->allocate(nullptr, true); + } + ModuleBase::timer::tick("ModuleIO", "init_from_adj"); +} + +void ModuleIO::set_rR_from_hR(const UnitCell& ucell, + const Grid_Driver& GridD, + const LCAO_Orbitals& orb, + const Parallel_Orbitals* pv, + cal_r_overlap_R& r_calculator, + const hamilt::HContainer>* hR, + ModuleBase::Vector3*>& rR) +{ + ModuleBase::TITLE("ModuleIO", "set_rR_from_hR"); + ModuleBase::timer::tick("ModuleIO", "set_rR_from_hR"); + + // init + std::vector adjs_all; + init_from_adj(ucell, GridD, orb, pv, adjs_all, rR); + + for (int iat1 = 0; iat1 < ucell.nat; iat1++) + { + auto tau1 = ucell.get_tau(iat1); + int T1, I1; + ucell.iat2iait(iat1, &I1, &T1); + AdjacentAtomInfo& adjs = adjs_all[iat1]; + for (int ad = 0; ad < adjs.adj_num + 1; ++ad) + { + const int T2 = adjs.ntype[ad]; + const int I2 = adjs.natom[ad]; + const int iat2 = ucell.itia2iat(T2, I2); + const ModuleBase::Vector3& r_index = adjs.box[ad]; + ModuleBase::Vector3 dtau = ucell.cal_dtau(iat1, iat2, r_index); + + Atom& atom1 = ucell.atoms[T1]; + Atom& atom2 = ucell.atoms[T2]; + const int npol = ucell.get_npol(); + + const int* iw2l1 = atom1.iw2l.data(); + const int* iw2n1 = atom1.iw2n.data(); + const int* iw2m1 = atom1.iw2m.data(); + const int* iw2l2 = atom2.iw2l.data(); + const int* iw2n2 = atom2.iw2n.data(); + const int* iw2m2 = atom2.iw2m.data(); + + auto row_indexes = pv->get_indexes_row(iat1); + auto col_indexes = pv->get_indexes_col(iat2); + + const ModuleBase::Vector3& tau1 = ucell.get_tau(iat1); + // std::cout << "tau1: " << tau1 << " tau2: " << GlobalC::ucell.get_tau(iat2) << " r_index: " << r_index + // << std::endl; + const ModuleBase::Vector3 tau2 = tau1 + dtau; + for (int iw1l = 0; iw1l < row_indexes.size(); iw1l += npol) + { + const int iw1 = row_indexes[iw1l] / npol; + const int L1 = iw2l1[iw1]; + const int N1 = iw2n1[iw1]; + const int m1 = iw2m1[iw1]; + + for (int iw2l = 0; iw2l < col_indexes.size(); iw2l += npol) + { + const int iw2 = col_indexes[iw2l] / npol; + const int L2 = iw2l2[iw2]; + const int N2 = iw2n2[iw2]; + const int m2 = iw2m2[iw2]; + + // std::cout<<"L1: "< tmp_r = r_calculator.get_psi_r_psi(tau1 * ucell.lat0, + T1, + L1, + m1, + N1, + tau2 * ucell.lat0, + T2, + L2, + m2, + N2); + for (size_t i_alpha = 0; i_alpha != 3; ++i_alpha) + { + hamilt::BaseMatrix* HlocR = rR[i_alpha]->find_matrix(iat1, iat2, r_index); + if(HlocR!=nullptr)HlocR->add_element(iw1, iw2, tmp_r[i_alpha]); + // if (i_alpha == 2) + // { + // std::cout << "iw1: " << iw1 << " iw2: " << iw2 << " i_alpha: " << i_alpha + // << " tmp_r: " << tmp_r[i_alpha] << std::endl; + // } + } + } + } + } + } + ModuleBase::TITLE("ModuleIO", "set_rR_from_sR"); +} +template +void ModuleIO::sum_HR( + const UnitCell& ucell, + const Parallel_Orbitals& pv, + const K_Vectors& kv, + const hamilt::HContainer* hR, + hamilt::HContainer>* full_hR, + const Exx_NAO>& exx_nao +) +{ + ModuleBase::TITLE("ModuleIO", "sum_HR"); + ModuleBase::timer::tick("ModuleIO", "sum_HR"); + + // init complex full_hR + init_from_hR(hR, full_hR); +#ifdef __EXX + const bool use_cell_nearest = (ModuleBase::Vector3(std::fmod(kv.get_koffset(0), 1.0), + std::fmod(kv.get_koffset(1), 1.0), + std::fmod(kv.get_koffset(2), 1.0)).norm()< 1e-10); + RI::Cell_Nearest cell_nearest; + // reallocate full_hR for BvK used in EXX + if (GlobalC::exx_info.info_global.cal_exx) + { + const std::array Rs_period = {kv.nmp[0], kv.nmp[1], kv.nmp[2]}; + if (use_cell_nearest) + { + // set cell_nearest + std::map> atoms_pos; + for (int iat = 0; iat < ucell.nat; ++iat) + { + atoms_pos[iat] = RI_Util::Vector3_to_array3( + ucell.atoms[ucell.iat2it[iat]].tau[ucell.iat2ia[iat]]); + } + const std::array, 3> latvec = {RI_Util::Vector3_to_array3(ucell.a1), + RI_Util::Vector3_to_array3(ucell.a2), + RI_Util::Vector3_to_array3(ucell.a3)}; + cell_nearest.init(atoms_pos, latvec, Rs_period); + hamilt::reallocate_hcontainer(ucell.nat, full_hR, Rs_period, &cell_nearest); + } + else + hamilt::reallocate_hcontainer(ucell.nat, full_hR, Rs_period); + } +#endif + // add other hR + add_HR(hR, full_hR); + // add velocity complex hR + if (PARAM.inp.td_stype==1) + { + if (TD_info::td_vel_op == nullptr) + { + ModuleBase::WARNING_QUIT("ModuleIO::write_current", "velocity gauge infos is null!"); + } + const hamilt::HContainer>* velocity_hR = TD_info::td_vel_op->get_velocity_HR_pointer(); + add_HR(velocity_hR, full_hR); + } +#ifdef __EXX + // add HexxR to complex full_hR + if (GlobalC::exx_info.info_global.cal_exx) + { + for (size_t is = 0; is != PARAM.inp.nspin; ++is) + { + if (use_cell_nearest) + RI_2D_Comm::add_HexxR(is, GlobalC::exx_info.info_global.hybrid_alpha, exx_nao.exc->get_Hexxs(), pv, PARAM.globalv.npol, *full_hR, &cell_nearest); + else + RI_2D_Comm::add_HexxR(is, GlobalC::exx_info.info_global.hybrid_alpha, exx_nao.exc->get_Hexxs(), pv, PARAM.globalv.npol, *full_hR, nullptr); + } + } +#endif + + ModuleBase::timer::tick("ModuleIO", "sum_HR"); +} + +template +void ModuleIO::add_HR(const hamilt::HContainer* hR, hamilt::HContainer* full_hR) +{ + ModuleBase::TITLE("ModuleIO", "add_HR"); + ModuleBase::timer::tick("ModuleIO", "add_HR"); + + for (int ipair = 0; ipair < hR->size_atom_pairs(); ++ipair) + { + hamilt::AtomPair atom_ij = hR->get_atom_pair(ipair); + const int iat1 = atom_ij.get_atom_i(); + const int iat2 = atom_ij.get_atom_j(); + // loop R-index + for (int iR = 0; iR < atom_ij.get_R_size(); iR++) + { + const ModuleBase::Vector3 r_index = atom_ij.get_R_index(iR); + hamilt::BaseMatrix* full_HlocR = full_hR->find_matrix(iat1, iat2, r_index.x, r_index.y, r_index.z); + const hamilt::BaseMatrix* HlocR = hR->find_matrix(iat1, iat2, r_index.x, r_index.y, r_index.z); + + if (full_HlocR == nullptr || HlocR == nullptr) + ModuleBase::WARNING_QUIT("ModuleIO::add_HR", "HR cannot be nullptr!"); + + for (int i = 0; i < atom_ij.get_row_size(); ++i) + { + for (int j = 0; j < atom_ij.get_col_size(); ++j) + { + Tadd v = HlocR->get_value(i, j); + full_HlocR->add_element(i, j, Tfull(v)); + } + } + } + } + + ModuleBase::timer::tick("ModuleIO", "add_HR"); +} + +// for molecule, if vacuum size is small, the number of R of Hs is smaller than SR +// which may lead to some errors +template +void ModuleIO::cal_velocity_basis_k(const UnitCell& ucell, + const LCAO_Orbitals& orb, + const Parallel_Orbitals* pv, + const K_Vectors& kv, + const ModuleBase::Vector3*>& rR, + const hamilt::HContainer& sR, + const hamilt::HContainer>& hR, + std::vector*>>& velocity_basis_k) +{ + ModuleBase::TITLE("ModuleIO", "cal_velocity_basis_k"); + ModuleBase::timer::tick("ModuleIO", "cal_velocity_basis_k"); +#ifdef __MPI + const int nlocal = PARAM.globalv.nlocal; + const char N_char = 'N'; + const std::complex one_imag = ModuleBase::IMAG_UNIT; + const std::complex neg_one_imag = ModuleBase::NEG_IMAG_UNIT; + const std::complex one_real = ModuleBase::ONE; + const std::complex neg_one_real = ModuleBase::NEG_ONE; + const std::complex zero_complex = ModuleBase::ZERO; + std::cout<nloc<* hk = new std::complex[pv->nloc]; + std::complex* sk = new std::complex[pv->nloc]; + std::complex* partial_hk = new std::complex[pv->nloc]; + std::complex* partial_sk = new std::complex[pv->nloc]; + std::complex* rk = new std::complex[pv->nloc]; + std::complex* h_is = new std::complex[pv->nloc]; + std::complex* h_is_r = new std::complex[pv->nloc]; + std::complex* r_is = new std::complex[pv->nloc]; + std::complex* r_is_h = new std::complex[pv->nloc]; + std::complex* h_is_ps = new std::complex[pv->nloc]; + + // for (size_t i_alpha = 0; i_alpha != 3; ++i_alpha) + // { + // for (int i = 0; i < hR.size_atom_pairs(); ++i) + // { + // hamilt::AtomPair& tmp = rR[i_alpha]->get_atom_pair(i); + // std::cout<<"cal_velocity_basis_k: "<size_atom_pairs()<<" R_size: + // "< r_index = tmp.get_R_index(ir); + // std::cout<<"r_index: "<nloc); + 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); + else + hamilt::folding_HR(hR, hk, kv.kvec_d[ik], nrow, 1); + // 1.2 set S(k) + 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); + else + hamilt::folding_HR(sR, sk, kv.kvec_d[ik], nrow, 1); + // for (int ir = 0; ir < pv->nrow; ir++) + // { + // const int iwt1 = pv->local2global_row(ir); + // const int iat1 = GlobalC::ucell.iwt2iat[iwt1]; + // for (int ic = 0; ic < pv->ncol; ic++) + // { + // const int iwt2 = pv->local2global_col(ic); + // const int iat2 = GlobalC::ucell.iwt2iat[iwt2]; + // const int irc = ic * pv->nrow + ir; + // std::cout << "ik: " << ik << " iat1:" << iat1 << " iat2:" << iat2 << " iwt1: " << iwt1 + // << " iwt2: " << iwt2 << " hk: " << hk[irc] << std::endl; + // } + // } + // 2. set inverse S(k) -> sk will be changed to sk_inv + int* ipiv = new int[pv->nloc]; + int info = 0; + // 2.1 compute ipiv + ScalapackConnector::getrf(nlocal, nlocal, sk, 1, 1, pv->desc, ipiv, &info); + int lwork = -1; + int liwotk = -1; + std::vector> work(1, 0); + std::vector iwork(1, 0); + // 2.2 compute work + ScalapackConnector::getri(nlocal, sk, 1, 1, pv->desc, ipiv, work.data(), &lwork, iwork.data(), &liwotk, &info); + lwork = work[0].real(); + work.resize(lwork, 0); + liwotk = iwork[0]; + iwork.resize(liwotk, 0); + // 2.3 compute inverse matrix of Sk + ScalapackConnector::getri(nlocal, + sk, // return sk^-1 + 1, + 1, + pv->desc, + ipiv, + work.data(), + &lwork, + iwork.data(), + &liwotk, + &info); + delete[] ipiv; + assert(0 == info); + for (size_t i_alpha = 0; i_alpha != 3; ++i_alpha) + { + // 3. set partial_H(k), partial_S(k) and r(k) + // 3.1 set partial_H(k) + ModuleBase::GlobalFunc::ZEROS(partial_hk, pv->nloc); + if (elecstate::H_TDDFT_pw::stype == 2) + module_rt::folding_partial_HR_td(ucell, hR, partial_hk, kv.kvec_d[ik], TD_info::cart_At, i_alpha, nrow, 1); + else + module_rt::folding_partial_HR(ucell, hR, partial_hk, kv.kvec_d[ik], i_alpha, nrow, 1); + // 3.2 set partial S(k) + ModuleBase::GlobalFunc::ZEROS(partial_sk, pv->nloc); + if (elecstate::H_TDDFT_pw::stype == 2) + module_rt::folding_partial_HR_td(ucell, sR, partial_sk, kv.kvec_d[ik], TD_info::cart_At, i_alpha, nrow, 1); + else + module_rt::folding_partial_HR(ucell, sR, partial_sk, kv.kvec_d[ik], i_alpha, nrow, 1); + // if(i_alpha == 2) + // { + // for(int ir=0;ir< pv->nrow; ir++) + // { + // const int iwt1 = pv->local2global_row(ir); + // const int iat1 = GlobalC::ucell.iwt2iat[iwt1]; + // for(int ic=0;ic< pv->ncol; ic++) + // { + // const int iwt2 = pv->local2global_col(ic); + // const int iat2 = GlobalC::ucell.iwt2iat[iwt2]; + // const int irc=ic*pv->nrow + ir; + // std::cout<<"ik: "<nloc); + // 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); + else + hamilt::folding_HR(*rR[i_alpha], rk, kv.kvec_d[ik], nrow, 1); // set r(k) + // if (i_alpha == 2) + // { + // std::cout << "ik: " << ik << " i_alpha: " << i_alpha << std::endl; + // for (int ir = 0; ir < pv->nrow; ir++) + // { + // const int iwt1 = pv->local2global_row(ir); + // const int iat1 = GlobalC::ucell.iwt2iat[iwt1]; + // for (int ic = 0; ic < pv->ncol; ic++) + // { + // const int iwt2 = pv->local2global_col(ic); + // const int iat2 = GlobalC::ucell.iwt2iat[iwt2]; + // const int irc = ic * pv->nrow + ir; + // std::cout << " iat1: " << iat1 << " iat2: " << iat2 << " iw1: " << + // GlobalC::ucell.iwt2iw[iwt1] + // << " iw2: " << GlobalC::ucell.iwt2iw[iwt2] << " rk: " << rk[irc] << std::endl; + // } + // } + // } + // 4. calculate <\vu,k|v_a|\mu,k> = partial_Hk + IMAG_UNIT * (Hk * Sk_inv * rk) - IMAG_UNIT * (rk * Sk_inv * + // Hk) - Hk * Sk_inv * partial_Sk + // 4.1.1 Hk * Sk_inv (note 2.) + ModuleBase::GlobalFunc::ZEROS(h_is, pv->nloc); + ScalapackConnector::gemm(N_char, + N_char, + nlocal, + nlocal, + nlocal, + one_real, + hk, + 1, + 1, + pv->desc, + sk, + 1, + 1, + pv->desc, + zero_complex, + h_is, + 1, + 1, + pv->desc); + // 4.1.2 (Hk * Sk_inv) * rk + ModuleBase::GlobalFunc::ZEROS(h_is_r, pv->nloc); + ScalapackConnector::gemm(N_char, + N_char, + nlocal, + nlocal, + nlocal, + one_real, + h_is, + 1, + 1, + pv->desc, + rk, + 1, + 1, + pv->desc, + zero_complex, + h_is_r, + 1, + 1, + pv->desc); + // 4.2.1 rk * Sk_inv (note 2.) + ModuleBase::GlobalFunc::ZEROS(r_is, pv->nloc); + ScalapackConnector::gemm(N_char, + N_char, + nlocal, + nlocal, + nlocal, + one_real, + rk, + 1, + 1, + pv->desc, + sk, + 1, + 1, + pv->desc, + zero_complex, + r_is, + 1, + 1, + pv->desc); + // 4.2.2 (rk * Sk_inv) * Hk + ModuleBase::GlobalFunc::ZEROS(r_is_h, pv->nloc); + ScalapackConnector::gemm(N_char, + N_char, + nlocal, + nlocal, + nlocal, + one_real, + r_is, + 1, + 1, + pv->desc, + hk, + 1, + 1, + pv->desc, + zero_complex, + r_is_h, + 1, + 1, + pv->desc); + // 4.3.1 (Hk * Sk_inv) * partial_Sk + ModuleBase::GlobalFunc::ZEROS(h_is_ps, pv->nloc); + ScalapackConnector::gemm(N_char, + N_char, + nlocal, + nlocal, + nlocal, + one_real, + h_is, + 1, + 1, + pv->desc, + partial_sk, + 1, + 1, + pv->desc, + zero_complex, + h_is_ps, + 1, + 1, + pv->desc); + // 4.4 h_is_r will be changed to partial_Hk + IMAG_UNIT * (Hk * Sk_inv * rk) + ScalapackConnector::geadd('N', + nlocal, + nlocal, + one_real, + partial_hk, + 1, + 1, + pv->desc, + one_imag, + h_is_r, + 1, + 1, + pv->desc); + // 4.5 r_is_h will be changed to h_is_r - IMAG_UNIT * (rk * Sk_inv * Hk) + ScalapackConnector::geadd('N', + nlocal, + nlocal, + one_real, + h_is_r, + 1, + 1, + pv->desc, + neg_one_imag, + r_is_h, + 1, + 1, + pv->desc); + // 4.6 h_is_ps will be changed to r_is_h - Hk * Sk_inv * partial_Sk + ScalapackConnector::geadd('N', + nlocal, + nlocal, + one_real, + r_is_h, + 1, + 1, + pv->desc, + neg_one_real, + h_is_ps, + 1, + 1, + pv->desc); + // 5. copy h_is_ps to velocity_basis_k[ik][i_alpha] + BlasConnector::copy(pv->nloc, h_is_ps, 1, velocity_basis_k[ik][i_alpha], 1); + // if(i_alpha == 2) + // { + // for(int ir=0;ir< pv->nrow; ir++) + // { + // const int iwt1 = pv->local2global_row(ir); + // const int iat1 = GlobalC::ucell.iwt2iat[iwt1]; + // for(int ic=0;ic< pv->ncol; ic++) + // { + // const int iwt2 = pv->local2global_col(ic); + // const int iat2 = GlobalC::ucell.iwt2iat[iwt2]; + // const int irc=ic*pv->nrow + ir; + // std::cout<<"ik: "<>* psi, + const Parallel_Orbitals* pv, + const K_Vectors& kv, + const std::vector*>>& velocity_basis_k, + std::vector>& velocity_k) +{ + ModuleBase::TITLE("ModuleIO", "cal_velocity_matrix"); + ModuleBase::timer::tick("ModuleIO", "cal_velocity_matrix"); +#ifdef __MPI + const char N_char = 'N'; + const char C_char = 'C'; + const std::complex one_real = ModuleBase::ONE; + const std::complex zero_complex = ModuleBase::ZERO; + const double zero_double = 0.0; + const int nlocal = PARAM.globalv.nlocal; + const int nbands = PARAM.inp.nbands; + std::complex* vk_c = new std::complex[pv->ncol_bands * pv->nrow_bands]; // local one + std::complex* v_c = new std::complex[pv->nloc_wfc]; + + for (int ik = 0; ik < kv.get_nks(); ik++) + { + // 1. set C + psi->fix_k(ik); + // 2. set <\Psi_{n,\mu}|v_{\mu,\nu}|\Psi_{m,\nu}> = C^\dagger_{n,\mu} * v_{\mu,\nu} * C_{\nu,m} + for (size_t i_alpha = 0; i_alpha != 3; ++i_alpha) + { + ModuleBase::GlobalFunc::ZEROS(vk_c, pv->ncol_bands * pv->nrow_bands); + ModuleBase::GlobalFunc::ZEROS(v_c, pv->nloc_wfc); + // v_c_{\mu,m} = v_{\mu,\nu} * C_{\nu,m} + ScalapackConnector::gemm(N_char, + N_char, + nlocal, + nbands, + nlocal, + one_real, + velocity_basis_k[ik][i_alpha], + 1, + 1, + pv->desc, + psi[0].get_pointer(), + 1, + 1, + pv->desc_wfc, + zero_complex, + v_c, + 1, + 1, + pv->desc_wfc); + // velocity_k_{n,m} = C^\dagger_{n,\mu} * v_c_{\mu,m} + ScalapackConnector::gemm(C_char, + N_char, + nbands, + nbands, + nlocal, + one_real, + psi[0].get_pointer(), + 1, + 1, + pv->desc_wfc, + v_c, + 1, + 1, + pv->desc_wfc, + zero_complex, + vk_c, + 1, + 1, + pv->desc_Eij); + + for (int ir = 0; ir < PARAM.inp.nbands; ++ir) + { + // const int iwt1 = pv->local2global_row(ir); + // const int iat1 = GlobalC::ucell.iwt2iat[iwt1]; + for (int ic = 0; ic < PARAM.inp.nbands; ++ic) + { + const int irc = ic * pv->nrow + ir; + if (pv->in_this_processor(ir, ic)) + { + // const int iwt2 = pv->local2global_col(ic); + // const int iat2 = GlobalC::ucell.iwt2iat[iwt2]; + velocity_k[ik][i_alpha](ir, ic) = vk_c[irc]; + // if (i_alpha == 0) + // { + // std::cout<<"ik: "< +void ModuleIO::cal_current_comm_k(const UnitCell& ucell, + const Grid_Driver& GridD, + const LCAO_Orbitals& orb, + const Parallel_Orbitals* pv, + const K_Vectors& kv, + cal_r_overlap_R& r_calculator, + const hamilt::HContainer& sR, + const hamilt::HContainer>& hR, + const psi::Psi>* psi, + const elecstate::ElecState* pelec, + std::vector>& current_k) +{ + ModuleBase::TITLE("ModuleIO", "cal_current_exx"); + ModuleBase::timer::tick("ModuleIO", "cal_current_exx"); + + const int nlocal = PARAM.globalv.nlocal; + const int nbands = PARAM.inp.nbands; + // init + ModuleBase::Vector3*> rR(nullptr, nullptr, nullptr); + std::vector*>> velocity_basis_k; + std::vector> velocity_k; + velocity_basis_k.resize(kv.get_nks()); + velocity_k.resize(kv.get_nks()); + for (size_t i_alpha = 0; i_alpha != 3; ++i_alpha) + { + rR[i_alpha] = new hamilt::HContainer(pv); + for (int ik = 0; ik < kv.get_nks(); ik++) + { + velocity_basis_k[ik][i_alpha] = new std::complex[pv->nloc]; + ModuleBase::GlobalFunc::ZEROS(velocity_basis_k[ik][i_alpha], pv->nloc); + velocity_k[ik][i_alpha].create(nbands, nbands); + } + } + // set rR + set_rR_from_hR(ucell, GridD, orb, pv, r_calculator, &hR, rR); + // set velocity_basis_k + cal_velocity_basis_k(ucell, orb, pv, kv, rR, sR, hR, velocity_basis_k); + // set velocity_k + cal_velocity_matrix(psi, pv, kv, velocity_basis_k, velocity_k); + + // sum n and m for current_k + for (size_t ik = 0; ik != kv.get_nks(); ++ik) + for (size_t i_alpha = 0; i_alpha != 3; ++i_alpha) + { + for (size_t ib = 0; ib != PARAM.inp.nbands; ++ib) + current_k[ik][i_alpha] -= pelec->wg(ik, ib) * velocity_k[ik][i_alpha](ib, ib).real() / 2.0; // for unit + } + for (size_t i_alpha = 0; i_alpha < 3; ++i_alpha) + { + delete rR[i_alpha]; + for (int ik = 0; ik < kv.get_nks(); ik++) + delete[] velocity_basis_k[ik][i_alpha]; + } + + ModuleBase::TITLE("ModuleIO", "cal_current_exx"); +} +template +void ModuleIO::write_current( + const UnitCell& ucell, + const Grid_Driver& GridD, + const int istep, + const psi::Psi>* psi, + const elecstate::ElecState* pelec, + const K_Vectors& kv, + const Parallel_Orbitals* pv, + const LCAO_Orbitals& orb, + cal_r_overlap_R& r_calculator, + const hamilt::HContainer* sR, + const hamilt::HContainer* hR, + const Exx_NAO>& exx_nao +) +{ + ModuleBase::TITLE("ModuleIO", "write_current"); + ModuleBase::timer::tick("ModuleIO", "write_current"); + double omega = ucell.omega; + + std::vector> current_k; + hamilt::HContainer>* full_hR; + 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); + delete full_hR; + + int nspin0 = 1; + if (PARAM.inp.nspin == 2) + { + nspin0 = 2; + } + for (int is = 0; is < nspin0; ++is) + { + for (int ik = 0; ik < kv.get_nks(); ik++) + { + if (is == kv.isk[ik]) + { + if (GlobalV::MY_RANK == 0 && TD_info::out_current_k) + { + std::string filename = PARAM.globalv.global_out_dir + "currents" + std::to_string(is) + + "k" + std::to_string(ik) + "comm.txt"; + std::ofstream fout; + fout.open(filename, std::ios::app); + fout << std::setprecision(16); + fout << std::scientific; + fout << istep << " " << current_k[ik][0] / omega << " " << current_k[ik][1] / omega << " " + << current_k[ik][2] / omega << std::endl; + fout.close(); + } + } + } + } + + ModuleBase::Vector3 current_total; + for (int dir = 0; dir < 3; dir++) + for (int ik = 0; ik < kv.get_nks(); ik++) + current_total[dir] += current_k[ik][dir]; + if (GlobalV::MY_RANK == 0) + { + std::string filename = PARAM.globalv.global_out_dir + "current_tot_comm.txt"; + std::ofstream fout; + fout.open(filename, std::ios::app); + fout << std::setprecision(16); + fout << std::scientific; + fout << istep << " " << current_total[0] / omega << " " << current_total[1] / omega << " " + << current_total[2] / omega << std::endl; + fout.close(); + } + + ModuleBase::timer::tick("ModuleIO", "write_current"); +} +template +void ModuleIO::write_current( + const UnitCell& ucell, + const Grid_Driver& GridD, + const int istep, + const psi::Psi>* psi, + const elecstate::ElecState* pelec, + const K_Vectors& kv, + const Parallel_Orbitals* pv, + const LCAO_Orbitals& orb, + cal_r_overlap_R& r_calculator, + const hamilt::HContainer* sR, + const hamilt::HContainer* hR, + const Exx_NAO>& exx_nao +); + +template +void ModuleIO::write_current>( + const UnitCell& ucell, + const Grid_Driver& GridD, + const int istep, + const psi::Psi>* psi, + const elecstate::ElecState* pelec, + const K_Vectors& kv, + const Parallel_Orbitals* pv, + const LCAO_Orbitals& orb, + cal_r_overlap_R& r_calculator, + const hamilt::HContainer>* sR, + const hamilt::HContainer>* hR, + const Exx_NAO>& exx_nao +); +#endif //__LCAO diff --git a/source/source_io/test/for_testing_input_conv.h b/source/source_io/test/for_testing_input_conv.h index 8cca07361f1..a9951642d92 100644 --- a/source/source_io/test/for_testing_input_conv.h +++ b/source/source_io/test/for_testing_input_conv.h @@ -30,7 +30,7 @@ #undef private bool berryphase::berry_phase_flag = false; -bool TD_Velocity::out_current; +int TD_Velocity::out_current; bool TD_Velocity::out_current_k; bool TD_Velocity::out_vecpot; bool TD_Velocity::init_vecpot_file; diff --git a/source/source_lcao/hamilt_lcao.cpp b/source/source_lcao/hamilt_lcao.cpp index d474db3eba8..d0a5cb62fdf 100644 --- a/source/source_lcao/hamilt_lcao.cpp +++ b/source/source_lcao/hamilt_lcao.cpp @@ -430,17 +430,35 @@ HamiltLCAO::HamiltLCAO(const UnitCell& ucell, // Peize Lin add 2016-12-03 // set xc type before the first cal of xc in pelec->init_scf // and calculate Cs, Vs - Operator* exx = new OperatorEXX>(this->hsk, - this->hR, - ucell, - *kv, - Hexxd, - Hexxc, - Add_Hexx_Type::R, - istep, - exx_two_level_step, - !GlobalC::restart.info_load.restart_exx - && GlobalC::restart.info_load.load_H); + Operator* exx; + if (PARAM.inp.esolver_type == "tddft") + { + exx = new OperatorEXX>(this->hsk, + this->hR, + ucell, + *this->kv, + Hexxd, + Hexxc, + Add_Hexx_Type::k, + istep, + exx_two_level_step, + !GlobalC::restart.info_load.restart_exx + && GlobalC::restart.info_load.load_H); + } + else + { + exx = new OperatorEXX>(this->hsk, + this->hR, + ucell, + *kv, + Hexxd, + Hexxc, + Add_Hexx_Type::R, + istep, + exx_two_level_step, + !GlobalC::restart.info_load.restart_exx + && GlobalC::restart.info_load.load_H); + } this->getOperator()->add(exx); } #endif diff --git a/source/source_lcao/module_deepks/test/CMakeLists.txt b/source/source_lcao/module_deepks/test/CMakeLists.txt index 36327d3f1da..ed7d1ffe8c6 100644 --- a/source/source_lcao/module_deepks/test/CMakeLists.txt +++ b/source/source_lcao/module_deepks/test/CMakeLists.txt @@ -39,9 +39,8 @@ add_executable( ../../module_operator_lcao/deepks_lcao.cpp ../../module_operator_lcao/operator_lcao.cpp ../../../source_hamilt/operator.cpp - ../../module_rt/td_info.cpp ../../module_rt/td_folding.cpp - ../../../source_estate/module_pot/H_TDDFT_pw.cpp + mock_tdinfo.cpp ) target_link_libraries( diff --git a/source/source_lcao/module_deepks/test/mock_tdinfo.cpp b/source/source_lcao/module_deepks/test/mock_tdinfo.cpp new file mode 100644 index 00000000000..5c1af5ed8f2 --- /dev/null +++ b/source/source_lcao/module_deepks/test/mock_tdinfo.cpp @@ -0,0 +1,17 @@ +#include "source_base/vector3.h" +#include "source_cell/unitcell.h" +// mock of TD_info +class TD_info { +public: + TD_info() {} + ~TD_info() {} + const UnitCell* get_ucell() + { + return nullptr; + } + static ModuleBase::Vector3 cart_At; + static TD_info* td_vel_op; +}; +TD_info td_info; +TD_info* TD_info::td_vel_op = &td_info; +ModuleBase::Vector3 TD_info::cart_At(0.0, 0.0, 0.0); \ No newline at end of file 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 31e219fbd71..1d7e61409fb 100644 --- a/source/source_lcao/module_operator_lcao/op_exx_lcao.hpp +++ b/source/source_lcao/module_operator_lcao/op_exx_lcao.hpp @@ -7,6 +7,7 @@ #include "source_lcao/module_ri/RI_2D_Comm.h" #include "source_hamilt/module_xc/xc_functional.h" #include "source_io/restart_exx_csr.h" +#include "source_lcao/module_rt/td_info.h" #include "source_io/restart.h" namespace hamilt @@ -339,7 +340,7 @@ void OperatorEXX>::contributeHR() template void OperatorEXX>::contributeHk(int ik) { - ModuleBase::TITLE("OperatorEXX", "constributeHR"); + ModuleBase::TITLE("OperatorEXX", "constributeHk"); // Peize Lin add 2016-12-03 if (PARAM.inp.calculation != "nscf" && this->two_level_step != nullptr && *this->two_level_step == 0 && !this->restart) { return; } //in the non-exx loop, do nothing @@ -369,26 +370,40 @@ void OperatorEXX>::contributeHk(int ik) } } // cal H(k) from H(R) normally - - if (GlobalC::exx_info.info_ri.real_number) { - RI_2D_Comm::add_Hexx( - ucell, - this->kv, - ik, - GlobalC::exx_info.info_global.hybrid_alpha, - *this->Hexxd, - *this->hR->get_paraV(), - this->hsk->get_hk()); - } else { - RI_2D_Comm::add_Hexx( + if(PARAM.inp.esolver_type == "tddft" && PARAM.inp.td_stype == 2) + { + RI_2D_Comm::add_Hexx_td( ucell, this->kv, ik, GlobalC::exx_info.info_global.hybrid_alpha, *this->Hexxc, *this->hR->get_paraV(), + TD_info::td_vel_op->cart_At, this->hsk->get_hk()); -} + } + else + { + if (GlobalC::exx_info.info_ri.real_number) { + RI_2D_Comm::add_Hexx( + ucell, + this->kv, + ik, + GlobalC::exx_info.info_global.hybrid_alpha, + *this->Hexxd, + *this->hR->get_paraV(), + this->hsk->get_hk()); + } else { + RI_2D_Comm::add_Hexx( + ucell, + this->kv, + ik, + GlobalC::exx_info.info_global.hybrid_alpha, + *this->Hexxc, + *this->hR->get_paraV(), + this->hsk->get_hk()); + } + } } } diff --git a/source/source_lcao/module_operator_lcao/operator_lcao.cpp b/source/source_lcao/module_operator_lcao/operator_lcao.cpp index 09c63faefb2..ec86d103e54 100644 --- a/source/source_lcao/module_operator_lcao/operator_lcao.cpp +++ b/source/source_lcao/module_operator_lcao/operator_lcao.cpp @@ -194,10 +194,14 @@ void OperatorLCAO::init(const int ik_in) { case calculation_type::lcao_exx: { //update HR first - if (!this->hr_done) + if (!this->hr_done && PARAM.inp.esolver_type != "tddft") { this->contributeHR(); } + else if(PARAM.inp.esolver_type == "tddft") + { + this->contributeHk(ik_in); + } //update HK next //in cal_type=lcao_exx, HK only need to update from one node @@ -271,7 +275,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(*this->hR, this->hsk->get_hk(), this->kvec_d[ik], nrow, 1, TD_info::td_vel_op->get_ucell(), TD_info::cart_At); + 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); } else { @@ -283,7 +287,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(*this->hR, this->hsk->get_hk(), this->kvec_d[ik], ncol, 0, TD_info::td_vel_op->get_ucell(), TD_info::cart_At); + 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); } else { diff --git a/source/source_lcao/module_operator_lcao/overlap_new.cpp b/source/source_lcao/module_operator_lcao/overlap_new.cpp index 1d842e56f28..09c2b07a346 100644 --- a/source/source_lcao/module_operator_lcao/overlap_new.cpp +++ b/source/source_lcao/module_operator_lcao/overlap_new.cpp @@ -236,7 +236,7 @@ void hamilt::OverlapNew>::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(*this->SR, this->hsk->get_sk(), this->kvec_d[ik], nrow, 1, ucell, TD_info::cart_At); + module_rt::folding_HR_td(*ucell, *this->SR, this->hsk->get_sk(), this->kvec_d[ik], TD_info::cart_At, nrow, 1); } else { @@ -248,7 +248,7 @@ void hamilt::OverlapNew>::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(*this->SR, this->hsk->get_sk(), this->kvec_d[ik], ncol, 0, ucell, TD_info::cart_At); + module_rt::folding_HR_td(*ucell, *this->SR, this->hsk->get_sk(), this->kvec_d[ik], TD_info::cart_At, ncol, 0); } else { diff --git a/source/source_lcao/module_operator_lcao/td_ekinetic_lcao.cpp b/source/source_lcao/module_operator_lcao/td_ekinetic_lcao.cpp index 9ef6bbe126f..c8e3f8d3071 100644 --- a/source/source_lcao/module_operator_lcao/td_ekinetic_lcao.cpp +++ b/source/source_lcao/module_operator_lcao/td_ekinetic_lcao.cpp @@ -94,7 +94,7 @@ void TDEkinetic>::calculate_HR() hamilt::BaseMatrix>* tmp = this->hR_tmp->find_matrix(iat1, iat2, R_index2); if (tmp != nullptr) { - if (TD_info::out_current) + if (TD_info::out_current == 1) { std::complex* tmp_c[3] = {nullptr, nullptr, nullptr}; for (int i = 0; i < 3; i++) @@ -341,10 +341,14 @@ void TDEkinetic>::contributeHR() static_cast*>(this->next_sub_op)->set_HR_fixed(this->hR_tmp); } // initialize current term if needed - if (TD_info::out_current) + if (TD_info::out_current == 1) { TD_info::td_vel_op->initialize_current_term(this->hR_tmp, paraV); } + if (TD_info::out_current==2) + { + TD_info::td_vel_op->set_velocity_HR(this->hR_tmp); + } // calculate the values in hR_tmp this->update_td(); this->hR_tmp->set_zero(); diff --git a/source/source_lcao/module_operator_lcao/td_nonlocal_lcao.cpp b/source/source_lcao/module_operator_lcao/td_nonlocal_lcao.cpp index 0bbcb6fd528..42254ca531a 100644 --- a/source/source_lcao/module_operator_lcao/td_nonlocal_lcao.cpp +++ b/source/source_lcao/module_operator_lcao/td_nonlocal_lcao.cpp @@ -154,7 +154,7 @@ void hamilt::TDNonlocal>::calculate_HR() const Parallel_Orbitals* paraV = this->hR_tmp->get_atom_pair(0).get_paraV(); const int npol = this->ucell->get_npol(); - const int nlm_dim = TD_info::out_current ? 4 : 1; + const int nlm_dim = TD_info::out_current==1 ? 4 : 1; // 1. calculate for each pair of atoms for (int iat0 = 0; iat0 < this->ucell->nat; iat0++) @@ -220,7 +220,7 @@ void hamilt::TDNonlocal>::calculate_HR() tau0 * this->ucell->lat0, T0, cart_At, - TD_info::out_current); + TD_info::out_current==1); for (int dir = 0; dir < nlm_dim; dir++) { nlm_tot[ad][dir].insert({all_indexes[iw1l], nlm[dir]}); @@ -287,7 +287,7 @@ void hamilt::TDNonlocal>::calculate_HR() // if not found , skip this pair of atoms if (tmp != nullptr) { - if (TD_info::out_current) + if (TD_info::out_current==1) { std::complex* tmp_c[3] = {nullptr, nullptr, nullptr}; for (int ii = 0; ii < 3; ii++) @@ -337,7 +337,7 @@ void hamilt::TDNonlocal>::cal_HR_IJR( std::complex* data_pointer, std::complex** data_pointer_c) { - const int nlm_dim = TD_info::out_current ? 4 : 1; + const int nlm_dim = TD_info::out_current==1 ? 4 : 1; // npol is the number of polarizations, // 1 for non-magnetic (one Hamiltonian matrix only has spin-up or spin-down), // 2 for magnetic (one Hamiltonian matrix has both spin-up and spin-down) 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 f0d68d01334..72f59490d76 100644 --- a/source/source_lcao/module_operator_lcao/td_pot_hybrid.cpp +++ b/source/source_lcao/module_operator_lcao/td_pot_hybrid.cpp @@ -6,10 +6,6 @@ #include "source_lcao/module_operator_lcao/operator_lcao.h" #include "source_lcao/module_hcontainer/hcontainer_funcs.h" -// Constructor -template -cal_r_overlap_R hamilt::TD_pot_hybrid>::r_calculator; - template hamilt::TD_pot_hybrid>::TD_pot_hybrid( HS_Matrix_K* hsk_in, @@ -206,7 +202,7 @@ void hamilt::TD_pot_hybrid>::cal_HR_IJR(const int& const int N2 = iw2n2[iw2]; const int m2 = iw2m2[iw2]; - ModuleBase::Vector3 tmp_r = r_calculator.get_psi_r_psi(tau1 * this->ucell->lat0, T1, L1, m1, N1, tau2 * this->ucell->lat0, T2, L2, m2, N2); + ModuleBase::Vector3 tmp_r = r_calculator->get_psi_r_psi(tau1 * this->ucell->lat0, T1, L1, m1, N1, tau2 * this->ucell->lat0, T2, L2, m2, N2); // convert m (0,1,...2l) to M (-l, -l+1, ..., l-1, l) int M2 = (m2 % 2 == 0) ? -m2 / 2 : (m2 + 1) / 2; @@ -227,12 +223,7 @@ template void hamilt::TD_pot_hybrid>::init_td() { // initialize the r_calculator - if(TD_info::td_vel_op->get_istep()==(TD_info::estep_shift-1)) - { - //std::cout << "init_r_overlap" <hR->get_paraV(), orb_); - } - //hk_hybrid.resize(this->hR->get_paraV()->nloc); + this->r_calculator = &TD_info::td_vel_op->r_calculator; } template void hamilt::TD_pot_hybrid>::update_td() 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 82f17b6fed3..5790dd53a45 100644 --- a/source/source_lcao/module_operator_lcao/td_pot_hybrid.h +++ b/source/source_lcao/module_operator_lcao/td_pot_hybrid.h @@ -86,7 +86,7 @@ class TD_pot_hybrid> : public OperatorLCAO bool HR_fixed_done = false; //tddft part - static cal_r_overlap_R r_calculator; + cal_r_overlap_R* r_calculator; //ETD //std::vector> hk_hybrid; //ETD diff --git a/source/source_lcao/module_operator_lcao/test/CMakeLists.txt b/source/source_lcao/module_operator_lcao/test/CMakeLists.txt index c6c7f379e3e..8ea1c830793 100644 --- a/source/source_lcao/module_operator_lcao/test/CMakeLists.txt +++ b/source/source_lcao/module_operator_lcao/test/CMakeLists.txt @@ -9,9 +9,7 @@ AddTest( ../../../source_basis/module_ao/parallel_orbitals.cpp ../../../source_basis/module_ao/ORB_atomic_lm.cpp tmp_mocks.cpp ../../../source_hamilt/operator.cpp - ../../module_rt/td_info.cpp ../../module_rt/td_folding.cpp - ../../../source_estate/module_pot/H_TDDFT_pw.cpp ) AddTest( @@ -22,9 +20,7 @@ AddTest( ../../../source_basis/module_ao/parallel_orbitals.cpp ../../../source_basis/module_ao/ORB_atomic_lm.cpp tmp_mocks.cpp ../../../source_hamilt/operator.cpp - ../../module_rt/td_info.cpp ../../module_rt/td_folding.cpp - ../../../source_estate/module_pot/H_TDDFT_pw.cpp ) AddTest( @@ -35,9 +31,7 @@ AddTest( ../../../source_basis/module_ao/parallel_orbitals.cpp ../../../source_basis/module_ao/ORB_atomic_lm.cpp tmp_mocks.cpp ../../../source_hamilt/operator.cpp - ../../module_rt/td_info.cpp ../../module_rt/td_folding.cpp - ../../../source_estate/module_pot/H_TDDFT_pw.cpp ) AddTest( 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 09afff3eb72..2c27539199a 100644 --- a/source/source_lcao/module_operator_lcao/test/tmp_mocks.cpp +++ b/source/source_lcao/module_operator_lcao/test/tmp_mocks.cpp @@ -209,3 +209,12 @@ void Numerical_Orbital::set_orbital_info(const int&, const int&, const int*, const int&) {} + +// mock of TD_info +class TD_info { +public: + TD_info() {} + ~TD_info() {} + static ModuleBase::Vector3 cart_At; +}; +ModuleBase::Vector3 TD_info::cart_At(0.0, 0.0, 0.0); \ No newline at end of file diff --git a/source/source_lcao/module_ri/RI_2D_Comm.h b/source/source_lcao/module_ri/RI_2D_Comm.h index 63cbbe8e68f..9b90c8a3ae9 100644 --- a/source/source_lcao/module_ri/RI_2D_Comm.h +++ b/source/source_lcao/module_ri/RI_2D_Comm.h @@ -51,6 +51,18 @@ namespace RI_2D_Comm const std::vector>>>& Hs, const Parallel_Orbitals& pv, TK* hk); + + + template + extern void add_Hexx_td( + const UnitCell& ucell, + const K_Vectors& kv, + const int ik, + const double alpha, + const std::vector>>>& Hs, + const Parallel_Orbitals& pv, + const ModuleBase::Vector3& At, + TK* hk); template extern void add_HexxR( diff --git a/source/source_lcao/module_ri/RI_2D_Comm.hpp b/source/source_lcao/module_ri/RI_2D_Comm.hpp index b9c34acd43d..dfd11621042 100644 --- a/source/source_lcao/module_ri/RI_2D_Comm.hpp +++ b/source/source_lcao/module_ri/RI_2D_Comm.hpp @@ -168,7 +168,70 @@ void RI_2D_Comm::add_Hexx( } ModuleBase::timer::tick("RI_2D_Comm", "add_Hexx"); } +template +void RI_2D_Comm::add_Hexx_td( + const UnitCell &ucell, + const K_Vectors& kv, + const int ik, + const double alpha, + const std::vector>>>& Hs, + const Parallel_Orbitals& pv, + const ModuleBase::Vector3& At, + TK* hk) +{ + ModuleBase::TITLE("RI_2D_Comm", "add_Hexx_td"); + ModuleBase::timer::tick("RI_2D_Comm", "add_Hexx_td"); + + const std::map> is_list = {{1, {0}}, {2, {kv.isk[ik]}}, {4, {0, 1, 2, 3}}}; + for (const int is_b: is_list.at(PARAM.inp.nspin)) + { + int is0_b, is1_b; + std::tie(is0_b, is1_b) = RI_2D_Comm::split_is_block(is_b); + for (const auto& Hs_tmpA: Hs[is_b]) + { + const TA& iat0 = Hs_tmpA.first; + for (const auto& Hs_tmpB: Hs_tmpA.second) + { + 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); + // cal tddft phase for hybrid gauge + 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)); + const RI::Tensor& H = Hs_tmpB.second; + for (size_t iw0_b = 0; iw0_b < H.shape[0]; ++iw0_b) + { + const int iwt0 = RI_2D_Comm::get_iwt(ucell, iat0, iw0_b, is0_b); + if (pv.global2local_row(iwt0) < 0) + { + continue; + } + for (size_t iw1_b = 0; iw1_b < H.shape[1]; ++iw1_b) + { + const int iwt1 = RI_2D_Comm::get_iwt(ucell, iat1, iw1_b, is1_b); + if (pv.global2local_col(iwt1) < 0) + { + continue; + } + LCAO_domain::set_mat2d(iwt0, + iwt1, + RI::Global_Func::convert(H(iw0_b, iw1_b)) + * RI::Global_Func::convert(frac), + pv, + hk); + } + } + } + } + } + ModuleBase::timer::tick("RI_2D_Comm", "add_Hexx_td"); +} std::tuple RI_2D_Comm::get_iat_iw_is_block(const UnitCell& ucell,const int& iwt) { diff --git a/source/source_lcao/module_rt/td_folding.cpp b/source/source_lcao/module_rt/td_folding.cpp index 13d358c53ba..e0b34cd31f4 100644 --- a/source/source_lcao/module_rt/td_folding.cpp +++ b/source/source_lcao/module_rt/td_folding.cpp @@ -2,13 +2,13 @@ #include "source_base/libm/libm.h" namespace module_rt{ template -void folding_HR_td(const hamilt::HContainer& hR, +void folding_HR_td(const UnitCell& ucell, + const hamilt::HContainer& hR, std::complex* hk, const ModuleBase::Vector3& kvec_d_in, + const ModuleBase::Vector3& cart_At, const int ncol, - const int hk_type, - const UnitCell* ucell, - const ModuleBase::Vector3& cart_At) + const int hk_type) { #ifdef _OPENMP #pragma omp parallel for @@ -24,8 +24,8 @@ void folding_HR_td(const hamilt::HContainer& hR, //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; + ModuleBase::Vector3 dtau = ucell.cal_dtau(iat1, iat2, r_index); + const double arg_td = cart_At * dtau * ucell.lat0; //new // cal k_phase @@ -41,20 +41,128 @@ void folding_HR_td(const hamilt::HContainer& hR, } } } + +template +void folding_partial_HR(const UnitCell& ucell, + const hamilt::HContainer& hR, + std::complex* hk, + const ModuleBase::Vector3& kvec_d_in, + const int ix, + const int ncol, + const int hk_type) +{ +#ifdef _OPENMP +#pragma omp parallel for +#endif + 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); + const ModuleBase::Vector3 dR(r_index.x, r_index.y, r_index.z); + const double arg = (kvec_d_in * dR) * ModuleBase::TWO_PI; + double sinp, cosp; + ModuleBase::libm::sincos(arg, &sinp, &cosp); + std::complex kphase = std::complex(cosp, sinp); + const ModuleBase::Vector3 dR_car = dR * ucell.latvec * ucell.lat0; + + tmp.find_R(r_index); + tmp.add_to_matrix(hk, ncol, kphase * ModuleBase::IMAG_UNIT * std::complex(dR_car[ix]), hk_type); + } + } +} + +template +void folding_partial_HR_td(const UnitCell& ucell, + const hamilt::HContainer& hR, + std::complex* hk, + const ModuleBase::Vector3& kvec_d_in, + const ModuleBase::Vector3& cart_At, + const int ix, + const int ncol, + const int hk_type) +{ +#ifdef _OPENMP +#pragma omp parallel for +#endif + 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); + + //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; + double sinp, cosp; + ModuleBase::libm::sincos(arg, &sinp, &cosp); + std::complex kphase = std::complex(cosp, sinp); + const ModuleBase::Vector3 dR_car = dR * ucell.latvec * ucell.lat0; + + tmp.find_R(r_index); + tmp.add_to_matrix(hk, ncol, kphase * ModuleBase::IMAG_UNIT * std::complex(dR_car[ix]), 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); +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); +template +void folding_partial_HR>(const UnitCell& ucell, + const hamilt::HContainer>& hR, + std::complex* hk, + const ModuleBase::Vector3& kvec_d_in, + const int ix, + const int ncol, + const int hk_type); +template +void folding_partial_HR(const UnitCell& ucell, + const hamilt::HContainer& hR, + std::complex* hk, + const ModuleBase::Vector3& kvec_d_in, + const int ix, + const int ncol, + const int hk_type); template -void folding_HR_td(const hamilt::HContainer& hR, +void folding_partial_HR_td>(const UnitCell& ucell, + const hamilt::HContainer>& hR, std::complex* hk, const ModuleBase::Vector3& kvec_d_in, + const ModuleBase::Vector3& cart_At, + const int ix, const int ncol, - const int hk_type, - const UnitCell* ucell, - const ModuleBase::Vector3& At); + const int hk_type); template -void folding_HR_td>(const hamilt::HContainer>& hR, +void folding_partial_HR_td(const UnitCell& ucell, + const hamilt::HContainer& hR, std::complex* hk, const ModuleBase::Vector3& kvec_d_in, + const ModuleBase::Vector3& cart_At, + const int ix, const int ncol, - const int hk_type, - const UnitCell* ucell, - const ModuleBase::Vector3& At); + const int hk_type); }// namespace module_rt \ No newline at end of file diff --git a/source/source_lcao/module_rt/td_folding.h b/source/source_lcao/module_rt/td_folding.h index 954d01d7765..dae13e21af5 100644 --- a/source/source_lcao/module_rt/td_folding.h +++ b/source/source_lcao/module_rt/td_folding.h @@ -6,13 +6,30 @@ namespace module_rt{ // folding HR to hk, for hybrid gauge template -void folding_HR_td(const hamilt::HContainer& hR, +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, - const UnitCell* ucell, - const ModuleBase::Vector3& At); + const int hk_type); +template +void folding_partial_HR(const UnitCell& ucell, + const hamilt::HContainer& hR, + std::complex* hk, + const ModuleBase::Vector3& kvec_d_in, + const int ix, + const int ncol, + const int hk_type); +template +void folding_partial_HR_td(const UnitCell& ucell, + const hamilt::HContainer& hR, + std::complex* hk, + const ModuleBase::Vector3& kvec_d_in, + const ModuleBase::Vector3& cart_At, + const int ix, + const int ncol, + const int hk_type); }// namespace module_rt #endif \ No newline at end of file diff --git a/source/source_lcao/module_rt/td_info.cpp b/source/source_lcao/module_rt/td_info.cpp index 259580d47d5..974c4116d3e 100644 --- a/source/source_lcao/module_rt/td_info.cpp +++ b/source/source_lcao/module_rt/td_info.cpp @@ -5,7 +5,7 @@ bool TD_info::out_mat_R = false; bool TD_info::out_vecpot = false; -bool TD_info::out_current = false; +int TD_info::out_current = 0; bool TD_info::out_current_k = false; bool TD_info::init_vecpot_file = false; bool TD_info::evolve_once = false; @@ -18,7 +18,7 @@ int TD_info::max_istep = -1; ModuleBase::Vector3 TD_info::cart_At; std::vector> TD_info::At_from_file; -TD_info::TD_info(const UnitCell* ucell_in) +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) @@ -39,6 +39,10 @@ TD_info::TD_info(const UnitCell* ucell_in) //std::cout<<"estep_shift"<istep += estep_shift; + if(out_current==2||elecstate::H_TDDFT_pw::stype == 2) + { + r_calculator.init(*ucell, pv, orb); + } return; } TD_info::~TD_info() diff --git a/source/source_lcao/module_rt/td_info.h b/source/source_lcao/module_rt/td_info.h index 6b2f3a839ee..57dccd7e9f6 100644 --- a/source/source_lcao/module_rt/td_info.h +++ b/source/source_lcao/module_rt/td_info.h @@ -3,17 +3,16 @@ #include "source_base/abfs-vector3_order.h" #include "source_base/timer.h" #include "source_lcao/module_hcontainer/hcontainer.h" +#include "source_io/cal_r_overlap_R.h" #include // Class to store TDDFT infos, mainly for periodic system. class TD_info { public: - TD_info(const UnitCell* ucell_in); + TD_info(const UnitCell* ucell_in,const Parallel_Orbitals& pv, const LCAO_Orbitals& orb); ~TD_info(); - void init(); - /// @brief switch to control the output of HR static bool out_mat_R; @@ -24,7 +23,7 @@ class TD_info static bool out_vecpot; /// @brief switch to control the output of current - static bool out_current; + static int out_current; /// @brief switch to control the format of the output current, in total or in each k-point static bool out_current_k; @@ -56,6 +55,15 @@ class TD_info { return this->current_term[i]; } + // set velocity HR. + void set_velocity_HR(hamilt::HContainer>* HR) + { + this->velocity_HR = HR; + } + hamilt::HContainer>* get_velocity_HR_pointer() const + { + return this->velocity_HR; + } int get_istep() { @@ -70,6 +78,9 @@ class TD_info // For TDDFT velocity gauge, to fix the output of HR std::map, std::map>>> HR_sparse_td_vel[2]; + //r_calculator + cal_r_overlap_R r_calculator; + private: /// @brief pointer to the unit cell const UnitCell* ucell = nullptr; @@ -94,6 +105,9 @@ class TD_info /// @brief part of Momentum operator, -iāˆ‡ - i[r,Vnl]. Used to calculate current. std::vector>*> current_term = {nullptr, nullptr, nullptr}; + + /// @brief store kinetic hamilton + hamilt::HContainer>* velocity_HR = nullptr; }; #endif