← Back to branch summary

~ubuntu-branches/ubuntu/karmic/psicode/karmic

« back to all changes in this revision

Viewing changes to src/bin/cctriples/ET_AAA.cc

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck, Michael Banck, Daniel Leidert
  • Date: 2009-02-23 00:12:02 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20090223001202-rutldoy3dimfpesc
Tags: 3.4.0-1
* New upstream release.

[ Michael Banck ]
* debian/patches/01_DESTDIR.dpatch: Refreshed.
* debian/patches/02_FHS.dpatch: Removed, applied upstream.
* debian/patches/03_debian_docdir: Likewise.
* debian/patches/04_man.dpatch: Likewise.
* debian/patches/06_466828_fix_gcc_43_ftbfs.dpatch: Likewise.
* debian/patches/07_464867_move_executables: Fixed and refreshed.
* debian/patches/00list: Adjusted.
* debian/control: Improved description.
* debian/patches-held: Removed.
* debian/rules (install/psi3): Do not ship the ruby bindings for now.

[ Daniel Leidert ]
* debian/rules: Fix txtdir via DEB_MAKE_INSTALL_TARGET.
* debian/patches/01_DESTDIR.dpatch: Refreshed.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/bin/cctriples/ET_AAA.cc
 
1
/*! \file
 
2
    \ingroup CCTRIPLES
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
 
 
6
/*! \defgroup CCTRIPLES cctriples: Evaluate triple excitations */
 
7
 
 
8
#include <cstdio>
 
9
#include <cmath>
 
10
#include <libdpd/dpd.h>
 
11
#include "MOInfo.h"
 
12
#include "Params.h"
 
13
#define EXTERN
 
14
#include "globals.h"
 
