3
#include <libdpd/dpd.h>
7
double ET_UHF_ABB_noddy(void)
11
int Gi, Gj, Gk, Ga, Gb, Gc, Ge, Gm;
12
int Gji, Gij, Gjk, Gik, Gbc, Gac, Gba, Gab;
13
int I, J, K, A, B, C, E, M;
14
int i, j, k, a, b, c, e, m;
15
int ij, ji, ik, ki, jk, kj;
16
int ab, ba, ac, ca, bc, cb;
17
int ae, be, eb, ce, ec, ie, je, ke;
18
int im, jm, mj, km, mk, ma, mb, mc;
19
int *aoccpi, *avirtpi, *aocc_off, *avir_off;
20
int *boccpi, *bvirtpi, *bocc_off, *bvir_off;
21
double value_c, value_d, denom, ET_ABB;
22
double t_ijae, t_ijeb, t_ijec, t_jkbe, t_jkce, t_ikae, t_ikeb, t_ikec;
23
double F_kebc, F_keca, F_keba, F_ieac, F_ieab, F_jebc, F_jeca, F_jeba;
24
double t_imac, t_imab, t_jmbc, t_mjac, t_mjab, t_kmbc, t_mkac, t_mkab;
25
double E_jkmb, E_jkmc, E_kima, E_ikmb, E_ikmc, E_jima, E_ijmb, E_ijmc;
26
double t_ia, t_jb, t_jc, t_kb, t_kc;
27
double D_jkbc, D_ikac, D_ikab, D_ijac, D_ijab;
29
dpdbuf4 FBBints, FABints, FBAints;
30
dpdbuf4 EBBints, EABints, EBAints;
31
dpdbuf4 DBBints, DABints;
32
dpdfile2 T1A, T1B, fIJ, fij, fAB, fab;
34
nirreps = moinfo.nirreps;
35
aoccpi = moinfo.aoccpi;
36
avirtpi = moinfo.avirtpi;
37
aocc_off = moinfo.aocc_off;
38
avir_off = moinfo.avir_off;
39
boccpi = moinfo.boccpi;
40
bvirtpi = moinfo.bvirtpi;
41
bocc_off = moinfo.bocc_off;
42
bvir_off = moinfo.bvir_off;
44
dpd_file2_init(&fIJ, CC_OEI, 0, 0, 0, "fIJ");
45
dpd_file2_init(&fij, CC_OEI, 0, 2, 2, "fij");
46
dpd_file2_init(&fAB, CC_OEI, 0, 1, 1, "fAB");
47
dpd_file2_init(&fab, CC_OEI, 0, 3, 3, "fab");
48
dpd_file2_mat_init(&fIJ);
49
dpd_file2_mat_init(&fij);
50
dpd_file2_mat_init(&fAB);
51
dpd_file2_mat_init(&fab);
52
dpd_file2_mat_rd(&fIJ);
53
dpd_file2_mat_rd(&fij);
54
dpd_file2_mat_rd(&fAB);
55
dpd_file2_mat_rd(&fab);
57
dpd_file2_init(&T1A, CC_OEI, 0, 0, 1, "tIA");
58
dpd_file2_mat_init(&T1A);
59
dpd_file2_mat_rd(&T1A);
60
dpd_file2_init(&T1B, CC_OEI, 0, 2, 3, "tia");
61
dpd_file2_mat_init(&T1B);
62
dpd_file2_mat_rd(&T1B);
64
dpd_buf4_init(&T2BB, CC_TAMPS, 0, 10, 15, 12, 17, 0, "tijab");
65
dpd_buf4_init(&T2AB, CC_TAMPS, 0, 22, 28, 22, 28, 0, "tIjAb");
67
dpd_buf4_init(&FBBints, CC_FINTS, 0, 30, 15, 30, 15, 1, "F <ia|bc>");
68
dpd_buf4_init(&FABints, CC_FINTS, 0, 24, 28, 24, 28, 0, "F <Ia|Bc>");
69
dpd_buf4_init(&FBAints, CC_FINTS, 0, 27, 29, 27, 29, 0, "F <iA|bC>");
71
dpd_buf4_init(&EBBints, CC_EINTS, 0, 10, 30, 12, 30, 0, "E <ij||ka> (i>j,ka)");
72
dpd_buf4_init(&EABints, CC_EINTS, 0, 22, 24, 22, 24, 0, "E <Ij|Ka>");
73
dpd_buf4_init(&EBAints, CC_EINTS, 0, 23, 27, 23, 27, 0, "E <iJ|kA>");
75
dpd_buf4_init(&DBBints, CC_DINTS, 0, 10, 15, 10, 15, 0, "D <ij||ab>");
76
dpd_buf4_init(&DABints, CC_DINTS, 0, 22, 28, 22, 28, 0, "D <Ij|Ab>");
78
for(h=0; h < nirreps; h++) {
79
dpd_buf4_mat_irrep_init(&T2BB, h);
80
dpd_buf4_mat_irrep_rd(&T2BB, h);
82
dpd_buf4_mat_irrep_init(&T2AB, h);
83
dpd_buf4_mat_irrep_rd(&T2AB, h);
85
dpd_buf4_mat_irrep_init(&FBBints, h);
86
dpd_buf4_mat_irrep_rd(&FBBints, h);
88
dpd_buf4_mat_irrep_init(&FABints, h);
89
dpd_buf4_mat_irrep_rd(&FABints, h);
91
dpd_buf4_mat_irrep_init(&FBAints, h);
92
dpd_buf4_mat_irrep_rd(&FBAints, h);
94
dpd_buf4_mat_irrep_init(&EBBints, h);
95
dpd_buf4_mat_irrep_rd(&EBBints, h);
97
dpd_buf4_mat_irrep_init(&EABints, h);
98
dpd_buf4_mat_irrep_rd(&EABints, h);
100
dpd_buf4_mat_irrep_init(&EBAints, h);
101
dpd_buf4_mat_irrep_rd(&EBAints, h);
103
dpd_buf4_mat_irrep_init(&DBBints, h);
104
dpd_buf4_mat_irrep_rd(&DBBints, h);
106
dpd_buf4_mat_irrep_init(&DABints, h);
107
dpd_buf4_mat_irrep_rd(&DABints, h);
113
for(Gi=0; Gi < nirreps; Gi++) {
114
for(Gj=0; Gj < nirreps; Gj++) {
115
for(Gk=0; Gk < nirreps; Gk++) {
120
for(i=0; i < aoccpi[Gi]; i++) {
121
I = aocc_off[Gi] + i;
122
for(j=0; j < boccpi[Gj]; j++) {
123
J = bocc_off[Gj] + j;
124
for(k=0; k < boccpi[Gk]; k++) {
125
K = bocc_off[Gk] + k;
127
ij = EABints.params->rowidx[I][J];
128
ji = EBAints.params->rowidx[J][I];
129
jk = EBBints.params->rowidx[J][K];
130
kj = EBBints.params->rowidx[K][J];
131
ik = EABints.params->rowidx[I][K];
132
ki = EBAints.params->rowidx[K][I];
136
for(Ga=0; Ga < nirreps; Ga++) {
137
for(Gb=0; Gb < nirreps; Gb++) {
138
Gc = Gi ^ Gj ^ Gk ^ Ga ^ Gb;
144
for(a=0; a < avirtpi[Ga]; a++) {
145
A = avir_off[Ga] + a;
146
for(b=0; b < bvirtpi[Gb]; b++) {
147
B = bvir_off[Gb] + b;
148
for(c=0; c < bvirtpi[Gc]; c++) {
149
C = bvir_off[Gc] + c;
151
ab = FABints.params->colidx[A][B];
152
ba = FBAints.params->colidx[B][A];
153
bc = FBBints.params->colidx[B][C];
154
cb = FBBints.params->colidx[C][B];
155
ac = FABints.params->colidx[A][C];
156
ca = FBAints.params->colidx[C][A];
160
/** <ov||vv> --> connected triples **/
162
/* +t_jkbe * F_IeAc */
164
for(e=0; e < bvirtpi[Ge]; e++) {
165
E = bvir_off[Ge] + e;
167
be = T2BB.params->colidx[B][E];
168
ie = FABints.params->rowidx[I][E];
170
t_jkbe = F_ieac = 0.0;
172
if(T2BB.params->rowtot[Gjk] && T2BB.params->coltot[Gjk])
173
t_jkbe = T2BB.matrix[Gjk][jk][be];
175
if(FABints.params->rowtot[Gac] && FABints.params->coltot[Gac])
176
F_ieac = FABints.matrix[Gac][ie][ac];
178
value_c += t_jkbe * F_ieac;
181
/* -t_jkce * F_IeAb */
183
for(e=0; e < bvirtpi[Ge]; e++) {
184
E = bvir_off[Ge] + e;
186
ce = T2BB.params->colidx[C][E];
187
ie = FABints.params->rowidx[I][E];
189
t_jkce = F_ieab = 0.0;
191
if(T2BB.params->rowtot[Gjk] && T2BB.params->coltot[Gjk])
192
t_jkce = T2BB.matrix[Gjk][jk][ce];
194
if(FABints.params->rowtot[Gba] && FABints.params->coltot[Gba])
195
F_ieab = FABints.matrix[Gba][ie][ab];
197
value_c -= t_jkce * F_ieab;
200
/* +t_IkAe * F_jebc */
202
for(e=0; e < bvirtpi[Ge]; e++) {
203
E = bvir_off[Ge] + e;
205
ae = T2AB.params->colidx[A][E];
206
je = FBBints.params->rowidx[J][E];
208
t_ikae = F_jebc = 0.0;
210
if(T2AB.params->rowtot[Gik] && T2AB.params->coltot[Gik])
211
t_ikae = T2AB.matrix[Gik][ik][ae];
213
if(FBBints.params->rowtot[Gbc] && FBBints.params->coltot[Gbc])
214
F_jebc = FBBints.matrix[Gbc][je][bc];
216
value_c += t_ikae * F_jebc;
219
/* -t_IkEb * F_jEcA */
221
for(e=0; e < avirtpi[Ge]; e++) {
222
E = avir_off[Ge] + e;
224
eb = T2AB.params->colidx[E][B];
225
je = FBAints.params->rowidx[J][E];
227
t_ikeb = F_jeca = 0.0;
229
if(T2AB.params->rowtot[Gik] && T2AB.params->coltot[Gik])
230
t_ikeb = T2AB.matrix[Gik][ik][eb];
232
if(FBAints.params->rowtot[Gac] && FBAints.params->coltot[Gac])
233
F_jeca = FBAints.matrix[Gac][je][ca];
235
value_c -= t_ikeb * F_jeca;
238
/* +t_IkEc * F_jEbA */
240
for(e=0; e < avirtpi[Ge]; e++) {
241
E = avir_off[Ge] + e;
243
ec = T2AB.params->colidx[E][C];
244
je = FBAints.params->rowidx[J][E];
246
t_ikec = F_jeba = 0.0;
248
if(T2AB.params->rowtot[Gik] && T2AB.params->coltot[Gik])
249
t_ikec = T2AB.matrix[Gik][ik][ec];
251
if(FBAints.params->rowtot[Gba] && FBAints.params->coltot[Gba])
252
F_jeba = FBAints.matrix[Gba][je][ba];
254
value_c += t_ikec * F_jeba;
257
/* -t_IjAe * F_kebc */
259
for(e=0; e < bvirtpi[Ge]; e++) {
260
E = bvir_off[Ge] + e;
262
ae = T2AB.params->colidx[A][E];
263
ke = FBBints.params->rowidx[K][E];
265
t_ijae = F_kebc = 0.0;
267
if(T2AB.params->rowtot[Gij] && T2AB.params->coltot[Gij])
268
t_ijae = T2AB.matrix[Gij][ij][ae];
270
if(FBBints.params->rowtot[Gbc] && FBBints.params->coltot[Gbc])
271
F_kebc = FBBints.matrix[Gbc][ke][bc];
273
value_c -= t_ijae * F_kebc;
276
/* +t_IjEb * F_kEcA */
278
for(e=0; e < avirtpi[Ge]; e++) {
279
E = avir_off[Ge] + e;
281
eb = T2AB.params->colidx[E][B];
282
ke = FBAints.params->rowidx[K][E];
284
t_ijeb = F_keca = 0.0;
286
if(T2AB.params->rowtot[Gij] && T2AB.params->coltot[Gij])
287
t_ijeb = T2AB.matrix[Gij][ij][eb];
289
if(FBAints.params->rowtot[Gac] && FBAints.params->coltot[Gac])
290
F_keca = FBAints.matrix[Gac][ke][ca];
292
value_c += t_ijeb * F_keca;
295
/* -t_IjEc * F_kEbA */
297
for(e=0; e < avirtpi[Ge]; e++) {
298
E = avir_off[Ge] + e;
300
ec = T2AB.params->colidx[E][C];
301
ke = FBAints.params->rowidx[K][E];
303
t_ijec = F_keba = 0.0;
305
if(T2AB.params->rowtot[Gij] && T2AB.params->coltot[Gij])
306
t_ijec = T2AB.matrix[Gij][ij][ec];
308
if(FBAints.params->rowtot[Gba] && FBAints.params->coltot[Gba])
309
F_keba = FBAints.matrix[Gba][ke][ba];
311
value_c -= t_ijec * F_keba;
314
/** <oo||ov> --> connected triples **/
316
/* +t_ImAc * E_jkmb */
318
for(m=0; m < boccpi[Gm]; m++) {
319
M = bocc_off[Gm] + m;
321
im = T2AB.params->rowidx[I][M];
322
mb = EBBints.params->colidx[M][B];
324
t_imac = E_jkmb = 0.0;
326
if(T2AB.params->rowtot[Gac] && T2AB.params->coltot[Gac])
327
t_imac = T2AB.matrix[Gac][im][ac];
329
if(EBBints.params->rowtot[Gjk] && EBBints.params->coltot[Gjk])
330
E_jkmb = EBBints.matrix[Gjk][jk][mb];
332
value_c += t_imac * E_jkmb;
335
/* -t_ImAb * E_jkmc */
337
for(m=0; m < boccpi[Gm]; m++) {
338
M = bocc_off[Gm] + m;
340
im = T2AB.params->rowidx[I][M];
341
mc = EBBints.params->colidx[M][C];
343
t_imab = E_jkmc = 0.0;
345
if(T2AB.params->rowtot[Gba] && T2AB.params->coltot[Gba])
346
t_imab = T2AB.matrix[Gba][im][ab];
348
if(EBBints.params->rowtot[Gjk] && EBBints.params->coltot[Gjk])
349
E_jkmc = EBBints.matrix[Gjk][jk][mc];
351
value_c -= t_imab * E_jkmc;
354
/* -t_jmbc * E_kImA */
356
for(m=0; m < boccpi[Gm]; m++) {
357
M = bocc_off[Gm] + m;
359
jm = T2BB.params->rowidx[J][M];
360
ma = EBAints.params->colidx[M][A];
362
t_jmbc = E_kima = 0.0;
364
if(T2BB.params->rowtot[Gbc] && T2BB.params->coltot[Gbc])
365
t_jmbc = T2BB.matrix[Gbc][jm][bc];
367
if(EBAints.params->rowtot[Gik] && EBAints.params->coltot[Gik])
368
E_kima = EBAints.matrix[Gik][ki][ma];
370
value_c -= t_jmbc * E_kima;
373
/* +t_MjAc * E_IkMb */
375
for(m=0; m < aoccpi[Gm]; m++) {
376
M = aocc_off[Gm] + m;
378
mj = T2AB.params->rowidx[M][J];
379
mb = EABints.params->colidx[M][B];
381
t_mjac = E_ikmb = 0.0;
383
if(T2AB.params->rowtot[Gac] && T2AB.params->coltot[Gac])
384
t_mjac = T2AB.matrix[Gac][mj][ac];
386
if(EABints.params->rowtot[Gik] && EABints.params->coltot[Gik])
387
E_ikmb = EABints.matrix[Gik][ik][mb];
389
value_c += t_mjac * E_ikmb;
392
/* -t_MjAb * E_IkMc */
394
for(m=0; m < aoccpi[Gm]; m++) {
395
M = aocc_off[Gm] + m;
397
mj = T2AB.params->rowidx[M][J];
398
mc = EABints.params->colidx[M][C];
400
t_mjab = E_ikmc = 0.0;
402
if(T2AB.params->rowtot[Gba] && T2AB.params->coltot[Gba])
403
t_mjab = T2AB.matrix[Gba][mj][ab];
405
if(EABints.params->rowtot[Gik] && EABints.params->coltot[Gik])
406
E_ikmc = EABints.matrix[Gik][ik][mc];
408
value_c -= t_mjab * E_ikmc;
411
/* +t_kmbc * E_jImA */
413
for(m=0; m < boccpi[Gm]; m++) {
414
M = bocc_off[Gm] + m;
416
km = T2BB.params->rowidx[K][M];
417
ma = EBAints.params->colidx[M][A];
419
t_kmbc = E_jima = 0.0;
421
if(T2BB.params->rowtot[Gbc] && T2BB.params->coltot[Gbc])
422
t_kmbc = T2BB.matrix[Gbc][km][bc];
424
if(EBAints.params->rowtot[Gji] && EBAints.params->coltot[Gji])
425
E_jima = EBAints.matrix[Gji][ji][ma];
427
value_c += t_kmbc * E_jima;
430
/* -t_MkAc * E_IjMb */
432
for(m=0; m < aoccpi[Gm]; m++) {
433
M = aocc_off[Gm] + m;
435
mk = T2AB.params->rowidx[M][K];
436
mb = EABints.params->colidx[M][B];
438
t_mkac = E_ijmb = 0.0;
440
if(T2AB.params->rowtot[Gac] && T2AB.params->coltot[Gac])
441
t_mkac = T2AB.matrix[Gac][mk][ac];
443
if(EABints.params->rowtot[Gji] && EABints.params->coltot[Gji])
444
E_ijmb = EABints.matrix[Gji][ij][mb];
446
value_c -= t_mkac * E_ijmb;
449
/* +t_MkAb * E_IjMc */
451
for(m=0; m < aoccpi[Gm]; m++) {
452
M = aocc_off[Gm] + m;
454
mk = T2AB.params->rowidx[M][K];
455
mc = EABints.params->colidx[M][C];
457
t_mkab = E_ijmc = 0.0;
459
if(T2AB.params->rowtot[Gba] && T2AB.params->coltot[Gba])
460
t_mkab = T2AB.matrix[Gba][mk][ab];
462
if(EABints.params->rowtot[Gji] && EABints.params->coltot[Gji])
463
E_ijmc = EABints.matrix[Gji][ij][mc];
465
value_c += t_mkab * E_ijmc;
468
/** disconnected triples **/
473
if(Gi == Ga && Gjk == Gbc) {
476
if(T1A.params->rowtot[Gi] && T1A.params->coltot[Gi])
477
t_ia = T1A.matrix[Gi][i][a];
479
if(DBBints.params->rowtot[Gjk] && DBBints.params->coltot[Gjk])
480
D_jkbc = DBBints.matrix[Gjk][jk][bc];
482
value_d += t_ia * D_jkbc;
486
if(Gj == Gb && Gik == Gac) {
489
if(T1B.params->rowtot[Gj] && T1B.params->coltot[Gj])
490
t_jb = T1B.matrix[Gj][j][b];
492
if(DABints.params->rowtot[Gik] && DABints.params->coltot[Gik])
493
D_ikac = DABints.matrix[Gik][ik][ac];
495
value_d += t_jb * D_ikac;
499
if(Gj == Gc && Gik == Gba) {
502
if(T1B.params->rowtot[Gj] && T1B.params->coltot[Gj])
503
t_jc = T1B.matrix[Gj][j][c];
505
if(DABints.params->rowtot[Gik] && DABints.params->coltot[Gik])
506
D_ikab = DABints.matrix[Gik][ik][ab];
508
value_d -= t_jc * D_ikab;
512
if(Gk == Gb && Gji == Gac) {
515
if(T1B.params->rowtot[Gk] && T1B.params->coltot[Gk])
516
t_kb = T1B.matrix[Gk][k][b];
518
if(DABints.params->rowtot[Gji] && DABints.params->coltot[Gji])
519
D_ijac = DABints.matrix[Gji][ij][ac];
521
value_d -= t_kb * D_ijac;
525
if(Gk == Gc && Gji == Gba) {
528
if(T1B.params->rowtot[Gk] && T1B.params->coltot[Gk])
529
t_kc = T1B.matrix[Gk][k][c];
531
if(DABints.params->rowtot[Gji] && DABints.params->coltot[Gji])
532
D_ijab = DABints.matrix[Gji][ij][ab];
534
value_d += t_kc * D_ijab;
538
if(fabs(value_c) > 1e-7) {
539
fprintf(outfile, "%d %d %d %d %d %d %20.15f\n", I, J, K, A, B, C, value_c);
543
/* Compute the Fock denominator */
545
if(fIJ.params->rowtot[Gi])
546
denom += fIJ.matrix[Gi][i][i];
547
if(fij.params->rowtot[Gj])
548
denom += fij.matrix[Gj][j][j];
549
if(fij.params->rowtot[Gk])
550
denom += fij.matrix[Gk][k][k];
551
if(fAB.params->rowtot[Ga])
552
denom -= fAB.matrix[Ga][a][a];
553
if(fab.params->rowtot[Gb])
554
denom -= fab.matrix[Gb][b][b];
555
if(fab.params->rowtot[Gc])
556
denom -= fab.matrix[Gc][c][c];
558
/* ET_ABB += (value_d + value_c) * value_c / denom; */
559
ET_ABB += (value_c + value_d) * value_c / denom;
578
/* fprintf(outfile, "cnt = %d\n", cnt); */
581
for(h=0; h < nirreps; h++) {
582
dpd_buf4_mat_irrep_close(&T2BB, h);
583
dpd_buf4_mat_irrep_close(&T2AB, h);
584
dpd_buf4_mat_irrep_close(&FBBints, h);
585
dpd_buf4_mat_irrep_close(&FABints, h);
586
dpd_buf4_mat_irrep_close(&FBAints, h);
587
dpd_buf4_mat_irrep_close(&EBBints, h);
588
dpd_buf4_mat_irrep_close(&EABints, h);
589
dpd_buf4_mat_irrep_close(&EBAints, h);
590
dpd_buf4_mat_irrep_close(&DBBints, h);
591
dpd_buf4_mat_irrep_close(&DABints, h);
594
dpd_buf4_close(&T2BB);
595
dpd_buf4_close(&T2AB);
596
dpd_buf4_close(&FBBints);
597
dpd_buf4_close(&FABints);
598
dpd_buf4_close(&FBAints);
599
dpd_buf4_close(&EBBints);
600
dpd_buf4_close(&EABints);
601
dpd_buf4_close(&EBAints);
602
dpd_buf4_close(&DBBints);
603
dpd_buf4_close(&DABints);
605
dpd_file2_mat_close(&T1A);
606
dpd_file2_close(&T1A);
607
dpd_file2_mat_close(&T1B);
608
dpd_file2_close(&T1B);
610
dpd_file2_mat_close(&fIJ);
611
dpd_file2_mat_close(&fij);
612
dpd_file2_mat_close(&fAB);
613
dpd_file2_mat_close(&fab);
614
dpd_file2_close(&fIJ);
615
dpd_file2_close(&fij);
616
dpd_file2_close(&fAB);
617
dpd_file2_close(&fab);