← Back to branch summary

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

« back to all changes in this revision

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