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

« back to all changes in this revision

Viewing changes to src/bin/cctriples/ET_UHF_ABB_noddy.c

  • 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:
1
 
#include <stdio.h>
2
 
#include <math.h>
3
 
#include <libdpd/dpd.h>
4
 
#define EXTERN
5
 
#include "globals.h"
6
 
 
7
 
double ET_UHF_ABB_noddy(void)
8
 
{
9
 
  int cnt;
10
 
  int h, nirreps;
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;
28
 
  dpdbuf4 T2AB, T2BB;
29
 
  dpdbuf4 FBBints, FABints, FBAints;
30
 
  dpdbuf4 EBBints, EABints, EBAints;
31
 
  dpdbuf4 DBBints, DABints;
32
 
  dpdfile2 T1A, T1B, fIJ, fij, fAB, fab;
33
 
 
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;
43
 
 
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);
56
 
 
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);
63
 
 
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");
66
 
 
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>");
70
 
 
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>");
74
 
 
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>");
77
 
 
78
 
  for(h=0; h < nirreps; h++) {
79
 
    dpd_buf4_mat_irrep_init(&T2BB, h);
80
 
    dpd_buf4_mat_irrep_rd(&T2BB, h);
81
 
 
82
 
    dpd_buf4_mat_irrep_init(&T2AB, h);
83
 
    dpd_buf4_mat_irrep_rd(&T2AB, h);
84
 
 
85
 
    dpd_buf4_mat_irrep_init(&FBBints, h);
86
 
    dpd_buf4_mat_irrep_rd(&FBBints, h);
87
 
 
88
 
    dpd_buf4_mat_irrep_init(&FABints, h);
89
 
    dpd_buf4_mat_irrep_rd(&FABints, h);
90
 
 
91
 
    dpd_buf4_mat_irrep_init(&FBAints, h);
92
 
    dpd_buf4_mat_irrep_rd(&FBAints, h);
93
 
 
94
 
    dpd_buf4_mat_irrep_init(&EBBints, h);
95
 
    dpd_buf4_mat_irrep_rd(&EBBints, h);
96
 
 
97
 
    dpd_buf4_mat_irrep_init(&EABints, h);
98
 
    dpd_buf4_mat_irrep_rd(&EABints, h);
99
 
 
100
 
    dpd_buf4_mat_irrep_init(&EBAints, h);
101
 
    dpd_buf4_mat_irrep_rd(&EBAints, h);
102
 
 
103
 
    dpd_buf4_mat_irrep_init(&DBBints, h);
104
 
    dpd_buf4_mat_irrep_rd(&DBBints, h);
105
 
 
106
 
    dpd_buf4_mat_irrep_init(&DABints, h);
107
 
    dpd_buf4_mat_irrep_rd(&DABints, h);
108
 
  }
109
 
 
110
 
  cnt = 0;
111
 
  ET_ABB = 0.0;
