Skip to content

Commit d11867d

Browse files
maki49maki49
committed
Fix some factor bugs in ceZc, K_ii terms (visible in DZP)
1 parent f183cab commit d11867d

5 files changed

Lines changed: 36 additions & 24 deletions

File tree

‎source/source_lcao/module_lr/Grad/CVCX/CVCX_parallel.cpp‎

Lines changed: 8 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -63,9 +63,9 @@ namespace LR
6363

6464
//AX_istate=[cX^T]^T[c^TV]^T (nvirt major)
6565
pdgemm_(&trans, &trans, &nvirt, &nocc, &naos,
66-
&one, cx.data<double>(), &i1, &i1, pcx.desc,
66+
&factor, cx.data<double>(), &i1, &i1, pcx.desc,
6767
cv.data<double>(), &i1, &i1, pcv.desc,
68-
add_on ? &factor : &zero, AX_istate + start, &i1, &i1, px.desc);
68+
add_on ? &one : &zero, AX_istate + start, &i1, &i1, px.desc);
6969
}
7070
}
7171

@@ -127,9 +127,9 @@ namespace LR
127127

128128
//AX_istate=[cX^T]^T[c^TV]^T (nvirt major)
129129
pzgemm_(&trans, &trans, &nvirt, &nocc, &naos,
130-
&one, cx.data<std::complex<double>>(), &i1, &i1, pcx.desc,
130+
&factor, cx.data<std::complex<double>>(), &i1, &i1, pcx.desc,
131131
cv.data<std::complex<double>>(), &i1, &i1, pcv.desc,
132-
add_on ? &factor : &zero, AX_istate + start, &i1, &i1, px.desc);
132+
add_on ? &one : &zero, AX_istate + start, &i1, &i1, px.desc);
133133
}
134134
}
135135

@@ -191,9 +191,9 @@ namespace LR
191191

192192
//AX_istate=[VC]^T[X^TC^T]^T (nvirt major)
193193
pdgemm_(&trans, &trans, &nvirt, &nocc, &naos,
194-
&one, cv.data<double>(), &i1, &i1, pcv.desc,
194+
&factor, cv.data<double>(), &i1, &i1, pcv.desc,
195195
cx.data<double>(), &i1, &i1, pcx.desc,
196-
add_on ? &factor : &zero, AX_istate + start, &i1, &i1, px.desc);
196+
add_on ? &one : &zero, AX_istate + start, &i1, &i1, px.desc);
197197
}
198198
}
199199

@@ -255,9 +255,9 @@ namespace LR
255255

256256
//AX_istate=[VC]^T[X^TC^T]^T (nvirt major)
257257
pzgemm_(&trans, &trans, &nvirt, &nocc, &naos,
258-
&one, cv.data<std::complex<double>>(), &i1, &i1, pcv.desc,
258+
&factor, cv.data<std::complex<double>>(), &i1, &i1, pcv.desc,
259259
cx.data<std::complex<double>>(), &i1, &i1, pcx.desc,
260-
add_on ? &factor : &zero, AX_istate + start, &i1, &i1, px.desc);
260+
add_on ? &one : &zero, AX_istate + start, &i1, &i1, px.desc);
261261
}
262262
}
263263
}

‎source/source_lcao/module_lr/Grad/CVCX/CVCX_serial.cpp‎

Lines changed: 8 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -99,8 +99,8 @@ namespace LR
9999
cx.data<double>(), &naos);
100100

101101
//AX_istate=[cX^T]^T[c^TV]^T (nvirt major)
102-
dgemm_(&trans, &trans, &nvirt, &nocc, &naos, &one,
103-
cx.data<double>(), &naos, cv.data<double>(), &nocc, add_on ? &factor : &zero,
102+
dgemm_(&trans, &trans, &nvirt, &nocc, &naos, &factor,
103+
cx.data<double>(), &naos, cv.data<double>(), &nocc, add_on ? &one : &zero,
104104
AX_istate + start, &nvirt);
105105
}
106106
}
@@ -145,8 +145,8 @@ namespace LR
145145
cx.data<std::complex<double>>(), &naos);
146146

