← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/cctriples/ET_UHF_AAB_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:
src/bin/cctriples/ET_UHF_AAB_noddy.c
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_AAB_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;
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, ec, ke, ie, je;
18
 
  int im, jm, 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_AAB;
22
 
  double t_ijae, t_ijbe, t_jkae, t_jkbe, t_jkec, t_ikae, t_ikbe, t_ikec;
23
 
  double F_kecb, F_keca, F_iebc, F_ieac, F_ieab, F_jebc, F_jeac, F_jeab;
24
 
  double t_imbc, t_imac, t_imba, t_jmbc, t_jmac, t_jmba, t_mkbc, t_mkac;
25
 
  double E_kjma, E_kjmb, E_jkmc, E_kima, E_kimb, E_ikmc, E_jima, E_jimb;
26
 
  double t_ia, t_ib, t_ja, t_jb, t_kc;
27
 
  double D_jkbc, D_jkac, D_ikbc, D_ikac, D_jiba;
28
 
  dpdbuf4 T2AB, T2AA;
29
 
  dpdbuf4 FAAints, FABints, FBAints;
30
 
  dpdbuf4 EAAints, EABints, EBAints;
31
 
  dpdbuf4 DAAints, 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(&T2AA, CC_TAMPS, 0, 0, 5, 2, 7, 0, "tIJAB");
65
 
  dpd_buf4_init(&T2AB, CC_TAMPS, 0, 22, 28, 22, 28, 0, "tIjAb");
