3
\brief Enter brief description of file here
9
#include <libciomr/libciomr.h>
10
#include <libdpd/dpd.h>
16
namespace psi { namespace cclambda {
18
/* WefabL2(): Computes the contribution of the Wefab HBAR matrix
19
** elements to the Lambda double de-excitation amplitude equations.
20
** These contributions are given in spin-orbitals as:
22
** L_ij^ab = 1/2 L_ij^ef Wefab
24
** where W_efab is defined in spin orbitals as:
26
** Wefab = <ef||ab> - P(ef) t_m^f <em||ab> + 1/2 tau_mn^ef <mn||ab>
28
** and tau_mn^ef = t_mn^ef + t_m^e t_n^f - t_m^f t_n^e.
30
** [cf. Gauss and Stanton, JCP 103, 3561-3577 (1995).]
32
** NB: Wefab is not symmetric, so one must be careful when defining
33
** intermediate quantities for efficient contractions. I use the
34
** following contraction steps for each spin case:
36
** Wefab term II spin cases:
38
** L_IJ^AB <-- 1/2 ( -t_M^F <EM||AB> L_IJ^EF + t_M^E <FM||AB> L_IJ^EF )
40
** Z(IJ,EM) = -t_M^F L_IJ^EF
42
** L_IJ^AB <-- Z(IJ,EM) <EM||AB>
44
** L_ij^ab <-- 1/2 ( -t_m^f <em||ab> L_ij^ef + t_m^e <fm||ab> L_ij^ef )
46
** Z(ij,em) = -t_m^f L_ij^ef
48
** L_ij^ab <-- Z(ij,em) <em||ab>
50
** L_Ij^Ab <-- -t_m^f <Em|Ab> L_Ij^Ef - t_M^E <Mf|Ab> L_Ij^Ef
52
** Z(Ij,Em) = -t_m^f L_Ij^Ef
53
** Z(Ij,Mf) = -t_M^E L_Ij^Ef
55
** L_Ij^Ab <-- Z(Ij,Em) <Em|Ab> + Z(Ij,Mf) <Mf|Ab>
59
** L_IJ^AB <-- 1/4 tau_MN^EF <MN||AB> L_IJ^EF
61
** Z(IJ,MN) = 1/2 tau_MN^EF L_IJ^EF
63
** L_IJ^AB <-- 1/2 Z(IJ,MN) <MN||AB>
65
** L_ij^ab <-- 1/4 tau_mn^ef <mn||ab> L_ij^ef
67
** Z(ij,mn) = 1/2 tau_mn^ef L_ij^ef
69
** L_ij^ab <-- 1/2 Z(ij,mn) <mn||ab>
71
** L_Ij^Ab <-- tau_Mn^Ef <Mn|Ab> L_Ij^Ef
73
** Z(Ij,Mn) = tau_Mn^Ef L_Ij^Ef
75
** L_Ij^Ab <-- Z(Ij,Mn) <Mn|Ab>
81
void WefabL2(int L_irr)
83
dpdbuf4 Lijab, LIJAB, LIjAb;
84
dpdbuf4 newLijab, newLIJAB, newLIjAb;
85
dpdbuf4 Tau, T2, Z, Z1, Z2, L, L2, B, D, F, Ltmp;
87
dpdbuf4 tau_a, tau_s, tau;
91
double **B_diag, **tau_diag;
93
int nbuckets, rows_per_bucket, rows_left, m, row_start, ab, cd, dc, d;
94
int nrows, ncols, nlinks;
97
/* RHS += Wefab*Lijef */
98
if(params.ref == 0) { /** RHF **/
100
if(!strcmp(params.abcd,"OLD")) {
101
dpd_buf4_init(&LIjAb, CC_LAMBDA, L_irr, 0, 5, 0, 5, 0, "LIjAb");
102
dpd_buf4_init(&Z, CC_TMP0, L_irr, 5, 0, 5, 0, 0, "ZAbIj");
103
dpd_buf4_init(&B, CC_BINTS, 0, 5, 5, 5, 5, 0, "B <ab|cd>");
104
dpd_contract444(&B, &LIjAb, &Z, 0, 0, 1, 0);
106
dpd_buf4_sort_axpy(&Z, CC_LAMBDA, rspq, 0, 5, "New LIjAb", 1);
108
dpd_buf4_close(&LIjAb);
110
else if(!strcmp(params.abcd,"NEW")) {
111
timer_on("ABCD:new");
113
/* L_a(-)(ij,ab) (i>j, a>b) = L(ij,ab) - L(ij,ba) */
114
dpd_buf4_init(&tau_a, CC_LAMBDA, L_irr, 4, 9, 0, 5, 1, "LIjAb");
115
dpd_buf4_copy(&tau_a, CC_LAMBDA, "L(-)(ij,ab)");
116
dpd_buf4_close(&tau_a);
118
/* L_s(+)(ij,ab) (i>=j, a>=b) = L(ij,ab) + L(ij,ba) */
119
dpd_buf4_init(&tau_a, CC_LAMBDA, L_irr, 0, 5, 0, 5, 0, "LIjAb");
120
dpd_buf4_copy(&tau_a, CC_TMP0, "L(+)(ij,ab)");
121
dpd_buf4_sort_axpy(&tau_a, CC_TMP0, pqsr, 0, 5, "L(+)(ij,ab)", 1);
122
dpd_buf4_close(&tau_a);
123
dpd_buf4_init(&tau_a, CC_TMP0, L_irr, 3, 8, 0, 5, 0, "L(+)(ij,ab)");
124
dpd_buf4_copy(&tau_a, CC_LAMBDA, "L(+)(ij,ab)");
125
dpd_buf4_close(&tau_a);
128
dpd_buf4_init(&tau_s, CC_LAMBDA, L_irr, 3, 8, 3, 8, 0, "L(+)(ij,ab)");
129
dpd_buf4_init(&B_s, CC_BINTS, 0, 8, 8, 8, 8, 0, "B(+) <ab|cd> + <ab|dc>");
130
dpd_buf4_init(&S, CC_TMP0, L_irr, 8, 3, 8, 3, 0, "S(ab,ij)");
131
dpd_contract444(&B_s, &tau_s, &S, 0, 0, 0.5, 0);
133
dpd_buf4_close(&B_s);
134
dpd_buf4_close(&tau_s);
137
/* L_diag(ij,c) = 2 * L(ij,cc)*/
139
/* NB: Gcc = 0, and B is totally symmetric, so Gab = 0 */
140
/* But Gij = L_irr ^ Gab = L_irr */
141
dpd_buf4_init(&tau, CC_LAMBDA, L_irr, 3, 8, 3, 8, 0, "L(+)(ij,ab)");
142
dpd_buf4_mat_irrep_init(&tau, L_irr);
143
dpd_buf4_mat_irrep_rd(&tau, L_irr);
144
tau_diag = dpd_block_matrix(tau.params->rowtot[L_irr], moinfo.nvirt);
145
for(ij=0; ij < tau.params->rowtot[L_irr]; ij++)
146
for(Gc=0; Gc < moinfo.nirreps; Gc++)
147
for(C=0; C < moinfo.virtpi[Gc]; C++) {
148
c = C + moinfo.vir_off[Gc];
149
cc = tau.params->colidx[c][c];
150
tau_diag[ij][c] = tau.matrix[L_irr][ij][cc];
152
dpd_buf4_mat_irrep_close(&tau, L_irr);
154
dpd_buf4_init(&B_s, CC_BINTS, 0, 8, 8, 8, 8, 0, "B(+) <ab|cd> + <ab|dc>");
155
dpd_buf4_init(&S, CC_TMP0, L_irr, 8, 3, 8, 3, 0, "S(ab,ij)");
156
dpd_buf4_mat_irrep_init(&S, 0);
157
dpd_buf4_mat_irrep_rd(&S, 0);
159
rows_per_bucket = dpd_memfree()/(B_s.params->coltot[0] + moinfo.nvirt);
160
if(rows_per_bucket > B_s.params->rowtot[0]) rows_per_bucket = B_s.params->rowtot[0];
161
nbuckets = (int) ceil((double) B_s.params->rowtot[0]/(double) rows_per_bucket);
162
rows_left = B_s.params->rowtot[0] % rows_per_bucket;
164
B_diag = dpd_block_matrix(rows_per_bucket, moinfo.nvirt);
166
ncols = tau.params->rowtot[L_irr];
167
nlinks = moinfo.nvirt;
168
for(m=0; m < (rows_left ? nbuckets-1:nbuckets); m++) {
169
row_start = m * rows_per_bucket;
170
nrows = rows_per_bucket;
171
if(nrows && ncols && nlinks) {
172
psio_read(CC_BINTS,"B(+) <ab|cc>",(char *) B_diag[0],nrows*nlinks*sizeof(double),next, &next);
173
C_DGEMM('n', 't', nrows, ncols, nlinks, -0.25, B_diag[0], nlinks,
174
tau_diag[0], nlinks, 1, S.matrix[0][row_start], ncols);
179
row_start = m * rows_per_bucket;
181
if(nrows && ncols && nlinks) {
182
psio_read(CC_BINTS,"B(+) <ab|cc>",(char *) B_diag[0],nrows*nlinks*sizeof(double),next, &next);
183
C_DGEMM('n', 't', nrows, ncols, nlinks, -0.25, B_diag[0], nlinks,
184
tau_diag[0], nlinks, 1, S.matrix[0][row_start], ncols);
187
dpd_buf4_mat_irrep_wrt(&S, 0);
188
dpd_buf4_mat_irrep_close(&S, 0);
190
dpd_buf4_close(&B_s);
191
dpd_free_block(B_diag, rows_per_bucket, moinfo.nvirt);
192
dpd_free_block(tau_diag, tau.params->rowtot[L_irr], moinfo.nvirt);
193
dpd_buf4_close(&tau);
196
dpd_buf4_init(&tau_a, CC_LAMBDA, L_irr, 4, 9, 4, 9, 0, "L(-)(ij,ab)");
197
dpd_buf4_init(&B_a, CC_BINTS, 0, 9, 9, 9, 9, 0, "B(-) <ab|cd> - <ab|dc>");
198
dpd_buf4_init(&A, CC_TMP0, L_irr, 9, 4, 9, 4, 0, "A(ab,ij)");
199
dpd_contract444(&B_a, &tau_a, &A, 0, 0, 0.5, 0);
201
dpd_buf4_close(&B_a);
202
dpd_buf4_close(&tau_a);
205
timer_on("ABCD:axpy");
206
dpd_buf4_init(&S, CC_TMP0, L_irr, 5, 0, 8, 3, 0, "S(ab,ij)");
207
dpd_buf4_sort_axpy(&S, CC_LAMBDA, rspq, 0, 5, "New LIjAb", 1);
209
dpd_buf4_init(&A, CC_TMP0, L_irr, 5, 0, 9, 4, 0, "A(ab,ij)");
210
dpd_buf4_sort_axpy(&A, CC_LAMBDA, rspq, 0, 5, "New LIjAb", 1);
212
timer_off("ABCD:axpy");
213
timer_off("ABCD:new");
216
dpd_buf4_init(&newLIjAb, CC_LAMBDA, L_irr, 0, 5, 0, 5, 0, "New LIjAb");
218
dpd_file2_init(&tIA, CC_OEI, 0, 0, 1, "tIA");
219
dpd_buf4_init(&LIjAb, CC_LAMBDA, L_irr, 0, 5, 0, 5, 0, "LIjAb");
220
dpd_buf4_init(&Z, CC_TMP0, L_irr, 10, 0, 10, 0, 0, "Z(Mf,Ij)");
221
dpd_contract244(&tIA, &LIjAb, &Z, 1, 2, 0, 1, 0);
223
dpd_buf4_close(&LIjAb);
224
dpd_file2_close(&tIA);
226
dpd_buf4_init(&F, CC_FINTS, 0, 10, 5, 10, 5, 0, "F <ia|bc>");
227
dpd_buf4_init(&Z, CC_TMP0, L_irr, 10, 0, 10, 0, 0, "Z(Mf,Ij)");
228
dpd_buf4_init(&Z1, CC_TMP0, L_irr, 5, 0, 5, 0, 0, "Z(Ab,Ij)");
229
dpd_contract444(&F, &Z, &Z1, 1, 1, -1, 0);
232
dpd_buf4_close(&newLIjAb);
233
dpd_buf4_sort_axpy(&Z1, CC_LAMBDA, srqp, 0, 5, "New LIjAb", 1);
234
dpd_buf4_sort_axpy(&Z1, CC_LAMBDA, rspq, 0, 5, "New LIjAb", 1);
235
dpd_buf4_init(&newLIjAb, CC_LAMBDA, L_irr, 0, 5, 0, 5, 0, "New LIjAb");
238
dpd_buf4_init(&Z, CC_TMP0, L_irr, 0, 0, 0, 0, 0, "Z(Ij,Mn)");
239
dpd_buf4_init(&Tau, CC_TAMPS, 0, 0, 5, 0, 5, 0, "tauIjAb");
240
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 0, 5, 0, 5, 0, "LIjAb");
241
dpd_contract444(&L2, &Tau, &Z, 0, 0, 1.0, 0.0);
243
dpd_buf4_close(&Tau);
244
dpd_buf4_init(&D, CC_DINTS, 0, 0, 5, 0, 5, 0, "D <ij|ab>");
245
dpd_contract444(&Z, &D, &newLIjAb, 0, 1, 1.0, 1.0);
249
dpd_buf4_close(&newLIjAb);
252
else if(params.ref == 1) { /** ROHF **/
254
dpd_file2_init(&tIA, CC_OEI, 0, 0, 1, "tIA");
255
dpd_file2_init(&tia, CC_OEI, 0, 0, 1, "tia");
257
dpd_buf4_init(&LIJAB, CC_LAMBDA, L_irr, 2, 7, 2, 7, 0, "LIJAB");
258
dpd_buf4_init(&Z, CC_TMP2, L_irr, 7, 2, 7, 2, 0, "ZABIJ");
259
dpd_buf4_init(&B, CC_BINTS, 0, 7, 7, 5, 5, 1, "B <ab|cd>");
260
dpd_contract444(&B, &LIJAB, &Z, 0, 0, 1, 0);
262
dpd_buf4_sort_axpy(&Z, CC_LAMBDA, rspq, 2, 7, "New LIJAB", 1);
263
dpd_buf4_close(&LIJAB);
266
dpd_buf4_init(&newLIJAB, CC_LAMBDA, L_irr, 2, 7, 2, 7, 0, "New LIJAB");
268
dpd_buf4_init(&LIJAB, CC_LAMBDA, L_irr, 2, 5, 2, 7, 0, "LIJAB");
269
dpd_buf4_init(&Ltmp, CC_TMP0, L_irr, 2, 10, 2, 10, 0, "Ltmp (I>J,MF)");
270
dpd_buf4_init(&F, CC_FINTS, 0, 10, 7, 10, 5, 1, "F <ia|bc>");
271
dpd_contract244(&tIA, &LIJAB, &Ltmp, 1, 2, 1, 1.0, 0.0);
272
dpd_contract444(&Ltmp, &F, &newLIJAB, 0, 1, -1.0, 1.0);
274
dpd_buf4_close(&Ltmp);
275
dpd_buf4_close(&LIJAB);
277
dpd_buf4_init(&Z, CC_TMP0, L_irr, 2, 2, 2, 2, 0, "Z(IJ,MN)");
278
dpd_buf4_init(&Tau, CC_TAMPS, 0, 2, 7, 2, 7, 0, "tauIJAB");
279
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 2, 7, 2, 7, 0, "LIJAB");
280
dpd_contract444(&L2, &Tau, &Z, 0, 0, 1.0, 0.0);
282
dpd_buf4_close(&Tau);
283
dpd_buf4_init(&D, CC_DINTS, 0, 2, 7, 2, 7, 0, "D <ij||ab> (i>j,a>b)");
284
dpd_contract444(&Z, &D, &newLIJAB, 0, 1, 1.0, 1.0);
287
dpd_buf4_close(&newLIJAB);
289
dpd_buf4_init(&Lijab, CC_LAMBDA, L_irr, 2, 7, 2, 7, 0, "Lijab");
290
dpd_buf4_init(&Z, CC_TMP2, L_irr, 7, 2, 7, 2, 0, "Zabij");
291
dpd_buf4_init(&B, CC_BINTS, 0, 7, 7, 5, 5, 1, "B <ab|cd>");
292
dpd_contract444(&B, &Lijab, &Z, 0, 0, 1, 0);
294
dpd_buf4_close(&Lijab);
295
dpd_buf4_sort_axpy(&Z, CC_LAMBDA, rspq, 2, 7, "New Lijab", 1);
298
dpd_buf4_init(&newLijab, CC_LAMBDA, L_irr, 2, 7, 2, 7, 0, "New Lijab");
300
dpd_buf4_init(&Lijab, CC_LAMBDA, L_irr, 2, 5, 2, 7, 0, "Lijab");
301
dpd_buf4_init(&Ltmp, CC_TMP0, L_irr, 2, 10, 2, 10, 0, "Ltmp (i>j,mf)");
302
dpd_buf4_init(&F, CC_FINTS, 0, 10, 7, 10, 5, 1, "F <ia|bc>");
303
dpd_contract244(&tia, &Lijab, &Ltmp, 1, 2, 1, 1.0, 0.0);
304
dpd_contract444(&Ltmp, &F, &newLijab, 0, 1, -1.0, 1.0);
306
dpd_buf4_close(&Ltmp);
307
dpd_buf4_close(&Lijab);
309
dpd_buf4_init(&Z, CC_TMP0, L_irr, 2, 2, 2, 2, 0, "Z(ij,mn)");
310
dpd_buf4_init(&Tau, CC_TAMPS, 0, 2, 7, 2, 7, 0, "tauijab");
311
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 2, 7, 2, 7, 0, "Lijab");
312
dpd_contract444(&L2, &Tau, &Z, 0, 0, 1.0, 0.0);
314
dpd_buf4_close(&Tau);
315
dpd_buf4_init(&D, CC_DINTS, 0, 2, 7, 2, 7, 0, "D <ij||ab> (i>j,a>b)");
316
dpd_contract444(&Z, &D, &newLijab, 0, 1, 1.0, 1.0);
319
dpd_buf4_close(&newLijab);
322
dpd_buf4_init(&LIjAb, CC_LAMBDA, L_irr, 0, 5, 0, 5, 0, "LIjAb");
323
dpd_buf4_init(&Z, CC_TMP2, L_irr, 5, 0, 5, 0, 0, "ZAbIj");
324
dpd_buf4_init(&B, CC_BINTS, 0, 5, 5, 5, 5, 0, "B <ab|cd>");
325
dpd_contract444(&B, &LIjAb, &Z, 0, 0, 1, 0);
327
dpd_buf4_sort_axpy(&Z, CC_LAMBDA, rspq, 0, 5, "New LIjAb", 1);
330
dpd_buf4_init(&newLIjAb, CC_LAMBDA, L_irr, 0, 5, 0, 5, 0, "New LIjAb");
332
dpd_buf4_init(&Ltmp, CC_TMP1, L_irr, 0, 11, 0, 11, 0, "Lt (Ij,Em)");
333
dpd_contract424(&LIjAb, &tia, &Ltmp, 3, 1, 0, 1.0, 0.0);
334
dpd_buf4_sort(&Ltmp, CC_TMP2, pqsr, 0, 10, "Lt (Ij,mE)");
335
dpd_buf4_close(&Ltmp);
336
dpd_buf4_init(&Ltmp, CC_TMP3, L_irr, 0, 10, 0, 10, 0, "Lt (Ij,Mf)");
337
dpd_contract244(&tIA, &LIjAb, &Ltmp, 1, 2, 1, 1.0, 0.0);
338
dpd_buf4_close(&Ltmp);
340
dpd_buf4_close(&LIjAb);
342
dpd_buf4_init(&F, CC_FINTS, 0, 10, 5, 10, 5, 0, "F <ia|bc>");
343
dpd_buf4_init(&Ltmp, CC_TMP3, L_irr, 0, 10, 0, 10, 0, "Lt (Ij,Mf)");
344
dpd_contract444(&Ltmp, &F, &newLIjAb, 0, 1, -1.0, 1.0);
345
dpd_buf4_close(&Ltmp);
346
dpd_buf4_sort(&F, CC_TMP0, pqsr, 10, 5, "<me|ab> (mE,Ab)");
349
dpd_buf4_init(&F, CC_TMP0, 0, 10, 5, 10, 5, 0, "<me|ab> (mE,Ab)");
350
dpd_buf4_init(&Ltmp, CC_TMP2, L_irr, 0, 10, 0, 10, 0, "Lt (Ij,mE)");
351
dpd_contract444(&Ltmp, &F, &newLIjAb, 0, 1, -1.0, 1.0);
352
dpd_buf4_close(&Ltmp);
355
dpd_buf4_init(&Z, CC_TMP0, L_irr, 0, 0, 0, 0, 0, "Z(Ij,Mn)");
356
dpd_buf4_init(&Tau, CC_TAMPS, 0, 0, 5, 0, 5, 0, "tauIjAb");
357
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 0, 5, 0, 5, 0, "LIjAb");
358
dpd_contract444(&L2, &Tau, &Z, 0, 0, 1.0, 0.0);
360
dpd_buf4_close(&Tau);
361
dpd_buf4_init(&D, CC_DINTS, 0, 0, 5, 0, 5, 0, "D <ij|ab>");
362
dpd_contract444(&Z, &D, &newLIjAb, 0, 1, 1.0, 1.0);
366
dpd_buf4_close(&newLIjAb);
368
dpd_file2_close(&tIA);
369
dpd_file2_close(&tia);
371
else if(params.ref == 2) { /** UHF **/
373
dpd_file2_init(&tIA, CC_OEI, 0, 0, 1, "tIA");
374
dpd_file2_init(&tia, CC_OEI, 0, 2, 3, "tia");
377
/** Z(AB,IJ) = <AB||CD> L(IJ,CD) **/
378
dpd_buf4_init(&Z, CC_TMP1, L_irr, 7, 2, 7, 2, 0, "Z(AB,IJ)");
379
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 2, 7, 2, 7, 0, "LIJAB");
380
dpd_buf4_init(&B, CC_BINTS, 0, 7, 7, 5, 5, 1, "B <AB|CD>");
381
dpd_contract444(&B, &L2, &Z, 0, 0, 1, 0);
384
/** Z(AB,IJ) --> New L(IJ,AB) **/
385
dpd_buf4_sort_axpy(&Z, CC_LAMBDA, rspq, 2, 7, "New LIJAB", 1);
388
/** Z(IJ,EM) = -L(IJ,EFf) t(M,F) **/
389
dpd_buf4_init(&Z, CC_TMP1, L_irr, 2, 21, 2, 21, 0, "Z(IJ,EM)");
390
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 2, 5, 2, 7, 0, "LIJAB");
391
dpd_contract424(&L2, &tIA, &Z, 3, 1, 0, -1, 0);
393
/** New L(IJ,AB) <-- Z(IJ,EM) <EM||AB> **/
394
dpd_buf4_init(&F, CC_FINTS, 0, 21, 7, 21, 5, 1, "F <AI|BC>");
395
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 2, 7, 2, 7, 0, "New LIJAB");
396
dpd_contract444(&Z, &F, &L2, 0, 1, 1, 1);
401
/** Z(IJ,MN) = 1/2 L(IJ,EF) tau_MN^EF **/
402
dpd_buf4_init(&Z, CC_TMP1, L_irr, 2, 2, 2, 2, 0, "Z(IJ,MN)");
403
dpd_buf4_init(&T2, CC_TAMPS, 0, 2, 7, 2, 7, 0, "tauIJAB");
404
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 2, 7, 2, 7, 0, "LIJAB");
405
dpd_contract444(&L2, &T2, &Z, 0, 0, 1, 0);
408
/** New L(IJ,AB) <-- 1/2 Z(IJ,MN) <MN||AB> **/
409
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 2, 7, 2, 7, 0, "New LIJAB");
410
dpd_buf4_init(&D, CC_DINTS, 0, 2, 7, 2, 7, 0, "D <IJ||AB> (I>J,A>B)");
411
dpd_contract444(&Z, &D, &L2, 0, 1, 1, 1);
417
/** Z(ab,ij) = <ab||cd> L(ij,cd) **/
418
dpd_buf4_init(&Z, CC_TMP1, L_irr, 17, 12, 17, 12, 0, "Z(ab,ij)");
419
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 12, 17, 12, 17, 0, "Lijab");
420
dpd_buf4_init(&B, CC_BINTS, 0, 17, 17, 15, 15, 1, "B <ab|cd>");
421
dpd_contract444(&B, &L2, &Z, 0, 0, 1, 0);
424
/** Z(ab,ij) --> New L(ij,ab) **/
425
dpd_buf4_sort_axpy(&Z, CC_LAMBDA, rspq, 12, 17, "New Lijab", 1);
428
/** Z(ij,em) = -L(ij,ef) t(m,f) **/
429
dpd_buf4_init(&Z, CC_TMP1, L_irr, 12, 31, 12, 31, 0, "Z(ij,em)");
430
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 12, 15, 12, 17, 0, "Lijab");
431
dpd_contract424(&L2, &tia, &Z, 3, 1, 0, -1, 0);
433
/** New L(ij,ab) <-- Z(ij,em) <em||ab> **/
434
dpd_buf4_init(&F, CC_FINTS, 0, 31, 17, 31, 15, 1, "F <ai|bc>");
435
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 12, 17, 12, 17, 0, "New Lijab");
436
dpd_contract444(&Z, &F, &L2, 0, 1, 1, 1);
441
/** Z(ij,mn) = 1/2 L(ij,ef) tau_mn^ef **/
442
dpd_buf4_init(&Z, CC_TMP1, L_irr, 12, 12, 12, 12, 0, "Z(ij,mn)");
443
dpd_buf4_init(&T2, CC_TAMPS, 0, 12, 17, 12, 17, 0, "tauijab");
444
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 12, 17, 12, 17, 0, "Lijab");
445
dpd_contract444(&L2, &T2, &Z, 0, 0, 1, 0);
448
/** New L(ij,ab) <-- 1/2 Z(ij,mn) <mn||ab> **/
449
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 12, 17, 12, 17, 0, "New Lijab");
450
dpd_buf4_init(&D, CC_DINTS, 0, 12, 17, 12, 17, 0, "D <ij||ab> (i>j,a>b)");
451
dpd_contract444(&Z, &D, &L2, 0, 1, 1, 1);
457
/** Z(Ab,Ij) = <Ab|Cd> L(Ij,Cd) **/
458
dpd_buf4_init(&Z, CC_TMP1, L_irr, 28, 22, 28, 22, 0, "Z(Ab,Ij)");
459
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 22, 28, 22, 28, 0, "LIjAb");
460
dpd_buf4_init(&B, CC_BINTS, 0, 28, 28, 28, 28, 0, "B <Ab|Cd>");
461
dpd_contract444(&B, &L2, &Z, 0, 0, 1, 0);
464
/** Z(Ab,Ij) --> New L(Ij,Ab) **/
465
dpd_buf4_sort_axpy(&Z, CC_LAMBDA, rspq, 22, 28, "New LIjAb", 1);
468
/** Z(Ij,Em) = -L(Ij,Ef) t(m,f) **/
469
dpd_buf4_init(&Z, CC_TMP1, L_irr, 22, 26, 22, 26, 0, "Z(Ij,Em)");
470
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 22, 28, 22, 28, 0, "LIjAb");
471
dpd_contract424(&L2, &tia, &Z, 3, 1, 0, -1, 0);
473
/** New L(Ij,Ab) <-- Z(Ij,Em) <Em|Ab> **/
474
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 22, 28, 22, 28, 0, "New LIjAb");
475
dpd_buf4_init(&F, CC_FINTS, 0, 26, 28, 26, 28, 0, "F <Ai|Bc>");
476
dpd_contract444(&Z, &F, &L2, 0, 1, 1, 1);
481
/** Z(Ij,Mf) = -t(M,E) L(Ij,Ef) **/
482
dpd_buf4_init(&Z, CC_TMP1, L_irr, 22, 24, 22, 24, 0, "Z(Ij,Mf)");
483
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 22, 28, 22, 28, 0, "LIjAb");
484
dpd_contract244(&tIA, &L2, &Z, 1, 2, 1, -1, 0);
486
/** New L(Ij,Ab) <-- Z(Ij,Mf) <Mf|Ab> **/
487
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 22, 28, 22, 28, 0, "New LIjAb");
488
dpd_buf4_init(&F, CC_FINTS, 0, 24, 28, 24, 28, 0, "F <Ia|Bc>");
489
dpd_contract444(&Z, &F, &L2, 0, 1, 1, 1);
494
/** Z(Ij,Mn) = L(Ij,Ef) tau(Mn,Ef) **/
495
dpd_buf4_init(&Z, CC_TMP1, L_irr, 22, 22, 22, 22, 0, "Z(Ij,Mn)");
496
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 22, 28, 22, 28, 0, "LIjAb");
497
dpd_buf4_init(&T2, CC_TAMPS, 0, 22, 28, 22, 28, 0, "tauIjAb");
498
dpd_contract444(&L2, &T2, &Z, 0, 0, 1, 0);
501
/** New L(Ij,Ab) <-- Z(Ij,Mn) <Mn|Ab> **/
502
dpd_buf4_init(&L2, CC_LAMBDA, L_irr, 22, 28, 22, 28, 0, "New LIjAb");
503
dpd_buf4_init(&D, CC_DINTS, 0, 22, 28, 22, 28, 0, "D <Ij|Ab>");
504
dpd_contract444(&Z, &D, &L2, 0, 1, 1, 1);
509
dpd_file2_close(&tIA);
510
dpd_file2_close(&tia);
516
}} // namespace psi::cclambda
added
removed
src/bin/cclambda/WefabL2.cc
