Skip to content

Commit 375762d

Browse files
maki49maki49
committed
calculate W multiplier (multi-k not finished)
1 parent 92aa141 commit 375762d

File tree

7 files changed

+169
-26
lines changed

7 files changed

+169
-26
lines changed

‎source/module_lr/Grad/multipliers/cal_multiplier_w_from_z.h‎

Lines changed: 132 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,132 @@
1+
#pragma once
2+
#include "module_hamilt_general/hamilt.h"
3+
#include "module_elecstate/module_dm/density_matrix.h"
4+
#include "module_lr/Grad/xc/pot_grad_xc.h"
5+
#include "module_lr/potentials/pot_hxc_lrtd.h"
6+
#include "module_lr/operator_casida/operator_lr_hxc.h"
7+
#include "module_basis/module_ao/parallel_orbitals.h"
8+
namespace LR
9+
{
10+
// $\sum_a (\Omega - \epsilon_a)\delta_{ij}$
11+
// !!! Not yet support multi-k
12+
template<typename T>
13+
void add_sum_eig_minus_evirt(T* const inout,
14+
const double eig_ext_istate,
15+
const double* const eig_ks,
16+
const int nocc,
17+
const int nvirt,
18+
const Parallel_2D& p_occ_occ)
19+
{
20+
for (int i = 0;i < nocc;++i)
21+
{
22+
if (p_occ_occ.in_this_processor(i, i)) // diagonal elements in this processor
23+
{
24+
const int lrow = p_occ_occ.global2local_row(i);
25+
const int lcol = p_occ_occ.global2local_col(i);
26+
for (int a = 0;a < nvirt;++i)
27+
{
28+
inout[lcol * p_occ_occ.get_row_size() + lrow] += (eig_ext_istate - eig_ks[/**ik, */nocc + i]);
29+
}
30+
}
31+
}
32+
}
33+
34+
template<typename T, typename TGint>
35+
void cal_W_from_Z(T* const W,
36+
const T* const Z,
37+
const T* const X,
38+
const double* const eig,
39+
const int& nspin,
40+
const int& naos,
41+
const std::vector<int>& nocc,
42+
const std::vector<int>& nvirt,
43+
const UnitCell& ucell_in,
44+
const std::vector<double>& orb_cutoff,
45+
const Grid_Driver& gd_in,
46+
const psi::Psi<T>& psi_ks_in,
47+
const ModuleBase::matrix& eig_ks,
48+
#ifdef __EXX
49+
std::weak_ptr<Exx_LRI<T>> exx_lri,
50+
const double& exx_alpha,
51+
#endif
52+
TGint* gint_in,
53+
std::weak_ptr<PotHxcLR> pot,
54+
const K_Vectors& kv,
55+
const std::vector<Parallel_2D>& px,
56+
const Parallel_2D& pc,
57+
const Parallel_2D& p_occ_occ, // < for W
58+
const Parallel_Orbitals& pmat,
59+
const std::string& spin_type = "singlet")
60+
{
61+
ModuleBase::TITLE("Z_vector_R", "Z_vector_R");
62+
using ATYPE = typename OperatorLRHxc<T>::MO_TO_AO_TYPE;
63+
const int nk = kv.get_nks() / nspin;
64+
// allocate memory for DMs
65+
elecstate::DensityMatrix<T, T> DM_trans(&pmat, 1, kv.kvec_d, nk); //DX
66+
elecstate::DensityMatrix<T, T> DM_diff_relaxed(&pmat, 1, kv.kvec_d, nk); //T+DZ
67+
68+
/// operators
69+
// 1. 0.5$H_{ia}[T]$, equals to $K_{ab}[T]$ when $T$ is symmetrized
70+
OperatorLRHxc<T> op_ht(nspin, naos, nocc, nvirt, psi_ks_in,
71+
DM_diff_relaxed, gint_in, pot, ucell_in, orb_cutoff, gd_in, kv, px, pc, pmat,
72+
{ 0 }, T(1.0), ATYPE::CC_oo);
73+
// 2. $2\sum_{jb,kc} g^{xc}_{ia, jb, kc}X_{jb}X_{kc}$
74+
PotGradXCLR pot_grad(pot.lock()->xc_kernel_components, pot.lock()->rho_basis, ucell_in, pot.lock()->nrxx);
75+
OperatorLRHxc<T> op_gxc(nspin, naos, nocc, nvirt, psi_ks_in,
76+
DM_trans, gint_in, pot_grad, ucell_in, orb_cutoff, gd_in, kv, px, pc, pmat,
77+
{ 0 }, T(-2.0), ATYPE::CC_oo);
78+
#ifdef __EXX
79+
// add EXX operators here
80+
#endif
81+
auto cal_dm_trans = [&](const int is, const T* const x_ptr)->void //DX
82+
{
83+
const auto psi_ks_is = LR_Util::get_psi_spin(psi_ks_in, is, nk);
84+
#ifdef __MPI
85+
std::vector<ct::Tensor> dm_trans_2d = cal_dm_trans_pblas(x_ptr, px[is], psi_ks_is, pc, naos, nocc[is], nvirt[is], pmat);
86+
for (auto& t : dm_trans_2d) LR_Util::matsym(t.data<T>(), naos, pmat);
87+
#else
88+
std::vector<ct::Tensor> dm_trans_2d = cal_dm_trans_blas(x_ptr, psi_ks_is, nocc[is], nvirt[is]);
89+
for (auto& t : dm_trans_2d) LR_Util::matsym(t.data<T>(), naos);
90+
#endif
91+
for (int ik = 0;ik < nk;++ik) { DM_trans->set_DMK_pointer(ik, dm_trans_2d[ik].data<T>()); }
92+
};
93+
94+
auto cal_dm_diff_relaxed = [&](const int& is, const T* const x_ptr, const T* const z_ptr)->void // T+DZ
95+
{
96+
const auto psi_ks_is = LR_Util::get_psi_spin(psi_ks_in, is, nk);
97+
#ifdef __MPI
98+
std::vector<ct::Tensor> z_2d = cal_dm_trans_pblas(z_ptr, px[is], psi_ks_is, pc, naos, nocc[is], nvirt[is], pmat);
99+
for (auto& t : dm_trans_2d) LR_Util::matsym(t.data<T>(), naos, pmat);
100+
#else
101+
std::vector<ct::Tensor> z_2d = cal_dm_trans_blas(z_ptr, psi_ks_is, nocc[is], nvirt[is]);
102+
for (auto& t : dm_trans_2d) LR_Util::matsym(t.data<T>(), naos);
103+
#endif
104+
105+
#ifdef __MPI
106+
std::vector<ct::Tensor> dm_diff_2d = cal_dm_diff_pblas(x_ptr, px[is], psi_ks_is, pc, naos, nocc[is], nvirt[is], pmat);
107+
#else
108+
std::vector<ct::Tensor> dm_diff_2d = cal_dm_diff_blas(x_ptr, psi_ks_is, naos, nocc[is], nvirt[is]);
109+
#endif
110+
for (int ik = 0;ik < nk;++ik)
111+
{
112+
dm_diff_2d[i] += z_2d[i];
113+
DM_diff->set_DMK_pointer(ik, dm_diff_2d[ik].data<T>());
114+
}
115+
};
116+
117+
// act the operators onto each state
118+
const int ld_vo = nk * px[0].get_local_size();
119+
const int ld_oo = nk * p_occ_occ.get_local_size();
120+
for (int ib = 0;ib < nstates;++ib)
121+
{
122+
const int offset_vo = ib * ld_vo;
123+
const int offset_oo = ib * ld_oo;
124+
cal_dm_trans(0, X + offset_vo, Z + offset_vo); // transition density matrix DX
125+
cal_dm_diff_relaxed(0, X + offset_vo, Z + offset_vo); // relaxed difference density matrix T+DZ
126+
// the 3 terms
127+
op_ht.act(/*nband=*/1, ld_vo, /*npol=*/1, X + offset_vo, W + offset_oo);
128+
op_gxc.act(/*nband=*/1, ld_vo, /*npol=*/1, X + offset_vo, W + offset_oo);
129+
add_sum_eig_minus_evirt(W + offset_oo, eig[ib], eig_ks.data(), nocc[0], nvirt[0], p_occ_occ);
130+
}
131+
}
132+
}

