diff --git a/source/source_esolver/esolver_ks_lcao_tddft.cpp b/source/source_esolver/esolver_ks_lcao_tddft.cpp index 29da11e602..483d732962 100644 --- a/source/source_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/source_esolver/esolver_ks_lcao_tddft.cpp @@ -1,11 +1,11 @@ #include "esolver_ks_lcao_tddft.h" +#include "source_estate/elecstate_tools.h" #include "source_io/cal_r_overlap_R.h" #include "source_io/dipole_io.h" #include "source_io/td_current_io.h" #include "source_io/write_HS.h" #include "source_io/write_HS_R.h" -#include "source_estate/elecstate_tools.h" //--------------temporary---------------------------- #include "source_base/blas_connector.h" @@ -17,17 +17,17 @@ #include "source_estate/module_dm/cal_edm_tddft.h" #include "source_estate/module_dm/density_matrix.h" #include "source_estate/occupy.h" +#include "source_io/print_info.h" #include "source_lcao/module_rt/evolve_elec.h" #include "source_lcao/module_rt/td_velocity.h" #include "source_pw/module_pwdft/global.h" -#include "source_io/print_info.h" //-----HSolver ElecState Hamilt-------- +#include "module_parameter/parameter.h" #include "source_estate/cal_ux.h" #include "source_estate/elecstate_lcao.h" -#include "source_lcao/hamilt_lcaodft/hamilt_lcao.h" #include "source_hsolver/hsolver_lcao.h" -#include "module_parameter/parameter.h" +#include "source_lcao/hamilt_lcaodft/hamilt_lcao.h" #include "source_psi/psi.h" //-----force& stress------------------- @@ -87,30 +87,30 @@ void ESolver_KS_LCAO_TDDFT::before_all_runners(UnitCell& ucell, const In template void ESolver_KS_LCAO_TDDFT::hamilt2rho_single(UnitCell& ucell, - const int istep, - const int iter, - const double ethr) + const int istep, + const int iter, + const double ethr) { if (PARAM.inp.init_wfc == "file") { if (istep >= 1) { module_rt::Evolve_elec::solve_psi(istep, - PARAM.inp.nbands, - PARAM.globalv.nlocal, - kv.get_nks(), - this->p_hamilt, - this->pv, - this->psi, - this->psi_laststep, - this->Hk_laststep, - this->Sk_laststep, - this->pelec->ekb, - GlobalV::ofs_running, - td_htype, - PARAM.inp.propagator, - use_tensor, - use_lapack); + PARAM.inp.nbands, + PARAM.globalv.nlocal, + kv.get_nks(), + this->p_hamilt, + this->pv, + this->psi, + this->psi_laststep, + this->Hk_laststep, + this->Sk_laststep, + this->pelec->ekb, + GlobalV::ofs_running, + td_htype, + PARAM.inp.propagator, + use_tensor, + use_lapack); this->weight_dm_rho(); } this->weight_dm_rho(); @@ -118,21 +118,21 @@ void ESolver_KS_LCAO_TDDFT::hamilt2rho_single(UnitCell& ucell, else if (istep >= 2) { module_rt::Evolve_elec::solve_psi(istep, - PARAM.inp.nbands, - PARAM.globalv.nlocal, - kv.get_nks(), - this->p_hamilt, - this->pv, - this->psi, - this->psi_laststep, - this->Hk_laststep, - this->Sk_laststep, - this->pelec->ekb, - GlobalV::ofs_running, - td_htype, - PARAM.inp.propagator, - use_tensor, - use_lapack); + PARAM.inp.nbands, + PARAM.globalv.nlocal, + kv.get_nks(), + this->p_hamilt, + this->pv, + this->psi, + this->psi_laststep, + this->Hk_laststep, + this->Sk_laststep, + this->pelec->ekb, + GlobalV::ofs_running, + td_htype, + PARAM.inp.propagator, + use_tensor, + use_lapack); this->weight_dm_rho(); } else @@ -163,18 +163,14 @@ void ESolver_KS_LCAO_TDDFT::hamilt2rho_single(UnitCell& ucell, } template -void ESolver_KS_LCAO_TDDFT::iter_finish( - UnitCell& ucell, - const int istep, - int& iter, - bool& conv_esolver) +void ESolver_KS_LCAO_TDDFT::iter_finish(UnitCell& ucell, const int istep, int& iter, bool& conv_esolver) { // print occupation of each band if (iter == 1 && istep <= 2) { - GlobalV::ofs_running << " ---------------------------------------------------------" - << std::endl; - GlobalV::ofs_running << " occupations of electrons" << std::endl; + GlobalV::ofs_running << " ----------------------------------------------------------" << std::endl; + GlobalV::ofs_running << " Occupations of electrons" << std::endl; + GlobalV::ofs_running << " ----------------------------------------------------------" << std::endl; GlobalV::ofs_running << " k-point state occupation" << std::endl; GlobalV::ofs_running << std::setiosflags(std::ios::showpoint); GlobalV::ofs_running << std::left; @@ -183,23 +179,21 @@ void ESolver_KS_LCAO_TDDFT::iter_finish( { for (int ib = 0; ib < PARAM.inp.nbands; ib++) { - GlobalV::ofs_running << " " << std::setw(9) - << ik+1 << std::setw(8) << ib + 1 - << std::setw(12) << this->pelec->wg(ik, ib) << std::endl; + GlobalV::ofs_running << " " << std::setw(9) << ik + 1 << std::setw(8) << ib + 1 << std::setw(12) + << this->pelec->wg(ik, ib) << std::endl; } } - GlobalV::ofs_running << " ---------------------------------------------------------" - << std::endl; + GlobalV::ofs_running << " ----------------------------------------------------------" << std::endl; } ESolver_KS_LCAO, double>::iter_finish(ucell, istep, iter, conv_esolver); } template -void ESolver_KS_LCAO_TDDFT::update_pot(UnitCell& ucell, - const int istep, - const int iter, - const bool conv_esolver) +void ESolver_KS_LCAO_TDDFT::update_pot(UnitCell& ucell, + const int istep, + const int iter, + const bool conv_esolver) { // Calculate new potential according to new Charge Density if (!conv_esolver) @@ -234,7 +228,6 @@ void ESolver_KS_LCAO_TDDFT::update_pot(UnitCell& ucell, nrow_tmp = nlocal; #endif this->psi_laststep = new psi::Psi>(kv.get_nks(), ncol_tmp, nrow_tmp, kv.ngk, true); - } // allocate memory for Hk_laststep and Sk_laststep @@ -282,8 +275,8 @@ void ESolver_KS_LCAO_TDDFT::update_pot(UnitCell& ucell, if (td_htype == 1) { this->p_hamilt->updateHk(ik); - hamilt::MatrixBlock > h_mat; - hamilt::MatrixBlock > s_mat; + hamilt::MatrixBlock> h_mat; + hamilt::MatrixBlock> s_mat; this->p_hamilt->matrix(h_mat, s_mat); if (use_tensor && use_lapack) @@ -323,31 +316,31 @@ void ESolver_KS_LCAO_TDDFT::update_pot(UnitCell& ucell, } // print "eigen value" for tddft -// it seems uncessary to print out E_ii because the band energies are printed -/* - if (conv_esolver) - { - GlobalV::ofs_running << "----------------------------------------------------------" - << std::endl; - GlobalV::ofs_running << " Print E= " << std::endl; - GlobalV::ofs_running << " k-point state energy (eV)" << std::endl; - GlobalV::ofs_running << "----------------------------------------------------------" - << std::endl; - GlobalV::ofs_running << std::setprecision(6); - GlobalV::ofs_running << std::setiosflags(std::ios::showpoint); - - for (int ik = 0; ik < kv.get_nks(); ik++) + // it seems unnecessary to print out E_ii because the band energies are printed + /* + if (conv_esolver) { - for (int ib = 0; ib < PARAM.inp.nbands; ib++) + GlobalV::ofs_running << "----------------------------------------------------------" + << std::endl; + GlobalV::ofs_running << " Print E= " << std::endl; + GlobalV::ofs_running << " k-point state energy (eV)" << std::endl; + GlobalV::ofs_running << "----------------------------------------------------------" + << std::endl; + GlobalV::ofs_running << std::setprecision(6); + GlobalV::ofs_running << std::setiosflags(std::ios::showpoint); + + for (int ik = 0; ik < kv.get_nks(); ik++) { - GlobalV::ofs_running << " " << std::setw(7) << ik + 1 - << std::setw(7) << ib + 1 - << std::setw(10) << this->pelec->ekb(ik, ib) * ModuleBase::Ry_to_eV - << std::endl; + for (int ib = 0; ib < PARAM.inp.nbands; ib++) + { + GlobalV::ofs_running << " " << std::setw(7) << ik + 1 + << std::setw(7) << ib + 1 + << std::setw(10) << this->pelec->ekb(ik, ib) * ModuleBase::Ry_to_eV + << std::endl; + } } } - } -*/ + */ } template @@ -365,16 +358,11 @@ void ESolver_KS_LCAO_TDDFT::after_scf(UnitCell& ucell, const int istep, { std::stringstream ss_dipole; ss_dipole << PARAM.globalv.global_out_dir << "SPIN" << is + 1 << "_DIPOLE"; - ModuleIO::write_dipole(ucell, - this->chr.rho_save[is], - this->chr.rhopw, - is, - istep, - ss_dipole.str()); + ModuleIO::write_dipole(ucell, this->chr.rho_save[is], this->chr.rhopw, is, istep, ss_dipole.str()); } } - // (2) write current information + // (2) write current information if (TD_Velocity::out_current == true) { elecstate::DensityMatrix, double>* tmp_DM @@ -392,7 +380,6 @@ void ESolver_KS_LCAO_TDDFT::after_scf(UnitCell& ucell, const int istep, this->RA); } - ModuleBase::timer::tick("ESolver_LCAO_TDDFT", "after_scf"); } @@ -410,7 +397,7 @@ void ESolver_KS_LCAO_TDDFT::weight_dm_rho() } // calculate Eband energy - elecstate::calEBand(this->pelec->ekb,this->pelec->wg,this->pelec->f_en); + elecstate::calEBand(this->pelec->ekb, this->pelec->wg, this->pelec->f_en); // calculate the density matrix ModuleBase::GlobalFunc::NOTE("Calculate the density matrix."); diff --git a/source/source_estate/module_dm/cal_edm_tddft.cpp b/source/source_estate/module_dm/cal_edm_tddft.cpp index 87e1658865..f30c88b8fa 100644 --- a/source/source_estate/module_dm/cal_edm_tddft.cpp +++ b/source/source_estate/module_dm/cal_edm_tddft.cpp @@ -4,7 +4,6 @@ #include "source_base/scalapack_connector.h" namespace elecstate { - // use the original formula (Hamiltonian matrix) to calculate energy density matrix void cal_edm_tddft(Parallel_Orbitals& pv, elecstate::ElecState* pelec, @@ -254,5 +253,5 @@ void cal_edm_tddft(Parallel_Orbitals& pv, #endif } return; -} -} // namespace ModuleESolver +} // cal_edm_tddft +} // namespace elecstate diff --git a/source/source_lcao/module_rt/evolve_elec.cpp b/source/source_lcao/module_rt/evolve_elec.cpp index a4e787b1ee..2d857d3784 100644 --- a/source/source_lcao/module_rt/evolve_elec.cpp +++ b/source/source_lcao/module_rt/evolve_elec.cpp @@ -83,8 +83,8 @@ void Evolve_elec::solve_psi(const int& istep, propagator, ofs_running, print_matrix); - // std::cout << "Print ekb: " << std::endl; - // ekb.print(std::cout); + // GlobalV::ofs_running << "Print ekb: " << std::endl; + // ekb.print(GlobalV::ofs_running); } else { @@ -118,7 +118,7 @@ void Evolve_elec::solve_psi(const int& istep, #ifdef __MPI // Access the rank of the calling process in the communicator int myid = 0; - int root_proc = 0; + const int root_proc = 0; MPI_Comm_rank(MPI_COMM_WORLD, &myid); // Gather psi to the root process @@ -203,8 +203,16 @@ void Evolve_elec::solve_psi(const int& istep, len_HS_laststep); syncmem_double_d2h_op()(&(ekb(ik, 0)), ekb_tensor.data(), nband); - // std::cout << "Print ekb tensor: " << std::endl; - // ekb.print(std::cout); +#ifdef __MPI + const int root_proc = 0; + if (use_lapack) + { + // Synchronize ekb to all MPI processes + MPI_Bcast(&(ekb(ik, 0)), nband, MPI_DOUBLE, root_proc, MPI_COMM_WORLD); + } +#endif + // GlobalV::ofs_running << "Print ekb: " << std::endl; + // ekb.print(GlobalV::ofs_running); } } else diff --git a/tests/15_rtTDDFT_GPU/CASES_GPU.txt b/tests/15_rtTDDFT_GPU/CASES_GPU.txt index bbee16ea6a..c930799511 100644 --- a/tests/15_rtTDDFT_GPU/CASES_GPU.txt +++ b/tests/15_rtTDDFT_GPU/CASES_GPU.txt @@ -4,7 +4,7 @@ 04_NO_CO_ocp_TDDFT_GPU 05_NO_cur_TDDFT_GPU 06_NO_dir_TDDFT_GPU -# 07_NO_EDM_TDDFT_GPU +07_NO_EDM_TDDFT_GPU 09_NO_HEAV_TDDFT_GPU 10_NO_HHG_TDDFT_GPU 11_NO_O3_TDDFT_GPU