← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/ccsort/fock.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/ccsort/fock.c
1
 
#include <stdio.h>
2
 
#include <stdlib.h>
3
 
#include <libpsio/psio.h>
4
 
#include <libciomr/libciomr.h>
5
 
#include <libdpd/dpd.h>
6
 
#include "Params.h"
7
 
#include "MOInfo.h"
8
 
#define EXTERN
9
 
#include "globals.h"
10
 
 
11
 
/* 
12
 
** fock(): Build the alpha and beta Fock matrices from the
13
 
** one-electron integrals/frozen-core operator and active two-electron
14
 
** integrals on disk.
15
 
**
16
 
** TDC, 1996
17
 
** Modified to include UHF references, TDC, June 2001.
18
 
**
19
 
** Notes:
20
 
**
21
 
** (1) This routine isn't absolutely necessary for RHF and UHF
22
 
** references, as one could simply read the Fock eigenvalues from
23
 
** PSIF_CHKPT and be happy with that. However, this code is useful as a
24
 
** partial check of the integral transformation and sorting routines.
25
 
**
26
 
** (2) An alternative but currently unused algorithm may be found in 
27
 
** fock_build.c.
28
 
*/
29
 
 
30
 
void fock_rhf(void);
31
 
void fock_uhf(void);
32
 
 
33
 
void fock(void)
34
 
{
35
 
  if(params.ref == 2) fock_uhf();
36
 
  else fock_rhf();
37
 
}
38
 
 
39
 
void fock_uhf(void)
40
 