147147
//AX_istate=[cX^T]^T[c^TV]^T (nvirt major)
148-
zgemm_(&trans, &trans, &nvirt, &nocc, &naos, &one,
149-
cx.data<std::complex<double>>(), &naos, cv.data<std::complex<double>>(), &nocc, add_on ? &factor : &zero,
148+
zgemm_(&trans, &trans, &nvirt, &nocc, &naos, &factor,
149+
cx.data<std::complex<double>>(), &naos, cv.data<std::complex<double>>(), &nocc, add_on ? &one : &zero,
150150
AX_istate + start, &nvirt);
151151
}
152152
}
@@ -247,8 +247,8 @@ namespace LR
247247
cx.data<double>(), &nocc);
248248

249249
//AX_istate=[VC]^T[X^TC^T]^T (nvirt major)
250-
dgemm_(&trans, &trans, &nvirt, &nocc, &naos, &one,
251-
cv.data<double>(), &naos, cx.data<double>(), &nocc, add_on ? &factor : &zero,
250+
dgemm_(&trans, &trans, &nvirt, &nocc, &naos, &factor,
251+
cv.data<double>(), &naos, cx.data<double>(), &nocc, add_on ? &one : &zero,
252252
AX_istate + start, &nvirt);
253253
}
254254
}
@@ -293,8 +293,8 @@ namespace LR
293293
cx.data<std::complex<double>>(), &nocc);
294294

295295
//AX_istate=[VC]^T[X^TC^T]^T (nvirt major)
296-
zgemm_(&trans, &trans, &nvirt, &nocc, &naos, &one,
297-
cv.data<std::complex<double>>(), &naos, cx.data<std::complex<double>>(), &nocc, add_on ? &factor : &zero,
296+
zgemm_(&trans, &trans, &nvirt, &nocc, &naos, &factor,
297+
cv.data<std::complex<double>>(), &naos, cx.data<std::complex<double>>(), &nocc, add_on ? &one : &zero,
298298
AX_istate+start, &nvirt);
299299
}
300300
}

‎source/source_lcao/module_lr/Grad/multipliers/cal_edm_from_multipliers.h‎

Lines changed: 12 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -103,14 +103,21 @@ namespace LR
103103
// 1. c * W * c
104104
const std::vector<ct::Tensor> cWc = cal_dm_trans_pblas(W, p_occ_occ, c, pc, naos, nocc, nvirt, pmat, (T)1., LR_Util::MO_TYPE::OO);
105105

106-
// 2. edm of Z : $\sum_i \sum_a c_a epsilon_i Z_{ai} c_i
107-
std::vector<T> epsi_Z(px.get_local_size());
108-
multiply_eig_onto_vec(Z, eig_ks, px, epsi_Z.data());
109-
std::vector<ct::Tensor> cZc = cal_dm_trans_pblas(Z, px, c, pc, naos, nocc, nvirt, pmat);
106+
// 2. edm of Z : $\sum_i \sum_a c_{\mu a} \epsilon_i Z_{ai} c_{\nu i}$
107+
std::vector<T> epsi_Z(px.get_local_size() * c.get_nk());
108+
for (int ik = 0;ik < c.get_nk();++ik)
109+
{
110+
multiply_eig_onto_vec(Z + ik * px.get_local_size(), eig_ks + ik * (nocc + nvirt),
111+
px, epsi_Z.data() + ik * px.get_local_size());
112+
}
113+
std::vector<ct::Tensor> cZc = cal_dm_trans_pblas(epsi_Z.data(), px, c, pc, naos, nocc, nvirt, pmat);
110114
std::for_each(cZc.begin(), cZc.end(), [&](ct::Tensor& s) { LR_Util::matsym(s.data<T>(), naos, pmat); });
111115

