Skip to content

Commit ceca2a8

Browse files
maki49maki49
committed
feature: output electron-hole density
fix: open-shell spectrum bug
1 parent a0ed499 commit ceca2a8

File tree

5 files changed

+145
-6
lines changed

5 files changed

+145
-6
lines changed

‎source/module_hamilt_pw/hamilt_pwdft/parallel_grid.h‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -40,6 +40,8 @@ class Parallel_Grid
4040
const int& ny = this->ncy;
4141
const int& nz = this->ncz;
4242

43+
int get_nrxx() const { return this->nrxx; }
44+
4345
private:
4446

4547
void z_distribution(void);

‎source/module_lr/Grad/multipliers/zeq_solver.hpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -125,6 +125,7 @@ namespace LR
125125
// 1. the right-hand side of Z-vector equation
126126
const int nloc_per_band = nk * px[0].get_local_size();
127127
container::Tensor R = LR_Util::newTensor<T>({ nstates, nloc_per_band });
128+
R.zero();
128129
Z_vector_R<T> ops_R(xc_kernel, nspin, naos, nocc, nvirt,
129130
ucell, orb_cutoff, gd, psi_ks, eig_ks,
130131
#ifdef __EXX

‎source/module_lr/esolver_lrtd_lcao.cpp‎

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@
66
#include "module_lr/potentials/pot_hxc_lrtd.h"
77
#include "module_lr/hsolver_lrtd.hpp"
88
#include "module_lr/lr_spectrum.h"
9+
#include "module_lr/lr_density.hpp"
910
#include <memory>
1011
#include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h"
1112
#include "module_io/read_wfc_nao.h"
@@ -577,6 +578,22 @@ void LR::ESolver_LR<T, TR>::after_all_runners(UnitCell& ucell)
577578
{
578579
ModuleBase::TITLE("ESolver_LR", "after_all_runners");
579580
if (input.ri_hartree_benchmark != "none") { return; } //no need to calculate the spectrum in the benchmark routine
581+
582+
// cal electron-hole density
583+
if (PARAM.inp.out_chg[0])
584+
{
585+
LR_Density<T> lr_density(gint_, ucell, kv, gd, *psi_ks, orb_cutoff_, Pgrid,
586+
nspin, nocc, nvirt, nbasis,
587+
paraX_, paraC_, paraMat_, openshell);
588+
589+
if (openshell)
590+
for (int is = 0;is < this->nspin;++is)
591+
lr_density.output_eh_density_all_states(this->X[0].template data<T>(), is, nstates);
592+
else
593+
for (int is = 0;is < this->X.size();++is)
594+
lr_density.output_eh_density_all_states(this->X[is].template data<T>(), is, nstates);
595+
}
596+
580597
//cal spectrum
581598
std::vector<double> freq(100);
582599
std::vector<double> abs_wavelen_range({ 20, 200 });//default range

‎source/module_lr/lr_density.hpp‎

Lines changed: 119 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,119 @@
1+
#pragma once
2+
#include "module_lr/utils/gint_template.h"
3+
#include "module_psi/psi.h"
4+
#include "module_lr/Grad/dm_diff/dm_diff.h"
5+
#include "module_lr/utils/lr_util_hcontainer.h"
6+
#include "module_io/cube_io.h"
7+
namespace LR
8+
{
9+
template<typename T>
10+
class LR_Density
11+
{
12+
typename TGint<T>::type* gint_;
13+
const UnitCell& ucell_;
14+
const K_Vectors& kv_;
15+
const Grid_Driver& gd_;
16+
const psi::Psi<T>& psi_ks_;
17+
const std::vector<double>& orb_cutoff_;
18+
const Parallel_Grid& pgrid_;
19+
const int nspin_;
20+
const std::vector<int>& nocc_;
21+
const std::vector<int>& nvirt_;
22+
const int nao_;
23+
const int nk_;
24+
const std::vector<Parallel_2D>& pX_;
25+
const Parallel_2D& pc_;
26+
const Parallel_Orbitals& pmat_;
27+
const bool openshell_;
28+
const std::vector<std::string> spintype_;
29+
30+
inline void dm_to_density(elecstate::DensityMatrix<double, double>& dm, double** density)
31+
{
32+
ModuleBase::TITLE("LR_Density", "dm_to_density");
33+
this->gint_->transfer_DM2DtoGrid(dm.get_DMR_vector());
34+
Gint_inout inout_rho(density, Gint_Tools::job_type::rho, 1, false);
35+
this->gint_->cal_gint(&inout_rho);
36+
}
37+
inline void dm_to_density(elecstate::DensityMatrix<complex<double>, complex<double>>& dm, double** density)
38+
{
39+
ModuleBase::TITLE("LR_Density", "dm_to_density");
40+
auto dm_to_density_real = [&](const char& part) -> void
41+
{
42+
elecstate::DensityMatrix<std::complex<double>, double> dm_real(&pmat_, 1, kv_.kvec_d, nk_);
43+
LR_Util::initialize_DMR<std::complex<double>, double>(dm_real, pmat_, ucell_, gd_, orb_cutoff_);
44+
LR_Util::get_DMR_real_imag_part(dm, dm_real, ucell_.nat, part);
45+
this->gint_->transfer_DM2DtoGrid(dm_real.get_DMR_vector());
46+
Gint_inout inout_rho(density, Gint_Tools::job_type::rho, 1, false);
47+
this->gint_->cal_gint(&inout_rho); // add-on
48+
};
49+
dm_to_density_real('R');
50+
dm_to_density_real('I');
51+
}
52+
53+
public:
54+
LR_Density(typename TGint<T>::type* gint,
55+
const UnitCell& ucell,
56+
const K_Vectors& kv,
57+
const Grid_Driver& gd,
58+
const psi::Psi<T>& psi_ks,
59+
const std::vector<double>& orb_cutoff,
60+
const Parallel_Grid& pgrid,
61+
const int& nspin,
62+
const std::vector<int>& nocc,
63+
const std::vector<int>& nvirt,
64+
const int& nao,
65+
const std::vector<Parallel_2D>& pX,
66+
const Parallel_2D& pc,
67+
const Parallel_Orbitals& pmat,
68+
const bool openshell = false) :
69+
gint_(gint), ucell_(ucell), kv_(kv), gd_(gd), psi_ks_(psi_ks),
70+
orb_cutoff_(orb_cutoff), pgrid_(pgrid), nspin_(nspin), nk_(kv.get_nks() / nspin),
71+
nocc_(nocc), nvirt_(nvirt), nao_(nao), pX_(pX), pc_(pc), pmat_(pmat),
72+
openshell_(openshell), spintype_(openshell ? std::vector<std::string>({ "up", "down" }) : std::vector<std::string>({ "singlet", "triplet" }))
73+
{
74+
}
75+
76+
/// @brief calculate the electron density from the density matrix in 2d-block distribution
77+
void cal_eh_density_single_state(const T* const X_istate, const int ispin, double** density)
78+
{
79+
ModuleBase::TITLE("LR_Density", "cal_eh_density_single_state");
80+
ModuleBase::GlobalFunc::ZEROS(density[0], this->pgrid_.get_nrxx());
81+
// 1. calculate the density matrix in AO basis
82+
auto c_spin = LR_Util::get_psi_spin(psi_ks_, ispin,nk_);
83+
const std::vector<ct::Tensor> dm_diff_k =
84+
cal_dm_diff_pblas<T>(X_istate, pX_[ispin], c_spin, pc_, nao_, nocc_[ispin], nvirt_[ispin], pmat_);
85+
// 2. calculate DM(R)
86+
elecstate::DensityMatrix<T, T> dm_diff=
87+
LR_Util::build_dm_from_dmk<T, T>(dm_diff_k,
88+
this->pmat_, this->nk_, this->kv_.kvec_d, this->ucell_, this->gd_, this->orb_cutoff_, /*symmetrize=*/false);
89+
// 3. calculate electron density from DM(R)
90+
this->dm_to_density(dm_diff, density);
91+
}
92+
93+
void write_density_single_state(const double* const* const density, const std::string& filepath)
94+
{
95+
ModuleIO::write_vdata_palgrid(pgrid_, density[0], 0, 1, 0, filepath, 0.0, &ucell_, PARAM.inp.out_chg[1], 0);
96+
}
97+
98+
void output_eh_density_all_states(const T* const X, const int ispin, const int nstate)
99+
{
100+
ModuleBase::TITLE("LR_Density", "cal_eh_density_all_states");
101+
const int offset_per_state = openshell_ ?
102+
this->nk_ * (this->pX_[0].get_local_size() + this->pX_[1].get_local_size())
103+
: this->nk_ * this->pX_[ispin].get_local_size();
104+
double** density;
105+
LR_Util::_allocate_2order_nested_ptr(density, 1, pgrid_.get_nrxx());
106+
for (int istate = 0;istate < nstate;++istate)
107+
{
108+
int offset = istate * offset_per_state;
109+
if (openshell_)
110+
offset += ispin * this->pX_[0].get_local_size();
111+
this->cal_eh_density_single_state(X + offset, ispin, density);
112+
const std::string filepath = PARAM.globalv.global_out_dir + "LR_e-h_density_" + spintype_[ispin] + "_" + std::to_string(istate + 1) + ".cube";
113+
this->write_density_single_state(density, filepath);
114+
}
115+
LR_Util::_deallocate_2order_nested_ptr(density, 1);
116+
}
117+
};
118+
119+
}

‎source/module_lr/lr_spectrum.cpp‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -100,10 +100,10 @@ ModuleBase::Vector3<std::complex<double>> LR::LR_Spectrum<std::complex<double>>:
100100
// 2. transition density
101101
double** rho_trans_real;
102102
double** rho_trans_imag;
103-
LR_Util::_allocate_2order_nested_ptr(rho_trans_real, 1, this->rho_basis.nrxx);
104-
LR_Util::_allocate_2order_nested_ptr(rho_trans_imag, 1, this->rho_basis.nrxx);
103+
LR_Util::_allocate_2order_nested_ptr(rho_trans_real, this->nspin_x, this->rho_basis.nrxx);
104+
LR_Util::_allocate_2order_nested_ptr(rho_trans_imag, this->nspin_x, this->rho_basis.nrxx);
105105

106-
elecstate::DensityMatrix<std::complex<double>, double> DM_trans_real_imag(&this->pmat, 1, this->kv.kvec_d, this->nk);
106+
elecstate::DensityMatrix<std::complex<double>, double> DM_trans_real_imag(&this->pmat, this->nspin_x, this->kv.kvec_d, this->nk);
107107
LR_Util::initialize_DMR(DM_trans_real_imag, this->pmat, this->ucell, this->gd_, this->orb_cutoff_);
108108

109109
// real part
@@ -128,10 +128,10 @@ ModuleBase::Vector3<std::complex<double>> LR::LR_Spectrum<std::complex<double>>:
128128
rd -= ModuleBase::Vector3<double>(0.5, 0.5, 0.5); //shift to the center of the grid (need ?)
129129
ModuleBase::Vector3<double> rc = rd * ucell.latvec * ucell.lat0; // real coordinate
130130
ModuleBase::Vector3<std::complex<double>> rc_complex(rc.x, rc.y, rc.z);
131-
trans_dipole += rc_complex * std::complex<double>(rho_trans_real[0][ir], rho_trans_imag[0][ir]);
131+
trans_dipole += rc_complex * std::complex<double>(rho_trans_real[is][ir], rho_trans_imag[is][ir]);
132132
}
133-
LR_Util::_deallocate_2order_nested_ptr(rho_trans_real, 1);
134-
LR_Util::_deallocate_2order_nested_ptr(rho_trans_imag, 1);
133+
LR_Util::_deallocate_2order_nested_ptr(rho_trans_real, this->nspin_x);
134+
LR_Util::_deallocate_2order_nested_ptr(rho_trans_imag, this->nspin_x);
135135
}
136136
trans_dipole *= (ucell.omega / static_cast<double>(gint->get_ncxyz())); // dv
137137
trans_dipole *= static_cast<double>(this->nk); // nk is divided inside DM_trans, now recover it

0 commit comments

Comments
 (0)