4
#include <libdpd/dpd.h>
9
void cc3_l3l2_RHF_AAA(void);
10
void cc3_l3l2_RHF_AAB(void);
12
void L3_AAA(double ***W1, int nirreps, int I, int Gi, int J, int Gj, int K, int Gk,
13
dpdbuf4 *T2, dpdbuf4 *F, dpdbuf4 *E, dpdfile2 *fIJ, dpdfile2 *fAB,
14
dpdbuf4 *D, dpdbuf4 *LIJAB, dpdfile2 *LIA, dpdfile2 *FME,
15
int *occpi, int *occ_off, int *virtpi, int *vir_off);
17
void L3_AAB(double ***W1, int nirreps, int I, int Gi, int J, int Gj, int K, int Gk,
18
dpdbuf4 *T2AA, dpdbuf4 *T2AB, dpdbuf4 *T2BA, dpdbuf4 *FAA, dpdbuf4 *FAB, dpdbuf4 *FBA,
19
dpdbuf4 *EAA, dpdbuf4 *EAB, dpdbuf4 *EBA, dpdfile2 *fIJ, dpdfile2 *fij,
20
dpdfile2 *fAB, dpdfile2 *fab, dpdbuf4 *DAA, dpdbuf4 *DAB, dpdbuf4 *LIJAB, dpdbuf4 *LIjAb,
21
dpdfile2 *LIA, dpdfile2 *Lia, dpdfile2 *FME, dpdfile2 *Fme,
22
int *aoccpi, int *aocc_off, int *boccpi, int *bocc_off,
23
int *avirtpi, int *avir_off, int *bvirtpi, int *bvir_off);
33
void cc3_l3l2_RHF_AAA(void)
37
int *vir_off, *virtpi;
54
dpdbuf4 L2new, L2, D2;
55
int Gjk, jk, Gid, id, Gik, ik;
61
int Gij, ij, Gmk, mk, am, Gbc, bc;
63
int nrows, ncols, nlinks;
66
nirreps = moinfo.nirreps;
68
occ_off = moinfo.occ_off;
69
virtpi = moinfo.virtpi;
70
vir_off = moinfo.vir_off;
72
dpd_buf4_init(&WMAFE, CC3_HET1, 0, 10, 5, 10, 7, 0, "CC3 WABEI (IE,B>A)");
73
dpd_buf4_init(&WMNIE, CC3_HET1, 0, 0, 10, 2, 10, 0, "CC3 WMBIJ (I>J,MB)");
74
for(h=0; h < nirreps; h++) {
75
dpd_buf4_mat_irrep_init(&WMNIE, h);
76
dpd_buf4_mat_irrep_rd(&WMNIE, h);
79
dpd_buf4_init(&L2new, CC3_MISC, 0, 0, 5, 0, 5, 0, "CC3 LIJAB");
80
for(h=0; h < nirreps; h++) dpd_buf4_mat_irrep_init(&L2new, h);
82
dpd_buf4_init(&ZIGDE, CC3_MISC, 0, 10, 5, 10, 5, 0, "CC3 ZIGDE");
83
dpd_buf4_scm(&ZIGDE, 0.0); /* must be cleared in each iteration */
85
dpd_buf4_init(&ZDMAE, CC3_MISC, 0, 10, 5, 10, 5, 0, "CC3 ZDMAE (MD,AE)");
86
dpd_buf4_scm(&ZDMAE, 0.0);
88
dpd_buf4_init(&ZLMAO, CC3_MISC, 0, 0, 11, 0, 11, 0, "CC3 ZLMAO");
89
for(h=0; h < nirreps; h++) dpd_buf4_mat_irrep_init(&ZLMAO, h);
91
dpd_buf4_init(&ZIMLE, CC3_MISC, 0, 0, 10, 0, 10, 0, "CC3 ZIMLE");
92
for(h=0; h < nirreps; h++) dpd_buf4_mat_irrep_init(&ZIMLE, h);
94
dpd_buf4_init(&T2, CC_TAMPS, 0, 0, 5, 2, 7, 0, "tIJAB");
95
for(h=0; h < nirreps; h++) {
96
dpd_buf4_mat_irrep_init(&T2, h);
97
dpd_buf4_mat_irrep_rd(&T2, h);
100
dpd_file2_init(&fIJ, CC_OEI, 0, 0, 0, "fIJ");
101
dpd_file2_init(&fAB, CC_OEI, 0, 1, 1, "fAB");
103
dpd_buf4_init(&L, CC_LAMBDA, 0, 0, 5, 2, 7, 0, "LIJAB");
104
dpd_buf4_init(&F, CC3_HET1, 0, 10, 5, 10, 7, 0, "CC3 WAMEF (MA,F>E)");
105
dpd_buf4_init(&E, CC3_HET1, 0, 0, 10, 2, 10, 0, "CC3 WMNIE (M>N,IE)");
107
dpd_buf4_init(&Dints, CC_DINTS, 0, 0, 5, 0, 5, 0, "D <ij||ab>");
108
dpd_buf4_init(&LIJAB, CC_LAMBDA, 0, 0, 5, 2, 7, 0, "LIJAB");
109
dpd_file2_init(&LIA, CC_LAMBDA, 0, 0, 1, "LIA");
110
dpd_file2_init(&FME, CC_OEI, 0, 0, 1, "FME");
112
W1 = (double ***) malloc(nirreps * sizeof(double **));
113
W2 = (double ***) malloc(nirreps * sizeof(double **));
115
for(Gi=0; Gi < nirreps; Gi++) {
116
for(Gj=0; Gj < nirreps; Gj++) {
118
for(Gk=0; Gk < nirreps; Gk++) {
124
for(Gab=0; Gab < nirreps; Gab++) {
125
Gc = Gab ^ Gijk; /* totally symmetric */
126
W1[Gab] = dpd_block_matrix(F.params->coltot[Gab], virtpi[Gc]);
128
for(Ga=0; Ga < nirreps; Ga++) {
129
Gbc = Ga ^ Gijk; /* totally symmetric */
130
W2[Ga] = dpd_block_matrix(virtpi[Ga], F.params->coltot[Gbc]);
133
for(i=0; i < occpi[Gi]; i++) {
135
for(j=0; j < occpi[Gj]; j++) {
137
for(k=0; k < occpi[Gk]; k++) {
140
L3_AAA(W1, nirreps, I, Gi, J, Gj, K, Gk, &L, &F, &E, &fIJ, &fAB,
141
&Dints, &LIJAB, &LIA, &FME, occpi, occ_off, virtpi, vir_off);
143
/* L_JKDC <-- +1/2 t_IJKABC W_ABID (IDAB) */
144
/* L_JKCD <-- -1/2 t_IJKABC W_ABID (IDAB) */
145
jk = L2new.params->rowidx[J][K];
146
for(Gd=0; Gd < nirreps; Gd++) {
147
Gab = Gid = Gi ^ Gd; /* assumes Wieab is totally symmetric */
148
Gc = Gab ^ Gijk; /* assumes T3 is totally symmetric */
150
id = WMAFE.row_offset[Gid][I];
152
Z = block_matrix(virtpi[Gc],virtpi[Gd]);
153
WMAFE.matrix[Gid] = dpd_block_matrix(virtpi[Gd], WMAFE.params->coltot[Gid]);
154
dpd_buf4_mat_irrep_rd_block(&WMAFE, Gid, id, virtpi[Gd]);
158
nlinks = WMAFE.params->coltot[Gid];
160
if(nrows && ncols && nlinks)
161
C_DGEMM('t', 't', nrows, ncols, nlinks, 0.5, W1[Gab][0], nrows,
162
WMAFE.matrix[Gid][0], nlinks, 0.0, Z[0], ncols);
164
for(c=0; c < virtpi[Gc]; c++) {
166
for(d=0; d < virtpi[Gd]; d++) {
167
DD = vir_off[Gd] + d;
168
cd = L2new.params->colidx[C][DD];
169
dc = L2new.params->colidx[DD][C];
170
L2new.matrix[Gjk][jk][dc] += Z[c][d];
171
L2new.matrix[Gjk][jk][cd] += -Z[c][d];
174
dpd_free_block(WMAFE.matrix[Gid], virtpi[Gd], WMAFE.params->coltot[Gid]);
178
/* t_MIAB <-- +1/2 t_IJKABC W_JKMC */
179
/* t_IMAB <-- -1/2 t_IJKABC W_JKMC */
180
jk = WMNIE.params->rowidx[J][K];
181
for(Gm=0; Gm < nirreps; Gm++) {
182
Gab = Gmi = Gm ^ Gi; /* assumes totally symmetric */
183
Gc = Gab ^ Gijk; /* assumes totally symmetric */
185
mc = WMNIE.col_offset[Gjk][Gm];
187
nrows = F.params->coltot[Gab];
191
Z = dpd_block_matrix(nrows, ncols);
193
if(nrows && ncols && nlinks)
194
C_DGEMM('n', 't', nrows, ncols, nlinks, 0.5, W1[Gab][0], nlinks,
195
&(WMNIE.matrix[Gjk][jk][mc]), nlinks, 0.0, Z[0], ncols);
197
for(m=0; m < ncols; m++) {
199
mi = L2new.params->rowidx[M][I];
200
im = L2new.params->rowidx[I][M];
201
for(ab=0; ab < nrows; ab++) {
202
L2new.matrix[Gmi][mi][ab] += Z[ab][m];
203
L2new.matrix[Gmi][im][ab] -= Z[ab][m];
207
dpd_free_block(Z, nrows, ncols);
210
/* Z_IDAB <-- 1/2 L_IJKABC t_JKDC */
212
jk = T2.params->rowidx[J][K];
213
for(Gab=0; Gab < nirreps; Gab++) {
214
Gid = Gab; /* totally symmetric */
215
Gc = Gab ^ Gijk; /* totally symmetric */
219
ncols = ZIGDE.params->coltot[Gid];
222
dc = T2.col_offset[Gjk][Gd];
223
id = ZIGDE.row_offset[Gid][I];
224
ZIGDE.matrix[Gid] = dpd_block_matrix(nrows, ncols);
225
dpd_buf4_mat_irrep_rd_block(&ZIGDE, Gid, id, nrows);
227
if(nrows && ncols && nlinks)
228
C_DGEMM('n', 't', nrows, ncols, nlinks, 0.5, &(T2.matrix[Gjk][jk][dc]), nlinks,
229
W1[Gab][0], nlinks, 1.0, ZIGDE.matrix[Gid][0], ncols);
231
dpd_buf4_mat_irrep_wrt_block(&ZIGDE, Gid, id, nrows);
232
dpd_free_block(ZIGDE.matrix[Gid], nrows, ncols);
235
/* Z_JDAB <-- 1/2 L_IJKABC t_IKDC */
236
ik = T2.params->rowidx[I][K];
237
for(Gab=0; Gab < nirreps; Gab++) {
238
Gjd = Gab; /* totally symmetric */
239
Gc = Gab ^ Gijk; /* totally symmetric */
243
ncols = ZDMAE.params->coltot[Gjd];
246
dc = T2.col_offset[Gik][Gd];
247
jd = ZDMAE.row_offset[Gjd][J];
248
ZDMAE.matrix[Gjd] = dpd_block_matrix(nrows, ncols);
249
dpd_buf4_mat_irrep_rd_block(&ZDMAE, Gjd, jd, nrows);
251
if(nrows && ncols && nlinks)
252
C_DGEMM('n', 't', nrows, ncols, nlinks, 0.5, &(T2.matrix[Gik][ik][dc]), nlinks,
253
W1[Gab][0], nlinks, 1.0, ZDMAE.matrix[Gjd][0], ncols);
255
dpd_buf4_mat_irrep_wrt_block(&ZDMAE, Gjd, jd, nrows);
256
dpd_free_block(ZDMAE.matrix[Gjd], nrows, ncols);
259
/* Z_IJAM <-- -1/2 L_IJKABC t_MKBC */
260
/* sort W(AB,C) to W(A,BC) */
261
for(Gab=0; Gab < nirreps; Gab++) {
263
for(ab=0; ab < F.params->coltot[Gab]; ab++) {
264
A = F.params->colorb[Gab][ab][0];
265
B = F.params->colorb[Gab][ab][1];
266
Ga = F.params->rsym[A];
268
for(c=0; c < virtpi[Gc]; c++) {
270
bc = F.params->colidx[B][C];
271
W2[Ga][a][bc] = W1[Gab][ab][c];
276
ij = ZLMAO.params->rowidx[I][J];
278
for(Gm=0; Gm < nirreps; Gm++) {
279
Gbc = Gmk = Gm ^ Gk; /* totally symmetric */
280
Ga = Gij ^ Gm; /* totally symmetric */
283
ncols = T2.params->coltot[Gmk];
285
for(m=0; m < occpi[Gm]; m++) {
287
mk = T2.params->rowidx[M][K];
288
am = ZLMAO.col_offset[Gij][Ga] + m;
291
C_DGEMV('n', nrows, ncols, -0.5, W2[Ga][0], ncols, T2.matrix[Gmk][mk], 1,
292
1.0, &(ZLMAO.matrix[Gij][ij][am]), occpi[Gm]);
296
/* Z_IJMB <-- -1/2 L_IJKABC t_MKAC */
297
/* sort W(AB,C) to W(B,AC) */
298
for(Gab=0; Gab < nirreps; Gab++) {
300
for(ab=0; ab < F.params->coltot[Gab]; ab++) {
301
A = F.params->colorb[Gab][ab][0];
302
B = F.params->colorb[Gab][ab][1];
303
Gb = F.params->ssym[B];
305
for(c=0; c < virtpi[Gc]; c++) {
307
ac = F.params->colidx[A][C];
308
W2[Gb][b][ac] = W1[Gab][ab][c];
313
ij = ZIMLE.params->rowidx[I][J];
315
for(Gm=0; Gm < nirreps; Gm++) {
316
Gb = Gm ^ Gij; /* totally symmetric */
320
ncols = T2.params->coltot[Gmk];
322
for(m=0; m < occpi[Gm]; m++) {
324
mk = T2.params->rowidx[M][K];
325
mb = ZIMLE.col_offset[Gij][Gm] + m * virtpi[Gb];
328
C_DGEMV('n', nrows, ncols, -0.5, W2[Gb][0], ncols, T2.matrix[Gmk][mk], 1,
329
1.0, &(ZIMLE.matrix[Gij][ij][mb]), 1);
339
for(Gab=0; Gab < nirreps; Gab++) {
340
Gc = Gab ^ Gijk; /* totally symmetric */
341
dpd_free_block(W1[Gab], F.params->coltot[Gab], virtpi[Gc]);
343
for(Ga=0; Ga < nirreps; Ga++) {
344
Gbc = Ga ^ Gijk; /* totally symmetric */
345
dpd_free_block(W2[Ga], virtpi[Ga], F.params->coltot[Gbc]);
359
dpd_file2_close(&fIJ);
360
dpd_file2_close(&fAB);
362
dpd_file2_close(&FME);
363
dpd_file2_close(&LIA);
364
dpd_buf4_close(&Dints);
365
dpd_buf4_close(&LIJAB);
367
dpd_buf4_close(&WMAFE);
368
for(h=0; h < nirreps; h++) dpd_buf4_mat_irrep_close(&WMNIE, h);
369
dpd_buf4_close(&WMNIE);
371
dpd_buf4_close(&ZIGDE);
372
dpd_buf4_close(&ZDMAE);
374
for(h=0; h < nirreps; h++) {
375
dpd_buf4_mat_irrep_wrt(&ZLMAO, h);
376
dpd_buf4_mat_irrep_close(&ZLMAO, h);
378
dpd_buf4_close(&ZLMAO);
380
for(h=0; h < nirreps; h++) {
381
dpd_buf4_mat_irrep_wrt(&ZIMLE, h);
382
dpd_buf4_mat_irrep_close(&ZIMLE, h);
384
dpd_buf4_close(&ZIMLE);
386
for(h=0; h < nirreps; h++) dpd_buf4_mat_irrep_close(&T2, h);
389
for(h=0; h < nirreps; h++) {
390
dpd_buf4_mat_irrep_wrt(&L2new, h);
391
dpd_buf4_mat_irrep_close(&L2new, h);
393
dpd_buf4_init(&D2, CC_DENOM, 0, 0, 5, 0, 5, 0, "dIjAb");
394
dpd_buf4_dirprd(&D2, &L2new);
396
dpd_buf4_init(&L2, CC_LAMBDA, 0, 0, 5, 2, 7, 0, "New LIJAB");
397
dpd_buf4_axpy(&L2new, &L2, 1);
399
dpd_buf4_close(&L2new);
403
void cc3_l3l2_RHF_AAB(void)
406
int *occ_off, *occpi;
407
int *vir_off, *virtpi;
408
int Gi, Gj, Gk, Gijk;
410
int i, j, k, I, J, K;
411
int a, b, c, A, B, C;
414
dpdbuf4 L2AA, L2AB, L2BA, EAA, EAB, EBA, FAA, FAB, FBA;
415
dpdfile2 fIJ, fAB, fij, fab;
416
dpdbuf4 DAAints, DABints, LIJAB, LIjAb;
417
dpdfile2 LIA, Lia, FME, Fme;
418
dpdbuf4 L2AAnew, L2ABnew, L2, D2;
419
dpdbuf4 WmAfE, WMnIe, WMAFE, WMaFe, WMNIE, WmNiE;
420
dpdbuf4 ZIGDE, T2AB, T2AA, ZIgDe;
421
dpdbuf4 ZDMAE, ZDmAe, ZdMAe;
422
dpdbuf4 ZLMAO, ZLmAo;
423
dpdbuf4 ZIMLE, ZImLe, ZImlE;
424
int nrows, ncols, nlinks;
426
int Gij, ij, Gji, ji, Gjk, jk, kj, Gkj;
427
int Gd, d, DD, ad, da, Gkd, kd;
428
int Gm, m, M, Gmi, mi, im, mc;
430
int Gac, ac, Gca, ca, bd, db;
431
int Gbc, bc, Gmk, mk, km, Gim, ma, am;
432
int Gik, ik, Gki, ki, Gjd, jd;
433
int Gmj, mj, cm, Gjm, jm;
437
nirreps = moinfo.nirreps;
438
occpi = moinfo.occpi;
439
occ_off = moinfo.occ_off;
440
virtpi = moinfo.virtpi;
441
vir_off = moinfo.vir_off;
443
dpd_buf4_init(&L2AAnew, CC3_MISC, 0, 0, 5, 0, 5, 0, "CC3 LIJAB");
444
for(h=0; h < nirreps; h++) dpd_buf4_mat_irrep_init(&L2AAnew, h);
446
dpd_buf4_init(&L2ABnew, CC3_MISC, 0, 0, 5, 0, 5, 0, "CC3 LIjAb");
447
for(h=0; h < nirreps; h++) dpd_buf4_mat_irrep_init(&L2ABnew, h);
449
dpd_buf4_init(&T2AB, CC_TAMPS, 0, 0, 5, 0, 5, 0, "tIjAb");
450
dpd_buf4_init(&T2AA, CC_TAMPS, 0, 0, 5, 2, 7, 0, "tIJAB");
451
for(h=0; h < nirreps; h++) {
452
dpd_buf4_mat_irrep_init(&T2AB, h);
453
dpd_buf4_mat_irrep_rd(&T2AB, h);
455
dpd_buf4_mat_irrep_init(&T2AA, h);
456
dpd_buf4_mat_irrep_rd(&T2AA, h);
459
dpd_buf4_init(&ZIGDE, CC3_MISC, 0, 10, 5, 10, 5, 0, "CC3 ZIGDE");
460
dpd_buf4_scm(&ZIGDE, 0.0); /* this must be cleared in each iteration */
461
dpd_buf4_init(&ZIgDe, CC3_MISC, 0, 10, 5, 10, 5, 0, "CC3 ZIgDe");
462
dpd_buf4_scm(&ZIgDe, 0.0); /* this must be cleared in each iteration */
464
dpd_buf4_init(&ZDMAE, CC3_MISC, 0, 10, 5, 10, 5, 0, "CC3 ZDMAE (MD,AE)");
465
dpd_buf4_scm(&ZDMAE, 0.0); /* must be cleared in each iteration */
466
dpd_buf4_init(&ZDmAe, CC3_MISC, 0, 10, 5, 10, 5, 0, "CC3 ZDmAe (mD,Ae)");
467
dpd_buf4_scm(&ZDmAe, 0.0); /* must be cleared in each iteration */
468
dpd_buf4_init(&ZdMAe, CC3_MISC, 0, 10, 5, 10, 5, 0, "CC3 ZdMAe (Md,Ae)");
469
dpd_buf4_scm(&ZdMAe, 0.0); /* must be cleared in each iteration */
471
dpd_buf4_init(&ZLMAO, CC3_MISC, 0, 0, 11, 0, 11, 0, "CC3 ZLMAO");
472
dpd_buf4_init(&ZLmAo, CC3_MISC, 0, 0, 11, 0, 11, 0, "CC3 ZLmAo");
473
for(h=0; h < nirreps; h++) {
474
dpd_buf4_mat_irrep_init(&ZLMAO, h);
475
dpd_buf4_mat_irrep_rd(&ZLMAO, h);
477
dpd_buf4_mat_irrep_init(&ZLmAo, h);
480
dpd_buf4_init(&ZIMLE, CC3_MISC, 0, 0, 10, 0, 10, 0, "CC3 ZIMLE");
481
dpd_buf4_init(&ZImLe, CC3_MISC, 0, 0, 10, 0, 10, 0, "CC3 ZImLe");
482
dpd_buf4_init(&ZImlE, CC3_MISC, 0, 0, 10, 0, 10, 0, "CC3 ZImlE");
483
for(h=0; h < nirreps; h++) {
484
dpd_buf4_mat_irrep_init(&ZIMLE, h);
485
dpd_buf4_mat_irrep_rd(&ZIMLE, h);
487
dpd_buf4_mat_irrep_init(&ZImLe, h);
488
dpd_buf4_mat_irrep_init(&ZImlE, h);
491
dpd_buf4_init(&WmAfE, CC3_HET1, 0, 10, 5, 10, 5, 0, "CC3 WAbEi (iE,bA)");
492
dpd_buf4_init(&WMAFE, CC3_HET1, 0, 10, 5, 10, 7, 0, "CC3 WABEI (IE,B>A)");
493
dpd_buf4_init(&WMaFe, CC3_HET1, 0, 10, 5, 10, 5, 0, "CC3 WaBeI (Ie,Ba)");
495
dpd_buf4_init(&WMnIe, CC3_HET1, 0, 0, 10, 0, 10, 0, "CC3 WMbIj (Ij,Mb)");
496
dpd_buf4_init(&WMNIE, CC3_HET1, 0, 0, 10, 2, 10, 0, "CC3 WMBIJ (I>J,MB)");
497
dpd_buf4_init(&WmNiE, CC3_HET1, 0, 0, 10, 0, 10, 0, "CC3 WmBiJ (iJ,mB)");
498
for(h=0; h < nirreps; h++) {
499
dpd_buf4_mat_irrep_init(&WMnIe, h);
500
dpd_buf4_mat_irrep_rd(&WMnIe, h);
501
dpd_buf4_mat_irrep_init(&WMNIE, h);
502
dpd_buf4_mat_irrep_rd(&WMNIE, h);
503
dpd_buf4_mat_irrep_init(&WmNiE, h);
504
dpd_buf4_mat_irrep_rd(&WmNiE, h);
507
dpd_file2_init(&fIJ, CC_OEI, 0, 0, 0, "fIJ");
508
dpd_file2_init(&fAB, CC_OEI, 0, 1, 1, "fAB");
509
dpd_file2_init(&fij, CC_OEI, 0, 0, 0, "fij");
510
dpd_file2_init(&fab, CC_OEI, 0, 1, 1, "fab");
512
dpd_buf4_init(&L2AA, CC_LAMBDA, 0, 0, 5, 2, 7, 0, "LIJAB");
513
dpd_buf4_init(&L2AB, CC_LAMBDA, 0, 0, 5, 0, 5, 0, "LIjAb");
514
dpd_buf4_init(&L2BA, CC_LAMBDA, 0, 0, 5, 0, 5, 0, "LiJaB");
515
dpd_buf4_init(&FAA, CC3_HET1, 0, 10, 5, 10, 7, 0, "CC3 WAMEF (MA,F>E)");
516
dpd_buf4_init(&FAB, CC3_HET1, 0, 10, 5, 10, 5, 0, "CC3 WaMeF (Ma,Fe)");
517
dpd_buf4_init(&FBA, CC3_HET1, 0, 10, 5, 10, 5, 0, "CC3 WAmEf (mA,fE)");
518
dpd_buf4_init(&EAA, CC3_HET1, 0, 0, 10, 2, 10, 0, "CC3 WMNIE (M>N,IE)");
519
dpd_buf4_init(&EAB, CC3_HET1, 0, 0, 10, 0, 10, 0, "CC3 WMnIe (Mn,Ie)");
520
dpd_buf4_init(&EBA, CC3_HET1, 0, 0, 10, 0, 10, 0, "CC3 WmNiE (mN,iE)");
522
dpd_buf4_init(&DAAints, CC_DINTS, 0, 0, 5, 0, 5, 0, "D <ij||ab>");
523
dpd_buf4_init(&DABints, CC_DINTS, 0, 0, 5, 0, 5, 0, "D <ij|ab>");
524
dpd_buf4_init(&LIJAB, CC_LAMBDA, 0, 0, 5, 2, 7, 0, "LIJAB");
525
dpd_buf4_init(&LIjAb, CC_LAMBDA, 0, 0, 5, 0, 5, 0, "LIjAb");
526
dpd_file2_init(&LIA, CC_LAMBDA, 0, 0, 1, "LIA");
527
dpd_file2_init(&Lia, CC_LAMBDA, 0, 0, 1, "Lia");
528
dpd_file2_init(&FME, CC_OEI, 0, 0, 1, "FME");
529
dpd_file2_init(&Fme, CC_OEI, 0, 0, 1, "Fme");
531
W1 = (double ***) malloc(nirreps * sizeof(double **));
532
W2 = (double ***) malloc(nirreps * sizeof(double **));
534
for(Gi=0; Gi < nirreps; Gi++) {
535
for(Gj=0; Gj < nirreps; Gj++) {
537
for(Gk=0; Gk < nirreps; Gk++) {
543
for(Gab=0; Gab < nirreps; Gab++) {
544
Gc = Gab ^ Gijk; /* totally symmetric */
545
W1[Gab] = dpd_block_matrix(FAA.params->coltot[Gab], virtpi[Gc]);
547
for(Ga=0; Ga < nirreps; Ga++) {
548
Gcb = Ga ^ Gijk; /* assumes totally symmetric */
549
W2[Ga] = dpd_block_matrix(virtpi[Ga], WmAfE.params->coltot[Gcb]); /* alpha-beta-alpha */
552
for(i=0; i < occpi[Gi]; i++) {
554
for(j=0; j < occpi[Gj]; j++) {
556
for(k=0; k < occpi[Gk]; k++) {
559
L3_AAB(W1, nirreps, I, Gi, J, Gj, K, Gk, &L2AA, &L2AB, &L2BA,
560
&FAA, &FAB, &FBA, &EAA, &EAB, &EBA, &fIJ, &fij, &fAB, &fab,
561
&DAAints, &DABints, &LIJAB, &LIjAb, &LIA, &Lia, &FME, &Fme,
562
occpi, occ_off, occpi, occ_off, virtpi, vir_off, virtpi, vir_off);
564
/* t_JIDA <-- t_IJkABc W_kDcB */
565
/* sort W1(AB,c) to W2(A,cB) */
566
for(Gab=0; Gab < nirreps; Gab++) {
568
for(ab=0; ab < FAA.params->coltot[Gab]; ab++) {
569
A = FAA.params->colorb[Gab][ab][0];
570
B = FAA.params->colorb[Gab][ab][1];
571
Ga = FAA.params->rsym[A];
573
for(c=0; c < virtpi[Gc]; c++) {
575
cb = WmAfE.params->colidx[C][B];
576
W2[Ga][a][cb] = W1[Gab][ab][c];
581
ji = L2AAnew.params->rowidx[J][I];
583
for(Gd=0; Gd < nirreps; Gd++) {
584
Gcb = Gkd = Gk ^ Gd; /* assumes totally symmetric */
585
Ga = Gd ^ Gij; /* assumes totally symmetric */
587
kd = WmAfE.row_offset[Gkd][K];
588
WmAfE.matrix[Gkd] = dpd_block_matrix(virtpi[Gd], WmAfE.params->coltot[Gkd]);
589
dpd_buf4_mat_irrep_rd_block(&WmAfE, Gkd, kd, virtpi[Gd]);
590
Z = block_matrix(virtpi[Ga], virtpi[Gd]);
594
nlinks = WmAfE.params->coltot[Gkd];
596
if(nrows && ncols && nlinks)
597
C_DGEMM('n', 't', nrows, ncols, nlinks, 1.0, W2[Ga][0], nlinks,
598
WmAfE.matrix[Gkd][0], nlinks, 0.0, Z[0], ncols);
600
for(a=0; a < virtpi[Ga]; a++) {
602
for(d=0; d < virtpi[Gd]; d++) {
603
DD = vir_off[Gd] + d;
604
ad = L2AAnew.params->colidx[A][DD];
605
da = L2AAnew.params->colidx[DD][A];
606
L2AAnew.matrix[Gij][ji][ad] += -Z[a][d];
607
L2AAnew.matrix[Gij][ji][da] += Z[a][d];
611
dpd_free_block(WmAfE.matrix[Gkd], virtpi[Gd], WmAfE.params->coltot[Gkd]);
615
/* t_MIAB <--- +t_IJkABc W_JkMc */
616
/* t_IMAB <--- -t_IJkABc W_JkMc */
618
jk = WMnIe.params->rowidx[J][K];
620
for(Gm=0; Gm < nirreps; Gm++) {
621
Gab = Gmi = Gm ^ Gi; /* assumes totally symmetric */
622
Gc = Gab ^ Gijk; /* assumes totally symmetric */
624
mc = WMnIe.col_offset[Gjk][Gm];
626
nrows = FAA.params->coltot[Gab];
630
Z = dpd_block_matrix(nrows, ncols);
632
if(nrows && ncols && nlinks)
633
C_DGEMM('n', 't', nrows, ncols, nlinks, 1.0, W1[Gab][0], nlinks,
634
&(WMnIe.matrix[Gjk][jk][mc]), nlinks, 0.0, Z[0], ncols);
636
for(m=0; m < ncols; m++) {
638
mi = L2AAnew.params->rowidx[M][I];
639
im = L2AAnew.params->rowidx[I][M];
640
for(ab=0; ab < nrows; ab++) {
641
L2AAnew.matrix[Gmi][mi][ab] += Z[ab][m];
642
L2AAnew.matrix[Gmi][im][ab] -= Z[ab][m];
646
dpd_free_block(Z, nrows, ncols);
649
/* t_JkDc <-- 1/2 t_IJkABc W_IDAB */
650
/* t_KjCd <-- 1/2 t_IJkABc W_IDAB */
652
jk = L2ABnew.params->rowidx[J][K];
653
kj = L2ABnew.params->rowidx[K][J];
655
for(Gd=0; Gd < nirreps; Gd++) {
656
Gab = Gid = Gi ^ Gd; /* assumes totally symmetric */
657
Gc = Gab ^ Gijk; /* assumes totally symmetric */
659
id = WMAFE.row_offset[Gid][I];
660
WMAFE.matrix[Gid] = dpd_block_matrix(virtpi[Gd], WMAFE.params->coltot[Gid]);
661
dpd_buf4_mat_irrep_rd_block(&WMAFE, Gid, id, virtpi[Gd]);
662
Z = block_matrix(virtpi[Gc], virtpi[Gd]);
666
nlinks = WMAFE.params->coltot[Gid];
668
if(nrows && ncols && nlinks)
669
C_DGEMM('t', 't', nrows, ncols, nlinks, 0.5, W1[Gab][0], nrows,
670
WMAFE.matrix[Gid][0], nlinks, 0.0, Z[0], ncols);
672
for(c=0; c < virtpi[Gc]; c++) {
674
for(d=0; d < virtpi[Gd]; d++) {
675
DD = vir_off[Gd] + d;
676
dc = L2ABnew.params->colidx[DD][C];
677
cd = L2ABnew.params->colidx[C][DD];
678
L2ABnew.matrix[Gjk][jk][dc] += Z[c][d];
679
L2ABnew.matrix[Gjk][kj][cd] += Z[c][d];
684
dpd_free_block(WMAFE.matrix[Gid], virtpi[Gd], WMAFE.params->coltot[Gid]);
687
/* t_JkBd <-- t_IJkABc W_IdAc */
688
/* t_KjBd <-- t_IJkABc W_IdAc */
689
/* sort W1(AB,c) to W2(B,Ac) */
690
for(Gab=0; Gab < nirreps; Gab++) {
692
for(ab=0; ab < FAA.params->coltot[Gab]; ab++) {
693
A = FAA.params->colorb[Gab][ab][0];
694
B = FAA.params->colorb[Gab][ab][1];
695
Gb = FAA.params->ssym[B];
697
for(c=0; c < virtpi[Gc]; c++) {
699
ac = WMaFe.params->colidx[A][C];
700
W2[Gb][b][ac] = W1[Gab][ab][c];
705
jk = L2ABnew.params->rowidx[J][K];
706
kj = L2ABnew.params->rowidx[K][J];
708
for(Gd=0; Gd < nirreps; Gd++) {
709
Gac = Gid = Gi ^ Gd; /* assumes totally symmetric */
710
Gb = Gac ^ Gijk; /* assumes totally symmetric */
712
id = WMaFe.row_offset[Gid][I];
713
WMaFe.matrix[Gid] = dpd_block_matrix(virtpi[Gd], WMaFe.params->coltot[Gid]);
714
dpd_buf4_mat_irrep_rd_block(&WMaFe, Gid, id, virtpi[Gd]);
715
Z = block_matrix(virtpi[Gb], virtpi[Gd]);
719
nlinks = WMaFe.params->coltot[Gid];
721
if(nrows && ncols && nlinks)
722
C_DGEMM('n', 't', nrows, ncols, nlinks, 1.0, W2[Gb][0], nlinks,
723
WMaFe.matrix[Gid][0], nlinks, 0.0, Z[0], ncols);
725
for(b=0; b < virtpi[Gb]; b++) {
727
for(d=0; d < virtpi[Gd]; d++) {
728
DD = vir_off[Gd] + d;
729
bd = L2ABnew.params->colidx[B][DD];
730
db = L2ABnew.params->colidx[DD][B];
731
L2ABnew.matrix[Gjk][jk][bd] += Z[b][d];
732
L2ABnew.matrix[Gjk][kj][db] += Z[b][d];
736
dpd_free_block(WMaFe.matrix[Gid], virtpi[Gd], WMaFe.params->coltot[Gid]);
740
/* t_MkBc <-- 1/2 t_IJkABc W_IJMA */
741
/* sort W(AB,c) to W(A,Bc) */
742
for(Gab=0; Gab < nirreps; Gab++) {
743
Gc = Gab ^ Gijk; /* assumes totally symmetric */
744
for(ab=0; ab < FAA.params->coltot[Gab]; ab++ ){
745
A = FAA.params->colorb[Gab][ab][0];
746
B = FAA.params->colorb[Gab][ab][1];
747
Ga = FAA.params->rsym[A];
749
for(c=0; c < virtpi[Gc]; c++) {
751
bc = L2ABnew.params->colidx[B][C];
752
W2[Ga][a][bc] = W1[Gab][ab][c];
757
ij = WMNIE.params->rowidx[I][J];
759
for(Gm=0; Gm < nirreps; Gm++) {
760
Gbc = Gmk = Gm ^ Gk; /* assumes totally symmetric */
761
Ga = Gbc ^ Gijk; /* assumes totally symmetric */
763
ma = WMNIE.col_offset[Gij][Gm];
765
nrows = L2ABnew.params->coltot[Gmk];
769
Z = dpd_block_matrix(nrows, ncols);
771
if(nrows && ncols && nlinks)
772
C_DGEMM('t', 't', nrows, ncols, nlinks, 0.5, W2[Ga][0], nrows,
773
&(WMNIE.matrix[Gij][ij][ma]), nlinks, 0.0, Z[0], ncols);
775
for(m=0; m < occpi[Gm]; m++) {
777
mk = L2ABnew.params->rowidx[M][K];
778
km = L2ABnew.params->rowidx[K][M];
779
for(Gb=0; Gb < nirreps; Gb++) {
781
for(b=0; b < virtpi[Gb]; b++) {
783
for(c=0; c < virtpi[Gc]; c++) {
785
bc = L2ABnew.params->colidx[B][C];
786
cb = L2ABnew.params->colidx[C][B];
787
L2ABnew.matrix[Gmk][mk][bc] += Z[bc][m];
788
L2ABnew.matrix[Gmk][km][cb] += Z[bc][m];
794
dpd_free_block(Z, nrows, ncols);
797
/* t_ImBc <-- t_IJkABc W_kJmA */
798
/* sort W(AB,c) to W(A,Bc) */
799
for(Gab=0; Gab < nirreps; Gab++) {
800
Gc = Gab ^ Gijk; /* assumes totally symmetric */
801
for(ab=0; ab < FAA.params->coltot[Gab]; ab++ ){
802
A = FAA.params->colorb[Gab][ab][0];
803
B = FAA.params->colorb[Gab][ab][1];
804
Ga = FAA.params->rsym[A];
806
for(c=0; c < virtpi[Gc]; c++) {
808
bc = L2ABnew.params->colidx[B][C];
809
W2[Ga][a][bc] = W1[Gab][ab][c];
814
kj = WmNiE.params->rowidx[K][J];
816
for(Gm=0; Gm < nirreps; Gm++) {
817
Gbc = Gim = Gi ^ Gm; /* assumes totally symmetric */
818
Ga = Gbc ^ Gijk; /* assumes totally symmetric */
820
ma = WmNiE.col_offset[Gjk][Gm];
822
nrows = L2ABnew.params->coltot[Gim];
826
Z = dpd_block_matrix(nrows, ncols);
828
if(nrows && ncols && nlinks)
829
C_DGEMM('t', 't', nrows, ncols, nlinks, 1.0, W2[Ga][0], nrows,
830
&(WmNiE.matrix[Gjk][kj][ma]), nlinks, 0.0, Z[0], ncols);
832
for(m=0; m < occpi[Gm]; m++) {
834
im = L2ABnew.params->rowidx[I][M];
835
mi = L2ABnew.params->rowidx[M][I];
836
for(Gb=0; Gb < nirreps; Gb++) {
838
for(b=0; b < virtpi[Gb]; b++) {
840
for(c=0; c < virtpi[Gc]; c++) {
842
bc = L2ABnew.params->colidx[B][C];
843
cb = L2ABnew.params->colidx[C][B];
844
L2ABnew.matrix[Gim][im][bc] += Z[bc][m];
845
L2ABnew.matrix[Gim][mi][cb] += Z[bc][m];
851
dpd_free_block(Z, nrows, ncols);
854
/* Z_IDAB <-- L_IJkABc t_JkDc */
856
jk = T2AB.params->rowidx[J][K];
857
for(Gab=0; Gab < nirreps; Gab++) {
858
Gid = Gab; /* totally symmetric */
859
Gc = Gab ^ Gijk; /* totally symmetric */
863
ncols = ZIGDE.params->coltot[Gid];
866
dc = T2AB.col_offset[Gjk][Gd];
867
id = ZIGDE.row_offset[Gid][I];
868
ZIGDE.matrix[Gid] = dpd_block_matrix(nrows, ncols);
870
if(nrows && ncols && nlinks) {
871
dpd_buf4_mat_irrep_rd_block(&ZIGDE, Gid, id, nrows);
873
C_DGEMM('n', 't', nrows, ncols, nlinks, 1.0, &(T2AB.matrix[Gjk][jk][dc]), nlinks,
874
W1[Gab][0], nlinks, 1.0, ZIGDE.matrix[Gid][0], ncols);
876
dpd_buf4_mat_irrep_wrt_block(&ZIGDE, Gid, id, nrows);
879
dpd_free_block(ZIGDE.matrix[Gid], nrows, ncols);
882
/* ZkDCa <-- 1/2 L_ijKabC t_ijdb */
884
ij = T2AA.params->rowidx[I][J];
886
/* sort W(ab,C) to W(b,Ca) */
887
for(Gab=0; Gab < nirreps; Gab++) {
888
Gc = Gab ^ Gijk; /* assumes totally symmetric */
889
for(ab=0; ab < FAA.params->coltot[Gab]; ab++ ){
890
A = FAA.params->colorb[Gab][ab][0];
891
B = FAA.params->colorb[Gab][ab][1];
892
Gb = FAA.params->ssym[B];
894
for(c=0; c < virtpi[Gc]; c++) {
896
ca = ZIgDe.params->colidx[C][A];
897
W2[Gb][b][ca] = W1[Gab][ab][c];
902
for(Gb=0; Gb < nirreps; Gb++) {
903
Gd = Gb ^ Gij; /* totally symmetric */
904
Gca = Gkd = Gk ^ Gd; /* totally symmetric */
907
ncols = ZIgDe.params->coltot[Gkd];
910
db = T2AA.col_offset[Gij][Gd];
911
kd = ZIgDe.row_offset[Gkd][K];
912
ZIgDe.matrix[Gkd] = dpd_block_matrix(nrows, ncols);
913
dpd_buf4_mat_irrep_rd_block(&ZIgDe, Gkd, kd, nrows);
915
if(nrows && ncols && nlinks)
916
C_DGEMM('n', 'n', nrows, ncols, nlinks, 0.5, &(T2AA.matrix[Gij][ij][db]), nlinks,
917
W2[Gb][0], ncols, 1.0, ZIgDe.matrix[Gkd][0], ncols);
919
dpd_buf4_mat_irrep_wrt_block(&ZIgDe, Gkd, kd, nrows);
920
dpd_free_block(ZIgDe.matrix[Gkd], nrows, ncols);
923
/* Z_IdAc <-- L_IJkABc t_JkBd */
925
jk = T2AB.params->rowidx[J][K];
927
/* sort W(AB,c) to W(B,Ac) */
928
for(Gab=0; Gab < nirreps; Gab++) {
929
Gc = Gab ^ Gijk; /* assumes totally symmetric */
930
for(ab=0; ab < FAA.params->coltot[Gab]; ab++ ){
931
A = FAA.params->colorb[Gab][ab][0];
932
B = FAA.params->colorb[Gab][ab][1];
933
Gb = FAA.params->ssym[B];
935
for(c=0; c < virtpi[Gc]; c++) {
937
ac = ZIgDe.params->colidx[A][C];
938
W2[Gb][b][ac] = W1[Gab][ab][c];
943
for(Gb=0; Gb < nirreps; Gb++) {
944
Gd = Gb ^ Gjk; /* totally symmetric */
945
Gac = Gid = Gi ^ Gd; /* totally symmetric */
948
ncols = ZIgDe.params->coltot[Gid];
951
bd = T2AB.col_offset[Gjk][Gb];
952
id = ZIgDe.row_offset[Gid][I];
953
ZIgDe.matrix[Gid] = dpd_block_matrix(nrows, ncols);
954
dpd_buf4_mat_irrep_rd_block(&ZIgDe, Gid, id, nrows);
956
if(nrows && ncols && nlinks)
957
C_DGEMM('t', 'n', nrows, ncols, nlinks, 1.0, &(T2AB.matrix[Gjk][jk][bd]), nrows,
958
W2[Gb][0], ncols, 1.0, ZIgDe.matrix[Gid][0], ncols);
960
dpd_buf4_mat_irrep_wrt_block(&ZIgDe, Gid, id, nrows);
961
dpd_free_block(ZIgDe.matrix[Gid], nrows, ncols);
964
/* Z_JDAB <-- 1/2 L_IJkABc t_IkDc */
965
ik = T2AB.params->rowidx[I][K];
966
for(Gab=0; Gab < nirreps; Gab++) {
967
Gjd = Gab; /* totally symmetric */
968
Gc = Gab ^ Gijk; /* totally symmetric */
972
ncols = ZDMAE.params->coltot[Gjd];
975
dc = T2AB.col_offset[Gik][Gd];
976
jd = ZDMAE.row_offset[Gjd][J];
977
ZDMAE.matrix[Gjd] = dpd_block_matrix(nrows, ncols);
978
dpd_buf4_mat_irrep_rd_block(&ZDMAE, Gjd, jd, nrows);
980
if(nrows && ncols && nlinks)
981
C_DGEMM('n', 't', nrows, ncols, nlinks, 1.0, &(T2AB.matrix[Gik][ik][dc]), nlinks,
982
W1[Gab][0], nlinks, 1.0, ZDMAE.matrix[Gjd][0], ncols);
984
dpd_buf4_mat_irrep_wrt_block(&ZDMAE, Gjd, jd, nrows);
985
dpd_free_block(ZDMAE.matrix[Gjd], nrows, ncols);
988
/* Z_kDAc <-- 1/2 L_IJkABc t_IJDB */
989
ij = T2AA.params->rowidx[I][J];
990
/* sort W(AB,c) to W(B,Ac) */
991
for(Gab=0; Gab < nirreps; Gab++) {
992
Gc = Gab ^ Gijk; /* assumes totally symmetric */
993
for(ab=0; ab < FAA.params->coltot[Gab]; ab++ ){
994
A = FAA.params->colorb[Gab][ab][0];
995
B = FAA.params->colorb[Gab][ab][1];
996
Gb = FAA.params->ssym[B];
998
for(c=0; c < virtpi[Gc]; c++) {
1000
ac = ZDmAe.params->colidx[A][C];
1001
W2[Gb][b][ac] = W1[Gab][ab][c];
1006
for(Gb=0; Gb < nirreps; Gb++) {
1007
Gd = Gb ^ Gij; /* totally symmetric */
1008
Gac = Gkd = Gk ^ Gd; /* totally symmetric */
1011
ncols = ZDmAe.params->coltot[Gkd];
1012
nlinks = virtpi[Gb];
1014
db = T2AA.col_offset[Gij][Gd];
1015
kd = ZDmAe.row_offset[Gkd][K];
1016
ZDmAe.matrix[Gkd] = dpd_block_matrix(nrows, ncols);
1017
dpd_buf4_mat_irrep_rd_block(&ZDmAe, Gkd, kd, nrows);
1019
if(nrows && ncols && nlinks)
1020
C_DGEMM('n', 'n', nrows, ncols, nlinks, 0.5, &(T2AA.matrix[Gij][ij][db]), nlinks,
1021
W2[Gb][0], ncols, 1.0, ZDmAe.matrix[Gkd][0], ncols);
1023
dpd_buf4_mat_irrep_wrt_block(&ZDmAe, Gkd, kd, nrows);
1024
dpd_free_block(ZDmAe.matrix[Gkd], nrows, ncols);
1027
/* Z_iDCa <-- L_ijKabC t_KjDb */
1028
kj = T2AB.params->rowidx[K][J];
1029
/* sort W(AB,c) to W(B,Ca) */
1030
for(Gab=0; Gab < nirreps; Gab++) {
1031
Gc = Gab ^ Gijk; /* assumes totally symmetric */
1032
for(ab=0; ab < FAA.params->coltot[Gab]; ab++ ){
1033
A = FAA.params->colorb[Gab][ab][0];
1034
B = FAA.params->colorb[Gab][ab][1];
1035
Gb = FAA.params->ssym[B];
1036
b = B - vir_off[Gb];
1037
for(c=0; c < virtpi[Gc]; c++) {
1038
C = vir_off[Gc] + c;
1039
ca = ZDmAe.params->colidx[C][A];
1040
W2[Gb][b][ca] = W1[Gab][ab][c];
1045
for(Gb=0; Gb < nirreps; Gb++) {
1046
Gd = Gb ^ Gkj; /* totally symmetric */
1047
Gca = Gid = Gi ^ Gd; /* totally symmetric */
1050
ncols = ZDmAe.params->coltot[Gid];
1051
nlinks = virtpi[Gb];
1053
db = T2AB.col_offset[Gkj][Gd];
1054
id = ZDmAe.row_offset[Gid][I];
1055
ZDmAe.matrix[Gid] = dpd_block_matrix(nrows, ncols);
1056
dpd_buf4_mat_irrep_rd_block(&ZDmAe, Gid, id, nrows);
1058
if(nrows && ncols && nlinks)
1059
C_DGEMM('n', 'n', nrows, ncols, nlinks, 1.0, &(T2AB.matrix[Gkj][kj][db]), nlinks,
1060
W2[Gb][0], ncols, 1.0, ZDmAe.matrix[Gid][0], ncols);
1062
dpd_buf4_mat_irrep_wrt_block(&ZDmAe, Gid, id, nrows);
1063
dpd_free_block(ZDmAe.matrix[Gid], nrows, ncols);
1066
/* Z_KdCa <-- -1/2 L_ijKabC t_ijdb */
1067
ij = T2AA.params->rowidx[I][J];
1068
/* sort W(AB,c) to W(B,Ca) */
1069
for(Gab=0; Gab < nirreps; Gab++) {
1070
Gc = Gab ^ Gijk; /* assumes totally symmetric */
1071
for(ab=0; ab < FAA.params->coltot[Gab]; ab++ ){
1072
A = FAA.params->colorb[Gab][ab][0];
1073
B = FAA.params->colorb[Gab][ab][1];
1074
Gb = FAA.params->ssym[B];
1075
b = B - vir_off[Gb];
1076
for(c=0; c < virtpi[Gc]; c++) {
1077
C = vir_off[Gc] + c;
1078
ca = ZdMAe.params->colidx[C][A];
1079
W2[Gb][b][ca] = W1[Gab][ab][c];
1084
for(Gb=0; Gb < nirreps; Gb++) {
1085
Gd = Gb ^ Gij; /* totally symmetric */
1086
Gca = Gkd = Gk ^ Gd; /* totally symmetric */
1089
ncols = ZdMAe.params->coltot[Gkd];
1090
nlinks = virtpi[Gb];
1092
db = T2AA.col_offset[Gij][Gd];
1093
kd = ZdMAe.row_offset[Gkd][K];
1094
ZdMAe.matrix[Gkd] = dpd_block_matrix(nrows, ncols);
1095
dpd_buf4_mat_irrep_rd_block(&ZdMAe, Gkd, kd, nrows);
1097
if(nrows && ncols && nlinks)
1098
C_DGEMM('n', 'n', nrows, ncols, nlinks, -0.5, &(T2AA.matrix[Gij][ij][db]), nlinks,
1099
W2[Gb][0], ncols, 1.0, ZdMAe.matrix[Gkd][0], ncols);
1101
dpd_buf4_mat_irrep_wrt_block(&ZdMAe, Gkd, kd, nrows);
1102
dpd_free_block(ZdMAe.matrix[Gkd], nrows, ncols);
1105
/* Z_JdAc <-- L_IJkABc t_IkBd */
1106
ik = T2AB.params->rowidx[I][K];
1108
/* sort W(AB,c) to W(B,Ca) */
1109
for(Gab=0; Gab < nirreps; Gab++) {
1110
Gc = Gab ^ Gijk; /* assumes totally symmetric */
1111
for(ab=0; ab < FAA.params->coltot[Gab]; ab++ ){
1112
A = FAA.params->colorb[Gab][ab][0];
1113
B = FAA.params->colorb[Gab][ab][1];
1114
Gb = FAA.params->ssym[B];
1115
b = B - vir_off[Gb];
1116
for(c=0; c < virtpi[Gc]; c++) {
1117
C = vir_off[Gc] + c;
1118
ac = ZdMAe.params->colidx[A][C];
1119
W2[Gb][b][ac] = W1[Gab][ab][c];
1124
for(Gb=0; Gb < nirreps; Gb++) {
1125
Gd = Gb ^ Gik; /* totally symmetric */
1126
Gca = Gjd = Gj ^ Gd; /* totally symmetric */
1129
ncols = ZdMAe.params->coltot[Gjd];
1130
nlinks = virtpi[Gb];
1132
bd = T2AB.col_offset[Gik][Gb];
1133
jd = ZdMAe.row_offset[Gjd][J];
1134
ZdMAe.matrix[Gjd] = dpd_block_matrix(nrows, ncols);
1135
dpd_buf4_mat_irrep_rd_block(&ZdMAe, Gjd, jd, nrows);
1137
if(nrows && ncols && nlinks)
1138
C_DGEMM('t', 'n', nrows, ncols, nlinks, 1.0, &(T2AB.matrix[Gik][ik][bd]), nrows,
1139
W2[Gb][0], ncols, 1.0, ZdMAe.matrix[Gjd][0], ncols);
1141
dpd_buf4_mat_irrep_wrt_block(&ZdMAe, Gjd, jd, nrows);
1142
dpd_free_block(ZdMAe.matrix[Gjd], nrows, ncols);
1145
/* Z_IJAM <-- -1/2 L_IJkABc t_MkBc */
1146
/* sort W(AB,C) to W(A,BC) */
1147
for(Gab=0; Gab < nirreps; Gab++) {
1149
for(ab=0; ab < FAA.params->coltot[Gab]; ab++) {
1150
A = FAA.params->colorb[Gab][ab][0];
1151
B = FAA.params->colorb[Gab][ab][1];
1152
Ga = FAA.params->rsym[A];
1153
a = A - vir_off[Ga];
1154
for(c=0; c < virtpi[Gc]; c++) {
1155
C = vir_off[Gc] + c;
1156
bc = FAA.params->colidx[B][C];
1157
W2[Ga][a][bc] = W1[Gab][ab][c];
1162
ij = ZLMAO.params->rowidx[I][J];
1164
for(Gm=0; Gm < nirreps; Gm++) {
1165
Gbc = Gmk = Gm ^ Gk; /* totally symmetric */
1166
Ga = Gij ^ Gm; /* totally symmetric */
1169
ncols = T2AB.params->coltot[Gmk];
1171
for(m=0; m < occpi[Gm]; m++) {
1172
M = occ_off[Gm] + m;
1173
mk = T2AB.params->rowidx[M][K];
1174
am = ZLMAO.col_offset[Gij][Ga] + m;
1177
C_DGEMV('n', nrows, ncols, -1.0, W2[Ga][0], ncols, T2AB.matrix[Gmk][mk], 1,
1178
1.0, &(ZLMAO.matrix[Gij][ij][am]), occpi[Gm]);
1182
/* Z_KiCm <-- -1/2 L_ijKabC t_mjab */
1183
ki = ZLmAo.params->rowidx[K][I];
1185
for(Gm=0; Gm < nirreps; Gm++) {
1186
Gab = Gmj = Gm ^ Gj; /* totally symmetric */
1187
Gc = Gm ^ Gki; /* totally symmetric */
1189
nrows = T2AA.params->coltot[Gmj];
1192
for(m=0; m < occpi[Gm]; m++) {
1193
M = occ_off[Gm] + m;
1194
mj = T2AA.params->rowidx[M][J];
1195
cm = ZLmAo.col_offset[Gki][Gc] + m;
1198
C_DGEMV('t', nrows, ncols, -0.5, W1[Gab][0], ncols, T2AA.matrix[Gmj][mj], 1,
1199
1.0, &(ZLmAo.matrix[Gki][ki][cm]), occpi[Gm]);
1203
/* Z_IkAm <-- - L_IJkABc t_mJcB */
1204
/* sort W(AB,C) to W(A,CB) */
1205
for(Gab=0; Gab < nirreps; Gab++) {
1207
for(ab=0; ab < FAA.params->coltot[Gab]; ab++) {
1208
A = FAA.params->colorb[Gab][ab][0];
1209
B = FAA.params->colorb[Gab][ab][1];
1210
Ga = FAA.params->rsym[A];
1211
a = A - vir_off[Ga];
1212
for(c=0; c < virtpi[Gc]; c++) {
1213
C = vir_off[Gc] + c;
1214
cb = FAA.params->colidx[C][B];
1215
W2[Ga][a][cb] = W1[Gab][ab][c];
1220
ik = ZLmAo.params->rowidx[I][K];
1222
for(Gm=0; Gm < nirreps; Gm++) {
1223
Gbc = Gmj = Gm ^ Gj; /* totally symmetric */
1224
Ga = Gm ^ Gik; /* totally symmetric */
1227
ncols = T2AB.params->coltot[Gmj];
1229
for(m=0; m < occpi[Gm]; m++) {
1230
M = occ_off[Gm] + m;
1231
mj = T2AB.params->rowidx[M][J];
1232
am = ZLmAo.col_offset[Gik][Ga] + m;
1235
C_DGEMV('n', nrows, ncols, -1.0, W2[Ga][0], ncols, T2AB.matrix[Gmj][mj], 1,
1236
1.0, &(ZLmAo.matrix[Gik][ik][am]), occpi[Gm]);
1240
/* Z_IJMB <-- - L_IJkABc t_MkAc */
1241
/* sort W(AB,C) to W(B,AC) */
1242
for(Gab=0; Gab < nirreps; Gab++) {
1244
for(ab=0; ab < FAA.params->coltot[Gab]; ab++) {
1245
A = FAA.params->colorb[Gab][ab][0];
1246
B = FAA.params->colorb[Gab][ab][1];
1247
Gb = FAA.params->ssym[B];
1248
b = B - vir_off[Gb];
1249
for(c=0; c < virtpi[Gc]; c++) {
1250
C = vir_off[Gc] + c;
1251
ac = FAA.params->colidx[A][C];
1252
W2[Gb][b][ac] = W1[Gab][ab][c];
1257
ij = ZIMLE.params->rowidx[I][J];
1259
for(Gm=0; Gm < nirreps; Gm++) {
1260
Gb = Gm ^ Gij; /* totally symmetric */
1264
ncols = T2AB.params->coltot[Gmk];
1266
for(m=0; m < occpi[Gm]; m++) {
1267
M = occ_off[Gm] + m;
1268
mk = T2AB.params->rowidx[M][K];
1269
mb = ZIMLE.col_offset[Gij][Gm] + m * virtpi[Gb];
1272
C_DGEMV('n', nrows, ncols, -1.0, W2[Gb][0], ncols, T2AB.matrix[Gmk][mk], 1,
1273
1.0, &(ZIMLE.matrix[Gij][ij][mb]), 1);
1278
/* Z_IkMc <-- -1/2 L_IJkABc t_MJAB */
1279
ik = ZImLe.params->rowidx[I][K];
1280
for(Gm=0; Gm < nirreps; Gm++) {
1281
Gc = Gm ^ Gik ; /* totally symmetric */
1282
Gab = Gmj = Gm ^ Gj; /* totally symmetric */
1284
nrows = T2AA.params->coltot[Gmj];
1287
for(m=0; m < occpi[Gm]; m++) {
1288
M = occ_off[Gm] + m;
1289
mj = T2AA.params->rowidx[M][J];
1290
mc = ZImLe.col_offset[Gik][Gm] + m * virtpi[Gc];
1293
C_DGEMV('t', nrows, ncols, -0.5, W1[Gab][0], ncols, T2AA.matrix[Gmj][mj], 1,
1294
1.0, &(ZImLe.matrix[Gik][ik][mc]), 1);
1298
/* Z_KiMa <-- - L_ijKabC t_MjCb */
1299
/* sort W(AB,C) to W(A,CB) */
1300
for(Gab=0; Gab < nirreps; Gab++) {
1302
for(ab=0; ab < FAA.params->coltot[Gab]; ab++) {
1303
A = FAA.params->colorb[Gab][ab][0];
1304
B = FAA.params->colorb[Gab][ab][1];
1305
Ga = FAA.params->rsym[A];
1306
a = A - vir_off[Ga];
1307
for(c=0; c < virtpi[Gc]; c++) {
1308
C = vir_off[Gc] + c;
1309
cb = FAA.params->colidx[C][B];
1310
W2[Ga][a][cb] = W1[Gab][ab][c];
1315
ki = ZImLe.params->rowidx[K][I];
1316
for(Gm=0; Gm < nirreps; Gm++) {
1317
Ga = Gm ^ Gki; /* totally symmetric */
1321
ncols = T2AB.params->coltot[Gmj];
1323
for(m=0; m < occpi[Gm]; m++) {
1324
M = occ_off[Gm] + m;
1325
mj = T2AB.params->rowidx[M][J];
1326
ma = ZImLe.col_offset[Gki][Gm] + m * virtpi[Ga];
1329
C_DGEMV('n', nrows, ncols, -1.0, W2[Ga][0], ncols, T2AB.matrix[Gmj][mj], 1,
1330
1.0, &(ZImLe.matrix[Gki][ki][ma]), 1);
1334
/* Z_KimC <-- 1/2 L_ijKabC t_mjab */
1335
ki = ZImlE.params->rowidx[K][I];
1336
for(Gm=0; Gm < nirreps; Gm++) {
1337
Gc = Gm ^ Gki; /* totally symmetric */
1338
Gab = Gmj = Gm ^ Gj;
1340
nrows = T2AA.params->coltot[Gmj];
1343
for(m=0; m < occpi[Gm]; m++) {
1344
M = occ_off[Gm] + m;
1345
mj = T2AA.params->rowidx[M][J];
1346
mc = ZImlE.col_offset[Gki][Gm] + m * virtpi[Gc];
1349
C_DGEMV('t', nrows, ncols, 0.5, W1[Gab][0], ncols, T2AA.matrix[Gmj][mj], 1,
1350
1.0, &(ZImlE.matrix[Gki][ki][mc]), 1);
1354
/* Z_IkmB <-- - l_IJkABc t_JmAc */
1355
ik = ZImlE.params->rowidx[I][K];
1356
/* sort W(AB,C) to W(B,AC) */
1357
for(Gab=0; Gab < nirreps; Gab++) {
1359
for(ab=0; ab < FAA.params->coltot[Gab]; ab++) {
1360
A = FAA.params->colorb[Gab][ab][0];
1361
B = FAA.params->colorb[Gab][ab][1];
1362
Gb = FAA.params->ssym[B];
1363
b = B - vir_off[Gb];
1364
for(c=0; c < virtpi[Gc]; c++) {
1365
C = vir_off[Gc] + c;
1366
ac = FAA.params->colidx[A][C];
1367
W2[Gb][b][ac] = W1[Gab][ab][c];
1372
for(Gm=0; Gm < nirreps; Gm++) {
1373
Gb = Gm ^ Gik; /* totally symmetric */
1377
ncols = T2AB.params->coltot[Gjm];
1379
for(m=0; m < occpi[Gm]; m++) {
1380
M = occ_off[Gm] + m;
1381
jm = T2AB.params->rowidx[J][M];
1382
mb = ZImlE.col_offset[Gki][Gm] + m * virtpi[Gb];
1385
C_DGEMV('n', nrows, ncols, -1.0, W2[Gb][0], ncols, T2AB.matrix[Gjm][jm], 1,
1386
1.0, &(ZImlE.matrix[Gik][ik][mb]), 1);
1394
for(Gab=0; Gab < nirreps; Gab++) {
1395
Gc = Gab ^ Gijk; /* totally symmetric */
1396
dpd_free_block(W1[Gab], FAA.params->coltot[Gab], virtpi[Gc]);
1398
for(Ga=0; Ga < nirreps; Ga++) {
1399
Gcb = Ga ^ Gijk; /* assumes totally symmetric */
1400
dpd_free_block(W2[Ga], virtpi[Ga], WmAfE.params->coltot[Gcb]);
1410
dpd_buf4_close(&EAA);
1411
dpd_buf4_close(&EAB);
1412
dpd_buf4_close(&EBA);
1413
dpd_buf4_close(&FAA);
1414
dpd_buf4_close(&FAB);
1415
dpd_buf4_close(&FBA);
1416
dpd_buf4_close(&L2AA);
1417
dpd_buf4_close(&L2AB);
1418
dpd_buf4_close(&L2BA);
1420
dpd_file2_close(&fIJ);
1421
dpd_file2_close(&fAB);
1422
dpd_file2_close(&fij);
1423
dpd_file2_close(&fab);
1425
dpd_file2_close(&FME);
1426
dpd_file2_close(&Fme);
1427
dpd_file2_close(&LIA);
1428
dpd_file2_close(&Lia);
1430
dpd_buf4_close(&DAAints);
1431
dpd_buf4_close(&DABints);
1432
dpd_buf4_close(&LIJAB);
1433
dpd_buf4_close(&LIjAb);
1435
dpd_buf4_close(&WmAfE);
1436
dpd_buf4_close(&WMAFE);
1437
dpd_buf4_close(&WMaFe);
1439
for(h=0; h < nirreps; h++) dpd_buf4_mat_irrep_close(&WMnIe, h);
1440
for(h=0; h < nirreps; h++) dpd_buf4_mat_irrep_close(&WMNIE, h);
1441
for(h=0; h < nirreps; h++) dpd_buf4_mat_irrep_close(&WmNiE, h);
1442
dpd_buf4_close(&WMnIe);
1443
dpd_buf4_close(&WMNIE);
1444
dpd_buf4_close(&WmNiE);
1446
dpd_buf4_close(&ZIgDe);
1447
dpd_buf4_close(&ZIGDE);
1449
dpd_buf4_close(&ZDMAE);
1450
dpd_buf4_close(&ZDmAe);
1451
dpd_buf4_close(&ZdMAe);
1453
for(h=0; h < nirreps; h++) {
1454
dpd_buf4_mat_irrep_wrt(&ZLMAO, h);
1455
dpd_buf4_mat_irrep_close(&ZLMAO, h);
1457
dpd_buf4_close(&ZLMAO);
1458
for(h=0; h < nirreps; h++) {
1459
dpd_buf4_mat_irrep_wrt(&ZLmAo, h);
1460
dpd_buf4_mat_irrep_close(&ZLmAo, h);
1462
dpd_buf4_close(&ZLmAo);
1464
for(h=0; h < nirreps; h++) {
1465
dpd_buf4_mat_irrep_wrt(&ZIMLE, h);
1466
dpd_buf4_mat_irrep_close(&ZIMLE, h);
1468
dpd_buf4_close(&ZIMLE);
1469
for(h=0; h < nirreps; h++) {
1470
dpd_buf4_mat_irrep_wrt(&ZImLe, h);
1471
dpd_buf4_mat_irrep_close(&ZImLe, h);
1473
dpd_buf4_close(&ZImLe);
1474
for(h=0; h < nirreps; h++) {
1475
dpd_buf4_mat_irrep_wrt(&ZImlE, h);
1476
dpd_buf4_mat_irrep_close(&ZImlE, h);
1478
dpd_buf4_close(&ZImlE);
1480
for(h=0; h < nirreps; h++) dpd_buf4_mat_irrep_close(&T2AB, h);
1481
for(h=0; h < nirreps; h++) dpd_buf4_mat_irrep_close(&T2AA, h);
1482
dpd_buf4_close(&T2AB);
1483
dpd_buf4_close(&T2AA);
1485
for(h=0; h < nirreps; h++) {
1486
dpd_buf4_mat_irrep_wrt(&L2AAnew, h);
1487
dpd_buf4_mat_irrep_close(&L2AAnew, h);
1489
dpd_buf4_init(&D2, CC_DENOM, 0, 0, 5, 0, 5, 0, "dIjAb");
1490
dpd_buf4_dirprd(&D2, &L2AAnew);
1491
dpd_buf4_close(&D2);
1492
dpd_buf4_init(&L2, CC_LAMBDA, 0, 0, 5, 2, 7, 0, "New LIJAB");
1493
dpd_buf4_axpy(&L2AAnew, &L2, 1);
1494
dpd_buf4_close(&L2);
1495
dpd_buf4_close(&L2AAnew);
1497
for(h=0; h < nirreps; h++) {
1498
dpd_buf4_mat_irrep_wrt(&L2ABnew, h);
1499
dpd_buf4_mat_irrep_close(&L2ABnew, h);
1501
dpd_buf4_init(&D2, CC_DENOM, 0, 0, 5, 0, 5, 0, "dIjAb");
1502
dpd_buf4_dirprd(&D2, &L2ABnew);
1503
dpd_buf4_close(&D2);
1504
dpd_buf4_init(&L2, CC_LAMBDA, 0, 0, 5, 0, 5, 0, "New LIjAb");
1505
dpd_buf4_axpy(&L2ABnew, &L2, 1);
1506
dpd_buf4_close(&L2);
1507
dpd_buf4_close(&L2ABnew);
1509
/* Spin adaptation will remove this. And yes, this means that all the above
1510
calculations for LIJAB were pointless... -TDC */
1511
dpd_buf4_init(&L2, CC_LAMBDA, 0, 2, 7, 0, 5, 1, "New LIjAb");
1512
dpd_buf4_copy(&L2, CC_LAMBDA, "New LIJAB");
1513
dpd_buf4_copy(&L2, CC_LAMBDA, "New Lijab");
1514
dpd_buf4_close(&L2);
added
removed
src/bin/cclambda/cc3_l3l2.c
