diff --git a/examples/lr-tddft/lcao_H2O/INPUT b/examples/lr-tddft/lcao_H2O/INPUT index 8e7164ae66..9a30e2de0c 100644 --- a/examples/lr-tddft/lcao_H2O/INPUT +++ b/examples/lr-tddft/lcao_H2O/INPUT @@ -6,6 +6,7 @@ orbital_dir ../../../tests/PP_ORB calculation scf nbands 23 symmetry -1 +nspin 2 #Parameters (2.Iteration) ecutwfc 60 ###Energy cutoff needs to be tested to ensure your calculation is reliable.[1] @@ -30,6 +31,7 @@ xc_kernel lda lr_solver dav lr_thr 1e-2 pw_diag_ndim 2 +# lr_unrestricted 1 ### use this to do TDUKS calculation for closeshell systems (openshell system will force TDUKS) esolver_type ks-lr out_alllog 1 @@ -39,6 +41,7 @@ out_alllog 1 nvirt 19 abs_wavelen_range 40 180 abs_broadening 0.01 +abs_gauge length ### [1] Energy cutoff determines the quality of numerical quadratures in your calculations. ### So it is strongly recommended to test whether your result (such as converged SCF energies) is diff --git a/examples/lr-tddft/lcao_Si2/INPUT b/examples/lr-tddft/lcao_Si2/INPUT index 11c9123140..62248494d4 100644 --- a/examples/lr-tddft/lcao_Si2/INPUT +++ b/examples/lr-tddft/lcao_Si2/INPUT @@ -5,7 +5,8 @@ pseudo_dir ../../../tests/PP_ORB orbital_dir ../../../tests/PP_ORB calculation scf nbands 23 -symmetry 0 +symmetry -1 +nspin 2 #Parameters (2.Iteration) ecutwfc 60 ###Energy cutoff needs to be tested to ensure your calculation is reliable.[1] @@ -37,6 +38,8 @@ out_alllog 1 nvirt 19 abs_wavelen_range 100 175 +abs_broadening 0.01 # in Ry +abs_gauge velocity ### velocity gauge is recommended for periodic systems ### [1] Energy cutoff determines the quality of numerical quadratures in your calculations. diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 7a0ff0aadf..e7dae3ebf1 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -731,6 +731,7 @@ OBJS_TENSOR=tensor.o\ xc_kernel.o\ pot_hxc_lrtd.o\ lr_spectrum.o\ + lr_spectrum_velocity.o\ hamilt_casida.o\ esolver_lrtd_lcao.o\ diff --git a/source/module_basis/module_nao/two_center_bundle.h b/source/module_basis/module_nao/two_center_bundle.h index bc4f3bcea5..0fcd888293 100644 --- a/source/module_basis/module_nao/two_center_bundle.h +++ b/source/module_basis/module_nao/two_center_bundle.h @@ -12,6 +12,7 @@ class TwoCenterBundle public: TwoCenterBundle() = default; ~TwoCenterBundle() = default; + TwoCenterBundle& operator=(TwoCenterBundle&&) = default; // NOTE: some variables might be set only on RANK-0 void build_orb(int ntype, const std::string* file_orb0); diff --git a/source/module_io/read_input_item_tddft.cpp b/source/module_io/read_input_item_tddft.cpp index bfd0542729..236b15545b 100644 --- a/source/module_io/read_input_item_tddft.cpp +++ b/source/module_io/read_input_item_tddft.cpp @@ -360,6 +360,12 @@ void ReadInput::item_lr_tddft() sync_doublevec(input.abs_wavelen_range, 2, 0.0); this->add_item(item); } + { + Input_Item item("abs_gauge"); + item.annotation = "whether to use length or velocity gauge to calculate the absorption spectrum in LR-TDDFT"; + read_sync_string(input.abs_gauge); + this->add_item(item); + } { Input_Item item("abs_broadening"); item.annotation = "the broadening (eta) for LR-TDDFT absorption spectrum"; diff --git a/source/module_io/test/read_input_ptest.cpp b/source/module_io/test/read_input_ptest.cpp index 11acb41f43..5e02408361 100644 --- a/source/module_io/test/read_input_ptest.cpp +++ b/source/module_io/test/read_input_ptest.cpp @@ -431,6 +431,7 @@ TEST_F(InputParaTest, ParaRead) EXPECT_EQ(param.inp.abs_wavelen_range.size(), 2); EXPECT_DOUBLE_EQ(param.inp.abs_wavelen_range[0], 0.0); EXPECT_DOUBLE_EQ(param.inp.abs_broadening, 0.01); + EXPECT_EQ(param.inp.abs_gauge, "length"); EXPECT_EQ(param.inp.rdmft, 0); EXPECT_DOUBLE_EQ(param.inp.rdmft_power_alpha, 0.656); } diff --git a/source/module_lr/CMakeLists.txt b/source/module_lr/CMakeLists.txt index a672e05c92..3459ca8d3a 100644 --- a/source/module_lr/CMakeLists.txt +++ b/source/module_lr/CMakeLists.txt @@ -16,6 +16,7 @@ if(ENABLE_LCAO) operator_casida/operator_lr_exx.cpp potentials/pot_hxc_lrtd.cpp lr_spectrum.cpp + lr_spectrum_velocity.cpp hamilt_casida.cpp esolver_lrtd_lcao.cpp potentials/xc_kernel.cpp) diff --git a/source/module_lr/esolver_lrtd_lcao.cpp b/source/module_lr/esolver_lrtd_lcao.cpp index 38fb300594..8c00ef2fdd 100644 --- a/source/module_lr/esolver_lrtd_lcao.cpp +++ b/source/module_lr/esolver_lrtd_lcao.cpp @@ -62,6 +62,21 @@ inline int cal_nupdown_form_occ(const ModuleBase::matrix& wg) return nupdown; } +inline void setup_2center_table(TwoCenterBundle& two_center_bundle, LCAO_Orbitals& orb, UnitCell& ucell) +{ + // set up 2-center table +#ifdef USE_NEW_TWO_CENTER + two_center_bundle.tabulate(); +#else + two_center_bundle.tabulate(inp.lcao_ecut, inp.lcao_dk, inp.lcao_dr, inp.lcao_rmax); +#endif + if (PARAM.inp.vnl_in_h) + { + ucell.infoNL.setupNonlocal(ucell.ntype, ucell.atoms, GlobalV::ofs_running, orb); + two_center_bundle.build_beta(ucell.ntype, ucell.infoNL.Beta); + } +} + template void LR::ESolver_LR::parameter_check()const { @@ -149,8 +164,7 @@ LR::ESolver_LR::ESolver_LR(ModuleESolver::ESolver_KS_LCAO&& ks_sol this->gd = std::move(ks_sol.gd); // xc kernel - this->xc_kernel = inp.xc_kernel; - std::transform(xc_kernel.begin(), xc_kernel.end(), xc_kernel.begin(), tolower); + this->xc_kernel = LR_Util::tolower(inp.xc_kernel); //kv this->kv = std::move(ks_sol.kv); @@ -218,8 +232,7 @@ LR::ESolver_LR::ESolver_LR(ModuleESolver::ESolver_KS_LCAO&& ks_sol if (xc_kernel == "hf" || xc_kernel == "hse") { // if the same kernel is calculated in the esolver_ks, move it - std::string dft_functional = input.dft_functional; - std::transform(dft_functional.begin(), dft_functional.end(), dft_functional.begin(), tolower); + std::string dft_functional = LR_Util::tolower(input.dft_functional); if (ks_sol.exx_lri_double && std::is_same::value && xc_kernel == dft_functional) { this->move_exx_lri(ks_sol.exx_lri_double); } else if (ks_sol.exx_lri_complex && std::is_same>::value && xc_kernel == dft_functional) { @@ -237,6 +250,10 @@ LR::ESolver_LR::ESolver_LR(ModuleESolver::ESolver_KS_LCAO&& ks_sol #endif this->pelec = new elecstate::ElecStateLCAO(); orb_cutoff_ = ks_sol.orb_.cutoffs(); + if (LR_Util::tolower(input.abs_gauge) == "velocity") + { + this->two_center_bundle_ = std::move(ks_sol.two_center_bundle_); + } } template @@ -247,8 +264,7 @@ LR::ESolver_LR::ESolver_LR(const Input_para& inp, UnitCell& ucell) : inpu { ModuleBase::TITLE("ESolver_LR", "ESolver_LR(from scratch)"); // xc kernel - this->xc_kernel = inp.xc_kernel; - std::transform(xc_kernel.begin(), xc_kernel.end(), xc_kernel.begin(), tolower); + this->xc_kernel = LR_Util::tolower(inp.xc_kernel); // necessary steps in ESolver_FP ESolver_FP::before_all_runners(ucell, inp); @@ -272,6 +288,10 @@ LR::ESolver_LR::ESolver_LR(const Input_para& inp, UnitCell& ucell) : inpu LCAO_Orbitals orb; two_center_bundle_.to_LCAO_Orbitals(orb, inp.lcao_ecut, inp.lcao_dk, inp.lcao_dr, inp.lcao_rmax); orb_cutoff_ = orb.cutoffs(); + if (LR_Util::tolower(input.abs_gauge) == "velocity") + { + setup_2center_table(this->two_center_bundle_, orb, ucell); + } this->set_dimension(); // setup 2d-block distribution for AO-matrix and KS wfc @@ -320,7 +340,6 @@ LR::ESolver_LR::ESolver_LR(const Input_para& inp, UnitCell& ucell) : inpu this->init_pot(chg_gs); // search adjacent atoms and init Gint - std::cout << "ucell.infoNL.get_rcutmax_Beta(): " << ucell.infoNL.get_rcutmax_Beta() << std::endl; double search_radius = -1.0; search_radius = atom_arrange::set_sr_NL(GlobalV::ofs_running, PARAM.inp.out_level, @@ -539,26 +558,21 @@ void LR::ESolver_LR::after_all_runners(UnitCell& ucell) auto spin_types = (nspin == 2 && !openshell) ? std::vector({ "singlet", "triplet" }) : std::vector({ "updown" }); for (int is = 0;is < this->X.size();++is) { - LR_Spectrum spectrum(nspin, - this->nbasis, - this->nocc, - this->nvirt, - this->gint_, - *this->pw_rho, - *this->psi_ks, - this->ucell, - this->kv, - this->gd, - this->orb_cutoff_, - this->paraX_, - this->paraC_, - this->paraMat_, - &this->pelec->ekb.c[is * nstates], - this->X[is].template data(), - nstates, - openshell); + LR_Spectrum spectrum(nspin, this->nbasis, this->nocc, this->nvirt, this->gint_, *this->pw_rho, *this->psi_ks, + this->ucell, this->kv, this->gd, this->orb_cutoff_, this->two_center_bundle_, + this->paraX_, this->paraC_, this->paraMat_, + &this->pelec->ekb.c[is * nstates], this->X[is].template data(), nstates, openshell, + LR_Util::tolower(input.abs_gauge)); spectrum.transition_analysis(spin_types[is]); - spectrum.optical_absorption(freq, input.abs_broadening, spin_types[is]); + if (spin_types[is] != "triplet") // triplets has no transition dipole and no contribution to the spectrum + { + spectrum.optical_absorption_method1(freq, input.abs_broadening); + // =============================================== for test ==================================================== + // spectrum.optical_absorption_method2(freq, input.abs_broadening, spin_types[is]); + // spectrum.test_transition_dipoles_velocity_ks(eig_ks.c); + // spectrum.write_transition_dipole(PARAM.globalv.global_out_dir + "dipole_velocity_ks.dat"); + // =============================================== for test ==================================================== + } } } diff --git a/source/module_lr/lr_spectrum.cpp b/source/module_lr/lr_spectrum.cpp index 3a919e9dc7..c6db92c67b 100644 --- a/source/module_lr/lr_spectrum.cpp +++ b/source/module_lr/lr_spectrum.cpp @@ -6,6 +6,34 @@ #include "module_lr/utils/lr_util.h" #include "module_lr/utils/lr_util_hcontainer.h" #include "module_lr/utils/lr_util_print.h" + +template +elecstate::DensityMatrix LR::LR_Spectrum::cal_transition_density_matrix(const int istate, const T* X_in, const bool need_R) +{ + const T* const X = X_in == nullptr ? this->X : X_in; + const int offset_b = istate * ldim; //start index of band istate + elecstate::DensityMatrix DM_trans(&this->pmat, this->nspin_x, this->kv.kvec_d, this->nk); + for (int is = 0;is < this->nspin_x; ++is) + { + const int offset_x = offset_b + is * nk * this->pX[0].get_local_size(); + //1. transition density +#ifdef __MPI + std::vector dm_trans_2d = cal_dm_trans_pblas(X + offset_x, this->pX[is], psi_ks[is], this->pc, this->naos, this->nocc[is], this->nvirt[is], this->pmat); + // if (this->tdm_sym) for (auto& t : dm_trans_2d) LR_Util::matsym(t.data(), naos, pmat); +#else + std::vector dm_trans_2d = cal_dm_trans_blas(X + offset_x, this->psi_ks[is], this->nocc[is], this->nvirt[is]); + // if (this->tdm_sym) for (auto& t : dm_trans_2d) LR_Util::matsym(t.data(), naos); +#endif + for (int ik = 0;ik < this->nk;++ik) { DM_trans.set_DMK_pointer(ik + is * nk, dm_trans_2d[ik].data()); } + } + if (need_R) + { + LR_Util::initialize_DMR(DM_trans, this->pmat, this->ucell, this->gd_, this->orb_cutoff_); + DM_trans.cal_DMR(); + } + return DM_trans; +} + template void LR::LR_Spectrum::cal_gint_rho(double** rho, const int& nrxx) { @@ -23,154 +51,158 @@ inline void check_sum_rule(const double& osc_tot) } } -template <> -void LR::LR_Spectrum::oscillator_strength(const Grid_Driver& gd, const std::vector& orb_cutoff) +template<> +ModuleBase::Vector3 LR::LR_Spectrum::cal_transition_dipole_istate_length(const int istate) { - ModuleBase::TITLE("LR::LR_Spectrum", "oscillator_strength"); - std::vector& osc = this->oscillator_strength_; // unit: Ry - osc.resize(nstate, 0.0); - double osc_tot = 0.0; - elecstate::DensityMatrix DM_trans(&this->pmat, 1, this->kv.kvec_d, this->nk); - LR_Util::initialize_DMR(DM_trans, this->pmat, this->ucell, gd, orb_cutoff); - this->transition_dipole_.resize(nstate, ModuleBase::Vector3(0.0, 0.0, 0.0)); - for (int istate = 0;istate < nstate;++istate) + ModuleBase::Vector3 trans_dipole(0.0, 0.0, 0.0); + // 1. transition density matrix + const elecstate::DensityMatrix& DM_trans = this->cal_transition_density_matrix(istate); + for (int is = 0;is < this->nspin_x;++is) { - const int offset_b = istate * ldim; //start index of band istate - for (int is = 0;is < this->nspin_x;++is) + this->gint->transfer_DM2DtoGrid({ DM_trans.get_DMR_vector().at(is) }); + + // 2. transition density + double** rho_trans; + LR_Util::_allocate_2order_nested_ptr(rho_trans, 1, this->rho_basis.nrxx); + this->cal_gint_rho(rho_trans, this->rho_basis.nrxx); + + // 3. transition dipole moment + for (int ir = 0; ir < rho_basis.nrxx; ++ir) { - const int offset_x = offset_b + is * nk * pX[0].get_local_size(); - //1. transition density -#ifdef __MPI - std::vector dm_trans_2d = cal_dm_trans_pblas(X + offset_x, this->pX[is], psi_ks[is], this->pc, this->naos, this->nocc[is], this->nvirt[is], this->pmat); - // if (this->tdm_sym) for (auto& t : dm_trans_2d) LR_Util::matsym(t.data(), naos, pmat); -#else - std::vector dm_trans_2d = cal_dm_trans_blas(X + offset_x, this->psi_ks[is], this->nocc[is], this->nvirt[is]); - // if (this->tdm_sym) for (auto& t : dm_trans_2d) LR_Util::matsym(t.data(), naos); -#endif - for (int ik = 0;ik < this->nk;++ik) { DM_trans.set_DMK_pointer(ik, dm_trans_2d[ik].data()); } - DM_trans.cal_DMR(); - this->gint->transfer_DM2DtoGrid(DM_trans.get_DMR_vector()); + int i = ir / (rho_basis.ny * rho_basis.nplane); + int j = ir / rho_basis.nplane - i * rho_basis.ny; + int k = ir % rho_basis.nplane + rho_basis.startz_current; + ModuleBase::Vector3 rd(static_cast(i) / rho_basis.nx, static_cast(j) / rho_basis.ny, static_cast(k) / rho_basis.nz); //+1/2 better? + rd -= ModuleBase::Vector3(0.5, 0.5, 0.5); //shift to the center of the grid (need ?) + ModuleBase::Vector3 rc = rd * ucell.latvec * ucell.lat0; // real coordinate + trans_dipole += rc * rho_trans[0][ir]; + } + LR_Util::_deallocate_2order_nested_ptr(rho_trans, 1); + } + trans_dipole *= (ucell.omega / static_cast(gint->get_ncxyz())); // dv + trans_dipole *= static_cast(this->nk); // nk is divided inside DM_trans, now recover it + if (this->nspin_x == 1) { trans_dipole *= sqrt(2.0); } // *2 for 2 spins, /sqrt(2) for the halfed dimension of X in the normalizaiton + Parallel_Reduce::reduce_all(trans_dipole.x); + Parallel_Reduce::reduce_all(trans_dipole.y); + Parallel_Reduce::reduce_all(trans_dipole.z); + return trans_dipole; +} - // 2. transition density - double** rho_trans; - LR_Util::_allocate_2order_nested_ptr(rho_trans, 1, this->rho_basis.nrxx); - this->cal_gint_rho(rho_trans, this->rho_basis.nrxx); +template<> +ModuleBase::Vector3> LR::LR_Spectrum>::cal_transition_dipole_istate_length(const int istate) +{ - // 3. transition dipole moment - for (int ir = 0; ir < rho_basis.nrxx; ++ir) - { - int i = ir / (rho_basis.ny * rho_basis.nplane); - int j = ir / rho_basis.nplane - i * rho_basis.ny; - int k = ir % rho_basis.nplane + rho_basis.startz_current; - ModuleBase::Vector3 rd(static_cast(i) / rho_basis.nx, static_cast(j) / rho_basis.ny, static_cast(k) / rho_basis.nz); //+1/2 better? - rd -= ModuleBase::Vector3(0.5, 0.5, 0.5); //shift to the center of the grid (need ?) - ModuleBase::Vector3 rc = rd * ucell.latvec * ucell.lat0; // real coordinate - transition_dipole_[istate] += rc * rho_trans[0][ir]; - } - LR_Util::_deallocate_2order_nested_ptr(rho_trans, 1); + //1. transition density matrix + ModuleBase::Vector3> trans_dipole(0.0, 0.0, 0.0); + const elecstate::DensityMatrix, std::complex>& DM_trans = this->cal_transition_density_matrix(istate); + for (int is = 0;is < this->nspin_x;++is) + { + // 2. transition density + double** rho_trans_real; + double** rho_trans_imag; + LR_Util::_allocate_2order_nested_ptr(rho_trans_real, 1, this->rho_basis.nrxx); + LR_Util::_allocate_2order_nested_ptr(rho_trans_imag, 1, this->rho_basis.nrxx); + + elecstate::DensityMatrix, double> DM_trans_real_imag(&this->pmat, 1, this->kv.kvec_d, this->nk); + LR_Util::initialize_DMR(DM_trans_real_imag, this->pmat, this->ucell, this->gd_, this->orb_cutoff_); + + // real part + LR_Util::get_DMR_real_imag_part(DM_trans, DM_trans_real_imag, ucell.nat, 'R'); + this->gint->transfer_DM2DtoGrid(DM_trans_real_imag.get_DMR_vector()); + this->cal_gint_rho(rho_trans_real, this->rho_basis.nrxx); + // LR_Util::print_grid_nonzero(rho_trans_real[0], this->rho_basis.nrxx, 10, "rho_trans"); + + // imag part + LR_Util::get_DMR_real_imag_part(DM_trans, DM_trans_real_imag, ucell.nat, 'I'); + this->gint->transfer_DM2DtoGrid(DM_trans_real_imag.get_DMR_vector()); + this->cal_gint_rho(rho_trans_imag, this->rho_basis.nrxx); + // LR_Util::print_grid_nonzero(rho_trans_imag[0], this->rho_basis.nrxx, 10, "rho_trans"); + + // 3. transition dipole moment + for (int ir = 0; ir < rho_basis.nrxx; ++ir) + { + int i = ir / (rho_basis.ny * rho_basis.nplane); + int j = ir / rho_basis.nplane - i * rho_basis.ny; + int k = ir % rho_basis.nplane + rho_basis.startz_current; + ModuleBase::Vector3 rd(static_cast(i) / rho_basis.nx, static_cast(j) / rho_basis.ny, static_cast(k) / rho_basis.nz); //+1/2 better? + rd -= ModuleBase::Vector3(0.5, 0.5, 0.5); //shift to the center of the grid (need ?) + ModuleBase::Vector3 rc = rd * ucell.latvec * ucell.lat0; // real coordinate + ModuleBase::Vector3> rc_complex(rc.x, rc.y, rc.z); + trans_dipole += rc_complex * std::complex(rho_trans_real[0][ir], rho_trans_imag[0][ir]); } - transition_dipole_[istate] *= (ucell.omega / static_cast(gint->get_ncxyz())); // dv - Parallel_Reduce::reduce_all(transition_dipole_[istate].x); - Parallel_Reduce::reduce_all(transition_dipole_[istate].y); - Parallel_Reduce::reduce_all(transition_dipole_[istate].z); - osc[istate] = transition_dipole_[istate].norm2() * eig[istate] * 2. / 3.; - osc_tot += osc[istate] / 2.; //Ry to Hartree (1/2) + LR_Util::_deallocate_2order_nested_ptr(rho_trans_real, 1); + LR_Util::_deallocate_2order_nested_ptr(rho_trans_imag, 1); } - check_sum_rule(osc_tot); + trans_dipole *= (ucell.omega / static_cast(gint->get_ncxyz())); // dv + trans_dipole *= static_cast(this->nk); // nk is divided inside DM_trans, now recover it + if (this->nspin_x == 1) { trans_dipole *= sqrt(2.0); } // *2 for 2 spins, /sqrt(2) for the halfed dimension of X in the normalizaiton + Parallel_Reduce::reduce_all(trans_dipole.x); + Parallel_Reduce::reduce_all(trans_dipole.y); + Parallel_Reduce::reduce_all(trans_dipole.z); + return trans_dipole; } -template <> -void LR::LR_Spectrum>::oscillator_strength(const Grid_Driver& gd, - const std::vector& orb_cutoff) +template<> double LR::LR_Spectrum::cal_mean_squared_dipole(ModuleBase::Vector3 dipole) +{ + return dipole.norm2() / 3.; +} +template<> double LR::LR_Spectrum>::cal_mean_squared_dipole(ModuleBase::Vector3> dipole) +{ + return dipole.norm2().real() / 3.; +} + +template +void LR::LR_Spectrum::cal_transition_dipoles_length() +{ + transition_dipole_.resize(nstate); + this->mean_squared_transition_dipole_.resize(nstate); + for (int istate = 0;istate < nstate;++istate) + { + transition_dipole_[istate] = cal_transition_dipole_istate_length(istate); + mean_squared_transition_dipole_[istate] = cal_mean_squared_dipole(transition_dipole_[istate]); + } +} + +template +void LR::LR_Spectrum::oscillator_strength() { ModuleBase::TITLE("LR::LR_Spectrum", "oscillator_strength"); std::vector& osc = this->oscillator_strength_; // unit: Ry osc.resize(nstate, 0.0); double osc_tot = 0.0; - elecstate::DensityMatrix, std::complex> DM_trans(&this->pmat, 1, this->kv.kvec_d, this->nk); - LR_Util::initialize_DMR(DM_trans, this->pmat, this->ucell, gd, orb_cutoff); - elecstate::DensityMatrix, double> DM_trans_real_imag(&this->pmat, 1, this->kv.kvec_d, this->nk); - LR_Util::initialize_DMR(DM_trans_real_imag, this->pmat, this->ucell, gd, orb_cutoff); - - this->transition_dipole_.resize(nstate, ModuleBase::Vector3>(0.0, 0.0, 0.0)); for (int istate = 0;istate < nstate;++istate) { - const int offset_b = istate * ldim; //start index of band istate - for (int is = 0;is < this->nspin_x;++is) - { - const int offset_x = offset_b + is * nk * pX[0].get_local_size(); - //1. transition density -#ifdef __MPI - std::vector dm_trans_2d = cal_dm_trans_pblas(X + offset_x, this->pX[is], psi_ks[is], this->pc, this->naos, this->nocc[is], this->nvirt[is], this->pmat, /*renorm_k=*/false, 1); - // if (this->tdm_sym) for (auto& t : dm_trans_2d) LR_Util::matsym(t.data(), naos, pmat); -#else - std::vector dm_trans_2d = cal_dm_trans_blas(X + offset_x, psi_ks[is], this->nocc[is], this->nvirt[is],/*renorm_k=*/false, 1); - // if (this->tdm_sym) for (auto& t : dm_trans_2d) LR_Util::matsym(t.data(), naos); -#endif - for (int ik = 0;ik < this->nk;++ik) { DM_trans.set_DMK_pointer(ik, dm_trans_2d[ik].data>()); } - // for (int ik = 0;ik < this->nk;++ik) - // LR_Util::print_tensor>(dm_trans_2d[ik], "1.DMK[ik=" + std::to_string(ik) + "]", dynamic_cast(&this->pmat)); - DM_trans.cal_DMR(); - - // 2. transition density - double** rho_trans_real; - double** rho_trans_imag; - LR_Util::_allocate_2order_nested_ptr(rho_trans_real, 1, this->rho_basis.nrxx); - LR_Util::_allocate_2order_nested_ptr(rho_trans_imag, 1, this->rho_basis.nrxx); - // real part - LR_Util::get_DMR_real_imag_part(DM_trans, DM_trans_real_imag, ucell.nat, 'R'); - this->gint->transfer_DM2DtoGrid(DM_trans_real_imag.get_DMR_vector()); - this->cal_gint_rho(rho_trans_real, this->rho_basis.nrxx); - // LR_Util::print_grid_nonzero(rho_trans_real[0], this->rho_basis.nrxx, 10, "rho_trans"); - - // imag part - LR_Util::get_DMR_real_imag_part(DM_trans, DM_trans_real_imag, ucell.nat, 'I'); - this->gint->transfer_DM2DtoGrid(DM_trans_real_imag.get_DMR_vector()); - this->cal_gint_rho(rho_trans_imag, this->rho_basis.nrxx); - // LR_Util::print_grid_nonzero(rho_trans_imag[0], this->rho_basis.nrxx, 10, "rho_trans"); - - // 3. transition dipole moment - for (int ir = 0; ir < rho_basis.nrxx; ++ir) - { - int i = ir / (rho_basis.ny * rho_basis.nplane); - int j = ir / rho_basis.nplane - i * rho_basis.ny; - int k = ir % rho_basis.nplane + rho_basis.startz_current; - ModuleBase::Vector3 rd(static_cast(i) / rho_basis.nx, static_cast(j) / rho_basis.ny, static_cast(k) / rho_basis.nz); //+1/2 better? - rd -= ModuleBase::Vector3(0.5, 0.5, 0.5); //shift to the center of the grid (need ?) - ModuleBase::Vector3 rc = rd * ucell.latvec * ucell.lat0; // real coordinate - ModuleBase::Vector3> rc_complex(rc.x, rc.y, rc.z); - transition_dipole_[istate] += rc_complex * std::complex(rho_trans_real[0][ir], rho_trans_imag[0][ir]); - } - LR_Util::_deallocate_2order_nested_ptr(rho_trans_real, 1); - LR_Util::_deallocate_2order_nested_ptr(rho_trans_imag, 1); - } - transition_dipole_[istate] *= (ucell.omega / static_cast(gint->get_ncxyz())); // dv - Parallel_Reduce::reduce_all(transition_dipole_[istate].x); - Parallel_Reduce::reduce_all(transition_dipole_[istate].y); - Parallel_Reduce::reduce_all(transition_dipole_[istate].z); - auto norm2 = [](const ModuleBase::Vector3>& v) -> double - { - return v.x.real() * v.x.real() + v.y.real() * v.y.real() + v.z.real() * v.z.real() - + v.x.imag() * v.x.imag() + v.y.imag() * v.y.imag() + v.z.imag() * v.z.imag(); - }; - osc[istate] = norm2(transition_dipole_[istate]) * eig[istate] * 2. / 3.; - osc_tot += osc[istate] / 2.; // Ry to Hartree (1/2) + osc[istate] = this->mean_squared_transition_dipole_[istate] * this->eig[istate] * 2.; + osc_tot += osc[istate] / 2.; //Ry to Hartree (1/2) } check_sum_rule(osc_tot); } + template -void LR::LR_Spectrum::optical_absorption(const std::vector& freq, const double eta, const std::string& spintype) +void LR::LR_Spectrum::optical_absorption_method1(const std::vector& freq, const double eta) { + // ============test dipole================ + // this->cal_transition_dipoles_length(); + // this->write_transition_dipole(PARAM.globalv.global_out_dir + "dipole_length.dat"); + // this->cal_transition_dipoles_velocity(); + // this->write_transition_dipole(PARAM.globalv.global_out_dir + "dipole_velocity.dat"); + // exit(0); + // ============test dipole================ ModuleBase::TITLE("LR::LR_Spectrum", "optical_absorption"); + // 4*pi^2/V * mean_squared_dipole *delta(w-Omega_S) + // = -8*pi*Omega_S/V * mean_squared_dipole * Im[1/[(w+i\eta)^2-\Omega_S^2]] + // = -4*pi/V * oscilator_strength * Im[1/[(w+i\eta)^2-\Omega_S^2]] std::vector& osc = this->oscillator_strength_; - std::ofstream ofs(PARAM.globalv.global_out_dir + "absorption_" + spintype + ".dat"); + std::ofstream ofs(PARAM.globalv.global_out_dir + "absorption.dat"); if (GlobalV::MY_RANK == 0) { ofs << "Frequency (eV) | wave length(nm) | Absorption (a.u.)" << std::endl; } double FourPI_div_c = ModuleBase::FOUR_PI / 137.036; + double fac = 4 * M_PI / ucell.omega * ModuleBase::e2 / this->nk; // e2 for Ry to Hartree in the denominator for (int f = 0;f < freq.size();++f) { std::complex f_complex = std::complex(freq[f], eta); double abs = 0.0; - for (int i = 0;i < osc.size();++i) { abs += (osc[i] / (f_complex * f_complex - eig[i] * eig[i])).imag() * freq[f] * FourPI_div_c; } + // for (int i = 0;i < osc.size();++i) { abs += (osc[i] / (f_complex * f_complex - eig[i] * eig[i])).imag() * freq[f] * FourPI_div_c; } + for (int i = 0;i < osc.size();++i) { abs += (osc[i] / (f_complex * f_complex - eig[i] * eig[i])).imag() * fac; } if (GlobalV::MY_RANK == 0) { ofs << freq[f] * ModuleBase::Ry_to_eV << "\t" << 91.126664 / freq[f] << "\t" << std::abs(abs) << std::endl; } } ofs.close(); @@ -247,5 +279,21 @@ std::map LR::LR_Spectrum::get_pair_info(const int i) return { {"ispin", ispin}, {"ik", ik}, {"iocc", iocc}, {"ivirt", ivirt} }; } +template +void LR::LR_Spectrum::write_transition_dipole(const std::string& filename) +{ + std::ofstream ofs(filename); + ofs << "Transition dipole moment (a.u.)" << std::endl; + ofs << std::setw(20) << "State" << std::setw(20) << "x" << std::setw(20) << "y" << std::setw(20) << "z" << std::setw(20) << "average" << std::endl; + for (int istate = 0;istate < nstate;++istate) + { + ofs << std::setw(20) << istate << std::setw(20) << transition_dipole_[istate].x << std::setw(20) + << transition_dipole_[istate].y << std::setw(20) + << transition_dipole_[istate].z << std::setw(20) + << mean_squared_transition_dipole_[istate] << std::endl; + } + ofs.close(); +} + template class LR::LR_Spectrum; template class LR::LR_Spectrum>; \ No newline at end of file diff --git a/source/module_lr/lr_spectrum.h b/source/module_lr/lr_spectrum.h index 80e1326343..1ef65320a6 100644 --- a/source/module_lr/lr_spectrum.h +++ b/source/module_lr/lr_spectrum.h @@ -1,79 +1,92 @@ +#pragma once #include "module_cell/klist.h" #include "module_lr/utils/gint_template.h" #include "module_psi/psi.h" #include "module_elecstate/module_dm/density_matrix.h" #include "module_lr/utils/lr_util.h" - +#include "module_basis/module_nao/two_center_bundle.h" +#include "module_hamilt_lcao/module_tddft/td_current.h" namespace LR { template class LR_Spectrum { public: - LR_Spectrum(const int& nspin_global, - const int& naos, - const std::vector& nocc, - const std::vector& nvirt, - typename TGint::type* gint, - const ModulePW::PW_Basis& rho_basis, - psi::Psi& psi_ks_in, - const UnitCell& ucell, - const K_Vectors& kv_in, - const Grid_Driver& gd, - const std::vector& orb_cutoff, - const std::vector& pX_in, - const Parallel_2D& pc_in, - const Parallel_Orbitals& pmat_in, - const double* eig, - const T* X, - const int& nstate, - const bool& openshell) - : nspin_x(openshell ? 2 : 1), naos(naos), nocc(nocc), nvirt(nvirt), nk(kv_in.get_nks() / nspin_global), - gint(gint), rho_basis(rho_basis), ucell(ucell), kv(kv_in), pX(pX_in), pc(pc_in), pmat(pmat_in), eig(eig), - X(X), nstate(nstate), - ldim(nk - * (nspin_x == 2 ? pX_in[0].get_local_size() + pX_in[1].get_local_size() : pX_in[0].get_local_size())), - gdim(nk * std::inner_product(nocc.begin(), nocc.end(), nvirt.begin(), 0)) - { - for (int is = 0; is < nspin_global; ++is) - { - psi_ks.emplace_back(LR_Util::get_psi_spin(psi_ks_in, is, nk)); - } - this->oscillator_strength(gd, orb_cutoff); + LR_Spectrum(const int& nspin_global, const int& naos, const std::vector& nocc, const std::vector& nvirt, + typename TGint::type* gint, const ModulePW::PW_Basis& rho_basis, psi::Psi& psi_ks_in, + const UnitCell& ucell, const K_Vectors& kv_in, const Grid_Driver& gd, const std::vector& orb_cutoff, + const TwoCenterBundle& two_center_bundle_, + const std::vector& pX_in, const Parallel_2D& pc_in, const Parallel_Orbitals& pmat_in, + const double* eig, const T* X, const int& nstate, const bool& openshell, + const std::string& gauge = "length") : + nspin_x(openshell ? 2 : 1), naos(naos), nocc(nocc), nvirt(nvirt), nk(kv_in.get_nks() / nspin_global), + gint(gint), rho_basis(rho_basis), ucell(ucell), kv(kv_in), gd_(gd), + orb_cutoff_(orb_cutoff), two_center_bundle_(two_center_bundle_), + pX(pX_in), pc(pc_in), pmat(pmat_in), + eig(eig), X(X), nstate(nstate), + ldim(nk* (nspin_x == 2 ? pX_in[0].get_local_size() + pX_in[1].get_local_size() : pX_in[0].get_local_size())), + gdim(nk* std::inner_product(nocc.begin(), nocc.end(), nvirt.begin(), 0)) + { + for (int is = 0;is < nspin_global;++is) { psi_ks.emplace_back(LR_Util::get_psi_spin(psi_ks_in, is, nk)); } + gauge == "velocity" ? this->cal_transition_dipoles_velocity() : this->cal_transition_dipoles_length(); + this->oscillator_strength(); }; - /// @brief calculate the optical absorption spectrum - void optical_absorption(const std::vector& freq, const double eta, const std::string& spintype); + /// @brief calculate the optical absorption spectrum with $Im[1/[(w+i\eta)^2-\Omega_S^2]]$ + void optical_absorption_method1(const std::vector& freq, const double eta); + /// @brief calculate the optical absorption spectrum with lorentzian delta function + void optical_absorption_method2(const std::vector& freq, const double eta); /// @brief print out the transition dipole moment and the main contributions to the transition amplitude void transition_analysis(const std::string& spintype); + + //========================================== test functions ============================================== + /// @brief write transition dipole + void write_transition_dipole(const std::string& filename); + /// @brief calculate transition dipole in velocity gauge using ks eigenvalues instead of excitation energies + void test_transition_dipoles_velocity_ks(const double* const ks_eig); + //====================================================================================================== private: /// $$2/3\Omega\sum_{ia\sigma} |\braket{\psi_{i}|\mathbf{r}|\psi_{a}} |^2\int \rho_{\alpha\beta}(\mathbf{r}) \mathbf{r} d\mathbf{r}$$ - void oscillator_strength(const Grid_Driver& gd, const std::vector& orb_cutoff); - const int nspin_x = 1; ///< 1 for singlet/triplet, 2 for updown(openshell) - const int naos = 1; - const std::vector& nocc; - const std::vector& nvirt; - const int nk = 1; - const int nstate = 1; - const int ldim = 1; ///< local leading dimension of X, or the data size of each state - const int gdim = 1; ///< global leading dimension of X - const double ana_thr = 0.3; ///< {abs(X) > thr} will appear in the transition analysis log - const double* eig; - const T* X; - const K_Vectors& kv; - std::vector> psi_ks; - const std::vector& pX; - const Parallel_2D& pc; - const Parallel_Orbitals& pmat; - typename TGint::type* gint = nullptr; - const ModulePW::PW_Basis& rho_basis; - const UnitCell& ucell; + void oscillator_strength(); + /// calculate the transition dipole of state S in length gauge: $\sum_{iak}X^S_{iak}$ + ModuleBase::Vector3 cal_transition_dipole_istate_length(const int istate); + /// calculate the transition dipole of all states in length gauge + void cal_transition_dipoles_length(); + /// calculate the transition dipole of state S in velocity gauge: $i(\sum_{iak}X^S_{iak})/\Omega_S$ + ModuleBase::Vector3 cal_transition_dipole_istate_velocity_R(const int istate, const TD_current& vR); + ModuleBase::Vector3 cal_transition_dipole_istate_velocity_k(const int istate, const TD_current& vR); + /// calculate the transition dipole of all states in velocity gauge + void cal_transition_dipoles_velocity(); + double cal_mean_squared_dipole(ModuleBase::Vector3 dipole); + /// calculate the transition density matrix + elecstate::DensityMatrix cal_transition_density_matrix(const int istate, const T* X_in = nullptr, const bool need_R = true); + const int nspin_x = 1; ///< 1 for singlet/triplet, 2 for updown(openshell) + const int naos = 1; + const std::vector& nocc; + const std::vector& nvirt; + const int nk = 1; + const int nstate = 1; + const int ldim = 1;///< local leading dimension of X, or the data size of each state + const int gdim = 1;///< global leading dimension of X + const double ana_thr = 0.3; ///< {abs(X) > thr} will appear in the transition analysis log + const double* eig; + const T* X; + const K_Vectors& kv; + std::vector> psi_ks; + const std::vector& pX; + const Parallel_2D& pc; + const Parallel_Orbitals& pmat; + typename TGint::type* gint = nullptr; + const ModulePW::PW_Basis& rho_basis; + const Grid_Driver& gd_; + const UnitCell& ucell; + const std::vector& orb_cutoff_; + const TwoCenterBundle& two_center_bundle_; - void cal_gint_rho(double** rho, const int& nrxx); - std::map get_pair_info( - const int i); ///< given the index in X, return its ispin, ik, iocc, ivirt + void cal_gint_rho(double** rho, const int& nrxx); + std::map get_pair_info(const int i); ///< given the index in X, return its ispin, ik, iocc, ivirt - std::vector> transition_dipole_; ///< $\braket{ \psi_{i} | \mathbf{r} | \psi_{a} }$ - std::vector - oscillator_strength_; ///< $2/3\Omega |\sum_{ia\sigma} \braket{\psi_{i}|\mathbf{r}|\psi_{a}} |^2$ + std::vector> transition_dipole_; ///< $\braket{ \psi_{i} | \mathbf{r} | \psi_{a} }$ + std::vector mean_squared_transition_dipole_; /// $|dipole|^2/3$, atomic unit (Hartree) + std::vector oscillator_strength_;///< $2/3\Omega |\sum_{ia\sigma} \braket{\psi_{i}|\mathbf{r}|\psi_{a}} |^2$, atomic unit (Hartree) }; } diff --git a/source/module_lr/lr_spectrum_velocity.cpp b/source/module_lr/lr_spectrum_velocity.cpp new file mode 100644 index 0000000000..5a6927f74d --- /dev/null +++ b/source/module_lr/lr_spectrum_velocity.cpp @@ -0,0 +1,181 @@ +#include "lr_spectrum.h" +#include "module_lr/dm_trans/dm_trans.h" +#include "module_lr/utils/lr_util_hcontainer.h" +#include "math.h" +#include "module_parameter/parameter.h" +#include "module_hamilt_lcao/module_hcontainer/hcontainer_funcs.h" +namespace LR +{ + /// get the velocity matrix v(R) + inline TD_current get_velocity_matrix_R(const UnitCell& ucell, + const Grid_Driver& gd, + const Parallel_Orbitals& pmat, + const TwoCenterBundle& two_center_bundle) + { + // convert the orbital object to the old class for TD_current + LCAO_Orbitals orb; + const auto& inp = PARAM.inp; + two_center_bundle.to_LCAO_Orbitals(orb, inp.lcao_ecut, inp.lcao_dk, inp.lcao_dr, inp.lcao_rmax); + // actually this class calculates the velocity matrix v(R) at A=0 + TD_current vR(&ucell, &gd, &pmat, orb, two_center_bundle.overlap_orb.get()); + vR.calculate_vcomm_r(); // $<\mu, 0|[Vnl, r]|\nu, R>$ + vR.calculate_grad_term(); // $<\mu, 0|\nabla|\nu, R>$ + return vR; + } + + inline double lorentz_delta(const double freq_diff, const double eta) + { + return eta / (freq_diff * freq_diff + eta * eta) / M_PI; + } + + template inline ModuleBase::Vector3 convert_vector_to_vector3(const std::vector>& vec); + template<> inline ModuleBase::Vector3 convert_vector_to_vector3(const std::vector>& vec) + { + assert(vec.size() == 3); + return ModuleBase::Vector3(vec[0].real(), vec[1].real(), vec[2].real()); + } + template<> inline ModuleBase::Vector3> convert_vector_to_vector3(const std::vector>& vec) + { + assert(vec.size() == 3); + return ModuleBase::Vector3>(vec[0], vec[1], vec[2]); + } + + /// this algorithm has bug in multi-k cases, just for test + template + ModuleBase::Vector3 LR::LR_Spectrum::cal_transition_dipole_istate_velocity_R(const int istate, const TD_current& vR) + { + // transition density matrix D(R) + const elecstate::DensityMatrix& DM_trans = this->cal_transition_density_matrix(istate); + + std::vector> trans_dipole(3, 0.0); // $=\sum_{uvR} v(R) D(R) = \sum_{iak}X_{iak}$ + const std::complex fac = ModuleBase::IMAG_UNIT / (eig[istate] / ModuleBase::e2); // eV to Hartree + for (int i = 0; i < 3; i++) + { + for (int is = 0;is < this->nspin_x; ++is) + { + trans_dipole[i] += LR_Util::dot_R_matrix(*vR.get_current_term_pointer(i), *DM_trans.get_DMR_pointer(is + 1), ucell.nat) * fac; + } // end for spin_x, only matter in open-shell system + trans_dipole[i] *= static_cast(this->nk); // nk is divided inside DM_trans, now recover it + if (this->nspin_x == 1) { trans_dipole[i] *= sqrt(2.0); } // *2 for 2 spins, /sqrt(2) for the halfed dimension of X in the normalizaiton + Parallel_Reduce::reduce_all(trans_dipole[i]); + } // end for direction + return convert_vector_to_vector3(trans_dipole); + } + + // this algorithm is actually in use + template + ModuleBase::Vector3 LR::LR_Spectrum::cal_transition_dipole_istate_velocity_k(const int istate, const TD_current& vR) + { + // transition density matrix D(R) + const elecstate::DensityMatrix& DM_trans = this->cal_transition_density_matrix(istate, this->X, false); + + std::vector> trans_dipole(3, 0.0); // $=\sum_{uvR} v(R) D(R) = \sum_{iak}X_{iak}$ + const std::complex fac = ModuleBase::IMAG_UNIT / (eig[istate] / ModuleBase::e2); // eV to Hartree + for (int i = 0; i < 3; i++) + { + for (int is = 0;is < this->nspin_x;++is) + { + for (int ik = 0;ik < nk;++ik) + { + std::vector> vk(pmat.get_local_size(), 0.0); + hamilt::folding_HR(*vR.get_current_term_pointer(i), vk.data(), kv.kvec_d[ik], pmat.get_row_size(), 1); + trans_dipole[i] += std::inner_product(vk.begin(), vk.end(), DM_trans.get_DMK_pointer(is * nk + ik), std::complex(0., 0.)) * fac; + } + } // end for spin_x, only matter in open-shell system + trans_dipole[i] *= static_cast(this->nk); // nk is divided inside DM_trans, now recover it + if (this->nspin_x == 1) { trans_dipole[i] *= sqrt(2.0); } // *2 for 2 spins, /sqrt(2) for the halfed dimension of X in the normalizaiton + Parallel_Reduce::reduce_all(trans_dipole[i]); + } // end for direction + return convert_vector_to_vector3(trans_dipole); + } + + template + void LR::LR_Spectrum::cal_transition_dipoles_velocity() + { + const TD_current& vR = get_velocity_matrix_R(ucell, gd_, pmat, two_center_bundle_); // velocity matrix v(R) + transition_dipole_.resize(nstate); + this->mean_squared_transition_dipole_.resize(nstate); + for (int istate = 0;istate < nstate;++istate) + { + transition_dipole_[istate] = cal_transition_dipole_istate_velocity_k(istate, vR); + mean_squared_transition_dipole_[istate] = cal_mean_squared_dipole(transition_dipole_[istate]); + } + } + + template + void LR::LR_Spectrum::optical_absorption_method2(const std::vector& freq, const double eta) + { + ModuleBase::TITLE("LR::LR_Spectrum", "optical_absorption_velocity"); + // 4*pi^2/V * mean_squared_dipole *delta(w-Omega_S) + std::ofstream ofs(PARAM.globalv.global_out_dir + "absorption.dat"); + if (GlobalV::MY_RANK == 0) { ofs << "Frequency (eV) | wave length(nm) | Absorption (a.u.)" << std::endl; } + const double fac = 4 * M_PI * M_PI / ucell.omega * ModuleBase::e2 / this->nk; // e2: Ry to Hartree in the denominator + for (int f = 0;f < freq.size();++f) + { + double abs_value = 0.0; + for (int i = 0;i < nstate;++i) + { + abs_value += this->mean_squared_transition_dipole_[i] * lorentz_delta((freq[f] - eig[i]), eta); + } + abs_value *= fac; + if (GlobalV::MY_RANK == 0) { ofs << freq[f] * ModuleBase::Ry_to_eV << "\t" << 91.126664 / freq[f] << "\t" << abs_value << std::endl; } + } + } + + inline void cal_eig_ks_diff(double* const eig_ks_diff, const double* const eig_ks, const Parallel_2D& px, const int nk, const int nocc, const int nvirt) + { + for (int ik = 0;ik < nk;++ik) + { + const int& start_k = ik * (nocc + nvirt); + for (int io = 0;io < px.get_col_size();++io) //nocc_local + { + for (int iv = 0;iv < px.get_row_size();++iv) //nvirt_local + { + int io_g = px.local2global_col(io); + int iv_g = px.local2global_row(iv); + eig_ks_diff[ik * px.get_local_size() + io * px.get_row_size() + iv] = (eig_ks[start_k + nocc + iv_g] - eig_ks[start_k + io_g]) / ModuleBase::e2; // eV to Hartree + } + } + } + } + + template + void LR::LR_Spectrum::test_transition_dipoles_velocity_ks(const double* const ks_eig) + { + // velocity matrix v(R) + const TD_current& vR = get_velocity_matrix_R(ucell, gd_, pmat, two_center_bundle_); + // (e_c-e_v) of KS eigenvalues + std::vector eig_ks_diff(this->ldim); + for (int is = 0;is < this->nspin_x;++is) + { + cal_eig_ks_diff(eig_ks_diff.data() + is * nk * pX[0].get_local_size(), ks_eig, pX[is], nk, nocc[is], nvirt[is]); + } + // X/(ec-ev) + std::vector X_div_ks_eig(nstate * this->ldim); + for (int istate = 0;istate < nstate;++istate) + { + const int st = istate * this->ldim; + std::transform(X + st, X + st + ldim, eig_ks_diff.begin(), X_div_ks_eig.data() + st, std::divides()); + } + + this->transition_dipole_.resize(nstate); + this->mean_squared_transition_dipole_.resize(nstate); + for (int istate = 0;istate < nstate;++istate) + { + // transition density matrix D(R) + const elecstate::DensityMatrix& DM_trans = this->cal_transition_density_matrix(istate, X_div_ks_eig.data()); + std::vector> tmp_trans_dipole(3, 0.0); + for (int i = 0; i < 3; i++) + { + for (int is = 0;is < this->nspin_x; ++is) + { + tmp_trans_dipole[i] += LR_Util::dot_R_matrix(*vR.get_current_term_pointer(i), *DM_trans.get_DMR_pointer(is + 1), ucell.nat) * ModuleBase::IMAG_UNIT; + } // end for spin_x, only matter in open-shell system + } // end for direction + this->transition_dipole_[istate] = convert_vector_to_vector3(tmp_trans_dipole); + this->mean_squared_transition_dipole_[istate] = cal_mean_squared_dipole(transition_dipole_[istate]); + } + } +} +template class LR::LR_Spectrum; +template class LR::LR_Spectrum>; \ No newline at end of file diff --git a/source/module_lr/utils/lr_util.cpp b/source/module_lr/utils/lr_util.cpp index 1418764126..ceb9501e0a 100644 --- a/source/module_lr/utils/lr_util.cpp +++ b/source/module_lr/utils/lr_util.cpp @@ -179,4 +179,17 @@ namespace LR_Util vl.data(), &ldvl, vr.data(), &ldvr, work2.data(), &lwork, rwork.data(), &info); if (info) { std::cout << "ERROR: Lapack solver zgeev, info=" << info << std::endl; } } + + std::string tolower(const std::string& str) + { + std::string str_lower = str; + std::transform(str_lower.begin(), str_lower.end(), str_lower.begin(), ::tolower); + return str_lower; + } + std::string toupper(const std::string& str) + { + std::string str_upper = str; + std::transform(str_upper.begin(), str_upper.end(), str_upper.begin(), ::toupper); + return str_upper; + } } \ No newline at end of file diff --git a/source/module_lr/utils/lr_util.h b/source/module_lr/utils/lr_util.h index 9e5ec32d91..3dafe2dd0a 100644 --- a/source/module_lr/utils/lr_util.h +++ b/source/module_lr/utils/lr_util.h @@ -103,5 +103,9 @@ namespace LR_Util /// @brief diagonalize a general matrix void diag_lapack_nh(const int& n, double* mat, std::complex* eig); void diag_lapack_nh(const int& n, std::complex* mat, std::complex* eig); + + ///=================string option==================== + std::string tolower(const std::string& str); + std::string toupper(const std::string& str); } #include "lr_util.hpp" diff --git a/source/module_lr/utils/lr_util_hcontainer.h b/source/module_lr/utils/lr_util_hcontainer.h index 715e0dad8b..8937d9d733 100644 --- a/source/module_lr/utils/lr_util_hcontainer.h +++ b/source/module_lr/utils/lr_util_hcontainer.h @@ -1,5 +1,7 @@ #pragma once #include "module_elecstate/module_dm/density_matrix.h" +#include +#include "module_base/parallel_reduce.h" namespace LR_Util { template @@ -86,4 +88,33 @@ namespace LR_Util initialize_HR(hR_tmp, ucell, gd, orb_cutoff); dm.init_DMR(hR_tmp); } + + /// $\sum_{uvR} H1_{uv}(R) H2_{uv}(R)$ + template + TR1 dot_R_matrix(const hamilt::HContainer& h1, const hamilt::HContainer& h2, const int& nat) + { + const auto& pmat = *h1.get_paraV(); + TR1 sum = 0; + // in case of the different order of atom pair and R-index in h1 and h2, we search by value instead of index + for (int iat1 = 0;iat1 < nat;++iat1) + { + for (int iat2 = 0;iat2 < nat;++iat2) + { + auto ap1 = h1.find_pair(iat1, iat2); + if (!ap1) { continue; } + auto ap2 = h2.find_pair(iat1, iat2); + assert(ap2); + for (int iR = 0;iR < ap1->get_R_size();++iR) + { + const ModuleBase::Vector3& R = ap1->get_R_index(iR); + auto mat1 = ap1->get_HR_values(R.x, R.y, R.z); + auto mat2 = ap2->get_HR_values(R.x, R.y, R.z); + sum += std::inner_product(mat1.get_pointer(), mat1.get_pointer() + mat1.get_memory_size(), mat2.get_pointer(), (TR1)0.0); + } + } + } + Parallel_Reduce::reduce_all(sum); + return sum; + } + } \ No newline at end of file diff --git a/source/module_parameter/input_parameter.h b/source/module_parameter/input_parameter.h index 9b2151c945..75b2647f76 100644 --- a/source/module_parameter/input_parameter.h +++ b/source/module_parameter/input_parameter.h @@ -311,6 +311,7 @@ struct Input_para bool lr_unrestricted = false; ///< whether to use the unrestricted construction for LR-TDDFT std::vector abs_wavelen_range = {}; ///< the range of wavelength(nm) to output the absorption spectrum double abs_broadening = 0.01; ///< the broadening (eta) for LR-TDDFT absorption spectrum + std::string abs_gauge = "length"; ///< whether to use length or velocity gauge to calculate the absorption spectrum in LR-TDDFT std::string ri_hartree_benchmark = "none"; ///< whether to use the RI approximation for the Hartree potential in LR-TDDFT for benchmark (with FHI-aims/ABACUS read-in style) std::vector aims_nbasis = {}; ///< the number of basis functions for each atom type used in FHI-aims (for benchmark) // ============== #Parameters (11.Output) ===========================