{
41
 
  int h, nirreps;
42
 
  int i, j, I, J, Gi, Gj, IM, JM, MI, MJ;
43
 
  int a, b, A, B, Ga, Gb, AM, BM, MA, MB;
44
 
  int m, M, Gm;
45
 
  int *aoccpi, *boccpi, *aocc_off, *bocc_off;
46
 
  int *avirtpi, *bvirtpi, *avir_off, *bvir_off;
47
 
  dpdfile2 fIJ, fij, fAB, fab, fIA, fia;
48
 
  dpdfile2 hIJ, hij, hAB, hab, hIA, hia;
49
 
  dpdbuf4 A_AA, A_BB, A_AB, C_AA, C_BB, C_AB, E_AA, E_BB, E_AB;
50
 
 
51
 
  nirreps = moinfo.nirreps;
52
 
  aoccpi = moinfo.aoccpi;
53
 
  boccpi = moinfo.boccpi;
54
 
  avirtpi = moinfo.avirtpi;
55
 
  bvirtpi = moinfo.bvirtpi;
56
 
  aocc_off = moinfo.aocc_off;
57
 
  bocc_off = moinfo.bocc_off;
58
 
  avir_off = moinfo.avir_off;
59
 
  bvir_off = moinfo.bvir_off;
60
 
 
61
 
  dpd_file2_init(&hIJ, CC_OEI, 0, 0, 0, "h(I,J)");
62
 
  dpd_file2_init(&hij, CC_OEI, 0, 2, 2, "h(i,j)");
63
 
  dpd_file2_init(&hAB, CC_OEI, 0, 1, 1, "h(A,B)");
64
 
  dpd_file2_init(&hab, CC_OEI, 0, 3, 3, "h(a,b)");
65
 
  dpd_file2_init(&hIA, CC_OEI, 0, 0, 1, "h(I,A)");
66
 
  dpd_file2_init(&hia, CC_OEI, 0, 2, 3, "h(i,a)");
67
 
 
68
 
  dpd_file2_mat_init(&hIJ);
69
 
  dpd_file2_mat_init(&hij);
70
 
  dpd_file2_mat_init(&hAB);
71
 
  dpd_file2_mat_init(&hab);
72
 
  dpd_file2_mat_init(&hIA);
73
 
  dpd_file2_mat_init(&hia);
74
 
 
75
 
  dpd_file2_mat_rd(&hIJ);
76
 
  dpd_file2_mat_rd(&hij);
77
 
  dpd_file2_mat_rd(&hAB);
78
 
  dpd_file2_mat_rd(&hab);
79
 
  dpd_file2_mat_rd(&hIA);
80
 
  dpd_file2_mat_rd(&hia);
81
 
 
82
 
  dpd_file2_init(&fIJ, CC_OEI, 0, 0, 0, "fIJ");
83
 
  dpd_file2_init(&fij, CC_OEI, 0, 2, 2, "fij");
84
 
 
85
 
  dpd_file2_mat_init(&fIJ);
86
 
  dpd_file2_mat_init(&fij);
87
 
 
88
 
  for(h=0; h < nirreps; h++) {
89
 
 
90
 
    for(i=0; i < aoccpi[h]; i++)
91
 
      for(j=0; j < aoccpi[h]; j++)
92
 
        fIJ.matrix[h][i][j] = hIJ.matrix[h][i][j];
93
 
 
94
 
    for(i=0; i < boccpi[h]; i++)
95
 
      for(j=0; j < boccpi[h]; j++)
96
 
        fij.matrix[h][i][j] = hij.matrix[h][i][j];
97
 
  }
98
 
 
99
 
  dpd_buf4_init(&A_AA, CC_AINTS, 0, 0, 0, 0, 0, 1, "A <IJ|KL>");
100
 
  dpd_buf4_init(&A_AB, CC_AINTS, 0, 22, 22, 22, 22, 0, "A <Ij|Kl>");
101
 
  for(h=0; h < nirreps; h++) {
102
 
    dpd_buf4_mat_irrep_init(&A_AA, h);
103
 
    dpd_buf4_mat_irrep_rd(&A_AA, h);
104
 
    dpd_buf4_mat_irrep_init(&A_AB, h);
105
 
    dpd_buf4_mat_irrep_rd(&A_AB, h);
106
 
    for(Gi=0; Gi < nirreps; Gi++) {
107
 
      Gj = Gi; Gm = Gi^h;
108
 
      for(i=0; i < aoccpi[Gi]; i++) {
109
 
        I = aocc_off[Gi] + i;
110
 
        for(j=0; j < aoccpi[Gj]; j++) {
111
 
          J = aocc_off[Gj] + j;
112
 
          for(m=0; m < aoccpi[Gm]; m++) {
113
 
            M = aocc_off[Gm] + m;
114
 
            IM = A_AA.params->rowidx[I][M];
115
 
            JM = A_AA.params->colidx[J][M];
116
 
            fIJ.matrix[Gi][i][j] += A_AA.matrix[h][IM][JM];
117
 
          }
118
 
          for(m=0; m < boccpi[Gm]; m++) {
119
 
            M = bocc_off[Gm] + m;
120
 
            IM = A_AB.params->rowidx[I][M];
121
 
            JM = A_AB.params->colidx[J][M];
122
 
            fIJ.matrix[Gi][i][j] += A_AB.matrix[h][IM][JM];
123
 
          }
124
 
        }
125
 
      }
126
 
    }
127
 
    dpd_buf4_mat_irrep_close(&A_AA, h);
128
 
    dpd_buf4_mat_irrep_close(&A_AB, h);
129
 
  }
130
 
  dpd_buf4_close(&A_AA);
131
 
  dpd_buf4_close(&A_AB);
132
 
 
133
 
  dpd_buf4_init(&A_BB, CC_AINTS, 0, 10, 10, 10, 10, 1, "A <ij|kl>");
134
 
  dpd_buf4_init(&A_AB, CC_AINTS, 0, 22, 22, 22, 22, 0, "A <Ij|Kl>");
135
 
  for(h=0; h < nirreps; h++) {
136
 
    dpd_buf4_mat_irrep_init(&A_BB, h);
137
 
    dpd_buf4_mat_irrep_rd(&A_BB, h);
138
 
    dpd_buf4_mat_irrep_init(&A_AB, h);
139
 
    dpd_buf4_mat_irrep_rd(&A_AB, h);
140
 
    for(Gi=0; Gi < nirreps; Gi++) {
141
 
      Gj = Gi; Gm = Gi^h;
142
 
      for(i=0; i < boccpi[Gi]; i++) {
143
 
        I = bocc_off[Gi] + i;
144
 
        for(j=0; j < boccpi[Gj]; j++) {
145
 
          J = bocc_off[Gj] + j;
146
 
          for(m=0; m < boccpi[Gm]; m++) {
147
 
            M = bocc_off[Gm] + m;
148
 
            IM = A_BB.params->rowidx[I][M];
149
 
            JM = A_BB.params->colidx[J][M];
150
 
            fij.matrix[Gi][i][j] += A_BB.matrix[h][IM][JM];
151
 
          }
152
 
          for(m=0; m < aoccpi[Gm]; m++) {
153
 
            M = aocc_off[Gm] + m;
154
 
            MI = A_AB.params->rowidx[M][I];
155
 
            MJ = A_AB.params->colidx[M][J];
156
 
            fij.matrix[Gi][i][j] += A_AB.matrix[h][MI][MJ];
157
 
          }
158
 
        }
159
 
      }
160
 
    }
161
 
    dpd_buf4_mat_irrep_close(&A_BB, h);
162
 
    dpd_buf4_mat_irrep_close(&A_AB, h);
163
 
  }
164
 
  dpd_buf4_close(&A_BB);
165
 
  dpd_buf4_close(&A_AB);
166
 
 
167
 
  dpd_file2_mat_wrt(&fIJ);
168
 
  dpd_file2_mat_wrt(&fij);
169
 
 
170
 
  dpd_file2_close(&fIJ);
171
 
  dpd_file2_close(&fij);
172
 
 
173
 
  dpd_file2_init(&fAB, CC_OEI, 0, 1, 1, "fAB");
174
 
  dpd_file2_init(&fab, CC_OEI, 0, 3, 3, "fab");
175
 
 
176
 
  dpd_file2_mat_init(&fAB);
177
 
  dpd_file2_mat_init(&fab);
178
 
 
179
 
  for(h=0; h < nirreps; h++) {
180
 
 
181
 
    for(a=0; a < avirtpi[h]; a++)
182
 
      for(b=0; b < avirtpi[h]; b++)
183
 
        fAB.matrix[h][a][b] = hAB.matrix[h][a][b];
184
 
 
185
 
    for(a=0; a < bvirtpi[h]; a++)
186
 
      for(b=0; b < bvirtpi[h]; b++)
187
 
        fab.matrix[h][a][b] = hab.matrix[h][a][b];
188
 
  }
189
 
 
190
 
  dpd_buf4_init(&C_AA, CC_CINTS, 0, 20, 20, 20, 20, 0, "C <IA||JB>");
191
 
  dpd_buf4_init(&C_AB, CC_CINTS, 0, 26, 26, 26, 26, 0, "C <Ai|Bj>");
192
 
  for(h=0; h < nirreps; h++) {
193
 
    dpd_buf4_mat_irrep_init(&C_AA, h);
194
 
    dpd_buf4_mat_irrep_rd(&C_AA, h);
195
 
    dpd_buf4_mat_irrep_init(&C_AB, h);
196
 
    dpd_buf4_mat_irrep_rd(&C_AB, h);
197
 
    for(Ga=0; Ga < nirreps; Ga++) {
198
 
      Gb = Ga; Gm = Ga^h;
199
 
      for(a=0; a < avirtpi[Ga]; a++) {
200
 
        A = avir_off[Ga] + a;
201
 
        for(b=0; b < avirtpi[Gb]; b++) {
202
 
          B = avir_off[Gb] + b;
203
 
          for(m=0; m < aoccpi[Gm]; m++) {
204
 
            M = aocc_off[Gm] + m;
205
 
            MA = C_AA.params->rowidx[M][A];
206
 
            MB = C_AA.params->colidx[M][B];
207
 
            fAB.matrix[Ga][a][b] += C_AA.matrix[h][MA][MB];
208
 
          }
209
 
          for(m=0; m < boccpi[Gm]; m++) {
210
 
            M = bocc_off[Gm] + m;
211
 
            AM = C_AB.params->rowidx[A][M];
212
 
            BM = C_AB.params->colidx[B][M];
213
 
            fAB.matrix[Ga][a][b] += C_AB.matrix[h][AM][BM];
214
 
          }
215
 
        }
216
 
      }
217
 
    }
218
 
    dpd_buf4_mat_irrep_close(&C_AA, h);
219
 
    dpd_buf4_mat_irrep_close(&C_AB, h);
220
 
  }
221
 
  dpd_buf4_close(&C_AA);
222
 
  dpd_buf4_close(&C_AB);
223
 
 
224
 
  dpd_buf4_init(&C_BB, CC_CINTS, 0, 30, 30, 30, 30, 0, "C <ia||jb>");
225
 
  dpd_buf4_init(&C_AB, CC_CINTS, 0, 24, 24, 24, 24, 0, "C <Ia|Jb>");
226
 
  for(h=0; h < nirreps; h++) {
227
 
    dpd_buf4_mat_irrep_init(&C_BB, h);
228
 
    dpd_buf4_mat_irrep_rd(&C_BB, h);
229
 
    dpd_buf4_mat_irrep_init(&C_AB, h);
230
 
    dpd_buf4_mat_irrep_rd(&C_AB, h);
231
 
    for(Ga=0; Ga < nirreps; Ga++) {
232
 
      Gb = Ga; Gm = Ga^h;
233
 
      for(a=0; a < bvirtpi[Ga]; a++) {
234
 
        A = bvir_off[Ga] + a;
235
 
        for(b=0; b < bvirtpi[Gb]; b++) {
236
 
          B = bvir_off[Gb] + b;
237
 
          for(m=0; m < boccpi[Gm]; m++) {
238
 
            M = bocc_off[Gm] + m;
239
 
            MA = C_BB.params->rowidx[M][A];
240
 
            MB = C_BB.params->colidx[M][B];
241
 
            fab.matrix[Ga][a][b] += C_BB.matrix[h][MA][MB];
242
 
          }
243
 
          for(m=0; m < aoccpi[Gm]; m++) {
244
 
            M = aocc_off[Gm] + m;
245
 
            MA = C_AB.params->rowidx[M][A];
246
 
            MB = C_AB.params->colidx[M][B];
247
 
            fab.matrix[Ga][a][b] += C_AB.matrix[h][MA][MB];
248
 
          }
249
 
        }
250
 
      }
251
 
    }
252
 
    dpd_buf4_mat_irrep_close(&C_BB, h);
253
 
    dpd_buf4_mat_irrep_close(&C_AB, h);
254
 
  }
255
 
  dpd_buf4_close(&C_BB);
256
 
  dpd_buf4_close(&C_AB);
257
 
 
258
 
  dpd_file2_mat_wrt(&fAB);
259
 
  dpd_file2_mat_wrt(&fab);
260
 
 
261
 
  dpd_file2_close(&fAB);
262
 
  dpd_file2_close(&fab);
263
 
 
264
 
  /* Prepare the alpha and beta occ-vir Fock matrix files */
265
 
  dpd_file2_init(&fIA, CC_OEI, 0, 0, 1, "fIA");
266
 
  dpd_file2_init(&fia, CC_OEI, 0, 2, 3, "fia");
267
 
  dpd_file2_mat_init(&fIA);
268
 
  dpd_file2_mat_init(&fia);
269
 
 
270
 
  /* One-electron (frozen-core) contributions */
271
 
  for(h=0; h < nirreps; h++) {
272
 
 
273
 
    for(i=0; i < aoccpi[h]; i++) 
274
 
      for(a=0; a < avirtpi[h]; a++) 
275
 
        fIA.matrix[h][i][a] = hIA.matrix[h][i][a];   
276
 
 
277
 
    for(i=0; i < boccpi[h]; i++) 
278
 
      for(a=0; a < bvirtpi[h]; a++) 
279
 
        fia.matrix[h][i][a] = hia.matrix[h][i][a];   
280
 
  }
281
 
 
282
 
  /* Two-electron contributions */
283
 
 
284
 
  /* Prepare the E integral buffers */
285
 
  dpd_buf4_init(&E_AA, CC_EINTS, 0, 21, 0, 21, 0, 1, "E <AI|JK>");
286
 
  dpd_buf4_init(&E_AB, CC_EINTS, 0, 26, 22, 26, 22, 0, "E <Ai|Jk>");
287
 
 
288
 
  for(h=0; h < nirreps; h++) {
289
 
 
290
 
    dpd_buf4_mat_irrep_init(&E_AA, h);
291
 
    dpd_buf4_mat_irrep_init(&E_AB, h);
292
 
    dpd_buf4_mat_irrep_rd(&E_AA, h);
293
 
    dpd_buf4_mat_irrep_rd(&E_AB, h);
294
 
 
295
 
    /* Loop over irreps of the target */
296
 
    for(Gi=0; Gi < nirreps; Gi++) {
297
 
      Ga = Gi; Gm = h^Gi;
298
 
 
299
 
      /* Loop over orbitals of the target */
300
 
      for(i=0; i < aoccpi[Gi]; i++) {
301
 
        I = aocc_off[Gi] + i;
302
 
        for(a=0; a < avirtpi[Ga]; a++) {
303
 
          A = avir_off[Ga] + a;
304
 
 
305
 
          for(m=0; m < aoccpi[Gm]; m++) {
306
 
            M = aocc_off[Gm] + m;
307
 
 
308
 
            AM = E_AA.params->rowidx[A][M];
309
 
            IM = E_AA.params->colidx[I][M];
310
 
 
311
 
            fIA.matrix[Gi][i][a] += E_AA.matrix[h][AM][IM];
312
 
 
313
 
          }
314
 
        }
315
 
      }
316
 
 
317
 
      /* Loop over orbitals of the target */
318
 
      for(i=0; i < aoccpi[Gi]; i++) {
319
 
        I = aocc_off[Gi] + i;
320
 
        for(a=0; a < avirtpi[Ga]; a++) {
321
 
          A = avir_off[Ga] + a;
322
 
 
323
 
          for(m=0; m < boccpi[Gm]; m++) {
324
 
            M = bocc_off[Gm] + m;
325
 
 
326
 
            AM = E_AB.params->rowidx[A][M];
327
 
            IM = E_AB.params->colidx[I][M];
328
 
 
329
 
            fIA.matrix[Gi][i][a] += E_AB.matrix[h][AM][IM];
330
 
 
331
 
          }
332
 
        }
333
 
      }
334
 
 
335
 
    }
336
 
 
337
 
    dpd_buf4_mat_irrep_close(&E_AA, h);
338
 
    dpd_buf4_mat_irrep_close(&E_AB, h);
339
 
  }
340
 
 
341
 
  dpd_buf4_close(&E_AA);
342
 
  dpd_buf4_close(&E_AB);
343
 
 
344
 
  dpd_buf4_init(&E_BB, CC_EINTS, 0, 31, 10, 31, 10, 1, "E <ai|jk>");
345
 
  dpd_buf4_init(&E_AB, CC_EINTS, 0, 22, 24, 22, 24, 0, "E <Ij|Ka>");
346
 
 
347
 
  for(h=0; h < nirreps; h++) {
348
 
 
349
 
    dpd_buf4_mat_irrep_init(&E_BB, h);
350
 
    dpd_buf4_mat_irrep_init(&E_AB, h);
351
 
    dpd_buf4_mat_irrep_rd(&E_BB, h);
352
 
    dpd_buf4_mat_irrep_rd(&E_AB, h);
353
 
 
354
 
    /* Loop over irreps of the target */
355
 
    for(Gi=0; Gi < nirreps; Gi++) {
356
 
      Ga = Gi; Gm = h^Gi;
357
 
 
358
 
      /* Loop over orbitals of the target */
359
 
      for(i=0; i < boccpi[Gi]; i++) {
360
 
        I = bocc_off[Gi] + i;
361
 
        for(a=0; a < bvirtpi[Ga]; a++) {
362
 
          A = bvir_off[Ga] + a;
363
 
 
364
 
          for(m=0; m < boccpi[Gm]; m++) {
365
 
            M = bocc_off[Gm] + m;
366
 
 
367
 
            AM = E_BB.params->rowidx[A][M];
368
 
            IM = E_BB.params->colidx[I][M];
369
 
 
370
 
            fia.matrix[Gi][i][a] += E_BB.matrix[h][AM][IM];
371
 
 
372
 
          }
373
 
        }
374
 
      }
375
 
   
376
 
      /* Loop over orbitals of the target */
377
 
      for(i=0; i < boccpi[Gi]; i++) {
378
 
        I = bocc_off[Gi] + i;
379
 
        for(a=0; a < bvirtpi[Ga]; a++) {
380
 
          A = bvir_off[Ga] + a;
381
 
 
382
 
          for(m=0; m < aoccpi[Gm]; m++) {
383
 
            M = aocc_off[Gm] + m;
384
 
 
385
 
            MI = E_AB.params->rowidx[M][I];
386
 
            MA = E_AB.params->colidx[M][A];
387
 
 
388
 
            fia.matrix[Gi][i][a] += E_AB.matrix[h][MI][MA];
389
 
 
390
 
          }
391
 
        }
392
 
      }
393
 
    }
394
 
 
395
 
    dpd_buf4_mat_irrep_close(&E_BB, h);
396
 
    dpd_buf4_mat_irrep_close(&E_AB, h);
397
 
 
398
 
  }
399
 
 
400
 
  /* Close the E integral buffers */
401
 
  dpd_buf4_close(&E_BB);
402
 
  dpd_buf4_close(&E_AB);
403
 
 
404
 
  /* Close the alpha and beta occ-vir Fock matrix files */
405
 
  dpd_file2_mat_wrt(&fIA);
406
 
  dpd_file2_mat_wrt(&fia);
407
 
  dpd_file2_mat_close(&fIA);
408
 
  dpd_file2_mat_close(&fia);
409
 
  dpd_file2_close(&fIA);
410
 
  dpd_file2_close(&fia);
411
 
 
412
 
  dpd_file2_mat_close(&hIJ);
413
 
  dpd_file2_mat_close(&hij);
414
 
  dpd_file2_mat_close(&hAB);
415
 
  dpd_file2_mat_close(&hab);
416
 
  dpd_file2_mat_close(&hIA);
417
 
  dpd_file2_mat_close(&hia);
418
 
 
419
 
  dpd_file2_close(&hIJ);
420
 
  dpd_file2_close(&hij);
421
 
  dpd_file2_close(&hAB);
422
 
  dpd_file2_close(&hab);
423
 
  dpd_file2_close(&hIA);
424
 
  dpd_file2_close(&hia);
425
 
}
426
 
 
427
 