15
 
 
16
namespace psi { namespace cctriples {
 
17
 
 
18
double ET_AAA(void)
 
19
{
 
20
  int cnt;
 
21
  int h, nirreps;
 
22
  int Gi, Gj, Gk, Ga, Gb, Gc, Ge, Gm;
 
23
  int Gji, Gij, Gjk, Gik, Gbc, Gac, Gba;
 
24
  int I, J, K, A, B, C, E, M;
 
25
  int i, j, k, a, b, c, e, m;
 
26
  int ij, ji, ik, jk, bc, ac, ba;
 
27
  int im, jm, km, ma, mb, mc;
 
28
  int ae, be, ce, ke, ie, je;
 
29
  int *occpi, *virtpi, *occ_off, *vir_off;
 
30
  double value_c, value_d, denom, ET_AAA;
 
31
  double t_ijae, t_ijbe, t_ijce, t_jkae, t_jkbe, t_jkce, t_ikae, t_ikbe, t_ikce;
 
32
  double F_kebc, F_keac, F_keba, F_iebc, F_ieac, F_ieba, F_jebc, F_jeac, F_jeba;
 
33
  double t_imbc, t_imac, t_imba, t_jmbc, t_jmac, t_jmba, t_kmbc, t_kmac, t_kmba;
 
34
  double E_jkma, E_jkmb, E_jkmc, E_ikma, E_ikmb, E_ikmc, E_jima, E_jimb, E_jimc;
 
35
  double t_ia, t_ib, t_ic, t_ja, t_jb, t_jc, t_ka, t_kb, t_kc;
 
36
  double D_jkbc, D_jkac, D_jkba, D_ikbc, D_ikac, D_ikba, D_jibc, D_jiac, D_jiba;
 
37
  dpdbuf4 T2, Fints, Eints, Dints;
 
38
  dpdfile2 fIJ, fAB, T1;
 
39
 
 
40
  nirreps = moinfo.nirreps;
 
41
  occpi = moinfo.occpi; virtpi = moinfo.virtpi;
 
42
  occ_off = moinfo.occ_off;
 
43
  vir_off = moinfo.vir_off;
 
44
 
 
45
  dpd_file2_init(&fIJ, CC_OEI, 0, 0, 0, "fIJ");
 
46
  dpd_file2_init(&fAB, CC_OEI, 0, 1, 1, "fAB");
 
47
  dpd_file2_mat_init(&fIJ);
 
48
  dpd_file2_mat_init(&fAB);
 
49
  dpd_file2_mat_rd(&fIJ);
 
50
  dpd_file2_mat_rd(&fAB);
 
51
 
 
52
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, "tIA");
 
53
  dpd_file2_mat_init(&T1);
 
54
  dpd_file2_mat_rd(&T1);
 
55
 
 
56
  dpd_buf4_init(&T2, CC_TAMPS, 0, 0, 5, 2, 7, 0, "tIJAB");
 
57
  dpd_buf4_init(&Fints, CC_FINTS, 0, 10, 5, 10, 5, 1, "F <ia|bc>");
 
58
  dpd_buf4_init(&Eints, CC_EINTS, 0, 0, 10, 2, 10, 0, "E <ij||ka> (i>j,ka)");
 
59
  dpd_buf4_init(&Dints, CC_DINTS, 0, 0, 5, 0, 5, 0, "D <ij||ab>");
 
60
  for(h=0; h < nirreps; h++) {
 
61
    dpd_buf4_mat_irrep_init(&T2, h);
 
62
    dpd_buf4_mat_irrep_rd(&T2, h);
 
63
 
 
64
    dpd_buf4_mat_irrep_init(&Fints, h);
 
65
    dpd_buf4_mat_irrep_rd(&Fints, h);
 
66
 
 
67
    dpd_buf4_mat_irrep_init(&Eints, h);
 
68
    dpd_buf4_mat_irrep_rd(&Eints, h);
 
69
 
 
70
    dpd_buf4_mat_irrep_init(&Dints, h);
 
71
    dpd_buf4_mat_irrep_rd(&Dints, h);
 
72
  }
 
73
 
 
74
  cnt = 0;
 
75
  ET_AAA = 0.0;
 
76
 
 
77
  for(Gi=0; Gi < nirreps; Gi++) {
 
78
    for(Gj=0; Gj < nirreps; Gj++) {
 
79
      Gij = Gji = Gi ^ Gj;
 
80
      for(Gk=0; Gk < nirreps; Gk++) {
 
81
        Gjk = Gj ^ Gk;
 
82
        Gik = Gi ^ Gk;
 
83
 
 
84
        for(Ga=0; Ga < nirreps; Ga++) {
 
85
          for(Gb=0; Gb < nirreps; Gb++) {
 
86
            Gc = Gi ^ Gj ^ Gk ^ Ga ^ Gb;
 
87
 
 
88
            Gbc = Gb^Gc;
 
89
            Gac = Ga^Gc;
 
90
            Gba = Gb^Ga;
 
91
 
 
92
            for(i=0; i < occpi[Gi]; i++) {
 
93
              I = occ_off[Gi] + i;
 
94
              for(j=0; j < occpi[Gj]; j++) {
 
95
                J = occ_off[Gj] + j;
 
96
                for(k=0; k < occpi[Gk]; k++) {
 
97
                  K = occ_off[Gk] + k;
 
98
 
 
99
                  ij = T2.params->rowidx[I][J];
 
100
                  ji = T2.params->rowidx[J][I];
 
101
                  jk = T2.params->rowidx[J][K];
 
102
                  ik = T2.params->rowidx[I][K];
 
103
 
 
104
                  for(a=0; a < virtpi[Ga]; a++) {
 
105
                    A = vir_off[Ga] + a;
 
106
                    for(b=0; b < virtpi[Gb]; b++) {
 
107
                      B = vir_off[Gb] + b;
 
108
                      for(c=0; c < virtpi[Gc]; c++) {
 
109
                        C = vir_off[Gc] + c;
 
110
 
 
111
                        bc = Fints.params->colidx[B][C];
 
112
                        ac = Fints.params->colidx[A][C];
 
113
                        ba = Fints.params->colidx[B][A];
 
114
 
 
115
                        value_c = 0.0;
 
116
 
 
117
                        /** <ov||vv> --> connected triples **/
 
118
 
 
119
                        /* -t_jkae * F_iebc */
 
120
                        Ge = Gj ^ Gk ^ Ga;
 
121
                        for(e=0; e < virtpi[Ge]; e++) {
 
122
                          E = vir_off[Ge] + e;
 
123
 
 
124
                          ae = T2.params->colidx[A][E];
 
125
                          ie = Fints.params->rowidx[I][E];
 
126
 
 
127
                          t_jkae = F_iebc = 0.0;
 
128
 
 
129
                          if(T2.params->rowtot[Gjk] && T2.params->coltot[Gjk])
 
130
                            t_jkae = T2.matrix[Gjk][jk][ae];
 
131
 
 
132
                          if(Fints.params->rowtot[Gbc] && Fints.params->coltot[Gbc])
 
133
                            F_iebc = Fints.matrix[Gbc][ie][bc];
 
134
 
 
135
                          value_c -= t_jkae * F_iebc;
 
136
                        }
 
137
 
 
138
                        /* +t_jkbe * F_ieac */
 
139
                        Ge = Gj ^ Gk ^ Gb;
 
140
                        for(e=0; e < virtpi[Ge]; e++) {
 
141
                          E = vir_off[Ge] + e;
 
142
 
 
143
                          be = T2.params->colidx[B][E];
 
144
                          ie = Fints.params->rowidx[I][E];
 
145
 
 
146
                          t_jkbe = F_ieac = 0.0;
 
147
 
 
148
                          if(T2.params->rowtot[Gjk] && T2.params->coltot[Gjk])
 
149
                            t_jkbe = T2.matrix[Gjk][jk][be];
 
150
 
 
151
                          if(Fints.params->rowtot[Gac] && Fints.params->coltot[Gac])
 
152
                            F_ieac = Fints.matrix[Gac][ie][ac];
 
153
 
 
154
                          value_c += t_jkbe * F_ieac;
 
155
                        }
 
156
 
 
157
                        /* +t_jkce * F_ieba */
 
158
                        Ge = Gj ^ Gk ^ Gc;
 
159
                        for(e=0; e < virtpi[Ge]; e++) {
 
160
                          E = vir_off[Ge] + e;
 
161
 
 
162
                          ce = T2.params->colidx[C][E];
 
163
                          ie = Fints.params->rowidx[I][E];
 
164
 
 
165
                          t_jkce = F_ieba = 0.0;
 
166
 
 
167
                          if(T2.params->rowtot[Gjk] && T2.params->coltot[Gjk])
 
168
                            t_jkce = T2.matrix[Gjk][jk][ce];
 
169
 
 
170
                          if(Fints.params->rowtot[Gba] && Fints.params->coltot[Gba])
 
171
                            F_ieba = Fints.matrix[Gba][ie][ba];
 
172
 
 
173
                          value_c += t_jkce * F_ieba;
 
174
                        }
 
175
 
 
176
                        /* +t_ikae * F_jebc */
 
177
                        Ge = Gi ^ Gk ^ Ga;
 
178
                        for(e=0; e < virtpi[Ge]; e++) {
 
179
                          E = vir_off[Ge] + e;
 
180
 
 
181
                          ae = T2.params->colidx[A][E];
 
182
                          je = Fints.params->rowidx[J][E];
 
183
 
 
184
                          t_ikae = F_jebc = 0.0;
 
185
 
 
186
                          if(T2.params->rowtot[Gik] && T2.params->coltot[Gik])
 
187
                            t_ikae = T2.matrix[Gik][ik][ae];
 
188
 
 
189
                          if(Fints.params->rowtot[Gbc] && Fints.params->coltot[Gbc]) 
 
190
                            F_jebc = Fints.matrix[Gbc][je][bc];
 
191
 
 
192
                          value_c += t_ikae * F_jebc;
 
193
                        }
 
194
 
 
195
                        /* -t_ikbe * F_jeac */
 
196
                        Ge = Gi ^ Gk ^ Gb;
 
197
                        for(e=0; e < virtpi[Ge]; e++) {
 
198
                          E = vir_off[Ge] + e;
 
199
 
 
200
                          be = T2.params->colidx[B][E];
 
201
                          je = Fints.params->rowidx[J][E];
 
202
 
 
203
                          t_ikbe = F_jeac = 0.0;
 
204
 
 
205
                          if(T2.params->rowtot[Gik] && T2.params->coltot[Gik])
 
206
                            t_ikbe = T2.matrix[Gik][ik][be];
 
207
 
 
208
                          if(Fints.params->rowtot[Gac] && Fints.params->coltot[Gac]) 
 
209
                            F_jeac = Fints.matrix[Gac][je][ac];
 
210
 
 
211
                          value_c -= t_ikbe * F_jeac;
 
212
                        }
 
213
 
 
214
                        /* -t_ikce * F_jeba */
 
215
                        Ge = Gi ^ Gk ^ Gc;
 
216
                        for(e=0; e < virtpi[Ge]; e++) {
 
217
                          E = vir_off[Ge] + e;
 
218
 
 
219
                          ce = T2.params->colidx[C][E];
 
220
                          je = Fints.params->rowidx[J][E];
 
221
 
 
222
                          t_ikce = F_jeba = 0.0;
 
223
 
 
224
                          if(T2.params->rowtot[Gik] && T2.params->coltot[Gik])
 
225
                            t_ikce = T2.matrix[Gik][ik][ce];
 
226
 
 
227
                          if(Fints.params->rowtot[Gba] && Fints.params->coltot[Gba])
 
228
                            F_jeba = Fints.matrix[Gba][je][ba];
 
229
 
 
230
                          value_c -= t_ikce * F_jeba;
 
231
                        }
 
232
 
 
233
                        /* -t_ijae * F_kebc */
 
234
                        Ge = Gi ^ Gj ^ Ga;
 
235
                        for(e=0; e < virtpi[Ge]; e++) {
 
236
                          E = vir_off[Ge] + e;
 
237
 
 
238
                          ae = T2.params->colidx[A][E];
 
239
                          ke = Fints.params->rowidx[K][E];
 
240
 
 
241
                          t_ijae = F_kebc = 0.0;
 
242
 
 
243
                          if(T2.params->rowtot[Gij] && T2.params->coltot[Gij])
 
244
                            t_ijae = T2.matrix[Gij][ij][ae];
 
245
 
 
246
                          if(Fints.params->rowtot[Gbc] && Fints.params->coltot[Gbc])
 
247
                            F_kebc = Fints.matrix[Gbc][ke][bc];
 
248
 
 
249
                          value_c -= t_ijae * F_kebc;
 
250
                        }
 
251
 
 
252
                        /* +t_ijbe * F_keac */
 
253
                        Ge = Gi ^ Gj ^ Gb;
 
254
                        for(e=0; e < virtpi[Ge]; e++) {
 
255
                          E = vir_off[Ge] + e;
 
256
 
 
257
                          be = T2.params->colidx[B][E];
 
258
                          ke = Fints.params->rowidx[K][E];
 
259
 
 
260
                          t_ijbe = F_keac = 0.0;
 
261
 
 
262
                          if(T2.params->rowtot[Gij] && T2.params->coltot[Gij])
 
263
                            t_ijbe = T2.matrix[Gij][ij][be];
 
264
 
 
265
                          if(Fints.params->rowtot[Gac] && Fints.params->coltot[Gac]) 
 
266
                            F_keac = Fints.matrix[Gac][ke][ac];
 
267
 
 
268
                          value_c += t_ijbe * F_keac;
 
269
                        }
 
270
 
 
271
                        /* +t_ijce * F_keba */
 
272
                        Ge = Gi ^ Gj ^ Gc;
 
273
                        for(e=0; e < virtpi[Ge]; e++) {
 
274
                          E = vir_off[Ge] + e;
 
275
 
 
276
                          ce = T2.params->colidx[C][E];
 
277
                          ke = Fints.params->rowidx[K][E];
 
278
 
 
279
                          t_ijce = F_keba = 0.0;
 
280
 
 
281
                          if(T2.params->rowtot[Gij] && T2.params->coltot[Gij])
 
282
                            t_ijce = T2.matrix[Gij][ij][ce];
 
283
 
 
284
                          if(Fints.params->rowtot[Gba] && Fints.params->coltot[Gba]) 
 
285
                            F_keba = Fints.matrix[Gba][ke][ba];
 
286
 
 
287
                          value_c += t_ijce * F_keba;
 
288
                        }
 
289
 
 
290
                        /** <oo||ov> --> connected triples **/
 
291
 
 
292
                        /* -t_imbc * E_jkma */
 
293
                        Gm = Gi ^ Gb ^ Gc;
 
294
                        for(m=0; m < occpi[Gm]; m++) {
 
295
                          M = occ_off[Gm] + m;
 
296
 
 
297
                          im = T2.params->rowidx[I][M];
 
298
                          ma = Eints.params->colidx[M][A];
 
299
 
 
300
                          t_imbc = E_jkma = 0.0;
 
301
 
 
302
                          if(T2.params->rowtot[Gbc] && T2.params->coltot[Gbc])
 
303
                            t_imbc = T2.matrix[Gbc][im][bc];
 
304
 
 
305
                          if(Eints.params->rowtot[Gjk] && Eints.params->coltot[Gjk])
 
306
                            E_jkma = Eints.matrix[Gjk][jk][ma];
 
307
 
 
308
                          value_c -= t_imbc * E_jkma;
 
309
                        }
 
310
 
 
311
                        /* +t_imac * E_jkmb */
 
312
                        Gm = Gi ^ Ga ^ Gc;
 
313
                        for(m=0; m < occpi[Gm]; m++) {
 
314
                          M = occ_off[Gm] + m;
 
315
 
 
316
                          im = T2.params->rowidx[I][M];
 
317
                          mb = Eints.params->colidx[M][B];
 
318
 
 
319
                          t_imac = E_jkmb = 0.0;
 
320
 
 
321
                          if(T2.params->rowtot[Gac] && T2.params->coltot[Gac])
 
322
                            t_imac = T2.matrix[Gac][im][ac];
 
323
 
 
324
                          if(Eints.params->rowtot[Gjk] && Eints.params->coltot[Gjk])
 
325
                            E_jkmb = Eints.matrix[Gjk][jk][mb];
 
326
 
 
327
                          value_c += t_imac * E_jkmb;
 
328
                        }
 
329
 
 
330
                        /* +t_imba * E_jkmc */
 
331
                        Gm = Gi ^ Gb ^ Ga;
 
332
                        for(m=0; m < occpi[Gm]; m++) {
 
333
                          M = occ_off[Gm] + m;
 
334
 
 
335
                          im = T2.params->rowidx[I][M];
 
336
                          mc = Eints.params->colidx[M][C];
 
337
 
 
338
                          t_imba = E_jkmc = 0.0;
 
339
 
 
340
                          if(T2.params->rowtot[Gba] && T2.params->coltot[Gba])
 
341
                            t_imba = T2.matrix[Gba][im][ba];
 
342
 
 
343
                          if(Eints.params->rowtot[Gjk] && Eints.params->coltot[Gjk])
 
344
                            E_jkmc = Eints.matrix[Gjk][jk][mc];
 
345
 
 
346
                          value_c += t_imba * E_jkmc;
 
347
                        }
 
348
 
 
349
                        /* +t_jmbc * E_ikma */
 
350
                        Gm = Gj ^ Gb ^ Gc;
 
351
                        for(m=0; m < occpi[Gm]; m++) {
 
352
                          M = occ_off[Gm] + m;
 
353
 
 
354
                          jm = T2.params->rowidx[J][M];
 
355
                          ma = Eints.params->colidx[M][A];
 
356
 
 
357
                          t_jmbc = E_ikma = 0.0;
 
358
 
 
359
                          if(T2.params->rowtot[Gbc] && T2.params->coltot[Gbc])
 
360
                            t_jmbc = T2.matrix[Gbc][jm][bc];
 
361
 
 
362
                          if(Eints.params->rowtot[Gik] && Eints.params->coltot[Gik])
 
363
                            E_ikma = Eints.matrix[Gik][ik][ma];
 
364
 
 
365
                          value_c += t_jmbc * E_ikma;
 
366
                        }
 
367
 
 
368
                        /* -t_jmac * E_ikmb */
 
369
                        Gm = Gj ^ Ga ^ Gc;
 
370
                        for(m=0; m < occpi[Gm]; m++) {
 
371
                          M = occ_off[Gm] + m;
 
372
 
 
373
                          jm = T2.params->rowidx[J][M];
 
374
                          mb = Eints.params->colidx[M][B];
 
375
 
 
376
                          t_jmac = E_ikmb = 0.0;
 
377
 
 
378
                          if(T2.params->rowtot[Gac] && T2.params->coltot[Gac])
 
379
                            t_jmac = T2.matrix[Gac][jm][ac];
 
380
 
 
381
                          if(Eints.params->rowtot[Gik] && Eints.params->coltot[Gik])
 
382
                            E_ikmb = Eints.matrix[Gik][ik][mb];
 
383
 
 
384
                          value_c -= t_jmac * E_ikmb;
 
385
                        }
 
386
 
 
387
                        /* -t_jmba * E_ikmc */
 
388
                        Gm = Gj ^ Gb ^ Ga;
 
389
                        for(m=0; m < occpi[Gm]; m++) {
 
390
                          M = occ_off[Gm] + m;
 
391
 
 
392
                          jm = T2.params->rowidx[J][M];
 
393
                          mc = Eints.params->colidx[M][C];
 
394
 
 
395
                          t_jmba = E_ikmc = 0.0;
 
396
 
 
397
                          if(T2.params->rowtot[Gba] && T2.params->coltot[Gba])
 
398
                            t_jmba = T2.matrix[Gba][jm][ba];
 
399
 
 
400
                          if(Eints.params->rowtot[Gik] && Eints.params->coltot[Gik])
 
401
                            E_ikmc = Eints.matrix[Gik][ik][mc];
 
402
 
 
403
                          value_c -= t_jmba * E_ikmc;
 
404
                        }
 
405
 
 
406
                        /* +t_kmbc * E_jima */
 
407
                        Gm = Gk ^ Gb ^ Gc;
 
408
                        for(m=0; m < occpi[Gm]; m++) {
 
409
                          M = occ_off[Gm] + m;
 
410
 
 
411
                          km = T2.params->rowidx[K][M];
 
412
                          ma = Eints.params->colidx[M][A];
 
413
 
 
414
                          t_kmbc = E_jima = 0.0;
 
415
 
 
416
                          if(T2.params->rowtot[Gbc] && T2.params->coltot[Gbc])
 
417
                            t_kmbc = T2.matrix[Gbc][km][bc];
 
418
 
 
419
                          if(Eints.params->rowtot[Gji] && Eints.params->coltot[Gji])
 
420
                            E_jima = Eints.matrix[Gji][ji][ma];
 
421
 
 
422
                          value_c += t_kmbc * E_jima;
 
423
                        }
 
424
 
 
425
                        /* -t_kmac * E_jimb */
 
426
                        Gm = Gk ^ Ga ^ Gc;
 
427
                        for(m=0; m < occpi[Gm]; m++) {
 
428
                          M = occ_off[Gm] + m;
 
429
 
 
430
                          km = T2.params->rowidx[K][M];
 
431
                          mb = Eints.params->colidx[M][B];
 
432
 
 
433
                          t_kmac = E_jimb = 0.0;
 
434
 
 
435
                          if(T2.params->rowtot[Gac] && T2.params->coltot[Gac])
 
436
                            t_kmac = T2.matrix[Gac][km][ac];
 
437
 
 
438
                          if(Eints.params->rowtot[Gji] && Eints.params->coltot[Gji])
 
439
                            E_jimb = Eints.matrix[Gji][ji][mb];
 
440
 
 
441
                          value_c -= t_kmac * E_jimb;
 
442
                        }
 
443
 
 
444
                        /* -t_kmba * E_jimc */
 
445
                        Gm = Gk ^ Gb ^ Ga;
 
446
                        for(m=0; m < occpi[Gm]; m++) {
 
447
                          M = occ_off[Gm] + m;
 
448
 
 
449
                          km = T2.params->rowidx[K][M];
 
450
                          mc = Eints.params->colidx[M][C];
 
451
 
 
452
                          t_kmba = E_jimc = 0.0;
 
453
 
 
454
                          if(T2.params->rowtot[Gba] && T2.params->coltot[Gba])
 
455
                            t_kmba = T2.matrix[Gba][km][ba];
 
456
 
 
457
                          if(Eints.params->rowtot[Gji] && Eints.params->coltot[Gji])
 
458
                            E_jimc = Eints.matrix[Gji][ji][mc];
 
459
 
 
460
                          value_c -= t_kmba * E_jimc;
 
461
                        }
 
462
 
 
463
                        /** disconnected triples **/
 
464
 
 
465
                        value_d = 0.0;
 
466
 
 
467
                        /* +t_ia * D_jkbc */
 
468
                        if(Gi == Ga && Gjk == Gbc) {
 
469
                          t_ia = D_jkbc = 0.0;
 
470
 
 
471
                          if(T1.params->rowtot[Gi] && T1.params->coltot[Gi])
 
472
                            t_ia = T1.matrix[Gi][i][a];
 
473
 
 
474
                          if(Dints.params->rowtot[Gjk] && Dints.params->coltot[Gjk])
 
475
                            D_jkbc = Dints.matrix[Gjk][jk][bc];
 
476
 
 
477
                          value_d += t_ia * D_jkbc;
 
478
                        }
 
479
 
 
480
                        /* -t_ib * D_jkac */
 
481
                        if(Gi == Gb && Gjk == Gac) {
 
482
                          t_ib = D_jkac = 0.0;
 
483
 
 
484
                          if(T1.params->rowtot[Gi] && T1.params->coltot[Gi])
 
485
                            t_ib = T1.matrix[Gi][i][b];
 
486
 
 
487
                          if(Dints.params->rowtot[Gjk] && Dints.params->coltot[Gjk])
 
488
                            D_jkac = Dints.matrix[Gjk][jk][ac];
 
489
 
 
490
                          value_d -= t_ib * D_jkac;
 
491
                        }
 
492
 
 
493
                        /* -t_ic * D_jkba */
 
494
                        if(Gi == Gc && Gjk == Gba) {
 
495
                          t_ic = D_jkba = 0.0;
 
496
 
 
497
                          if(T1.params->rowtot[Gi] && T1.params->coltot[Gi])
 
498
                            t_ic = T1.matrix[Gi][i][c];
 
499
 
 
500
                          if(Dints.params->rowtot[Gjk] && Dints.params->coltot[Gjk])
 
501
                            D_jkba = Dints.matrix[Gjk][jk][ba];
 
502
 
 
503
                          value_d -= t_ic * D_jkba;
 
504
                        }
 
505
 
 
506
                        /* -t_ja * D_ikbc */
 
507
                        if(Gj == Ga && Gik == Gbc) {
 
508
                          t_ja = D_ikbc = 0.0;
 
509
 
 
510
                          if(T1.params->rowtot[Gj] && T1.params->coltot[Gj])
 
511
                            t_ja = T1.matrix[Gj][j][a];
 
512
 
 
513
                          if(Dints.params->rowtot[Gik] && Dints.params->coltot[Gik])
 
514
                            D_ikbc = Dints.matrix[Gik][ik][bc];
 
515
 
 
516
                          value_d -= t_ja * D_ikbc;
 
517
                        }
 
518
 
 
519
                        /* +t_jb * D_ikac */
 
520
                        if(Gj == Gb && Gik == Gac) {
 
521
                          t_jb = D_ikac = 0.0;
 
522
 
 
523
                          if(T1.params->rowtot[Gj] && T1.params->coltot[Gj])
 
524
                            t_jb = T1.matrix[Gj][j][b];
 
525
 
 
526
                          if(Dints.params->rowtot[Gik] && Dints.params->coltot[Gik])
 
527
                            D_ikac = Dints.matrix[Gik][ik][ac];
 
528
 
 
529
                          value_d += t_jb * D_ikac;
 
530
                        }
 
531
 
 
532
                        /* +t_jc * D_ikba */
 
533
                        if(Gj == Gc && Gik == Gba) {
 
534
                          t_jc = D_ikba = 0.0;
 
535
 
 
536
                          if(T1.params->rowtot[Gj] && T1.params->coltot[Gj])
 
537
                            t_jc = T1.matrix[Gj][j][c];
 
538
 
 
539
                          if(Dints.params->rowtot[Gik] && Dints.params->coltot[Gik])
 
540
                            D_ikba = Dints.matrix[Gik][ik][ba];
 
541
 
 
542
                          value_d += t_jc * D_ikba;
 
543
                        }
 
544
 
 
545
                        /* -t_ka * D_jibc */
 
546
                        if(Gk == Ga && Gji == Gbc) {
 
547
                          t_ka = D_jibc = 0.0;
 
548
 
 
549
                          if(T1.params->rowtot[Gk] && T1.params->coltot[Gk])
 
550
                            t_ka = T1.matrix[Gk][k][a];
 
551
 
 
552
                          if(Dints.params->rowtot[Gji] && Dints.params->coltot[Gji])
 
553
                            D_jibc = Dints.matrix[Gji][ji][bc];
 
554
 
 
555
                          value_d -= t_ka * D_jibc;
 
556
                        }
 
557
 
 
558
                        /* +t_kb * D_jiac */
 
559
                        if(Gk == Gb && Gji == Gac) {
 
560
                          t_kb = D_jiac = 0.0;
 
561
 
 
562
                          if(T1.params->rowtot[Gk] && T1.params->coltot[Gk])
 
563
                            t_kb = T1.matrix[Gk][k][b];
 
564
 
 
565
                          if(Dints.params->rowtot[Gji] && Dints.params->coltot[Gji])
 
566
                            D_jiac = Dints.matrix[Gji][ji][ac];
 
567
 
 
568
                          value_d += t_kb * D_jiac;
 
569
                        }
 
570
 
 
571
                        /* +t_kc * D_jiba */
 
572
                        if(Gk == Gc && Gji == Gba) {
 
573
                          t_kc = D_jiba = 0.0;
 
574
 
 
575
                          if(T1.params->rowtot[Gk] && T1.params->coltot[Gk])
 
576
                            t_kc = T1.matrix[Gk][k][c];
 
577
 
 
578
                          if(Dints.params->rowtot[Gji] && Dints.params->coltot[Gji])
 
579
                            D_jiba = Dints.matrix[Gji][ji][ba];
 
580
 
 
581
                          value_d += t_kc * D_jiba;
 
582
                        }
 
583
 
 
584
                        /*
 
585
                        if(fabs(value_c) > 1e-7) {
 
586
                          cnt++;
 
587
                          fprintf(outfile, "%d %d %d %d %d %d %20.14f\n", I, J, K, A, B, C, value_c);
 
588
                        }
 
589
                        */
 
590
 
 
591
                        /* Compute the Fock denominator */
 
592
                        denom = 0.0;
 
593
                        if(fIJ.params->rowtot[Gi])
 
594
                          denom += fIJ.matrix[Gi][i][i];
 
595
                        if(fIJ.params->rowtot[Gj])
 
596
                          denom += fIJ.matrix[Gj][j][j];
 
597
                        if(fIJ.params->rowtot[Gk])
 
598
                          denom += fIJ.matrix[Gk][k][k];
 
599
                        if(fAB.params->rowtot[Ga])
 
600
                          denom -= fAB.matrix[Ga][a][a];
 
601
                        if(fAB.params->rowtot[Gb])
 
602
                          denom -= fAB.matrix[Gb][b][b];
 
603
                        if(fAB.params->rowtot[Gc])
 
604
                          denom -= fAB.matrix[Gc][c][c];
 
605
 
 
606
                        ET_AAA += (value_d + value_c) * value_c / denom;
 
607
 
 
608
 
 
609
                      } /* c */
 
610
                    } /* b */
 
611
                  } /* a */
 
612
 
 
613
                } /* k */
 
614
              } /* j */
 
615
            } /* i */
 
616
 
 
617
          } /* Gb */
 
618
        } /* Ga */
 
619
 
 
620
      } /* Gk */
 
621
    } /* Gj */
 
622
  } /* Gi */
 
623
 
 
624
  /*  fprintf(outfile, "cnt = %d\n", cnt); */
 
625
  ET_AAA /= 36.0;
 
626
  /*  fprintf(outfile, "ET_AAA = %20.14f\n", ET_AAA); */
 
627
 
 
628
  for(h=0; h < nirreps; h++) {
 
629
    dpd_buf4_mat_irrep_close(&T2, h);
 
630
    dpd_buf4_mat_irrep_close(&Fints, h);
 
631
    dpd_buf4_mat_irrep_close(&Eints, h);
 
632
    dpd_buf4_mat_irrep_close(&Dints, h);
 
633
  }
 
634
 
 
635
  dpd_buf4_close(&T2);
 
636
  dpd_buf4_close(&Fints);
 
637
  dpd_buf4_close(&Eints);
 
638
  dpd_buf4_close(&Dints);
 
639
 
 
640
  dpd_file2_mat_close(&T1);
 
641
  dpd_file2_close(&T1);
 
642
 
 
643
  dpd_file2_mat_close(&fIJ);
 
644
  dpd_file2_mat_close(&fAB);
 
645
  dpd_file2_close(&fIJ);
 
646
  dpd_file2_close(&fAB);
 
647
 
 
648
  return ET_AAA;
 
649
}
 
650
 
 
651
}} // namespace psi::cctriples