66
 
 
67
 
  dpd_buf4_init(&FAAints, CC_FINTS, 0, 20, 5, 20, 5, 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(&EAAints, CC_EINTS, 0, 0, 20, 2, 20, 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(&DAAints, CC_DINTS, 0, 0, 5, 0, 5, 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(&T2AA, h);
80
 
    dpd_buf4_mat_irrep_rd(&T2AA, 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(&FAAints, h);
86
 
    dpd_buf4_mat_irrep_rd(&FAAints, 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(&EAAints, h);
95
 
    dpd_buf4_mat_irrep_rd(&EAAints, 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(&DAAints, h);
104
 
    dpd_buf4_mat_irrep_rd(&DAAints, 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_AAB = 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
 
 
117
 
        Gij = Gji = Gi ^ Gj;
118
 
        Gjk = Gj ^ Gk;
119
 
        Gik = Gi ^ Gk;
120
 
 
121
 
        for(i=0; i < aoccpi[Gi]; i++) {
122
 
          I = aocc_off[Gi] + i;
123
 
          for(j=0; j < aoccpi[Gj]; j++) {
124
 
            J = aocc_off[Gj] + j;
125
 
            for(k=0; k < boccpi[Gk]; k++) {
126
 
              K = bocc_off[Gk] + k;
127
 
 
128
 
              if(I >= J) {
129
 
 
130
 
                ij = EAAints.params->rowidx[I][J];
131
 
                ji = EAAints.params->rowidx[J][I];
132
 
                jk = EABints.params->rowidx[J][K];
133
 
                kj = EBAints.params->rowidx[K][J];
134
 
                ik = EABints.params->rowidx[I][K];
135
 
                ki = EBAints.params->rowidx[K][I];
136
 
 
137
 
                for(Ga=0; Ga < nirreps; Ga++) {
138
 
                  for(Gb=0; Gb < nirreps; Gb++) {
139
 
                    Gc = Gi ^ Gj ^ Gk ^ Ga ^ Gb;
140
 
 
141
 
                    Gbc = Gb^Gc;
142
 
                    Gac = Ga^Gc;
143
 
                    Gba = Gb^Ga;
144
 
 
145
 
                    for(a=0; a < avirtpi[Ga]; a++) {
146
 
                      A = avir_off[Ga] + a;
147
 
                      for(b=0; b < avirtpi[Gb]; b++) {
148
 
                        B = avir_off[Gb] + b;
149
 
                        for(c=0; c < bvirtpi[Gc]; c++) {
150
 
                          C = bvir_off[Gc] + c;
151
 
 
152
 
                          /*                      if(A >= B) { */
153
 
 
154
 
                            ab = FAAints.params->colidx[A][B];
155
 
                            ba = FAAints.params->colidx[B][A];
156
 
                            bc = FABints.params->colidx[B][C];
157
 
                            cb = FBAints.params->colidx[C][B];
158
 
                            ac = FABints.params->colidx[A][C];
159
 
                            ca = FBAints.params->colidx[C][A];
160
 
 
161
 
                            value_c = 0.0;
162
 
 
163
 
                            /** <ov||vv> --> connected triples **/
164
 
 
165
 
                            /* -t_JkAe * F_IeBc */
166
 
                            Ge = Gj ^ Gk ^ Ga;
167
 
                            for(e=0; e < bvirtpi[Ge]; e++) {
168
 
                              E = bvir_off[Ge] + e;
169
 
 
170
 
                              ae = T2AB.params->colidx[A][E];
171
 
                              ie = FABints.params->rowidx[I][E];
172
 
 
173
 
                              t_jkae = F_iebc = 0.0;
174
 
 
175
 
                              if(T2AB.params->rowtot[Gjk] && T2AB.params->coltot[Gjk])
176
 
                                t_jkae = T2AB.matrix[Gjk][jk][ae];
177
 
 
178
 
                              if(FABints.params->rowtot[Gbc] && FABints.params->coltot[Gbc])
179
 
                                F_iebc = FABints.matrix[Gbc][ie][bc];
180
 
 
181
 
                              value_c -= t_jkae * F_iebc;
182
 
                            }
183
 
 
184
 
                            /* +t_JkBe * F_IeAc */
185
 
                            Ge = Gj ^ Gk ^ Gb;
186
 
                            for(e=0; e < bvirtpi[Ge]; e++) {
187
 
                              E = bvir_off[Ge] + e;
188
 
 
189
 
                              be = T2AB.params->colidx[B][E];
190
 
                              ie = FABints.params->rowidx[I][E];
191
 
 
192
 
                              t_jkbe = F_ieac = 0.0;
193
 
 
194
 
                              if(T2AB.params->rowtot[Gjk] && T2AB.params->coltot[Gjk])
195
 
                                t_jkbe = T2AB.matrix[Gjk][jk][be];
196
 
 
197
 
                              if(FABints.params->rowtot[Gac] && FABints.params->coltot[Gac])
198
 
                                F_ieac = FABints.matrix[Gac][ie][ac];
199
 
 
200
 
                              value_c += t_jkbe * F_ieac;
201
 
                            }
202
 
 
203
 
                            /* +t_JkEc * F_IEAB */
204
 
                            Ge = Gj ^ Gk ^ Gc;
205
 
                            for(e=0; e < avirtpi[Ge]; e++) {
206
 
                              E = avir_off[Ge] + e;
207
 
 
208
 
                              ec = T2AB.params->colidx[E][C];
209
 
                              ie = FAAints.params->rowidx[I][E];
210
 
 
211
 
                              t_jkec = F_ieab = 0.0;
212
 
 
213
 
                              if(T2AB.params->rowtot[Gjk] && T2AB.params->coltot[Gjk])
214
 
                                t_jkec = T2AB.matrix[Gjk][jk][ec];
215
 
 
216
 
                              if(FAAints.params->rowtot[Gba] && FAAints.params->coltot[Gba])
217
 
                                F_ieab = FAAints.matrix[Gba][ie][ab];
218
 
 
219
 
                              value_c += t_jkec * F_ieab;
220
 
                            }
221
 
 
222
 
                            /* +t_IkAe * F_JeBc */
223
 
                            Ge = Gi ^ Gk ^ Ga;
224
 
                            for(e=0; e < bvirtpi[Ge]; e++) {
225
 
                              E = bvir_off[Ge] + e;
226
 
 
227
 
                              ae = T2AB.params->colidx[A][E];
228
 
                              je = FABints.params->rowidx[J][E];
229
 
 
230
 
                              t_ikae = F_jebc = 0.0;
231
 
 
232
 
                              if(T2AB.params->rowtot[Gik] && T2AB.params->coltot[Gik])
233
 
                                t_ikae = T2AB.matrix[Gik][ik][ae];
234
 
 
235
 
                              if(FABints.params->rowtot[Gbc] && FABints.params->coltot[Gbc])
236
 
                                F_jebc = FABints.matrix[Gbc][je][bc];
237
 
 
238
 
                              value_c += t_ikae * F_jebc;
239
 
                            }
240
 
 
241
 
                            /* -t_IkBe * F_JeAc */
242
 
                            Ge = Gi ^ Gk ^ Gb;
243
 
                            for(e=0; e < bvirtpi[Ge]; e++) {
244
 
                              E = bvir_off[Ge] + e;
245
 
 
246
 
                              be = T2AB.params->colidx[B][E];
247
 
                              je = FABints.params->rowidx[J][E];
248
 
 
249
 
                              t_ikbe = F_jeac = 0.0;
250
 
 
251
 
                              if(T2AB.params->rowtot[Gik] && T2AB.params->coltot[Gik])
252
 
                                t_ikbe = T2AB.matrix[Gik][ik][be];
253
 
 
254
 
                              if(FABints.params->rowtot[Gac] && FABints.params->coltot[Gac])
255
 
                                F_jeac = FABints.matrix[Gac][je][ac];
256
 
 
257
 
                              value_c -= t_ikbe * F_jeac;
258
 
                            }
259
 
 
260
 
                            /* -t_IkEc * F_JEAB */
261
 
                            Ge = Gi ^ Gk ^ Gc;
262
 
                            for(e=0; e < avirtpi[Ge]; e++) {
263
 
                              E = avir_off[Ge] + e;
264
 
 
265
 
                              ec = T2AB.params->colidx[E][C];
266
 
                              je = FAAints.params->rowidx[J][E];
267
 
 
268
 
                              t_ikec = F_jeab = 0.0;
269
 
 
270
 
                              if(T2AB.params->rowtot[Gik] && T2AB.params->coltot[Gik])
271
 
                                t_ikec = T2AB.matrix[Gik][ik][ec];
272
 
 
273
 
                              if(FAAints.params->rowtot[Gba] && FAAints.params->coltot[Gba])
274
 
                                F_jeab = FAAints.matrix[Gba][je][ab];
275
 
 
276
 
                              value_c -= t_ikec * F_jeab;
277
 
                            }
278
 
 
279
 
                            /* +t_IJAE * F_kEcB */
280
 
                            Ge = Gi ^ Gj ^ Ga;
281
 
                            for(e=0; e < avirtpi[Ge]; e++) {
282
 
                              E = avir_off[Ge] + e;
283
 
 
284
 
                              ae = T2AA.params->colidx[A][E];
285
 
                              ke = FBAints.params->rowidx[K][E];
286
 
 
287
 
                              t_ijae = F_kecb = 0.0;
288
 
 
289
 
                              if(T2AA.params->rowtot[Gij] && T2AA.params->coltot[Gij])
290
 
                                t_ijae = T2AA.matrix[Gij][ij][ae];
291
 
 
292
 
                              if(FBAints.params->rowtot[Gbc] && FBAints.params->coltot[Gbc])
293
 
                                F_kecb = FBAints.matrix[Gbc][ke][cb];
294
 
 
295
 
                              value_c += t_ijae * F_kecb;
296
 
                            }
297
 
 
298
 
                            /* -t_IJBE * F_kEcA */
299
 
                            Ge = Gi ^ Gj ^ Gb;
300
 
                            for(e=0; e < avirtpi[Ge]; e++) {
301
 
                              E = avir_off[Ge] + e;
302
 
 
303
 
                              be = T2AA.params->colidx[B][E];
304
 
                              ke = FBAints.params->rowidx[K][E];
305
 
 
306
 
                              t_ijbe = F_keca = 0.0;
307
 
 
308
 
                              if(T2AA.params->rowtot[Gij] && T2AA.params->coltot[Gij])
309
 
                                t_ijbe = T2AA.matrix[Gij][ij][be];
310
 
 
311
 
                              if(FBAints.params->rowtot[Gac] && FBAints.params->coltot[Gac])
312
 
                                F_keca = FBAints.matrix[Gac][ke][ca];
313
 
 
314
 
                              value_c -= t_ijbe * F_keca;
315
 
                            }
316
 
 
317
 
                            /** <oo||ov> --> connected triples **/
318
 
 
319
 
                            /* +t_ImBc * E_kJmA */
320
 
                            Gm = Gi ^ Gb ^ Gc;
321
 
                            for(m=0; m < boccpi[Gm]; m++) {
322
 
                              M = bocc_off[Gm] + m;
323
 
 
324
 
                              im = T2AB.params->rowidx[I][M];
325
 
                              ma = EBAints.params->colidx[M][A];
326
 
 
327
 
                              t_imbc = E_kjma = 0.0;
328
 
 
329
 
                              if(T2AB.params->rowtot[Gbc] && T2AB.params->coltot[Gbc])
330
 
                                t_imbc = T2AB.matrix[Gbc][im][bc];
331
 
 
332
 
                              if(EBAints.params->rowtot[Gjk] && EBAints.params->coltot[Gjk])
333
 
                                E_kjma = EBAints.matrix[Gjk][kj][ma];
334
 
 
335
 
                              value_c += t_imbc * E_kjma;
336
 
                            }
337
 
 
338
 
                            /* -t_ImAc * E_kJmB */
339
 
                            Gm = Gi ^ Ga ^ Gc;
340
 
                            for(m=0; m < boccpi[Gm]; m++) {
341
 
                              M = bocc_off[Gm] + m;
342
 
 
343
 
                              im = T2AB.params->rowidx[I][M];
344
 
                              mb = EBAints.params->colidx[M][B];
345
 
 
346
 
                              t_imac = E_kjmb = 0.0;
347
 
 
348
 
                              if(T2AB.params->rowtot[Gac] && T2AB.params->coltot[Gac])
349
 
                                t_imac = T2AB.matrix[Gac][im][ac];
350
 
 
351
 
                              if(EBAints.params->rowtot[Gjk] && EBAints.params->coltot[Gjk])
352
 
                                E_kjmb = EBAints.matrix[Gjk][kj][mb];
353
 
 
354
 
                              value_c -= t_imac * E_kjmb;
355
 
                            }
356
 
 
357
 
                            /* +t_IMBA * E_JkMc */
358
 
                            Gm = Gi ^ Gb ^ Ga;
359
 
                            for(m=0; m < aoccpi[Gm]; m++) {
360
 
                              M = aocc_off[Gm] + m;
361
 
 
362
 
                              im = T2AA.params->rowidx[I][M];
363
 
                              mc = EABints.params->colidx[M][C];
364
 
 
365
 
                              t_imba = E_jkmc = 0.0;
366
 
 
367
 
                              if(T2AA.params->rowtot[Gba] && T2AA.params->coltot[Gba])
368
 
                                t_imba = T2AA.matrix[Gba][im][ba];
369
 
 
370
 
                              if(EABints.params->rowtot[Gjk] && EABints.params->coltot[Gjk])
371
 
                                E_jkmc = EABints.matrix[Gjk][jk][mc];
372
 
 
373
 
                              value_c += t_imba * E_jkmc;
374
 
                            }
375
 
 
376
 
                            /* -t_JmBc * E_kImA */
377
 
                            Gm = Gj ^ Gb ^ Gc;
378
 
                            for(m=0; m < boccpi[Gm]; m++) {
379
 
                              M = bocc_off[Gm] + m;
380
 
 
381
 
                              jm = T2AB.params->rowidx[J][M];
382
 
                              ma = EBAints.params->colidx[M][A];
383
 
 
384
 
                              t_jmbc = E_kima = 0.0;
385
 
 
386
 
                              if(T2AB.params->rowtot[Gbc] && T2AB.params->coltot[Gbc])
387
 
                                t_jmbc = T2AB.matrix[Gbc][jm][bc];
388
 
 
389
 
                              if(EBAints.params->rowtot[Gik] && EBAints.params->coltot[Gik])
390
 
                                E_kima = EBAints.matrix[Gik][ki][ma];
391
 
 
392
 
                              value_c -= t_jmbc * E_kima;
393
 
                            }
394
 
 
395
 
                            /* +t_JmAc * E_kImB */
396
 
                            Gm = Gj ^ Ga ^ Gc;
397
 
                            for(m=0; m < boccpi[Gm]; m++) {
398
 
                              M = bocc_off[Gm] + m;
399
 
 
400
 
                              jm = T2AB.params->rowidx[J][M];
401
 
                              mb = EBAints.params->colidx[M][B];
402
 
 
403
 
                              t_jmac = E_kimb = 0.0;
404
 
 
405
 
                              if(T2AB.params->rowtot[Gac] && T2AB.params->coltot[Gac])
406
 
                                t_jmac = T2AB.matrix[Gac][jm][ac];
407
 
 
408
 
                              if(EBAints.params->rowtot[Gik] && EBAints.params->coltot[Gik])
409
 
                                E_kimb = EBAints.matrix[Gik][ki][mb];
410
 
 
411
 
                              value_c += t_jmac * E_kimb;
412
 
                            }
413
 
 
414
 
                            /* -t_JMBA * E_IkMc */
415
 
                            Gm = Gj ^ Gb ^ Ga;
416
 
                            for(m=0; m < aoccpi[Gm]; m++) {
417
 
                              M = aocc_off[Gm] + m;
418
 
 
419
 
                              jm = T2AA.params->rowidx[J][M];
420
 
                              mc = EABints.params->colidx[M][C];
421
 
 
422
 
                              t_jmba = E_ikmc = 0.0;
423
 
 
424
 
                              if(T2AA.params->rowtot[Gba] && T2AA.params->coltot[Gba])
425
 
                                t_jmba = T2AA.matrix[Gba][jm][ba];
426
 
 
427
 
                              if(EABints.params->rowtot[Gik] && EABints.params->coltot[Gik])
428
 
                                E_ikmc = EABints.matrix[Gik][ik][mc];
429
 
 
430
 
                              value_c -= t_jmba * E_ikmc;
431
 
                            }
432
 
 
433
 
                            /* -t_MkBc * E_JIMA */
434
 
                            Gm = Gk ^ Gb ^ Gc;
435
 
                            for(m=0; m < aoccpi[Gm]; m++) {
436
 
                              M = aocc_off[Gm] + m;
437
 
 
438
 
                              mk = T2AB.params->rowidx[M][K];
439
 
                              ma = EAAints.params->colidx[M][A];
440
 
 
441
 
                              t_mkbc = E_jima = 0.0;
442
 
 
443
 
                              if(T2AB.params->rowtot[Gbc] && T2AB.params->coltot[Gbc])
444
 
                                t_mkbc = T2AB.matrix[Gbc][mk][bc];
445
 
 
446
 
                              if(EAAints.params->rowtot[Gji] && EAAints.params->coltot[Gji])
447
 
                                E_jima = EAAints.matrix[Gji][ji][ma];
448
 
 
449
 
                              value_c -= t_mkbc * E_jima;
450
 
                            }
451
 
 
452
 
                            /* +t_MkAc * E_JIMB */
453
 
                            Gm = Gk ^ Ga ^ Gc;
454
 
                            for(m=0; m < aoccpi[Gm]; m++) {
455
 
                              M = aocc_off[Gm] + m;
456
 
 
457
 
                              mk = T2AB.params->rowidx[M][K];
458
 
                              mb = EAAints.params->colidx[M][B];
459
 
 
460
 
                              t_mkac = E_jimb = 0.0;
461
 
 
462
 
                              if(T2AB.params->rowtot[Gac] && T2AB.params->coltot[Gac])
463
 
                                t_mkac = T2AB.matrix[Gac][mk][ac];
464
 
 
465
 
                              if(EAAints.params->rowtot[Gji] && EAAints.params->coltot[Gji])
466
 
                                E_jimb = EAAints.matrix[Gji][ji][mb];
467
 
 
468
 
                              value_c += t_mkac * E_jimb;
469
 
                            }
470
 
 
471
 
                            /** disconnected triples **/
472
 
 
473
 
                            value_d = 0.0;
474
 
 
475
 
                            /* +t_IA * D_JkBc */
476
 
                            if(Gi == Ga && Gjk == Gbc) {
477
 
                              t_ia = D_jkbc = 0.0;
478
 
 
479
 
                              if(T1A.params->rowtot[Gi] && T1A.params->coltot[Gi])
480
 
                                t_ia = T1A.matrix[Gi][i][a];
481
 
 
482
 
                              if(DABints.params->rowtot[Gjk] && DABints.params->coltot[Gjk])
483
 
                                D_jkbc = DABints.matrix[Gjk][jk][bc];
484
 
 
485
 
                              value_d += t_ia * D_jkbc;
486
 
                            }
487
 
 
488
 
                            /* -t_IB * D_JkAc */
489
 
                            if(Gi == Gb && Gjk == Gac) {
490
 
                              t_ib = D_jkac = 0.0;
491
 
 
492
 
                              if(T1A.params->rowtot[Gi] && T1A.params->coltot[Gi])
493
 
                                t_ib = T1A.matrix[Gi][i][b];
494
 
 
495
 
                              if(DABints.params->rowtot[Gjk] && DABints.params->coltot[Gjk])
496
 
                                D_jkac = DABints.matrix[Gjk][jk][ac];
497
 
 
498
 
                              value_d -= t_ib * D_jkac;
499
 
                            }
500
 
 
501
 
                            /* -t_JA * D_IkBc */
502
 
                            if(Gj == Ga && Gik == Gbc) {
503
 
                              t_ja = D_ikbc = 0.0;
504
 
 
505
 
                              if(T1A.params->rowtot[Gj] && T1A.params->coltot[Gj])
506
 
                                t_ja = T1A.matrix[Gj][j][a];
507
 
 
508
 
                              if(DABints.params->rowtot[Gik] && DABints.params->coltot[Gik])
509
 
                                D_ikbc = DABints.matrix[Gik][ik][bc];
510
 
 
511
 
                              value_d -= t_ja * D_ikbc;
512
 
                            }
513
 
 
514
 
                            /* +t_JB * D_IkAc */
515
 
                            if(Gj == Gb && Gik == Gac) {
516
 
                              t_jb = D_ikac = 0.0;
517
 
 
518
 
                              if(T1A.params->rowtot[Gj] && T1A.params->coltot[Gj])
519
 
                                t_jb = T1A.matrix[Gj][j][b];
520
 
 
521
 
                              if(DABints.params->rowtot[Gik] && DABints.params->coltot[Gik])
522
 
                                D_ikac = DABints.matrix[Gik][ik][ac];
523
 
 
524
 
                              value_d += t_jb * D_ikac;
525
 
                            }
526
 
 
527
 
                            /* +t_kc * D_JIBA */
528
 
                            if(Gk == Gc && Gji == Gba) {
529
 
                              t_kc = D_jiba = 0.0;
530
 
 
531
 
                              if(T1B.params->rowtot[Gk] && T1B.params->coltot[Gk])
532
 
                                t_kc = T1B.matrix[Gk][k][c];
533
 
 
534
 
                              if(DAAints.params->rowtot[Gji] && DAAints.params->coltot[Gji])
535
 
                                D_jiba = DAAints.matrix[Gji][ji][ba];
536
 
 
537
 
                              value_d += t_kc * D_jiba;
538
 
                            }
539
 
 
540
 
                            /*
541
 
                              if(fabs(value_c) > 1e-7) {
542
 
                              cnt++;
543
 
                              fprintf(outfile, "%d %d %d %d %d %d %20.14f\n", I, J, K, A, B, C, value_c);
544
 
                              }
545
 
                            */
546
 
 
547
 
                            /* Compute the Fock denominator */
548
 
                            denom = 0.0;
549
 
                            if(fIJ.params->rowtot[Gi])
550
 
                              denom += fIJ.matrix[Gi][i][i];
551
 
                            if(fIJ.params->rowtot[Gj])
552
 
                              denom += fIJ.matrix[Gj][j][j];
553
 
                            if(fij.params->rowtot[Gk])
554
 
                              denom += fij.matrix[Gk][k][k];
555
 
                            if(fAB.params->rowtot[Ga])
556
 
                              denom -= fAB.matrix[Ga][a][a];
557
 
                            if(fAB.params->rowtot[Gb])
558
 
                              denom -= fAB.matrix[Gb][b][b];
559
 
                            if(fab.params->rowtot[Gc])
560
 
                              denom -= fab.matrix[Gc][c][c];
561
 
 
562
 
                            /*
563
 
                            if(fabs(value_c) > 1e-7)
564
 
                              fprintf(outfile, "%d %d %d %d %d %d %20.15f\n", I,J,K,A,B,C,value_c);
565
 
                            */
566
 
 
567
 
                            ET_AAB += (value_d + value_c) * value_c / denom;
568
 
 
569
 
                            /* } */ /* A >= B */
570
 
 
571
 
                        } /* c */
572
 
                      } /* b */
573
 
                    } /* a */
574
 
 
575
 
                  } /* Gb */
576
 
                } /* Ga */
577
 
              
578
 
              } /* I >= J */
579
 
 
580
 
            } /* k */
581
 
          } /* j */
582
 
        } /* i */
583
 
 
584
 
      } /* Gk */
585
 
    } /* Gj */
586
 
  } /* Gi */
587
 
 
588
 
  ET_AAB /= 2;
589
 
 
590
 
  /*  fprintf(outfile, "cnt = %d\n", cnt); */
591
 
  /* fprintf(outfile, "ET_AAB = %20.14f\n", ET_AAB); */
592
 
 
593
 
  for(h=0; h < nirreps; h++) {
594
 
    dpd_buf4_mat_irrep_close(&T2AA, h);
595
 
    dpd_buf4_mat_irrep_close(&T2AB, h);
596
 
    dpd_buf4_mat_irrep_close(&FAAints, h);
597
 
    dpd_buf4_mat_irrep_close(&FABints, h);
598
 
    dpd_buf4_mat_irrep_close(&FBAints, h);
599
 
    dpd_buf4_mat_irrep_close(&EAAints, h);
600
 
    dpd_buf4_mat_irrep_close(&EABints, h);
601
 
    dpd_buf4_mat_irrep_close(&EBAints, h);
602
 
    dpd_buf4_mat_irrep_close(&DAAints, h);
603
 
    dpd_buf4_mat_irrep_close(&DABints, h);
604
 
  }
605
 
 
606
 
  dpd_buf4_close(&T2AA);
607
 
  dpd_buf4_close(&T2AB);
608
 
  dpd_buf4_close(&FAAints);
609
 
  dpd_buf4_close(&FABints);
610
 
  dpd_buf4_close(&FBAints);
611
 
  dpd_buf4_close(&EAAints);
612
 
  dpd_buf4_close(&EABints);
613
 
  dpd_buf4_close(&EBAints);
614
 
  dpd_buf4_close(&DAAints);
615
 
  dpd_buf4_close(&DABints);
616
 
 
617
 
  dpd_file2_mat_close(&T1A);
618
 
  dpd_file2_close(&T1A);
619
 
  dpd_file2_mat_close(&T1B);
620
 
  dpd_file2_close(&T1B);
621
 
 
622
 
  dpd_file2_mat_close(&fIJ);
623
 
  dpd_file2_mat_close(&fij);
624
 
  dpd_file2_mat_close(&fAB);
625
 
  dpd_file2_mat_close(&fab);
626
 
  dpd_file2_close(&fIJ);
627
 
  dpd_file2_close(&fij);
628
 
  dpd_file2_close(&fAB);
629
 
  dpd_file2_close(&fab);
630
 
 
631
 
  return ET_AAB;
632
 
}