From 6af22e1f6d175eddbc36e7897a6ffc04d31738b7 Mon Sep 17 00:00:00 2001 From: Fisherd99 Date: Mon, 1 Dec 2025 12:25:20 +0800 Subject: [PATCH] enable out_velocity for RPA --- source/module_io/ctrl_output_lcao.cpp | 1 + .../module_lr/ao_to_mo_transformer/ao_to_mo.h | 33 ++- .../ao_to_mo_parallel.cpp | 26 +-- .../ao_to_mo_transformer/ao_to_mo_serial.cpp | 30 +-- source/module_lr/dm_trans/dm_trans.h | 10 +- source/module_lr/utils/spectrum_mo.hpp | 195 ++++++++++++++++++ source/module_ri/RPA_LRI.h | 6 + source/module_ri/RPA_LRI.hpp | 39 ++++ 8 files changed, 312 insertions(+), 28 deletions(-) create mode 100644 source/module_lr/utils/spectrum_mo.hpp diff --git a/source/module_io/ctrl_output_lcao.cpp b/source/module_io/ctrl_output_lcao.cpp index 31b9e589c3..30fff13ccc 100644 --- a/source/module_io/ctrl_output_lcao.cpp +++ b/source/module_io/ctrl_output_lcao.cpp @@ -340,6 +340,7 @@ void ctrl_output_lcao(UnitCell& ucell, orb); rpa_lri_double.init(MPI_COMM_WORLD, kv, orb.cutoffs()); rpa_lri_double.out_for_RPA(ucell, pv, *psi, pelec, kv, orb); + rpa_lri_double.out_velocity(ucell, gd, two_center_bundle, pv, *psi, pelec); } #endif diff --git a/source/module_lr/ao_to_mo_transformer/ao_to_mo.h b/source/module_lr/ao_to_mo_transformer/ao_to_mo.h index 23d42271d8..21e2f523d8 100644 --- a/source/module_lr/ao_to_mo_transformer/ao_to_mo.h +++ b/source/module_lr/ao_to_mo_transformer/ao_to_mo.h @@ -9,8 +9,39 @@ namespace LR { #ifndef MO_TYPE_H #define MO_TYPE_H - enum MO_TYPE { OO, VO, VV }; + enum MO_TYPE { OO, VO, VV, ALL }; #endif +/* +MO_TYPE: OO VO VV ALL +nmo1 nocc nocc nvirt nocc+nvirt +nmo2 nocc nvirt nvirt nocc+nvirt +imo1 0 0 nocc 0 +imo2 0 nocc nocc 0 +*/ + inline void set_dim(const MO_TYPE type, const int& nocc, const int& nvirt, + int& nmo1, int& nmo2, int& imo1, int& imo2) + { + switch(type) + { + case MO_TYPE::OO: + nmo1 = nocc; nmo2 = nocc; imo1 = 0; imo2 = 0; + break; + case MO_TYPE::VO: + nmo1 = nocc; nmo2 = nvirt; imo1 = 0; imo2 = nocc; + break; + case MO_TYPE::VV: + nmo1 = nvirt; nmo2 = nvirt; imo1 = nocc; imo2 = nocc; + break; + case MO_TYPE::ALL: + nmo1 = nocc + nvirt; + nmo2 = nocc + nvirt; + imo1 = 0; + imo2 = 0; + break; + default: + throw std::runtime_error("Error in LR::set_dim: unknown MO_TYPE"); + } + } template void ao_to_mo_forloop_serial( const std::vector& mat_ao, diff --git a/source/module_lr/ao_to_mo_transformer/ao_to_mo_parallel.cpp b/source/module_lr/ao_to_mo_transformer/ao_to_mo_parallel.cpp index 2cadd4ac0c..86f0f061a9 100644 --- a/source/module_lr/ao_to_mo_transformer/ao_to_mo_parallel.cpp +++ b/source/module_lr/ao_to_mo_transformer/ao_to_mo_parallel.cpp @@ -23,18 +23,19 @@ namespace LR const bool add_on, const MO_TYPE type) { - ModuleBase::TITLE("hamilt_lrtd", "ao_to_mo_pblas"); + ModuleBase::TITLE("LR", "ao_to_mo_pblas"); assert(pmat_ao.comm() == pcoeff.comm() && pmat_ao.comm() == pmat_mo.comm()); assert(pmat_ao.blacs_ctxt == pcoeff.blacs_ctxt && pmat_ao.blacs_ctxt == pmat_mo.blacs_ctxt); assert(pmat_mo.get_local_size() > 0); const int nks = mat_ao.size(); + int nmo1_set, nmo2_set, imo1_set, imo2_set; + set_dim(type, nocc, nvirt, nmo1_set, nmo2_set, imo1_set, imo2_set); + const int nmo1 = nmo1_set; + const int nmo2 = nmo2_set; + const int imo1 = imo1_set + 1; + const int imo2 = imo2_set + 1; const int i1 = 1; - const int ivirt = nocc + 1; - const int nmo1 = type == MO_TYPE::VV ? nvirt : nocc; - const int nmo2 = type == MO_TYPE::OO ? nocc : nvirt; - const int imo1 = type == MO_TYPE::VV ? ivirt : i1; - const int imo2 = type == MO_TYPE::OO ? i1 : ivirt; Parallel_2D pVc; // for intermediate Vc LR_Util::setup_2d_division(pVc, pmat_ao.get_block_size(), naos, nmo1, pmat_ao.blacs_ctxt); @@ -81,18 +82,19 @@ namespace LR const bool add_on, const MO_TYPE type) { - ModuleBase::TITLE("hamilt_lrtd", "cal_AX_plas"); + ModuleBase::TITLE("LR", "ao_to_mo_pblas"); assert(pmat_ao.comm() == pcoeff.comm() && pmat_ao.comm() == pmat_mo.comm()); assert(pmat_ao.blacs_ctxt == pcoeff.blacs_ctxt && pmat_ao.blacs_ctxt == pmat_mo.blacs_ctxt); assert(pmat_mo.get_local_size() > 0); const int nks = mat_ao.size(); + int nmo1_set, nmo2_set, imo1_set, imo2_set; + set_dim(type, nocc, nvirt, nmo1_set, nmo2_set, imo1_set, imo2_set); + const int nmo1 = nmo1_set; + const int nmo2 = nmo2_set; + const int imo1 = imo1_set + 1; + const int imo2 = imo2_set + 1; const int i1 = 1; - const int ivirt = nocc + 1; - const int nmo1 = type == MO_TYPE::VV ? nvirt : nocc; - const int nmo2 = type == MO_TYPE::OO ? nocc : nvirt; - const int imo1 = type == MO_TYPE::VV ? ivirt : i1; - const int imo2 = type == MO_TYPE::OO ? i1 : ivirt; Parallel_2D pVc; // for intermediate Vc LR_Util::setup_2d_division(pVc, pmat_ao.get_block_size(), naos, nmo1, pmat_ao.blacs_ctxt); diff --git a/source/module_lr/ao_to_mo_transformer/ao_to_mo_serial.cpp b/source/module_lr/ao_to_mo_transformer/ao_to_mo_serial.cpp index 04573381a7..1d8c09575c 100644 --- a/source/module_lr/ao_to_mo_transformer/ao_to_mo_serial.cpp +++ b/source/module_lr/ao_to_mo_transformer/ao_to_mo_serial.cpp @@ -13,7 +13,7 @@ namespace LR double* mat_mo, const MO_TYPE type) { - ModuleBase::TITLE("hamilt_lrtd", "ao_to_mo_forloop_serial"); + ModuleBase::TITLE("LR", "ao_to_mo_forloop_serial"); const int nks = mat_ao.size(); const int naos = coeff.get_nbasis(); const int nmo1 = type == MO_TYPE::VV ? nvirt : nocc; @@ -51,7 +51,7 @@ namespace LR std::complex* const mat_mo, const MO_TYPE type) { - ModuleBase::TITLE("hamilt_lrtd", "ao_to_mo_forloop_serial"); + ModuleBase::TITLE("LR", "ao_to_mo_forloop_serial"); const int nks = mat_ao.size(); const int naos = coeff.get_nbasis(); const int nmo1 = type == MO_TYPE::VV ? nvirt : nocc; @@ -90,13 +90,15 @@ namespace LR const bool add_on, const MO_TYPE type) { - ModuleBase::TITLE("hamilt_lrtd", "ao_to_mo_blas"); + ModuleBase::TITLE("LR", "ao_to_mo_blas"); const int nks = mat_ao.size(); - const int naos = coeff.get_nbasis(); - const int nmo1 = type == MO_TYPE::VV ? nvirt : nocc; - const int nmo2 = type == MO_TYPE::OO ? nocc : nvirt; - const int imo1 = type == MO_TYPE::VV ? nocc : 0; - const int imo2 = type == MO_TYPE::OO ? 0 : nocc; + const int naos = coeff.get_nbasis(); + int nmo1_set, nmo2_set, imo1_set, imo2_set; + set_dim(type, nocc, nvirt, nmo1_set, nmo2_set, imo1_set, imo2_set); + const int nmo1 = nmo1_set; + const int nmo2 = nmo2_set; + const int imo1 = imo1_set; + const int imo2 = imo2_set; for (int isk = 0;isk < nks;++isk) { @@ -131,13 +133,15 @@ namespace LR const bool add_on, const MO_TYPE type) { - ModuleBase::TITLE("hamilt_lrtd", "ao_to_mo_blas"); + ModuleBase::TITLE("LR", "ao_to_mo_blas"); const int nks = mat_ao.size(); const int naos = coeff.get_nbasis(); - const int nmo1 = type == MO_TYPE::VV ? nvirt : nocc; - const int nmo2 = type == MO_TYPE::OO ? nocc : nvirt; - const int imo1 = type == MO_TYPE::VV ? nocc : 0; - const int imo2 = type == MO_TYPE::OO ? 0 : nocc; + int nmo1_set, nmo2_set, imo1_set, imo2_set; + set_dim(type, nocc, nvirt, nmo1_set, nmo2_set, imo1_set, imo2_set); + const int nmo1 = nmo1_set; + const int nmo2 = nmo2_set; + const int imo1 = imo1_set; + const int imo2 = imo2_set; for (int isk = 0;isk < nks;++isk) { diff --git a/source/module_lr/dm_trans/dm_trans.h b/source/module_lr/dm_trans/dm_trans.h index a02f2999d0..109982b45b 100644 --- a/source/module_lr/dm_trans/dm_trans.h +++ b/source/module_lr/dm_trans/dm_trans.h @@ -11,9 +11,15 @@ namespace LR #ifndef MO_TYPE_H #define MO_TYPE_H - enum MO_TYPE { OO, VO, VV }; + enum MO_TYPE { OO, VO, VV, ALL }; #endif - +/* +MO_TYPE: OO VO VV ALL(not used in dm_trans) +nmo1 nocc nocc nvirt nocc+nvirt +nmo2 nocc nvirt nvirt nocc+nvirt +imo1 0 0 nocc 0 +imo2 0 nocc nocc 0 +*/ #ifdef __MPI /// @brief calculate the 2d-block transition density matrix in AO basis using p?gemm /// \f[ \tilde{\rho}_{\mu_j\mu_b}=\sum_{jb}c_{j,\mu_j}X_{jb}c^*_{b,\mu_b} \f] diff --git a/source/module_lr/utils/spectrum_mo.hpp b/source/module_lr/utils/spectrum_mo.hpp new file mode 100644 index 0000000000..21e9882b38 --- /dev/null +++ b/source/module_lr/utils/spectrum_mo.hpp @@ -0,0 +1,195 @@ +#pragma once +#include "module_parameter/parameter.h" +#include "module_basis/module_nao/two_center_bundle.h" +#include "module_cell/klist.h" +#include "module_hamilt_lcao/module_hcontainer/hcontainer_funcs.h" +#include "module_lr/ao_to_mo_transformer/ao_to_mo.h" +#include "module_lr/utils/lr_util.h" +#include "module_lr/utils/lr_util_print.h" +#include "module_hamilt_lcao/module_tddft/td_current.h" +#include "module_psi/psi.h" + +namespace LR_Util +{ +template +std::vector> cal_velocity_mo(const UnitCell& ucell, + const Grid_Driver& gd, + const TwoCenterBundle& two_center_bundle, + const Parallel_Orbitals& pmat,/*nbasis×nbasis*/ + const Parallel_2D& pc,/*nbasis×nbands*/ + const K_Vectors& kv, + const psi::Psi& psi_ks, + const int nk, + const int nspin_tmp, + const int nbasis, + const std::vector nocc, + const std::vector nvirt) +{ + ModuleBase::TITLE("LR::LR_Util", "cal_velocity_mo"); + ModuleBase::timer::tick("LR::LR_Util", "cal_velocity_mo"); + std::cout<<"Calculating velocity matrix in KS presentation..."<$ + vR.calculate_vcomm_r(); // $<\mu, 0|i[Vnl, r]|\nu, R>$ + + int nks = kv.get_nks(); // include spin + assert(nks == nk * nspin_tmp); + assert(psi_ks.get_nk() == nks); + int KS_num = nocc[0] + nvirt[0]; + std::vector> velocity_mo(nspin_tmp * 3 * nk * KS_num * KS_num, 0.0); + Parallel_2D pmo; + LR_Util::setup_2d_division(pmo, pmat.get_block_size(), KS_num, KS_num +#ifdef __MPI + , pc.blacs_ctxt +#endif + ); + + //1. psi_ks to c_psi_ks>, ensure complex for dipole calculation + psi::Psi> c_psi_ks(nks, + pc.get_col_size(), + pc.get_row_size(), + kv.ngk, + true); + for(int iks = 0; iks < nks; ++iks) + { + for(int ic = 0; ic < pc.get_col_size(); ++ic)//band + { + for(int ir = 0; ir < pc.get_row_size(); ++ir)//basis + { + c_psi_ks(iks, ic, ir) = std::complex(psi_ks(iks, ic, ir)); + } + } + } + + //2. calculate v_mo = c^\dagger v c + std::vector vk(nks, LR_Util::newTensor>({ pmat.get_col_size(), pmat.get_row_size() })); + for (int id = 0; id < 3; ++id) + { + for (auto& v : vk) v.zero(); + + std::vector> v_mo(nks * pmo.get_local_size(), 0.0); + + for (int is = 0; is < nspin_tmp; ++is) + { + assert(KS_num == nocc[is] + nvirt[is]); + + for (int ik = is*nk; ik < nk; ++ik) + { + hamilt::folding_HR(*vR.get_current_term_pointer(id), vk[ik].data>(), kv.kvec_d[ik], pmat.get_row_size(), 1/*column-major*/); + } + } +#ifdef __MPI + ao_to_mo_pblas(vk, pmat, c_psi_ks, pc, nbasis, + nocc[0], nvirt[0], pmo, v_mo.data(), + false, // add_on + LR::MO_TYPE::ALL); +#else + ao_to_mo_blas(vk, c_psi_ks, + nocc[0], nvirt[0], v_mo.data(), + false , //add_on + LR::MO_TYPE::ALL); +#endif + // gather local vk to global velocity_mo + for (int is = 0; is < nspin_tmp; ++is) + { + for (int ik = 0; ik < nk; ++ik) + { + LR_Util::gather_2d_to_full(pmo, v_mo.data() + ik * pmo.get_local_size(), + &velocity_mo[(is * 3 * nk + id * nk + ik) * KS_num * KS_num ], + false/*col_first*/, KS_num, KS_num); + + //std::cout<< "is" << is << "id: " << id << " ik: " << ik << " v_mo: " << std::endl; + //LR_Util::print_value(velocity_mo.data()+(is * 3 * nk + id * nk + ik) * KS_num * KS_num, KS_num, KS_num); + } + } + }//id + ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "Finish velocity matrix in KS presentation."); + ModuleBase::timer::tick("LR::LR_Util", "cal_velocity_mo"); + return velocity_mo; +} + +/// @brief output the velocity matrix in KS presentation in human-read friendly format +inline void output_spectrum_mo(const std::vector>& out_spectrum_mo, + const std::string& filename, + const double* const eig_ks, + const int nk, + const int nspin_tmp, + const int KS_num, + const K_Vectors& kv) +{ + assert(out_spectrum_mo.size() == nspin_tmp * 3 * nk * KS_num * KS_num); + std::ofstream ofs(PARAM.globalv.global_out_dir + filename +".dat"); + ofs << "Data Unit (Hartree * Bohr) " << filename << std::endl; + ofs << "NOTICE: KS_index are restricted in nocc and nvirt" << std::endl; + for (int is = 0; is < nspin_tmp; ++is) + { + ofs << "ispin: " << is << std::endl; + for (int ik = 0;ik < nk;++ik) + { + ofs << "k-point: " << ik << " " << kv.kvec_d[ik] << std::endl; + ofs << std::setw(4) << "KS1" << std::setw(10) << "E1(eV)" << std::setw(6) << "KS2" << std::setw(10) << "E2(eV)" + << std::setw(16) << "x" << std::setw(23) << "|x|^2" << std::setw(18) << "y" << std::setw(23) <<"|y|^2" + << std::setw(18) << "z" << std::setw(23) <<"|z|^2" << std::setw(13) << "average" << std::endl; + + for (int i = 0; i < KS_num; ++i) + { + for (int j = i; j < KS_num ; ++j) + { + int ipair = (is * 3 * nk + ik) * KS_num * KS_num + i * KS_num + j; + int step = nk * KS_num * KS_num; + double average = (std::norm(out_spectrum_mo[ipair]) + std::norm(out_spectrum_mo[ipair + step]) + std::norm(out_spectrum_mo[ipair + 2 * step])) / 3.0; + ofs << std::setw(3) << i << std::setw(12) << std::setprecision(6) << eig_ks[i + ik*KS_num] * ModuleBase::Ry_to_eV + << std::setw(4) << j << std::setw(12) << eig_ks[j + ik*KS_num] * ModuleBase::Ry_to_eV + << std::setw(28) << out_spectrum_mo[ipair] << std::setw(13) << std::norm(out_spectrum_mo[ipair]) + << std::setw(28) << out_spectrum_mo[ipair + step] << std::setw(13) << std::norm(out_spectrum_mo[ipair + step]) + << std::setw(28) << out_spectrum_mo[ipair + 2 * step] << std::setw(13) << std::norm(out_spectrum_mo[ipair + 2 * step] ) + << std::setw(13) << average << std::endl; + } + } + } + } + ofs.close(); +} + +/// @brief output the velocity matrix in KS presentation in LibRPA format +inline void output_spectrum_mo_librpa(const std::vector>& out_spectrum_mo, + const std::string& filename, + const int nk, + const int nspin_tmp, + const int KS_num, + const K_Vectors& kv) +{ + assert(out_spectrum_mo.size() == nspin_tmp * 3 * nk * KS_num * KS_num); + std::ofstream ofs(filename); + ofs << std::scientific << nk << std::endl; + ofs << nspin_tmp << std::endl; + ofs << PARAM.inp.nbands << std::endl; + ofs << PARAM.globalv.nlocal << std::endl; + double HaBohrToEvAng = ModuleBase::Hartree_to_eV * ModuleBase::BOHR_TO_A; // 27.211396 * 0.5291770 + for (int is = 0; is < nspin_tmp; ++is) + { + for (int ik = 0; ik < nk; ++ik) + { + for (int id = 0; id < 3; ++id) + { + ofs <<" " << id + 1 << " " << ik + 1 << " " << is + 1 << std::endl; + for (int ib1 = 0; ib1 < KS_num; ++ib1) + { + for (int ib2 = 0; ib2 < KS_num; ++ib2) + { + int ipair = (is * 3 * nk + id * nk + ik) * KS_num * KS_num + ib1 * KS_num + ib2; + ofs << std::setw(26) << std::fixed << std::setprecision(16) << out_spectrum_mo[ipair].real() * HaBohrToEvAng + << std::setw(26) << out_spectrum_mo[ipair].imag() * HaBohrToEvAng << std::endl; + } + } + } + + } + } +} +} \ No newline at end of file diff --git a/source/module_ri/RPA_LRI.h b/source/module_ri/RPA_LRI.h index 62829c4e0d..1eb1ac35a9 100644 --- a/source/module_ri/RPA_LRI.h +++ b/source/module_ri/RPA_LRI.h @@ -69,6 +69,12 @@ class RPA_LRI void out_pure_ri_tensor(const std::string fn, RI::Tensor>& olp, const double threshold); void out_pure_ri_tensor(const std::string fn, RI::Tensor& olp, const double threshold); void out_bands(const elecstate::ElecState* pelec); + void out_velocity(const UnitCell& ucell, + const Grid_Driver& gd, + const TwoCenterBundle& two_center_bundle, + const Parallel_Orbitals& parav, /*nbasis×nbasis*/ + const psi::Psi& psi, + const elecstate::ElecState* pelec); void out_Cs(const UnitCell& ucell, std::map>>& Cs_in, diff --git a/source/module_ri/RPA_LRI.hpp b/source/module_ri/RPA_LRI.hpp index 2e5ea448b6..a143608f84 100644 --- a/source/module_ri/RPA_LRI.hpp +++ b/source/module_ri/RPA_LRI.hpp @@ -8,6 +8,7 @@ #include "RPA_LRI.h" #include "module_parameter/parameter.h" #include "module_ri/module_exx_symmetry/symmetry_rotation.h" +#include "module_lr/utils/spectrum_mo.hpp" #include #include @@ -1156,6 +1157,44 @@ void RPA_LRI::out_coulomb_k(const UnitCell& ucell, ofs.close(); } +template +void RPA_LRI::out_velocity(const UnitCell &ucell, + const Grid_Driver &gd, + const TwoCenterBundle &two_center_bundle, + const Parallel_Orbitals ¶v,/*nbasis×nbasis*/ + const psi::Psi &psi, + const elecstate::ElecState* pelec) +{ + ModuleBase::TITLE("DFT_RPA_interface", "out_velocity"); + ModuleBase::timer::tick("RPA_LRI", "out_velocity"); + + Parallel_2D parac;/*nbasis×nbands*/ + LR_Util::setup_2d_division(parac, parav.get_block_size(), PARAM.globalv.nlocal, PARAM.inp.nbands + #ifdef __MPI + , parav.blacs_ctxt + #endif + ); + + const int nk = PARAM.inp.nspin == 2 ? p_kv->get_nks() / 2 : p_kv->get_nks(); + const int nspin_tmp = PARAM.inp.nspin == 2 ? 2 : 1; + + // nocc and nvirt dosen't matter, their sum is actually used + std::vector nocc(2, PARAM.inp.nbands); + std::vector nvirt(2, 0); + + std::vector> velocity_mo = LR_Util::cal_velocity_mo(ucell, gd, two_center_bundle, + parav, parac, *this->p_kv, psi, nk, nspin_tmp, PARAM.globalv.nlocal, nocc, nvirt); + if (GlobalV::MY_RANK == 0){ + // for librpa readable + LR_Util::output_spectrum_mo_librpa(velocity_mo, "velocity_matrix", + nk, nspin_tmp, PARAM.inp.nbands, *this->p_kv); + // for human readable + LR_Util::output_spectrum_mo(velocity_mo, "velocity_matrix_rpa", pelec->ekb.c, + nk, nspin_tmp, PARAM.inp.nbands, *this->p_kv); + } + ModuleBase::timer::tick("RPA_LRI", "out_velocity"); +} + // template // void RPA_LRI::init(const MPI_Comm &mpi_comm_in) // {