112116
//3. c * K_cvcx * c
113-
std::vector<ct::Tensor> cKc = cal_dm_trans_pblas(K_cvcx, px, c, pc, naos, nocc, nvirt, pmat, (T)2.0);
117+
// $\sum_{kl}K_{kl}[D^X](c_{\kappa k}X_{\lambda l}+X_{\kappa k}c_{\lambda l})$.
118+
// `K_cvcx` already carries the factor 2 of $W^X_{ij}=2K_{ij}[D^X]$ (see `op_K_cvcx` above),
119+
// `matsym` then supplies the 1/2 that turns $2\,X_\kappa K c_\lambda$ into the symmetric pair above.
120+
std::vector<ct::Tensor> cKc = cal_dm_trans_pblas(K_cvcx, px, c, pc, naos, nocc, nvirt, pmat, (T)1.0);
114121
std::for_each(cKc.begin(), cKc.end(), [&](ct::Tensor& s) { LR_Util::matsym(s.data<T>(), naos, pmat); });
115122

116123
// 4. $\sum_i (\Omega + \epsilon_i) \sum_{ab} C_{\mu a} X_{ia} C_{\nu b} X_{ib}$

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

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -50,10 +50,13 @@ namespace LR
5050
// Hessian (A+B) with GS XC kernel
5151
// 1. diag term in A
5252
this->ops = new OperatorLRDiag<T>(eig_ks.c, pX[0], kv.get_nks() / nspin, nocc[0], nvirt[0]);
53-
// 2. $H_{ia}[D^Z]$, equals to $2K_{ab}[D^Z]$ when $D^Z$ is symmetrized
53+
// 2. $H_{ia}[D^Z]$, equals to $2K_{ab}[D^Z]$ when $D^Z$ is symmetrized.
54+
// Factor 4 (not 2): the singlet kernel is $K^S_\text{Hxc}=2$`pot_hxc_gs` (not doubled),
55+
// while `pot` is the already-doubled singlet potential, so $H^S=2K^S$ here needs 4.
56+
// The EXX line below is already $2\alpha$ and is consistent.
5457
hamilt::Operator<T>* op_hz = new OperatorLRHxc<T>(nspin, naos, nocc, nvirt, psi_ks,
5558
*this->DM_trans, pot_hxc_gs, ucell, orb_cutoff, gd, kv, pX, pc, pmat,
56-
{ 0 }, 2.0, ATYPE::CC_vo);
59+
{ 0 }, 4.0, ATYPE::CC_vo);
5760
this->ops->add(op_hz);
5861
#ifdef __EXX
5962
if (exx_kernel_list().count(PARAM.inp.dft_functional))

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

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -70,9 +70,11 @@ namespace LR
7070
#endif
7171
// 2. $H_{ia}[T]$, equals to $2K_{ab}[T]$ when $T$ is symmetrized
7272
// kernel: ground state
73+
// Factor -4 (not -2), for the same reason as in `hamilt_zeq_left.h`:
74+
// $K^S_\text{Hxc}=2$`pot_hxc_gs`, so $H^S=2K^S$ needs 4.
7375
hamilt::Operator<T>* op_ht = new OperatorLRHxc<T>(nspin, naos, nocc, nvirt, psi_ks,
7476
*this->DM_diff, pot_hxc_gs, ucell, orb_cutoff, gd, kv, pX, pc, pmat,
75-
{ 0 }, T(-2.0), ATYPE::CC_vo, hamilt::calculation_type::lr_dmdiff_hxc);
77+
{ 0 }, T(-4.0), ATYPE::CC_vo, hamilt::calculation_type::lr_dmdiff_hxc);
7678
this->ops->add(op_ht);
7779
#ifdef __EXX
7880
if (exx_kernel_list().count(PARAM.inp.dft_functional))

0 commit comments

Comments
 (0)