‎source/module_lr/Grad/multipliers/hamilt_zeq_left.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -10,7 +10,7 @@ namespace LR
1010
template<typename T>
1111
class Z_vector_L : public HamiltLR<T>
1212
{
13-
using ATYPE = typename OperatorLRHxc<T>::AX_TYPE;
13+
using ATYPE = typename OperatorLRHxc<T>::MO_TO_AO_TYPE;
1414
public:
1515
template<typename TGint>
1616
Z_vector_L(const std::string& xc_kernel,
@@ -47,7 +47,7 @@ namespace LR
4747
// 2. $H_{ia}[T]$, equals to $2K_{ab}[T]$ when $T$ is symmetrized
4848
OperatorLRHxc<T>* op_hz = new OperatorLRHxc<T>(nspin, naos, nocc, nvirt, psi_ks_in,
4949
this->DM_trans, gint_in, pot_in, ucell_in, orb_cutoff, gd_in, kv_in, pX_in, pc_in, pmat_in,
50-
{ 0 }, 2.0, ATYPE::CC);
50+
{ 0 }, 2.0, ATYPE::CC_vo);
5151
this->ops->add(op_hz);
5252
#ifdef __EXX
5353
// add EXX operators here

‎source/module_lr/Grad/multipliers/hamilt_zeq_right.h‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -10,7 +10,7 @@ namespace LR
1010
template<typename T>
1111
class Z_vector_R : public HamiltLR<T>
1212
{
13-
using ATYPE = typename OperatorLRHxc<T>::AX_TYPE;
13+
using ATYPE = typename OperatorLRHxc<T>::MO_TO_AO_TYPE;
1414
public:
1515
template<typename TGint>
1616
Z_vector_R(const std::string& xc_kernel,
@@ -50,13 +50,13 @@ namespace LR
5050
// 2. $H_{ia}[T]$, equals to $2K_{ab}[T]$ when $T$ is symmetrized
5151
OperatorLRHxc<T>* op_ht = new OperatorLRHxc<T>(nspin, naos, nocc, nvirt, psi_ks_in,
5252
this->DM_diff, gint_in, pot_in, ucell_in, orb_cutoff, gd_in, kv_in, pX_in, pc_in, pmat_in,
53-
{ 0 }, T(-2.0), ATYPE::CC);
53+
{ 0 }, T(-2.0), ATYPE::CC_vo);
5454
this->ops->add(op_ht);
5555
// 3. $2\sum_{jb,kc} g^{xc}_{ia, jb, kc}X_{jb}X_{kc}$
5656
this->pot_grad = std::make_shared<PotGradXCLR>(pot_in.lock()->xc_kernel_components, pot_in.lock()->rho_basis, ucell_in, pot_in.lock()->nrxx);
5757
OperatorLRHxc<T>* op_gxc = new OperatorLRHxc<T>(nspin, naos, nocc, nvirt, psi_ks_in,
5858
this->DM_trans, gint_in, this->pot_grad, ucell_in, orb_cutoff, gd_in, kv_in, pX_in, pc_in, pmat_in,
59-
{ 0 }, T(-2.0), ATYPE::CC);
59+
{ 0 }, T(-2.0), ATYPE::CC_vo);
6060
this->ops->add(op_gxc);
6161
#ifdef __EXX
6262
// add EXX operators here

‎source/module_lr/Grad/operator_lr_hxc_common/operator_lr_hxc_common.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -105,7 +105,7 @@ namespace LR
105105
// 5. [AX]^{Hxc}_{ai}=\sum_{\mu,\nu}c^*_{a,\mu,}V^{Hxc}_{\mu,\nu}c_{\nu,i}
106106
switch (this->dm_pq)
107107
{
108-
case AX_TYPE::CXC:
108+
case MO_TO_AO_TYPE::CXC:
109109
#ifdef __MPI
110110
CVCX_virt_pblas(v_hxc_2d, this->pmat, *this->psi_ks, this->pc, psi_in_bfirst, this->pX,
111111
this->naos, this->nocc, this->nvirt, psi_out_bfirst, /*add_on=*/true, factor);
@@ -116,7 +116,7 @@ namespace LR
116116
CVCX_occ_blas(v_hxc_2d, *this->psi_ks, psi_in_bfirst, this->naos, this->nocc, this->nvirt, psi_out_bfirst, /*add_on=*/true, -factor);
117117
#endif
118118
break;
119-
case AX_TYPE::CC: // C_onebase_ai
119+
case MO_TO_AO_TYPE::CC_vo: // C_onebase_ai
120120
#ifdef __MPI
121121
cal_AX_pblas(v_hxc_2d, this->pmat, *this->psi_ks, this->pc, this->naos, this->nocc, this->nvirt, this->pX, psi_out_bfirst);
122122
#else

‎source/module_lr/Grad/operator_lr_hxc_common/operator_lr_hxc_common.h‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -17,7 +17,7 @@ namespace LR
1717
/// @brief AX type:
1818
/// CC: C_v^* C_o^T;
1919
/// CXC: C_v^* X^* C_v^T- C_o^* X^* C_o^T
20-
enum class AX_TYPE { CC, CXC };
20+
enum class MO_TO_AO_TYPE { CC, CXC };
2121

2222
//when nspin=2, nks is 2 times of real number of k-points. else (nspin=1 or 4), nks is the real number of k-points
2323
OperatorLRHxcCommon(const int& nspin,
@@ -36,7 +36,7 @@ namespace LR
3636
Parallel_2D* pc_in,
3737
Parallel_Orbitals* pmat_in,
3838
DM_TYPE dm_rs_in, ///< rs: the labels to be contracted with DM
39-
AX_TYPE dm_pq_in, ///< pq: the labels of the final result
39+
MO_TO_AO_TYPE dm_pq_in, ///< pq: the labels of the final result
4040
const double factor_in)
4141
: dm_pq(dm_pq_in), dm_rs(dm_rs_in), factor(factor_in),
4242
OperatorLRHxc<T, Device>(nspin, naos, nocc, nvirt, psi_ks_in, DM_trans_in,
@@ -55,7 +55,7 @@ namespace LR
5555
void act_to_bfirst(const psi::Psi<T>& psi_in_bfirst, psi::Psi<T>& psi_out_bfirst, const int nbands) const;
5656

5757
DM_TYPE dm_rs = DM_TYPE::X;
58-
AX_TYPE dm_pq = AX_TYPE::CC;
58+
MO_TO_AO_TYPE dm_pq = MO_TO_AO_TYPE::CC_vo;
5959
const T factor = (T)1.0;
6060
};
6161
}

‎source/module_lr/operator_casida/operator_lr_hxc.cpp‎

Lines changed: 19 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -42,16 +42,24 @@ namespace LR
4242
// for (int ik = 0;ik < nk;++ik)
4343
// LR_Util::print_tensor<T>(v_hxc_2d[ik], "4.V(k)[ik=" + std::to_string(ik) + "]", &this->pmat);
4444

45+
// 5. AO to MO transformation
4546
switch (this->dm_pq_)
4647
{
47-
case AX_TYPE::CC: // C_onebase_ai
48+
case MO_TO_AO_TYPE::CC_vo: //[AX]^{Hxc}_{ai}=\sum_{\mu,\nu}c^*_{a,\mu,}V^{Hxc}_{\mu,\nu}c_{\nu,i}
4849
#ifdef __MPI
49-
cal_AX_pblas(v_hxc_2d, this->pmat, psil_ks, this->pc, this->naos, this->nocc[sl], this->nvirt[sl], this->pX[sl], hpsi);
50+
ao_to_mo_pblas(v_hxc_2d, this->pmat, psil_ks, this->pc, this->naos, this->nocc[sl], this->nvirt[sl], this->pX[sl], hpsi);
5051
#else
51-
cal_AX_blas(v_hxc_2d, psil_ks, this->naos, this->nocc[sl], this->nvirt[sl], hpsi);
52+
ao_to_mo_blas(v_hxc_2d, psil_ks, this->nocc[sl], this->nvirt[sl], hpsi);
5253
#endif
5354
break;
54-
case AX_TYPE::CXC:
55+
case MO_TO_AO_TYPE::CC_oo: //[AX]^{Hxc}_{ai}=\sum_{\mu,\nu}c^*_{a,\mu,}V^{Hxc}_{\mu,\nu}c_{\nu,i}
56+
#ifdef __MPI
57+
ao_to_mo_pblas(v_hxc_2d, this->pmat, psil_ks, this->pc, this->naos, this->nocc[sl], this->nvirt[sl], this->pX[sl], hpsi, MO_TYPE::OO);
58+
#else
59+
ao_to_mo_blas(v_hxc_2d, psil_ks, this->nocc[sl], this->nvirt[sl], hpsi, MO_TYPE::OO);
60+
#endif
61+
break;
62+
case MO_TO_AO_TYPE::CXC:
5563
#ifdef __MPI
5664
CVCX_virt_pblas(v_hxc_2d, this->pmat, psil_ks, this->pc, psi_in, this->pX[sl],
5765
this->naos, this->nocc[sl], this->nvirt[sl], hpsi, /*add_on=*/true, this->factor_);
@@ -62,16 +70,17 @@ namespace LR
6270
CVCX_occ_blas(v_hxc_2d, *this->psi_ks, psi_in_bfirst, this->naos, this->nocc, this->nvirt, hpsi, /*add_on=*/true, -this->factor_);
6371
#endif
6472
break;
73+
case MO_TO_AO_TYPE::CXC_o:
74+
#ifdef __MPI
75+
CVCX_occ_pblas(v_hxc_2d, this->pmat, psil_ks, this->pc, psi_in, this->pX[sl],
76+
this->naos, this->nocc[sl], this->nvirt[sl], hpsi, /*add_on=*/true, this->factor_);
77+
#else
78+
CVCX_occ_blas(v_hxc_2d, *this->psi_ks, psi_in_bfirst, this->naos, this->nocc, this->nvirt, hpsi, /*add_on=*/true, this->factor_);
79+
#endif
6580
default:
6681
throw std::runtime_error("Unknown DM_TYPE");
6782
break;
6883
}
69-
// 5. [AX]^{Hxc}_{ai}=\sum_{\mu,\nu}c^*_{a,\mu,}V^{Hxc}_{\mu,\nu}c_{\nu,i}
70-
#ifdef __MPI
71-
ao_to_mo_pblas(v_hxc_2d, this->pmat, psil_ks, this->pc, naos, nocc[sl], nvirt[sl], this->pX[sl], hpsi);
72-
#else
73-
ao_to_mo_blas(v_hxc_2d, psil_ks, nocc[sl], nvirt[sl], hpsi);
74-
#endif
7584
}
7685

