Skip to content

enable out_velocity for RPA #1

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
AroundPeking merged 1 commit into from
Dec 1, 2025
Merged

enable out_velocity for RPA #1

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions source/module_io/ctrl_output_lcao.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
33 changes: 32 additions & 1 deletion source/module_lr/ao_to_mo_transformer/ao_to_mo.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<typename T>
void ao_to_mo_forloop_serial(
const std::vector<container::Tensor>& mat_ao,
Expand Down
26 changes: 14 additions & 12 deletions source/module_lr/ao_to_mo_transformer/ao_to_mo_parallel.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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);
Expand Down
30 changes: 17 additions & 13 deletions source/module_lr/ao_to_mo_transformer/ao_to_mo_serial.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -51,7 +51,7 @@ namespace LR
std::complex<double>* 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;
Expand Down Expand Up @@ -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)
{
Expand Down Expand Up @@ -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)
{
Expand Down
10 changes: 8 additions & 2 deletions source/module_lr/dm_trans/dm_trans.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down
195 changes: 195 additions & 0 deletions source/module_lr/utils/spectrum_mo.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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<typename T>
std::vector<std::complex<double>> 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<T>& psi_ks,
const int nk,
const int nspin_tmp,
const int nbasis,
const std::vector<int> nocc,
const std::vector<int> 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..."<<std::endl;
// get_velocity_matrix_R(ucell, gd_, pmat, two_center_bundle_);
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);

TD_current vR(&ucell, &gd, &pmat, orb, two_center_bundle.overlap_orb.get());
vR.calculate_grad_term(); // $<\mu, 0|-i∇r|\nu, R>$
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<std::complex<double>> 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<T> to c_psi_ks<complex<double>>, ensure complex<double> for dipole calculation
psi::Psi<std::complex<double>> 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<double>(psi_ks(iks, ic, ir));
}
}
}

//2. calculate v_mo = c^\dagger v c
std::vector<ct::Tensor> vk(nks, LR_Util::newTensor<std::complex<double>>({ pmat.get_col_size(), pmat.get_row_size() }));
for (int id = 0; id < 3; ++id)
{
for (auto& v : vk) v.zero();

std::vector<std::complex<double>> 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<std::complex<double>>(), 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<std::complex<double>>& 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<std::complex<double>>& 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;
}
}
}

}
}
}
}
6 changes: 6 additions & 0 deletions source/module_ri/RPA_LRI.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,6 +69,12 @@ class RPA_LRI
void out_pure_ri_tensor(const std::string fn, RI::Tensor<std::complex<double>>& olp, const double threshold);
void out_pure_ri_tensor(const std::string fn, RI::Tensor<double>& 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<T>& psi,
const elecstate::ElecState* pelec);

void out_Cs(const UnitCell& ucell,
std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>& Cs_in,
Expand Down
Loading