3
\brief Enter brief description of file here
6
#include <libdpd/dpd.h>
12
namespace psi { namespace cchbar {
14
/* WAbEi_UHF(): Computes all contributions to the AbEi spin case of
15
** the Wabei HBAR matrix elements. The final product is stored in
16
** (Ei,Ab) ordering and is referred to on disk as "WAbEi".
18
** The spin-orbital expression for the Wabei elements is:
20
** Wabei = <ab||ei> - Fme t_mi^ab + t_i^f <ab||ef>
21
** - P(ab) t_m^b <am||ef>t_i^f + 1/2 tau_mn^ab <mn||ef> t_i^f
22
** + 1/2 <mn||ei> tau_mn^ab - P(ab) <mb||ef> t_mi^af
23
** - P(ab) t_m^a { <mb||ei> - t_ni^bf <mn||ef> }
25
** (cf. Gauss and Stanton, JCP 103, 3561-3577 (1995).)
27
** For the AbEi spin case, we evaluate these contractions with two
28
** target orderings, (Ab,Ei) and (Ei,Ab), depending on the term.
29
** After all terms have been evaluated, the (Ab,Ei) terms are sorted
30
** into (Ei,Ab) ordering and both groups arer added together.
38
dpdbuf4 F, W, T2, B, Z, Z1, Z2, D, T, E, C;
42
/** W(Ei,Ab) <--- <Ei|Ab> **/
43
dpd_buf4_init(&F, CC_FINTS, 0, 26, 28, 26, 28, 0, "F <Ai|Bc>");
44
dpd_buf4_copy(&F, CC_HBAR, "WEiAb");
49
/** W(Ei,Ab) <--- - F_ME t_Mi^Ab **/
50
dpd_buf4_init(&T2, CC_TAMPS, 0, 22, 28, 22, 28, 0, "tIjAb");
51
dpd_file2_init(&Fme, CC_OEI, 0, 0, 1, "FME");
52
dpd_buf4_init(&W, CC_HBAR, 0, 26, 28, 26, 28, 0, "WEiAb");
53
dpd_contract244(&Fme, &T2, &W, 0, 0, 0, -1, 1);
55
dpd_file2_close(&Fme);
61
dpd_buf4_init(&W, CC_TMP0, 0, 28, 26, 28, 26, 0, "W'(Ab,Ei)");
62
dpd_buf4_init(&B, CC_BINTS, 0, 28, 28, 28, 28, 0, "B <Ab|Cd>");
63
dpd_file2_init(&T1, CC_OEI, 0, 2, 3, "tia");
64
dpd_contract424(&B, &T1, &W, 3, 1, 0, 1, 0);
71
/** WAbEi <-- - t_m^b <mA|fE> t_i^f - t_M^A <Mb|Ef> t_i^f
72
Evaluate in three steps:
73
(1) Z_mAEi = - <Am|Ef> t_i^f [stored (Am,Ei)]
74
(2) Z_MbEi = <Mb|Ef> t_i^f [stored (Mb,Ei)]
75
(3) WAbEi <-- t_m^b Z_mAEi - t_M^A Z_MbEi
76
Store targets in: W(Ei,Ab) and W'(Ab,Ei)
79
/** Z(Am,Ei) <-- - <Am|Ef> t_i^f **/
80
dpd_buf4_init(&Z, CC_TMP0, 0, 26, 26, 26, 26, 0, "Z(Am,Ei)");
81
dpd_buf4_init(&F, CC_FINTS, 0, 26, 28, 26, 28, 0, "F <Ai|Bc>");
82
dpd_file2_init(&T1, CC_OEI, 0, 2, 3, "tia");
83
dpd_contract424(&F, &T1, &Z, 3, 1, 0, -1, 0);
88
/** t_m^b Z(Am,Ei) --> W(Ei,Ab) **/
89
dpd_buf4_init(&W, CC_HBAR, 0, 26, 28, 26, 28, 0, "WEiAb");
90
dpd_buf4_init(&Z, CC_TMP0, 0, 26, 26, 26, 26, 0, "Z(Am,Ei)");
91
dpd_file2_init(&T1, CC_OEI, 0, 2, 3, "tia");
92
dpd_contract424(&Z, &T1, &W, 1, 0, 0, 1, 1);
97
/** Z(Mb,Ei) <-- <Mb|Ef> t_i^f **/
98
dpd_buf4_init(&Z, CC_TMP0, 0, 24, 26, 24, 26, 0, "Z(Mb,Ei)");
99
dpd_buf4_init(&F, CC_FINTS, 0, 24, 28, 24, 28, 0, "F <Ia|Bc>");
100
dpd_file2_init(&T1, CC_OEI, 0, 2, 3, "tia");
101
dpd_contract424(&F, &T1, &Z, 3, 1, 0, 1, 0);
102
dpd_file2_close(&T1);
106
/** - t_M^A Z(Mb,Ei) --> W'(Ab,Ei) **/
107
dpd_buf4_init(&W, CC_TMP0, 0, 28, 26, 28, 26, 0, "W'(Ab,Ei)");
108
dpd_buf4_init(&Z, CC_TMP0, 0, 24, 26, 24, 26, 0, "Z(Mb,Ei)");
109
dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
110
dpd_contract244(&T1, &Z, &W, 0, 0, 0, -1, 1);
111
dpd_file2_close(&T1);
117
/** WAbEi <-- tau_Mn^Ab <Mn|Ef> t_i^f
118
Evaluate in two steps:
119
(1) Z_MnEi = <Mn|Ef> t_i^f
120
(2) WAbEi <-- tau_Mn^Ab Z_MnEi
121
Store target in W'(Ab,Ei)
124
/** Z(Mn,Ei) <-- <Mn|Ef> t_i^f **/
125
dpd_buf4_init(&Z, CC_TMP0, 0, 22, 26, 22, 26, 0, "Z(Mn,Ei)");
126
dpd_buf4_init(&D, CC_DINTS, 0, 22, 28, 22, 28, 0, "D <Ij|Ab>");
127
dpd_file2_init(&T1, CC_OEI, 0, 2, 3, "tia");
128
dpd_contract424(&D, &T1, &Z, 3, 1, 0, 1, 0);
129
dpd_file2_close(&T1);
132
/** tau_Mn^Ab Z1(Mn,Ei) --> W'(Ab,Ei) **/
133
dpd_buf4_init(&Z, CC_TMP0, 0, 22, 26, 22, 26, 0, "Z(Mn,Ei)");
134
dpd_buf4_init(&W, CC_TMP0, 0, 28, 26, 28, 26, 0, "W'(Ab,Ei)");
135
dpd_buf4_init(&T2, CC_TAMPS, 0, 22, 28, 22, 28, 0, "tauIjAb");
136
dpd_contract444(&T2, &Z, &W, 1, 1, 1, 1);
143
/** tau_Mn^Ab <Mn|Ei> --> Z(Ab,Ei) **/
144
dpd_buf4_init(&W, CC_TMP0, 0, 28, 26, 28, 26, 0, "W'(Ab,Ei)");
145
dpd_buf4_init(&T2, CC_TAMPS, 0, 22, 28, 22, 28, 0, "tauIjAb");
146
dpd_buf4_init(&E, CC_EINTS, 0, 22, 26, 22, 26, 0, "E <Ij|Ak>");
147
dpd_contract444(&T2, &E, &W, 1, 1, 1, 1);
154
/** WAbEi <-- <bM|fE> t_iM^fA - <AM||EF> t_iM^bF - <Am|Ef> t_im^bf
155
Evaluate in five steps:
156
(1) Sort <bM|fE> to F(bE,Mf) ordering. (** Note that we assume a sort has already
157
been done for <AM||EF> and <Am|Ef> in the WABEI and Wabei terms. **)
158
(2) Z'(bE,iA) = F(bE,Mf) T(iA,Mf)
159
(3) Sort Z'(bE,iA) --> W(Ei,Ab)
160
(4) Z''(AE,ib) = - F(AE,MF) T(ib,MF) - F(AE,mf) T(ib,mf)
161
(5) Sort Z''(AE,ib) --> W(Ei,Ab)
163
NB: The storage for the sorts is expensive and will eventually require out-of-core
167
/** <bM|fE> --> F(bE,Mf) **/
168
dpd_buf4_init(&F, CC_FINTS, 0, 25, 29, 25, 29, 0, "F <aI|bC>");
169
dpd_buf4_sort(&F, CC_FINTS, psqr, 29, 24, "F <aI|bC> (aC,Ib)");
172
/** <bM|fE> t_iM^fA --> Z(bE,iA) **/
173
dpd_buf4_init(&Z, CC_TMP0, 0, 29, 27, 29, 27, 0, "Z(bE,iA)");
174
dpd_buf4_init(&F, CC_FINTS, 0, 29, 24, 29, 24, 0, "F <aI|bC> (aC,Ib)");
175
dpd_buf4_init(&T2, CC_TAMPS, 0, 27, 24, 27, 24, 0, "tiBJa");
176
dpd_contract444(&F, &T2, &Z, 0, 0, -1, 0);
181
/** Z(bE,iA) --> W(Ei,Ab) **/
182
dpd_buf4_init(&Z, CC_TMP0, 0, 29, 27, 29, 27, 0, "Z(bE,iA)");
183
dpd_buf4_sort_axpy(&Z, CC_HBAR, qrsp, 26, 28, "WEiAb", 1);
186
/** Z''(AE,ib) <-- - <AM||EF> t_iM^bF **/
187
dpd_buf4_init(&Z, CC_TMP0, 0, 5, 30, 5, 30, 0, "Z(AE,ib)");
188
dpd_buf4_init(&F, CC_FINTS, 0, 5, 20, 5, 20, 0, "F <AI||BC> (AB,IC)");
189
dpd_buf4_init(&T2, CC_TAMPS, 0, 30, 20, 30, 20, 0, "tiaJB");
190
dpd_contract444(&F, &T2, &Z, 0, 0, 1, 0);
195
/** Z''(AE,ib) <-- -<Am|Ef> t_im^bf **/
196
dpd_buf4_init(&Z, CC_TMP0, 0, 5, 30, 5, 30, 0, "Z(AE,ib)");
197
dpd_buf4_init(&F, CC_FINTS, 0, 5, 30, 5, 30, 0, "F <Ai|Bc> (AB,ic)");
198
dpd_buf4_init(&T2, CC_TAMPS, 0, 30, 30, 30, 30, 0, "tiajb");
199
dpd_contract444(&F, &T2, &Z, 0, 0, 1, 1);
204
/** Z''(AE,ib) --> W(Ei,Ab) **/
205
dpd_buf4_init(&Z, CC_TMP0, 0, 5, 30, 5, 30, 0, "Z(AE,ib)");
206
dpd_buf4_sort_axpy(&Z, CC_HBAR, qrps, 26, 28, "WEiAb", 1);
209
/**** Terms VIII and IX ****/
211
/** WAbEi <-- - t_M^A { <Mb|Ei> + t_iN^bF <MN||EF> + t_in^bf <Mn|Ef> }
212
+ t_m^b {-<mA|iE> + t_iN^fA <mN|fE> }
213
Evaluate in three steps:
214
(1) Z_MbEi = <Mb|Ei> + t_iN^bF <MN||EF> + tin^bf <Mn|Ef> [stored (Mb,Ei)]
215
(2) Z_mAEi = -<mA|iE> + t_iN^fA <mN|fE> [stored (Am,Ei)]
216
(3) WAbEi <-- - t_M^A Z_MbEi + t_m^b Z_mAEi
217
Store targets in W'(Ab,Ei) and W(Ei,AB)
220
/** Z(Mb,Ei) <-- <Mb|Ei> **/
221
dpd_buf4_init(&D, CC_DINTS, 0, 24, 26, 24, 26, 0, "D <Ij|Ab> (Ib,Aj)");
222
dpd_buf4_copy(&D, CC_TMP0, "Z(Mb,Ei)");
225
/** <MN||EF> t_iN^bF --> Z(ME,ib) **/
226
dpd_buf4_init(&Z, CC_TMP0, 0, 20, 30, 20, 30, 0, "Z(ME,ib)");
227
dpd_buf4_init(&D, CC_DINTS, 0, 20, 20, 20, 20, 0, "D <IJ||AB> (IA,JB)");
228
dpd_buf4_init(&T2, CC_TAMPS, 0, 30, 20, 30, 20, 0, "tiaJB");
229
dpd_contract444(&D, &T2, &Z, 0, 0, 1, 0);
234
/** <Mn|Ef> t_in^bf --> Z(ME,ib) **/
235
dpd_buf4_init(&Z, CC_TMP0, 0, 20, 30, 20, 30, 0, "Z(ME,ib)");
236
dpd_buf4_init(&D, CC_DINTS, 0, 20, 30, 20, 30, 0, "D <Ij|Ab> (IA,jb)");
237
dpd_buf4_init(&T2, CC_TAMPS, 0, 30, 30, 30, 30, 0, "tiajb");
238
dpd_contract444(&D, &T2, &Z, 0, 0, 1, 1);
243
/** Z(ME,ib) --> Z(Mb,Ei) **/
244
dpd_buf4_init(&Z, CC_TMP0, 0, 20, 30, 20, 30, 0, "Z(ME,ib)");
245
dpd_buf4_sort_axpy(&Z, CC_TMP0, psqr, 24, 26, "Z(Mb,Ei)", 1);
248
/** W'(Ab,Ei) <-- - t_M^A Z(Mb,Ei) **/
249
dpd_buf4_init(&W, CC_TMP0, 0, 28, 26, 28, 26, 0, "W'(Ab,Ei)");
250
dpd_buf4_init(&Z, CC_TMP0, 0, 24, 26, 24, 26, 0, "Z(Mb,Ei)");
251
dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
252
dpd_contract244(&T1, &Z, &W, 0, 0, 0, -1, 1);
253
dpd_file2_close(&T1);
257
/** Z(Am,Ei) <-- - <mA|iE> **/
258
dpd_buf4_init(&C, CC_CINTS, 0, 26, 26, 26, 26, 0, "C <Ai|Bj>");
259
dpd_buf4_copy(&C, CC_TMP0, "Z(Am,Ei)");
261
dpd_buf4_init(&Z, CC_TMP0, 0, 26, 26, 26, 26, 0, "Z(Am,Ei)");
262
dpd_buf4_scm(&Z, -1);
265
/** Z(mE,iA) <-- t_iN^fA <mN|fE> **/
266
dpd_buf4_init(&Z, CC_TMP0, 0, 27, 27, 27, 27, 0, "Z(mE,iA)");
267
dpd_buf4_init(&D, CC_DINTS, 0, 27, 24, 27, 24, 0, "D <iJ|aB> (iB,Ja)");
268
dpd_buf4_init(&T2, CC_TAMPS, 0, 27, 24, 27, 24, 0, "tiBJa");
269
dpd_contract444(&D, &T2, &Z, 0, 0, 1, 0);
274
/** Z(mE,iA) --> Z(Am,Ei) **/
275
dpd_buf4_init(&Z, CC_TMP0, 0, 27, 27, 27, 27, 0, "Z(mE,iA)");
276
dpd_buf4_sort_axpy(&Z, CC_TMP0, spqr, 26, 26, "Z(Am,Ei)", 1);
279
/** W(Ei,AB) <-- t_m^b Z_mAEi **/
280
dpd_buf4_init(&W, CC_HBAR, 0, 26, 28, 26, 28, 0, "WEiAb");
281
dpd_buf4_init(&Z, CC_TMP0, 0, 26, 26, 26, 26, 0, "Z(Am,Ei)");
282
dpd_file2_init(&T1, CC_OEI, 0, 2, 3, "tia");
283
dpd_contract424(&Z, &T1, &W, 1, 0, 0, 1, 1);
284
dpd_file2_close(&T1);
288
/**** Combine accumulated W'(Ab,Ei) and W(Ei,Ab) terms into WEiAb ****/
289
dpd_buf4_init(&W, CC_TMP0, 0, 28, 26, 28, 26, 0, "W'(Ab,Ei)");
290
dpd_buf4_sort_axpy(&W, CC_HBAR, rspq, 26, 28, "WEiAb", 1);
294
}} // namespace psi::cchbar