7786

‎source/module_lr/operator_casida/operator_lr_hxc.h‎

Lines changed: 8 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -14,10 +14,12 @@ namespace LR
1414
class OperatorLRHxc : public hamilt::Operator<T, Device>
1515
{
1616
public:
17-
/// @brief AX type:
18-
/// CC: C_v^* C_o^T;
19-
/// CXC: C_v^* X^* C_v^T- C_o^* X^* C_o^T
20-
enum class AX_TYPE { CC, CXC };
17+
/// @brief type of molecular orbital to atomic orbital transformation:
18+
/// CC_vo: MO = C_v^* AO C_o;
19+
/// CC_oo: MO = C_o^* AO C_o;
20+
/// CXC: MO = C_v^* AO X C_v- C_o^* X^* AO C_o
21+
/// CXC_o: MO = C_o^* X^* AO C_o
22+
enum class MO_TO_AO_TYPE { CC_vo, CC_oo, CXC, CXC_o };
2123

2224
//when nspin=2, nks is 2 times of real number of k-points. else (nspin=1 or 4), nks is the real number of k-points
2325
OperatorLRHxc(const int& nspin,
@@ -37,7 +39,7 @@ namespace LR
3739
const Parallel_Orbitals& pmat_in,
3840
const std::vector<int>& ispin_ks = { 0 },
3941
const T factor_in = (T)1.0,
40-
const AX_TYPE dm_pq_in = AX_TYPE::CC)
42+
const MO_TO_AO_TYPE dm_pq_in = MO_TO_AO_TYPE::CC_vo)
4143
: nspin(nspin), naos(naos), nocc(nocc), nvirt(nvirt), nk(kv_in.get_nks() / nspin), psi_ks(psi_ks_in),
4244
DM_trans(DM_trans_in), gint(gint_in), pot(pot_in), ucell(ucell_in), orb_cutoff_(orb_cutoff), gd(gd_in),
4345
kv(kv_in), pX(pX_in), pc(pc_in), pmat(pmat_in), ispin_ks(ispin_ks),
@@ -96,7 +98,7 @@ namespace LR
9698
std::vector<double> orb_cutoff_;
9799
const Grid_Driver& gd;
98100

99-
AX_TYPE dm_pq_ = AX_TYPE::CC;
101+
MO_TO_AO_TYPE dm_pq_ = MO_TO_AO_TYPE::CC_vo;
100102
const T factor_ = (T)1.0;
101103

102104
/// test

0 commit comments

Comments
 (0)