112
 
 
113
 
  for(Gi=0; Gi < nirreps; Gi++) {
114
 
    for(Gj=0; Gj < nirreps; Gj++) {
115
 
      for(Gk=0; Gk < nirreps; Gk++) {
116
 
        Gij = Gji = Gi ^ Gj;
117
 
        Gjk = Gj ^ Gk;
118
 
        Gik = Gi ^ Gk;
119
 
 
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;
126
 
 
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];
133
 
 
134
 
              if(J >= K) {
135
 
 
136
 
                for(Ga=0; Ga < nirreps; Ga++) {
137
 
                  for(Gb=0; Gb < nirreps; Gb++) {
138
 
                    Gc = Gi ^ Gj ^ Gk ^ Ga ^ Gb;
139
 
 
140
 
                    Gbc = Gb^Gc;
141
 
                    Gac = Ga^Gc;
142
 
                    Gba = Gb^Ga;
143
 
 
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;
150
 
 
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];
157
 
 
158
 
                          value_c = 0.0;
159
 
 
160
 
                          /** <ov||vv> --> connected triples **/
161
 
 
162
 
                          /* +t_jkbe * F_IeAc */
163
 
                          Ge = Gj ^ Gk ^ Gb;
164
 
                          for(e=0; e < bvirtpi[Ge]; e++) {
165
 
                            E = bvir_off[Ge] + e;
166
 
 
167
 
                            be = T2BB.params->colidx[B][E];
168
 
                            ie = FABints.params->rowidx[I][E];
169
 
 
170
 
                            t_jkbe = F_ieac = 0.0;
171
 
 
172
 
                            if(T2BB.params->rowtot[Gjk] && T2BB.params->coltot[Gjk])
173
 
                              t_jkbe = T2BB.matrix[Gjk][jk][be];
174
 
 
175
 
                            if(FABints.params->rowtot[Gac] && FABints.params->coltot[Gac])
176
 
                              F_ieac = FABints.matrix[Gac][ie][ac];
177
 
 
178
 
                            value_c += t_jkbe * F_ieac;
179
 
                          }
180
 
 
181
 
                          /* -t_jkce * F_IeAb */
182
 
                          Ge = Gj ^ Gk ^ Gc;
183
 
                          for(e=0; e < bvirtpi[Ge]; e++) {
184
 
                            E = bvir_off[Ge] + e;
185
 
 
186
 
                            ce = T2BB.params->colidx[C][E];
187
 
                            ie = FABints.params->rowidx[I][E];
188
 
 
189
 
                            t_jkce = F_ieab = 0.0;
190
 
 
191
 
                            if(T2BB.params->rowtot[Gjk] && T2BB.params->coltot[Gjk])
192
 
                              t_jkce = T2BB.matrix[Gjk][jk][ce];
193
 
 
194
 
                            if(FABints.params->rowtot[Gba] && FABints.params->coltot[Gba])
195
 
                              F_ieab = FABints.matrix[Gba][ie][ab];
196
 
 
197
 
                            value_c -= t_jkce * F_ieab;
198
 
                          }
199
 
 
200
 
                          /* +t_IkAe * F_jebc */
201
 
                          Ge = Gi ^ Gk ^ Ga;
202
 
                          for(e=0; e < bvirtpi[Ge]; e++) {
203
 
                            E = bvir_off[Ge] + e;
204
 
 
205
 
                            ae = T2AB.params->colidx[A][E];
206
 
                            je = FBBints.params->rowidx[J][E];
207
 
 
208
 
                            t_ikae = F_jebc = 0.0;
209
 
 
210
 
                            if(T2AB.params->rowtot[Gik] && T2AB.params->coltot[Gik])
211
 
                              t_ikae = T2AB.matrix[Gik][ik][ae];
212
 
 
213
 
                            if(FBBints.params->rowtot[Gbc] && FBBints.params->coltot[Gbc])
214
 
                              F_jebc = FBBints.matrix[Gbc][je][bc];
215
 
 
216
 
                            value_c += t_ikae * F_jebc;
217
 
                          }
218
 
 
219
 
                          /* -t_IkEb * F_jEcA */
220
 
                          Ge = Gi ^ Gk ^ Gb;
221
 
                          for(e=0; e < avirtpi[Ge]; e++) {
222
 
                            E = avir_off[Ge] + e;
223
 
 
224
 
                            eb = T2AB.params->colidx[E][B];
225
 
                            je = FBAints.params->rowidx[J][E];
226
 
 
227
 
                            t_ikeb = F_jeca = 0.0;
228
 
 
229
 
                            if(T2AB.params->rowtot[Gik] && T2AB.params->coltot[Gik])
230
 
                              t_ikeb = T2AB.matrix[Gik][ik][eb];
231
 
 
232
 
                            if(FBAints.params->rowtot[Gac] && FBAints.params->coltot[Gac])
233
 
                              F_jeca = FBAints.matrix[Gac][je][ca];
234
 
 
235
 
                            value_c -= t_ikeb * F_jeca;
236
 
                          }
237
 
 
238
 
                          /* +t_IkEc * F_jEbA */
239
 
                          Ge = Gi ^ Gk ^ Gc;
240
 
                          for(e=0; e < avirtpi[Ge]; e++) {
241
 
                            E = avir_off[Ge] + e;
242
 
 
243
 
                            ec = T2AB.params->colidx[E][C];
244
 
                            je = FBAints.params->rowidx[J][E];
245
 
 
246
 
                            t_ikec = F_jeba = 0.0;
247
 
 
248
 
                            if(T2AB.params->rowtot[Gik] && T2AB.params->coltot[Gik])
249
 
                              t_ikec = T2AB.matrix[Gik][ik][ec];
250
 
 
251
 
                            if(FBAints.params->rowtot[Gba] && FBAints.params->coltot[Gba])
252
 
                              F_jeba = FBAints.matrix[Gba][je][ba];
253
 
 
254
 
                            value_c += t_ikec * F_jeba;
255
 
                          }
256
 
 
257
 
                          /* -t_IjAe * F_kebc */
258
 
                          Ge = Gi ^ Gj ^ Ga;
259
 
                          for(e=0; e < bvirtpi[Ge]; e++) {
260
 
                            E = bvir_off[Ge] + e;
261
 
 
262
 
                            ae = T2AB.params->colidx[A][E];
263
 
                            ke = FBBints.params->rowidx[K][E];
264
 
 
265
 
                            t_ijae = F_kebc = 0.0;
266
 
 
267
 
                            if(T2AB.params->rowtot[Gij] && T2AB.params->coltot[Gij])
268
 
                              t_ijae = T2AB.matrix[Gij][ij][ae];
269
 
 
270
 
                            if(FBBints.params->rowtot[Gbc] && FBBints.params->coltot[Gbc])
271
 
                              F_kebc = FBBints.matrix[Gbc][ke][bc];
272
 
 
273
 
                            value_c -= t_ijae * F_kebc;
274
 
                          }
275
 
 
276
 
                          /* +t_IjEb * F_kEcA */
277
 
                          Ge = Gi ^ Gj ^ Gb;
278
 
                          for(e=0; e < avirtpi[Ge]; e++) {
279
 
                            E = avir_off[Ge] + e;
280
 
 
281
 
                            eb = T2AB.params->colidx[E][B];
282
 
                            ke = FBAints.params->rowidx[K][E];
283
 
 
284
 
                            t_ijeb = F_keca = 0.0;
285
 
 
286
 
                            if(T2AB.params->rowtot[Gij] && T2AB.params->coltot[Gij])
287
 
                              t_ijeb = T2AB.matrix[Gij][ij][eb];
288
 
 
289
 
                            if(FBAints.params->rowtot[Gac] && FBAints.params->coltot[Gac])
290
 
                              F_keca = FBAints.matrix[Gac][ke][ca];
291
 
 
292
 
                            value_c += t_ijeb * F_keca;
293
 
                          }
294
 
 
295
 
                          /* -t_IjEc * F_kEbA */
296
 
                          Ge = Gi ^ Gj ^ Gc;
297
 
                          for(e=0; e < avirtpi[Ge]; e++) {
298
 
                            E = avir_off[Ge] + e;
299
 
 
300
 
                            ec = T2AB.params->colidx[E][C];
301
 
                            ke = FBAints.params->rowidx[K][E];
302
 
 
303
 
                            t_ijec = F_keba = 0.0;
304
 
 
305
 
                            if(T2AB.params->rowtot[Gij] && T2AB.params->coltot[Gij])
306
 
                              t_ijec = T2AB.matrix[Gij][ij][ec];
307
 
 
308
 
                            if(FBAints.params->rowtot[Gba] && FBAints.params->coltot[Gba])
309
 
                              F_keba = FBAints.matrix[Gba][ke][ba];
310
 
 
311
 
                            value_c -= t_ijec * F_keba;
312
 
                          }
313
 
 
314
 
                          /** <oo||ov> --> connected triples **/
315
 
 
316
 
                          /* +t_ImAc * E_jkmb */
317
 
                          Gm = Gi ^ Ga ^ Gc;
318
 
                          for(m=0; m < boccpi[Gm]; m++) {
319
 
                            M = bocc_off[Gm] + m;
320
 
 
321
 
                            im = T2AB.params->rowidx[I][M];
322
 
                            mb = EBBints.params->colidx[M][B];
323
 
 
324
 
                            t_imac = E_jkmb = 0.0;
325
 
 
326
 
                            if(T2AB.params->rowtot[Gac] && T2AB.params->coltot[Gac])
327
 
                              t_imac = T2AB.matrix[Gac][im][ac];
328
 
 
329
 
                            if(EBBints.params->rowtot[Gjk] && EBBints.params->coltot[Gjk])
330
 
                              E_jkmb = EBBints.matrix[Gjk][jk][mb];
331
 
 
332
 
                            value_c += t_imac * E_jkmb;
333
 
                          }
334
 
 
335
 
                          /* -t_ImAb * E_jkmc */
336
 
                          Gm = Gi ^ Gb ^ Ga;
337
 
                          for(m=0; m < boccpi[Gm]; m++) {
338
 
                            M = bocc_off[Gm] + m;
339
 
 
340
 
                            im = T2AB.params->rowidx[I][M];
341
 
                            mc = EBBints.params->colidx[M][C];
342
 
 
343
 
                            t_imab = E_jkmc = 0.0;
344
 
 
345
 
                            if(T2AB.params->rowtot[Gba] && T2AB.params->coltot[Gba])
346
 
                              t_imab = T2AB.matrix[Gba][im][ab];
347
 
 
348
 
                            if(EBBints.params->rowtot[Gjk] && EBBints.params->coltot[Gjk])
349
 
                              E_jkmc = EBBints.matrix[Gjk][jk][mc];
350
 
 
351
 
                            value_c -= t_imab * E_jkmc;
352
 
                          }
353
 
 
354
 
                          /* -t_jmbc * E_kImA */
355
 
                          Gm = Gj ^ Gb ^ Gc;
356
 
                          for(m=0; m < boccpi[Gm]; m++) {
357
 
                            M = bocc_off[Gm] + m;
358
 
 
359
 
                            jm = T2BB.params->rowidx[J][M];
360
 
                            ma = EBAints.params->colidx[M][A];
361
 
 
362
 
                            t_jmbc = E_kima = 0.0;
363
 
 
364
 
                            if(T2BB.params->rowtot[Gbc] && T2BB.params->coltot[Gbc])
365
 
                              t_jmbc = T2BB.matrix[Gbc][jm][bc];
366
 
 
367
 
                            if(EBAints.params->rowtot[Gik] && EBAints.params->coltot[Gik])
368
 
                              E_kima = EBAints.matrix[Gik][ki][ma];
369
 
 
370
 
                            value_c -= t_jmbc * E_kima;
371
 
                          }
372
 
 
373
 
                          /* +t_MjAc * E_IkMb */
374
 
                          Gm = Gj ^ Ga ^ Gc;
375
 
                          for(m=0; m < aoccpi[Gm]; m++) {
376
 
                            M = aocc_off[Gm] + m;
377
 
 
378
 
                            mj = T2AB.params->rowidx[M][J];
379
 
                            mb = EABints.params->colidx[M][B];
380
 
 
381
 
                            t_mjac = E_ikmb = 0.0;
382
 
 
383
 
                            if(T2AB.params->rowtot[Gac] && T2AB.params->coltot[Gac])
384
 
                              t_mjac = T2AB.matrix[Gac][mj][ac];
385
 
 
386
 
                            if(EABints.params->rowtot[Gik] && EABints.params->coltot[Gik])
387
 
                              E_ikmb = EABints.matrix[Gik][ik][mb];
388
 
 
389
 
                            value_c += t_mjac * E_ikmb;
390
 
                          }
391
 
 
392
 
                          /* -t_MjAb * E_IkMc */
393
 
                          Gm = Gj ^ Gb ^ Ga;
394
 
                          for(m=0; m < aoccpi[Gm]; m++) {
395
 
                            M = aocc_off[Gm] + m;
396
 
 
397
 
                            mj = T2AB.params->rowidx[M][J];
398
 
                            mc = EABints.params->colidx[M][C];
399
 
 
400
 
                            t_mjab = E_ikmc = 0.0;
401
 
 
402
 
                            if(T2AB.params->rowtot[Gba] && T2AB.params->coltot[Gba])
403
 
                              t_mjab = T2AB.matrix[Gba][mj][ab];
404
 
 
405
 
                            if(EABints.params->rowtot[Gik] && EABints.params->coltot[Gik])
406
 
                              E_ikmc = EABints.matrix[Gik][ik][mc];
407
 
 
408
 
                            value_c -= t_mjab * E_ikmc;
409
 
                          }
410
 
 
411
 
                          /* +t_kmbc * E_jImA */
412
 
                          Gm = Gk ^ Gb ^ Gc;
413
 
                          for(m=0; m < boccpi[Gm]; m++) {
414
 
                            M = bocc_off[Gm] + m;
415
 
 
416
 
                            km = T2BB.params->rowidx[K][M];
417
 
                            ma = EBAints.params->colidx[M][A];
418
 
 
419
 
                            t_kmbc = E_jima = 0.0;
420
 
 
421
 
                            if(T2BB.params->rowtot[Gbc] && T2BB.params->coltot[Gbc])
422
 
                              t_kmbc = T2BB.matrix[Gbc][km][bc];
423
 
 
424
 
                            if(EBAints.params->rowtot[Gji] && EBAints.params->coltot[Gji])
425
 
                              E_jima = EBAints.matrix[Gji][ji][ma];
426
 
 
427
 
                            value_c += t_kmbc * E_jima;
428
 
                          }
429
 
 
430
 
                          /* -t_MkAc * E_IjMb */
431
 
                          Gm = Gk ^ Ga ^ Gc;
432
 
                          for(m=0; m < aoccpi[Gm]; m++) {
433
 
                            M = aocc_off[Gm] + m;
434
 
 
435
 
                            mk = T2AB.params->rowidx[M][K];
436
 
                            mb = EABints.params->colidx[M][B];
437
 
 
438
 
                            t_mkac = E_ijmb = 0.0;
439
 
 
440
 
                            if(T2AB.params->rowtot[Gac] && T2AB.params->coltot[Gac])
441
 
                              t_mkac = T2AB.matrix[Gac][mk][ac];
442
 
 
443
 
                            if(EABints.params->rowtot[Gji] && EABints.params->coltot[Gji])
444
 
                              E_ijmb = EABints.matrix[Gji][ij][mb];
445
 
 
446
 
                            value_c -= t_mkac * E_ijmb;
447
 
                          }
448
 
 
449
 
                          /* +t_MkAb * E_IjMc */
450
 
                          Gm = Gk ^ Gb ^ Ga;
451
 
                          for(m=0; m < aoccpi[Gm]; m++) {
452
 
                            M = aocc_off[Gm] + m;
453
 
 
454
 
                            mk = T2AB.params->rowidx[M][K];
455
 
                            mc = EABints.params->colidx[M][C];
456
 
 
457
 
                            t_mkab = E_ijmc = 0.0;
458
 
 
459
 
                            if(T2AB.params->rowtot[Gba] && T2AB.params->coltot[Gba])
460
 
                              t_mkab = T2AB.matrix[Gba][mk][ab];
461
 
 
462
 
                            if(EABints.params->rowtot[Gji] && EABints.params->coltot[Gji])
463
 
                              E_ijmc = EABints.matrix[Gji][ij][mc];
464
 
 
465
 
                            value_c += t_mkab * E_ijmc;
466
 
                          }
467
 
 
468
 
                          /** disconnected triples **/
469
 
 
470
 
                          value_d = 0.0;
471
 
 
472
 
                          /* +t_IA * D_jkbc */
473
 
                          if(Gi == Ga && Gjk == Gbc) {
474
 
                            t_ia = D_jkbc = 0.0;
475
 
 
476
 
                            if(T1A.params->rowtot[Gi] && T1A.params->coltot[Gi])
477
 
                              t_ia = T1A.matrix[Gi][i][a];
478
 
 
479
 
                            if(DBBints.params->rowtot[Gjk] && DBBints.params->coltot[Gjk])
480
 
                              D_jkbc = DBBints.matrix[Gjk][jk][bc];
481
 
 
482
 
                            value_d += t_ia * D_jkbc;
483
 
                          }
484
 
 
485
 
                          /* +t_jb * D_IkAc */
486
 
                          if(Gj == Gb && Gik == Gac) {
487
 
                            t_jb = D_ikac = 0.0;
488
 
 
489
 
                            if(T1B.params->rowtot[Gj] && T1B.params->coltot[Gj])
490
 
                              t_jb = T1B.matrix[Gj][j][b];
491
 
 
492
 
                            if(DABints.params->rowtot[Gik] && DABints.params->coltot[Gik])
493
 
                              D_ikac = DABints.matrix[Gik][ik][ac];
494
 
 
495
 
                            value_d += t_jb * D_ikac;
496
 
                          }
497
 
 
498
 
                          /* -t_jc * D_IkAb */
499
 
                          if(Gj == Gc && Gik == Gba) {
500
 
                            t_jc = D_ikab = 0.0;
501
 
 
502
 
                            if(T1B.params->rowtot[Gj] && T1B.params->coltot[Gj])
503
 
                              t_jc = T1B.matrix[Gj][j][c];
504
 
 
505
 
                            if(DABints.params->rowtot[Gik] && DABints.params->coltot[Gik])
506
 
                              D_ikab = DABints.matrix[Gik][ik][ab];
507
 
 
508
 
                            value_d -= t_jc * D_ikab;
509
 
                          }
510
 
 
511
 
                          /* -t_kb * D_IjAc */
512
 
                          if(Gk == Gb && Gji == Gac) {
513
 
                            t_kb = D_ijac = 0.0;
514
 
 
515
 
                            if(T1B.params->rowtot[Gk] && T1B.params->coltot[Gk])
516
 
                              t_kb = T1B.matrix[Gk][k][b];
517
 
 
518
 
                            if(DABints.params->rowtot[Gji] && DABints.params->coltot[Gji])
519
 
                              D_ijac = DABints.matrix[Gji][ij][ac];
520
 
 
521
 
                            value_d -= t_kb * D_ijac;
522
 
                          }
523
 
 
524
 
                          /* +t_kc * D_IjAb */
525
 
                          if(Gk == Gc && Gji == Gba) {
526
 
                            t_kc = D_ijab = 0.0;
527
 
 
528
 
                            if(T1B.params->rowtot[Gk] && T1B.params->coltot[Gk])
529
 
                              t_kc = T1B.matrix[Gk][k][c];
530
 
 
531
 
                            if(DABints.params->rowtot[Gji] && DABints.params->coltot[Gji])
532
 
                              D_ijab = DABints.matrix[Gji][ij][ab];
533
 
 
534
 
                            value_d += t_kc * D_ijab;
535
 
                          }
536
 
 
537
 
                          /*
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);
540
 
                          }
541
 
                          */
542
 
 
543
 
                          /* Compute the Fock denominator */
544
 
                          denom = 0.0;
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];
557
 
 
558
 
                          /*                    ET_ABB += (value_d + value_c) * value_c / denom; */
559
 
                          ET_ABB += (value_c + value_d) * value_c / denom;
560
 
 
561
 
                        } /* c */
562
 
                      } /* ab */
563
 
                    } /* Gab */
564
 
 
565
 
                  } /* J >= K */
566
 
 
567
 
                } /* k */
568
 
              } /* j */
569
 
            } /* i */
570
 
 
571
 
          } /* Gb */
572
 
        } /* Ga */
573
 
 
574
 
      } /* Gk */
575
 
    } /* Gj */
576
 
  } /* Gi */
577
 
 
578
 
  /*  fprintf(outfile, "cnt = %d\n", cnt); */
579
 
  ET_ABB /= 2.0;
580
 
 
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);
592
 
  }
593
 
 
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);
604
 
 
605
 
  dpd_file2_mat_close(&T1A);
606
 
  dpd_file2_close(&T1A);
607
 
  dpd_file2_mat_close(&T1B);
608
 
  dpd_file2_close(&T1B);
609
 
 
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);
618
 
 
619
 
  return ET_ABB;
620
 
}