void fock_rhf(void)
428
 
{
429
 
  int h, Gi, Gj, Ga, Gb, Gm;
430
 
  int i,j,a,b,m;
431
 
  int I, J, A, B, M;
432
 
  int IM, JM, MA, MB, AM;
433
 
  int nirreps;
434
 
  int *occpi, *virtpi;
435
 
  int *occ_off, *vir_off;
436
 
  int *occ_sym, *vir_sym;
437
 
  int *openpi;
438
 
  dpdfile2 fIJ, fij, fAB, fab, fIA, fia, Hoo, Hvv, Hov;
439
 
  dpdbuf4 AInts_anti, AInts, CInts, CInts_anti, EInts_anti, EInts;
440
 
 
441
 
  nirreps = moinfo.nirreps;
442
 
  occpi = moinfo.occpi; virtpi = moinfo.virtpi;
443
 
  occ_off = moinfo.occ_off; vir_off = moinfo.vir_off;
444
 
  occ_sym = moinfo.occ_sym; vir_sym = moinfo.vir_sym;
445
 
  openpi = moinfo.openpi;
446
 
 
447
 
  dpd_file2_init(&Hoo, CC_OEI, 0, 0, 0, "h(i,j)");
448
 
  dpd_file2_init(&Hvv, CC_OEI, 0, 1, 1, "h(a,b)");
449
 
  dpd_file2_init(&Hov, CC_OEI, 0, 0, 1, "h(i,a)");
450
 
  dpd_file2_mat_init(&Hoo);
451
 
  dpd_file2_mat_init(&Hvv);
452
 
  dpd_file2_mat_init(&Hov);
453
 
  dpd_file2_mat_rd(&Hoo);
454
 
  dpd_file2_mat_rd(&Hvv);
455
 
  dpd_file2_mat_rd(&Hov);
456
 
  
457
 
  /* Prepare the alpha and beta occ-occ Fock matrix files */
458
 
  dpd_file2_init(&fIJ, CC_OEI, 0, 0, 0, "fIJ");
459
 
  dpd_file2_init(&fij, CC_OEI, 0, 0, 0, "fij");
460
 
  dpd_file2_mat_init(&fIJ);
461
 
  dpd_file2_mat_init(&fij);
462
 
 
463
 
  /* One-electron (frozen-core) contributions */
464
 
  for(h=0; h < nirreps; h++) {
465
 
 
466
 
      for(i=0; i < occpi[h]; i++) 
467
 
          for(j=0; j < occpi[h]; j++) 
468
 
              fIJ.matrix[h][i][j] = Hoo.matrix[h][i][j];   
469
 
 
470
 
      for(i=0; i < (occpi[h]-openpi[h]); i++) 
471
 
          for(j=0; j < (occpi[h]-openpi[h]); j++) 
472
 
              fij.matrix[h][i][j] = Hoo.matrix[h][i][j];   
473
 
    }
474
 
 
475
 
  /* Two-electron contributions */
476
 
 
477
 
  /* Prepare the A integral buffers */
478
 
  dpd_buf4_init(&AInts_anti, CC_AINTS, 0, 0, 0, 0, 0, 1, "A <ij|kl>");
479
 
  dpd_buf4_init(&AInts, CC_AINTS, 0, 0, 0, 0, 0, 0, "A <ij|kl>");
480
 
 
481
 
  for(h=0; h < nirreps; h++) {
482
 
 
483
 
      dpd_buf4_mat_irrep_init(&AInts_anti, h);
484
 
      dpd_buf4_mat_irrep_rd(&AInts_anti, h);
485
 
 
486
 
      /* Loop over irreps of the target */
487
 
      for(Gi=0; Gi < nirreps; Gi++) {
488
 
          Gj=Gi; Gm=h^Gi;
489
 
 
490
 
          /* Loop over the orbitals of the target */
491
 
          for(i=0; i < occpi[Gi]; i++) {
492
 
              I = occ_off[Gi] + i;
493
 
              for(j=0; j < occpi[Gj]; j++) {
494
 
                  J = occ_off[Gj] + j;
495
 
                  for(m=0; m < occpi[Gm]; m++) {
496
 
                      M = occ_off[Gm] + m;
497
 
 
498
 
                      IM = AInts_anti.params->rowidx[I][M];
499
 
                      JM = AInts_anti.params->colidx[J][M];
500
 
 
501
 
                      fIJ.matrix[Gi][i][j] += AInts_anti.matrix[h][IM][JM];
502
 
 
503
 
                    }
504
 
                }
505
 
            }
506
 
 
507
 
          /* Loop over the orbitals of the target */
508
 
          for(i=0; i < (occpi[Gi] - openpi[Gi]); i++) {
509
 
              I = occ_off[Gi] + i;
510
 
              for(j=0; j < (occpi[Gj] - openpi[Gj]); j++) {
511
 
                  J = occ_off[Gj] + j;
512
 
                  for(m=0; m < (occpi[Gm] - openpi[Gm]); m++) {
513
 
                      M = occ_off[Gm] + m;
514
 
 
515
 
                      IM = AInts_anti.params->rowidx[I][M];
516
 
                      JM = AInts_anti.params->colidx[J][M];
517
 
 
518
 
                      fij.matrix[Gi][i][j] += AInts_anti.matrix[h][IM][JM];
519
 
 
520
 
                    }
521
 
                }
522
 
            }
523
 
 
524
 
        }
525
 
      
526
 
      dpd_buf4_mat_irrep_close(&AInts_anti, h);
527
 
 
528
 
      dpd_buf4_mat_irrep_init(&AInts, h);
529
 
      dpd_buf4_mat_irrep_rd(&AInts, h);
530
 
 
531
 
      /* Loop over irreps of the target */
532
 
      for(Gi=0; Gi < nirreps; Gi++) {
533
 
          Gj=Gi; Gm=h^Gi;
534
 
 
535
 
          /* Loop over the orbitals of the target */
536
 
          for(i=0; i < occpi[Gi]; i++) {
537
 
              I = occ_off[Gi] + i;
538
 
              for(j=0; j < occpi[Gj]; j++) {
539
 
                  J = occ_off[Gj] + j;
540
 
                  for(m=0; m < (occpi[Gm] - openpi[Gm]); m++) {
541
 
                      M = occ_off[Gm] + m;
542
 
 
543
 
                      IM = AInts.params->rowidx[I][M];
544
 
                      JM = AInts.params->colidx[J][M];
545
 
 
546
 
                      fIJ.matrix[Gi][i][j] += AInts.matrix[h][IM][JM];
547
 
 
548
 
                    }
549
 
                }
550
 
            }
551
 
 
552
 
          /* Loop over the orbitals of the target */
553
 
          for(i=0; i < (occpi[Gi] - openpi[Gi]); i++) {
554
 
              I = occ_off[Gi] + i;
555
 
              for(j=0; j < (occpi[Gj] - openpi[Gj]); j++) {
556
 
                  J = occ_off[Gj] + j;
557
 
                  for(m=0; m < occpi[Gm]; m++) {
558
 
                      M = occ_off[Gm] + m;
559
 
 
560
 
                      IM = AInts.params->rowidx[I][M];
561
 
                      JM = AInts.params->colidx[J][M];
562
 
 
563
 
                      fij.matrix[Gi][i][j] += AInts.matrix[h][IM][JM];
564
 
 
565
 
                    }
566
 
                }
567
 
            }
568
 
 
569
 
        }
570
 
      
571
 
      dpd_buf4_mat_irrep_close(&AInts, h);
572
 
 
573
 
    }
574
 
 
575
 
  /* Close the A Integral buffers */
576
 
  dpd_buf4_close(&AInts_anti);
577
 
  dpd_buf4_close(&AInts);
578
 
  
579
 
  /* Close the alpha and beta occ-occ Fock matrix files */
580
 
  dpd_file2_mat_wrt(&fIJ);
581
 
  dpd_file2_mat_wrt(&fij);
582
 
  dpd_file2_mat_close(&fIJ);
583
 
  dpd_file2_mat_close(&fij);
584
 
  dpd_file2_close(&fIJ);
585
 
  dpd_file2_close(&fij);
586
 
  
587
 
  /* Prepare the alpha and beta vir-vir Fock matrix files */
588
 
  dpd_file2_init(&fAB, CC_OEI, 0, 1, 1, "fAB");
589
 
  dpd_file2_init(&fab, CC_OEI, 0, 1, 1, "fab");
590
 
  dpd_file2_mat_init(&fAB);
591
 
  dpd_file2_mat_init(&fab);
592
 
 
593
 
  /* One-electron (frozen-core) contributions */
594
 
  for(h=0; h < nirreps; h++) {
595
 
 
596
 
      for(a=0; a < (virtpi[h] - openpi[h]); a++) 
597
 
          for(b=0; b < (virtpi[h] - openpi[h]); b++) 
598
 
              fAB.matrix[h][a][b] = Hvv.matrix[h][a][b];   
599
 
 
600
 
      for(a=0; a < virtpi[h]; a++) 
601
 
          for(b=0; b < virtpi[h]; b++) 
602
 
              fab.matrix[h][a][b] = Hvv.matrix[h][a][b];   
603
 
    }
604
 
 
605
 
  /* Two-electron contributions */
606
 
 
607
 
  /* Prepare the C integral buffers */
608
 
  dpd_buf4_init(&CInts_anti, CC_CINTS, 0, 10, 10, 10, 10, 0, "C <ia||jb>");
609
 
  dpd_buf4_init(&CInts, CC_CINTS, 0, 10, 10, 10, 10, 0, "C <ia|jb>");
610
 
 
611
 
  for(h=0; h < nirreps; h++) {
612
 
 
613
 
      dpd_buf4_mat_irrep_init(&CInts_anti, h);
614
 
      dpd_buf4_mat_irrep_rd(&CInts_anti, h);
615
 
 
616
 
      /* Loop over irreps of the target */
617
 
      for(Ga=0; Ga < nirreps; Ga++) {
618
 
          Gb = Ga; Gm = h^Ga;
619
 
 
620
 
          /* Loop over orbitals of the target */
621
 
          for(a=0; a < (virtpi[Ga] - openpi[Ga]); a++) {
622
 
              A = vir_off[Ga] + a;
623
 
              for(b=0; b < (virtpi[Gb] - openpi[Gb]); b++) {
624
 
                  B = vir_off[Gb] + b;
625
 
 
626
 
                  for(m=0; m < occpi[Gm]; m++) {
627
 
                      M = occ_off[Gm] + m;
628
 
 
629
 
                      MA = CInts_anti.params->rowidx[M][A];
630
 
                      MB = CInts_anti.params->colidx[M][B];
631
 
 
632
 
                      fAB.matrix[Ga][a][b] += CInts_anti.matrix[h][MA][MB];
633
 
 
634
 
                    }
635
 
                }
636
 
            }
637
 
 
638
 
          /* Loop over orbitals of the target */
639
 
          for(a=0; a < virtpi[Ga]; a++) {
640
 
              A = vir_off[Ga] + a;
641
 
              for(b=0; b < virtpi[Gb]; b++) {
642
 
                  B = vir_off[Gb] + b;
643
 
 
644
 
                  for(m=0; m < (occpi[Gm] - openpi[Gm]); m++) {
645
 
                      M = occ_off[Gm] + m;
646
 
 
647
 
                      MA = CInts_anti.params->rowidx[M][A];
648
 
                      MB = CInts_anti.params->colidx[M][B];
649
 
 
650
 
                      fab.matrix[Ga][a][b] += CInts_anti.matrix[h][MA][MB];
651
 
                    }
652
 
                }
653
 
            }
654
 
        }
655
 
 
656
 
      dpd_buf4_mat_irrep_close(&CInts_anti, h);
657
 
 
658
 
      dpd_buf4_mat_irrep_init(&CInts, h);
659
 
      dpd_buf4_mat_irrep_rd(&CInts, h);
660
 
 
661
 
      /* Loop over irreps of the target */
662
 
      for(Ga=0; Ga < nirreps; Ga++) {
663
 
          Gb = Ga; Gm = h^Ga;
664
 
 
665
 
          /* Loop over orbitals of the target */
666
 
          for(a=0; a < (virtpi[Ga] - openpi[Ga]); a++) {
667
 
              A = vir_off[Ga] + a;
668
 
              for(b=0; b < (virtpi[Gb] - openpi[Gb]); b++) {
669
 
                  B = vir_off[Gb] + b;
670
 
 
671
 
                  for(m=0; m < (occpi[Gm] - openpi[Gm]); m++) {
672
 
                      M = occ_off[Gm] + m;
673
 
 
674
 
                      MA = CInts.params->rowidx[M][A];
675
 
                      MB = CInts.params->colidx[M][B];
676
 
 
677
 
                      fAB.matrix[Ga][a][b] += CInts.matrix[h][MA][MB];
678
 
 
679
 
                    }
680
 
                }
681
 
            }
682
 
 
683
 
          /* Loop over orbitals of the target */
684
 
          for(a=0; a < virtpi[Ga]; a++) {
685
 
              A = vir_off[Ga] + a;
686
 
              for(b=0; b < virtpi[Gb]; b++) {
687
 
                  B = vir_off[Gb] + b;
688
 
 
689
 
                  for(m=0; m < occpi[Gm]; m++) {
690
 
                      M = occ_off[Gm] + m;
691
 
 
692
 
                      MA = CInts.params->rowidx[M][A];
693
 
                      MB = CInts.params->colidx[M][B];
694
 
 
695
 
                      fab.matrix[Ga][a][b] += CInts.matrix[h][MA][MB];
696
 
                    }
697
 
                }
698
 
            }
699
 
        }
700
 
 
701
 
      dpd_buf4_mat_irrep_close(&CInts, h);
702
 
    }
703
 
 
704
 
  /* Close the C integral buffers */
705
 
  dpd_buf4_close(&CInts_anti);
706
 
  dpd_buf4_close(&CInts);
707
 
 
708
 
  /* Close the alpha and beta vir-vir Fock matrix files */
709
 
  dpd_file2_mat_wrt(&fAB);
710
 
  dpd_file2_mat_wrt(&fab);
711
 
  dpd_file2_mat_close(&fAB);
712
 
  dpd_file2_mat_close(&fab);
713
 
  dpd_file2_close(&fAB);
714
 
  dpd_file2_close(&fab);
715
 
  
716
 
  /* Prepare the alpha and beta occ-vir Fock matrix files */
717
 
  dpd_file2_init(&fIA, CC_OEI, 0, 0, 1, "fIA");
718
 
  dpd_file2_init(&fia, CC_OEI, 0, 0, 1, "fia");
719
 
  dpd_file2_mat_init(&fIA);
720
 
  dpd_file2_mat_init(&fia);
721
 
 
722
 
  /* One-electron (frozen-core) contributions */
723
 
  for(h=0; h < nirreps; h++) {
724
 
 
725
 
      for(i=0; i < occpi[h]; i++) 
726
 
          for(a=0; a < (virtpi[h] - openpi[h]); a++) 
727
 
              fIA.matrix[h][i][a] = Hov.matrix[h][i][a];   
728
 
 
729
 
      for(i=0; i < (occpi[h] - openpi[h]); i++) 
730
 
          for(a=0; a < virtpi[h]; a++) 
731
 
              fia.matrix[h][i][a] = Hov.matrix[h][i][a];   
732
 
    }
733
 
 
734
 
  /* Two-electron contributions */
735
 
 
736
 
  /* Prepare the E integral buffers */
737
 
  dpd_buf4_init(&EInts_anti, CC_EINTS, 0, 11, 0, 11, 0, 1, "E <ai|jk>");
738
 
  dpd_buf4_init(&EInts, CC_EINTS, 0, 11, 0, 11, 0, 0, "E <ai|jk>");
739
 
 
740
 
  for(h=0; h < nirreps; h++) {
741
 
 
742
 
      dpd_buf4_mat_irrep_init(&EInts_anti, h);
743
 
      dpd_buf4_mat_irrep_rd(&EInts_anti, h);
744
 
 
745
 
      /* Loop over irreps of the target */
746
 
      for(Gi=0; Gi < nirreps; Gi++) {
747
 
          Ga = Gi; Gm = h^Gi;
748
 
 
749
 
          /* Loop over orbitals of the target */
750
 
          for(i=0; i < occpi[Gi]; i++) {
751
 
              I = occ_off[Gi] + i;
752
 
              for(a=0; a < (virtpi[Ga] - openpi[Ga]); a++) {
753
 
                  A = vir_off[Ga] + a;
754
 
 
755
 
                  for(m=0; m < occpi[Gm]; m++) {
756
 
                      M = occ_off[Gm] + m;
757
 
 
758
 
                      AM = EInts_anti.params->rowidx[A][M];
759
 
                      IM = EInts_anti.params->colidx[I][M];
760
 
 
761
 
                      fIA.matrix[Gi][i][a] += EInts_anti.matrix[h][AM][IM];
762
 
 
763
 
                    }
764
 
                }
765
 
            }
766
 
 
767
 
          /* Loop over orbitals of the target */
768
 
          for(i=0; i < (occpi[Gi] - openpi[Gi]); i++) {
769
 
              I = occ_off[Gi] + i;
770
 
              for(a=0; a < virtpi[Ga]; a++) {
771
 
                  A = vir_off[Ga] + a;
772
 
 
773
 
                  for(m=0; m < (occpi[Gm] - openpi[Gm]); m++) {
774
 
                      M = occ_off[Gm] + m;
775
 
 
776
 
                      AM = EInts_anti.params->rowidx[A][M];
777
 
                      IM = EInts_anti.params->colidx[I][M];
778
 
 
779
 
                      fia.matrix[Gi][i][a] += EInts_anti.matrix[h][AM][IM];
780
 
 
781
 
                    }
782
 
                }
783
 
            }
784
 
        }
785
 
 
786
 
      dpd_buf4_mat_irrep_close(&EInts_anti, h);
787
 
 
788
 
      dpd_buf4_mat_irrep_init(&EInts, h);
789
 
      dpd_buf4_mat_irrep_rd(&EInts, h);
790
 
 
791
 
      /* Loop over irreps of the target */
792
 
      for(Gi=0; Gi < nirreps; Gi++) {
793
 
          Ga = Gi; Gm = h^Gi;
794
 
 
795
 
          /* Loop over orbitals of the target */
796
 
          for(i=0; i < occpi[Gi]; i++) {
797
 
              I = occ_off[Gi] + i;
798
 
              for(a=0; a < (virtpi[Ga] - openpi[Ga]); a++) {
799
 
                  A = vir_off[Ga] + a;
800
 
 
801
 
                  for(m=0; m < (occpi[Gm] - openpi[Gm]); m++) {
802
 
                      M = occ_off[Gm] + m;
803
 
 
804
 
                      AM = EInts.params->rowidx[A][M];
805
 
                      IM = EInts.params->colidx[I][M];
806
 
 
807
 
                      fIA.matrix[Gi][i][a] += EInts.matrix[h][AM][IM];
808
 
 
809
 
                    }
810
 
                }
811
 
            }
812
 
 
813
 
          /* Loop over orbitals of the target */
814
 
          for(i=0; i < (occpi[Gi] - openpi[Gi]); i++) {
815
 
              I = occ_off[Gi] + i;
816
 
              for(a=0; a < virtpi[Ga]; a++) {
817
 
                  A = vir_off[Ga] + a;
818
 
 
819
 
                  for(m=0; m < occpi[Gm]; m++) {
820
 
                      M = occ_off[Gm] + m;
821
 
 
822
 
                      AM = EInts.params->rowidx[A][M];
823
 
                      IM = EInts.params->colidx[I][M];
824
 
 
825
 
                      fia.matrix[Gi][i][a] += EInts.matrix[h][AM][IM];
826
 
 
827
 
                    }
828
 
                }
829
 
            }
830
 
        }
831
 
 
832
 
      dpd_buf4_mat_irrep_close(&EInts, h);
833
 
 
834
 
    }
835
 
 
836
 
  /* Close the E integral buffers */
837
 
  dpd_buf4_close(&EInts_anti);
838
 
  dpd_buf4_close(&EInts);
839
 
 
840
 
  /* Close the alpha and beta occ-vir Fock matrix files */
841
 
  dpd_file2_mat_wrt(&fIA);
842
 
  dpd_file2_mat_wrt(&fia);
843
 
  dpd_file2_mat_close(&fIA);
844
 
  dpd_file2_mat_close(&fia);
845
 
  dpd_file2_close(&fIA);
846
 
  dpd_file2_close(&fia);
847
 
 
848
 
  dpd_file2_mat_close(&Hoo);
849
 
  dpd_file2_mat_close(&Hvv);
850
 
  dpd_file2_mat_close(&Hov);
851
 
  dpd_file2_close(&Hoo);
852
 
  dpd_file2_close(&Hvv);
853
 
  dpd_file2_close(&Hov);
854
 
 
855
 
}