← Back to branch summary

~ubuntu-branches/ubuntu/quantal/psicode/quantal

« back to all changes in this revision

Viewing changes to src/bin/mp2/fold.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/mp2/fold.c
1
 
#include <libdpd/dpd.h>
2
 
#define EXTERN
3
 
#include "globals.h"
4
 
 
5
 
/* FOLD(): Fold the Fock matrix contributions to the energy (or energy
6
 
** derivative) into the two-particle density matrix.  Here we are
7
 
** trying to convert from an energy expression of the form:
8
 
**
9
 
** E = sum_pq Dpq fpq + 1/4 sum_pqrs Gpqrs <pq||rs>
10
 
**
11
 
** to the form:
12
 
**
13
 
** E = sum_pq Dpq hpq + 1/4 sum_pqrs Gpqrs <pq||rs>
14
 
**
15
 
** We do this by shifting some one-particle density matrix components
16
 
** into appropriate two-particle density matrix components:
17
 
**
18
 
** G'pmrm = Dpr + 4 * Gpmrm
19
 
**
20
 
** One problem is that we need to make sure the resulting density,
21
 
** G'pqrs, is still antisymmetric to permutation of p and q or r and
22
 
** s.  So, for example, for the Gimkm component we compute:
23
 
**
24
 
** G'pmrm = Dpr + Gpmrm
25
 
** G'mprm = Dpr - Gmprm
26
 
** G'pmmr = Dpr - Gpmmr
27
 
** G'mpmr = Dpr + Gmpmr
28
 
** */
29
 
 
30
 
void rhf_sf_fold(void);
31
 
void uhf_sf_fold(void);
32
 
 
33
 
void fold(void)
34
 
{
35
 
  if(params.ref == 0) rhf_sf_fold();
36
 
  else if(params.ref == 2) uhf_sf_fold();
37
 
}
38
 
 
39
 
void rhf_sf_fold(void)
40
 
{
41
 
  int h, nirreps;
42
 
  int i, j, k, l, m, a, b;
43
 
  int I, J, K, L, M, A, B;
44
 
  int IM, JM, MI, MJ, MK, ML, MA, MB;
45
 
  int Gi, Gj, Gk, Gl, Gm, Ga, Gb;
46
 
  int *occpi, *virtpi;
47
 
  int *occ_off, *vir_off;
48
 
  int *occ_sym, *vir_sym;
49
 
  dpdfile2 D, D1, D2, F;
50
 
  dpdbuf4 G, Aints, E, C, DInts, FInts, BInts, G1, G2;
51
 
  double one_energy=0.0, two_energy=0.0, total_two_energy=0.0;
52
 
  double test_energy = 0.0, tmp;
53
 
  double this_energy;
54
 
 
55
 
  nirreps = mo.nirreps;
56
 
  occpi = mo.occpi; 
57
 
  virtpi = mo.virtpi;
58
 
  occ_off = mo.occ_off; 
59
 
  vir_off = mo.vir_off;
60
 
  occ_sym = mo.occ_sym; 
61
 
  vir_sym = mo.vir_sym;
62
 
 
63
 
  fprintf(outfile, "\n\tEnergies re-computed from Fock-adjusted CC density:\n");
64
 
  fprintf(outfile,   "\t---------------------------------------------------\n");
65
 
 
66
 
  dpd_file2_init(&D, CC_OEI, 0, 0, 0, "DIJ");
67
 
  dpd_file2_init(&F, CC_OEI, 0, 0, 0, "h(i,j)");
68
 
  this_energy = dpd_file2_dot(&D, &F);
69
 
  dpd_file2_close(&F);
70
 
  dpd_file2_close(&D);
71
 
 
72
 
  /*    fprintf(outfile, "\tDIJ = %20.15f\n", this_energy); */
73
 
  one_energy += this_energy;
74
 
 
75
 
  dpd_file2_init(&D, CC_OEI, 0, 0, 0, "Dij");
76
 
  dpd_file2_init(&F, CC_OEI, 0, 0, 0, "h(i,j)");
77
 
  this_energy = dpd_file2_dot(&D, &F);
78
 
  dpd_file2_close(&F);
79
 
  dpd_file2_close(&D);
80
 
 
81
 
  /*    fprintf(outfile, "\tDij = %20.15f\n", this_energy); */
82
 
  one_energy += this_energy;
83
 
 
84
 
  dpd_file2_init(&D, CC_OEI, 0, 1, 1, "DAB");
85
 
  dpd_file2_init(&F, CC_OEI, 0, 1, 1, "h(a,b)");
86
 
  this_energy = dpd_file2_dot(&D, &F);
87
 
  dpd_file2_close(&F);
88
 
  dpd_file2_close(&D);
89
 
 
90
 
  /*    fprintf(outfile, "\tDAB = %20.15f\n", this_energy); */
91
 
  one_energy += this_energy;
92
 
 
93
 
  dpd_file2_init(&D, CC_OEI, 0, 1, 1, "Dab");
94
 
  dpd_file2_init(&F, CC_OEI, 0, 1, 1, "h(a,b)");
95
 
  this_energy = dpd_file2_dot(&D, &F);
96
 
  dpd_file2_close(&F);
97
 
  dpd_file2_close(&D);
98
 
 
99
 
  /*    fprintf(outfile, "\tDab = %20.15f\n", this_energy); */
100
 
  one_energy += this_energy;
101
 
 
102
 
  dpd_file2_init(&D, CC_OEI, 0, 0, 1, "DIA");
103
 
  dpd_file2_init(&F, CC_OEI, 0, 0, 1, "h(i,a)");
104
 
  this_energy = dpd_file2_dot(&D, &F);
105
 
  dpd_file2_close(&F);
106
 
  dpd_file2_close(&D);
107
 
 
108
 
  /*    fprintf(outfile, "\tDIA = %20.15f\n", this_energy); */
109
 
  one_energy += this_energy;
110
 
 
111
 
  dpd_file2_init(&D, CC_OEI, 0, 0, 1, "Dia");
112
 
  dpd_file2_init(&F, CC_OEI, 0, 0, 1, "h(i,a)");
113
 
  this_energy = dpd_file2_dot(&D, &F);
114
 
  dpd_file2_close(&F);
115
 
  dpd_file2_close(&D);
116
 
 
117
 
  /*    fprintf(outfile, "\tDia = %20.15f\n", this_energy); */
118
 
  one_energy += this_energy;
119
 
 
120
 
  dpd_file2_init(&D, CC_OEI, 0, 0, 1, "DAI");
121
 
  dpd_file2_init(&F, CC_OEI, 0, 0, 1, "h(i,a)");
122
 
  this_energy = dpd_file2_dot(&D, &F);
123
 
  dpd_file2_close(&F);
124
 
    dpd_file2_close(&D);
125
 
 
126
 
  /*    fprintf(outfile, "\tDAI = %20.15f\n", this_energy); */
127
 
  one_energy += this_energy;
128
 
 
129
 
  dpd_file2_init(&D, CC_OEI, 0, 0, 1, "Dai");
130
 
  dpd_file2_init(&F, CC_OEI, 0, 0, 1, "h(i,a)");
131
 
  this_energy = dpd_file2_dot(&D, &F);
132
 
  dpd_file2_close(&F);
133
 
  dpd_file2_close(&D);
134
 
 
135
 
  /*    fprintf(outfile, "\tDai = %20.15f\n", this_energy); */
136
 
  one_energy += this_energy;
137
 
 
138
 
  fprintf(outfile, "\tOne-electron energy        = %20.15f\n", one_energy);
139
 
  fflush(outfile);
140
 
 
141
 
  dpd_file2_init(&D, CC_OEI, 0, 0, 0, "DIJ");
142
 
  dpd_file2_mat_init(&D);
143
 
  dpd_file2_mat_rd(&D);
144
 
 
145
 
  dpd_buf4_init(&G, CC_GAMMA, 0, 0, 0, 2, 2, 0, "GIJKL");
146
 
  dpd_buf4_scm(&G,0.0);
147
 
  dpd_buf4_close(&G);
148
 
 
149
 
  dpd_buf4_init(&G, CC_GAMMA, 0, 0, 0, 2, 2, 0, "GIJKL");
150
 
  for(h=0; h < nirreps; h++) {
151
 
    dpd_buf4_mat_irrep_init(&G, h);
152
 
    dpd_buf4_mat_irrep_rd(&G, h);
153
 
 
154
 
    for(Gm=0; Gm < nirreps; Gm++) {
155
 
      Gi = Gj = h^Gm;
156
 
 
157
 
      for(i=0; i < occpi[Gi]; i++) {
158
 
        I = occ_off[Gi] + i;
159
 
        for(j=0; j < occpi[Gj]; j++) {
160
 
          J = occ_off[Gj] + j;
161
 
          for(m=0; m < occpi[Gm]; m++) {
162
 
            M = occ_off[Gm] + m;
163
 
 
164
 
            IM = G.params->rowidx[I][M];
165
 
            JM = G.params->colidx[J][M];
166
 
            MI = G.params->rowidx[M][I];
167
 
            MJ = G.params->colidx[M][J];
168
 
 
169
 
            G.matrix[h][IM][JM] += D.matrix[Gi][i][j];
170
 
            G.matrix[h][IM][MJ] -= D.matrix[Gi][i][j];
171
 
            G.matrix[h][MI][MJ] += D.matrix[Gi][i][j];
172
 
            G.matrix[h][MI][JM] -= D.matrix[Gi][i][j];
173
 
          }
174
 
        }
175
 
      }
176
 
    }
177
 
 
178
 
    dpd_buf4_mat_irrep_wrt(&G, h);
179
 
    dpd_buf4_mat_irrep_close(&G, h);
180
 
  }
181
 
      
182
 
  two_energy = 0.0;
183
 
  dpd_buf4_init(&Aints, CC_AINTS, 0, 0, 0, 0, 0, 1, "A <ij|kl>");
184
 
  two_energy += 0.25 * dpd_buf4_dot(&Aints, &G);
185
 
  dpd_buf4_close(&Aints);
186
 
 
187
 
  dpd_buf4_close(&G);
188
 
 
189
 
  dpd_file2_mat_close(&D);
190
 
  dpd_file2_close(&D);
191
 
 
192
 
  dpd_file2_init(&D, CC_OEI, 0, 0, 0, "Dij");
193
 
  dpd_file2_mat_init(&D);
194
 
  dpd_file2_mat_rd(&D);
195
 
 
196
 
  dpd_buf4_init(&G, CC_GAMMA, 0, 0, 0, 2, 2, 0, "Gijkl");
197
 
  dpd_buf4_scm(&G,0.0);
198
 
  dpd_buf4_close(&G);
199
 
 
200
 
  dpd_buf4_init(&G, CC_GAMMA, 0, 0, 0, 2, 2, 0, "Gijkl");
201
 
  for(h=0; h < nirreps; h++) {
202
 
    dpd_buf4_mat_irrep_init(&G, h);
203
 
    dpd_buf4_mat_irrep_rd(&G, h);
204
 
 
205
 
    for(Gm=0; Gm < nirreps; Gm++) {
206
 
      Gi = Gj = h^Gm;
207
 
 
208
 
      for(i=0; i < occpi[Gi]; i++) {
209
 
        I = occ_off[Gi] + i;
210
 
        for(j=0; j < occpi[Gj]; j++) {
211
 
          J = occ_off[Gj] + j;
212
 
          for(m=0; m < occpi[Gm]; m++) {
213
 
            M = occ_off[Gm] + m;
214
 
 
215
 
            IM = G.params->rowidx[I][M];
216
 
            JM = G.params->colidx[J][M];
217
 
            MI = G.params->rowidx[M][I];
218
 
            MJ = G.params->colidx[M][J];
219
 
 
220
 
            G.matrix[h][IM][JM] += D.matrix[Gi][i][j];
221
 
            G.matrix[h][MI][JM] -= D.matrix[Gi][i][j];
222
 
            G.matrix[h][MI][MJ] += D.matrix[Gi][i][j];
223
 
            G.matrix[h][IM][MJ] -= D.matrix[Gi][i][j];
224
 
          }
225
 
        }
226
 
      }
227
 
    }
228
 
 
229
 
    dpd_buf4_mat_irrep_wrt(&G, h);
230
 
    dpd_buf4_mat_irrep_close(&G, h);
231
 
  }
232
 
 
233
 
  dpd_buf4_init(&Aints, CC_AINTS, 0, 0, 0, 0, 0, 1, "A <ij|kl>");
234
 
  two_energy += 0.25 * dpd_buf4_dot(&Aints, &G);
235
 
  dpd_buf4_close(&Aints);
236
 
 
237
 
  dpd_buf4_close(&G);
238
 
 
239
 
  dpd_file2_mat_close(&D);
240
 
  dpd_file2_close(&D);
241
 
 
242
 
  dpd_file2_init(&D, CC_OEI, 0, 0, 0, "DIJ");
243
 
  dpd_file2_mat_init(&D);
244
 
  dpd_file2_mat_rd(&D);
245
 
 
246
 
  dpd_buf4_init(&G, CC_GAMMA, 0, 0, 0, 0, 0, 0, "GIjKl");
247
 
  dpd_buf4_scm(&G,0.0);
248
 
  dpd_buf4_close(&G);
249
 
 
250
 
  dpd_buf4_init(&G, CC_GAMMA, 0, 0, 0, 0, 0, 0, "GIjKl");
251
 
  for(h=0; h < nirreps; h++) {
252
 
    dpd_buf4_mat_irrep_init(&G, h);
253
 
    dpd_buf4_mat_irrep_rd(&G, h);
254
 
 
255
 
    for(Gm=0; Gm < nirreps; Gm++) {
256
 
      Gi = Gj = h^Gm;
257
 
 
258
 
      for(i=0; i < occpi[Gi]; i++) {
259
 
        I = occ_off[Gi] + i;
260
 
        for(j=0; j < occpi[Gj]; j++) {
261
 
          J = occ_off[Gj] + j;
262
 
          for(m=0; m < occpi[Gm]; m++) {
263
 
            M = occ_off[Gm] + m;
264
 
 
265
 
            IM = G.params->rowidx[I][M];
266
 
            JM = G.params->colidx[J][M];
267
 
 
268
 
            G.matrix[h][IM][JM] += D.matrix[Gi][i][j];
269
 
          }
270
 
        }
271
 
      }
272
 
    }
273
 
 
274
 
    dpd_buf4_mat_irrep_wrt(&G, h); 
275
 
    dpd_buf4_mat_irrep_close(&G, h);
276
 
  }
277
 
  dpd_buf4_close(&G);
278
 
 
279
 
  dpd_file2_mat_close(&D);
280
 
  dpd_file2_close(&D);
281
 
 
282
 
 
283
 
  dpd_file2_init(&D, CC_OEI, 0, 0, 0, "Dij");
284
 
  dpd_file2_mat_init(&D);
285
 
  dpd_file2_mat_rd(&D);
286
 
 
287
 
  dpd_buf4_init(&G, CC_GAMMA, 0, 0, 0, 0, 0, 0, "GIjKl");
288
 
  for(h=0; h < nirreps; h++) {
289
 
    dpd_buf4_mat_irrep_init(&G, h);
290
 
    dpd_buf4_mat_irrep_rd(&G, h);
291
 
 
292
 
    for(Gm=0; Gm < nirreps; Gm++) {
293
 
      Gk = Gl = h^Gm;
294
 
 
295
 
      for(k=0; k < occpi[Gk]; k++) {
296
 
        K = occ_off[Gk] + k;
297
 
        for(l=0; l < occpi[Gl]; l++) {
298
 
          L = occ_off[Gl] + l;
299
 
          for(m=0; m < occpi[Gm]; m++) {
300
 
            M = occ_off[Gm] + m;
301
 
 
302
 
            MK = G.params->rowidx[M][K];
303
 
            ML = G.params->colidx[M][L];
304
 
 
305
 
            G.matrix[h][MK][ML] += D.matrix[Gk][k][l];
306
 
          }
307
 
        }
308
 
      }
309
 
    }
310
 
 
311
 
    dpd_buf4_mat_irrep_wrt(&G, h);
312
 
    dpd_buf4_mat_irrep_close(&G, h);
313
 
  }
314
 
 
315
 
  dpd_buf4_init(&Aints, CC_AINTS, 0, 0, 0, 0, 0, 0, "A <ij|kl>");
316
 
  two_energy += dpd_buf4_dot(&Aints, &G);
317
 
  dpd_buf4_close(&Aints);
318
 
 
319
 
  dpd_buf4_close(&G);
320
 
 
321
 
  total_two_energy += two_energy;
322
 
  fprintf(outfile, "\tIJKL energy                = %20.15f\n", two_energy);
323
 
  fflush(outfile);
324
 
 
325
 
  dpd_file2_mat_close(&D);
326
 
  dpd_file2_close(&D);
327
 
 
328
 
  dpd_file2_init(&D1, CC_OEI, 0, 0, 1, "DIA");
329
 
  dpd_file2_mat_init(&D1);
330
 
  dpd_file2_mat_rd(&D1);
331
 
  dpd_file2_init(&D2, CC_OEI, 0, 0, 1, "DAI");
332
 
  dpd_file2_mat_init(&D2);
333
 
  dpd_file2_mat_rd(&D2);
334
 
 
335
 
  dpd_buf4_init(&G, CC_GAMMA, 0, 0, 10, 2, 10, 0, "GIJKA");
336
 
  dpd_buf4_scm(&G,0.0);
337
 
  dpd_buf4_close(&G);
338
 
 
339
 
  dpd_buf4_init(&G, CC_GAMMA, 0, 0, 10, 2, 10, 0, "GIJKA");
340
 
  for(h=0; h < nirreps; h++) {
341
 
    dpd_buf4_mat_irrep_init(&G, h);
342
 
    dpd_buf4_mat_irrep_rd(&G, h);
343
 
 
344
 
    for(Gm=0; Gm < nirreps; Gm++) {
345
 
      Gi = Ga = h^Gm;
346
 
 
347
 
      for(i=0; i < occpi[Gi]; i++) {
348
 
        I = occ_off[Gi] + i;
349
 
        for(a=0; a < virtpi[Ga]; a++) {
350
 
          A = vir_off[Ga] + a;
351
 
          for(m=0; m < occpi[Gm]; m++) {
352
 
            M = occ_off[Gm] + m;
353
 
 
354
 
            MI = G.params->rowidx[M][I];
355
 
            IM = G.params->rowidx[I][M];
356
 
            MA = G.params->colidx[M][A];
357
 
 
358
 
            G.matrix[h][MI][MA] += 0.5 * (D1.matrix[Gi][i][a] +
359
 
                                          D2.matrix[Gi][i][a]);
360
 
            G.matrix[h][IM][MA] -= 0.5 * (D1.matrix[Gi][i][a] +
361
 
                                          D2.matrix[Gi][i][a]);
362
 
          }
363
 
        }
364
 
      }
365
 
    }
366
 
 
367
 
    dpd_buf4_mat_irrep_wrt(&G, h);
368
 
    dpd_buf4_mat_irrep_close(&G, h);
369
 
  }
370
 
 
371
 
  two_energy = 0.0;
372
 
  dpd_buf4_init(&E, CC_EINTS, 0, 0, 10, 2, 10, 0, "E <ij||ka> (i>j,ka)");
373
 
  two_energy += dpd_buf4_dot(&E, &G);
374
 
  dpd_buf4_close(&E);
375
 
 
376
 
  dpd_buf4_close(&G);
377
 
 
378
 
  dpd_file2_mat_close(&D1);
379
 
  dpd_file2_close(&D1);
380
 
  dpd_file2_mat_close(&D2);
381
 
  dpd_file2_close(&D2);
382
 
 
383
 
  dpd_file2_init(&D1, CC_OEI, 0, 0, 1, "Dia");
384
 
  dpd_file2_mat_init(&D1);
385
 
  dpd_file2_mat_rd(&D1);
386
 
  dpd_file2_init(&D2, CC_OEI, 0, 0, 1, "Dai");
387
 
  dpd_file2_mat_init(&D2);
388
 
  dpd_file2_mat_rd(&D2);
389
 
 
390
 
  dpd_buf4_init(&G, CC_GAMMA, 0, 0, 10, 2, 10, 0, "Gijka");
391
 
  dpd_buf4_scm(&G,0.0);
392
 
  dpd_buf4_close(&G);
393
 
 
394
 
  dpd_buf4_init(&G, CC_GAMMA, 0, 0, 10, 2, 10, 0, "Gijka");
395
 
  for(h=0; h < nirreps; h++) {
396
 
    dpd_buf4_mat_irrep_init(&G, h);
397
 
    dpd_buf4_mat_irrep_rd(&G, h);
398
 
 
399
 
    for(Gm=0; Gm < nirreps; Gm++) {
400
 
      Gi = Ga = h^Gm;
401
 
 
402
 
      for(i=0; i < occpi[Gi]; i++) {
403
 
        I = occ_off[Gi] + i;
404
 
        for(a=0; a < virtpi[Ga]; a++) {
405
 
          A = vir_off[Ga] + a;
406
 
          for(m=0; m < occpi[Gm]; m++) {
407
 
            M = occ_off[Gm] + m;
408
 
 
409
 
            MI = G.params->rowidx[M][I];
410
 
            IM = G.params->rowidx[I][M];
411
 
            MA = G.params->colidx[M][A];
412
 
 
413
 
            G.matrix[h][MI][MA] += 0.5 * (D1.matrix[Gi][i][a] +
414
 
                                          D2.matrix[Gi][i][a]);
415
 
            G.matrix[h][IM][MA] -= 0.5 * (D1.matrix[Gi][i][a] +
416
 
                                          D2.matrix[Gi][i][a]);
417
 
          }
418
 
        }
419
 
      }
420
 
    }
421
 
 
422
 
    dpd_buf4_mat_irrep_wrt(&G, h);
423
 
    dpd_buf4_mat_irrep_close(&G, h);
424
 
  }
425
 
 
426
 
  dpd_buf4_init(&E, CC_EINTS, 0, 0, 10, 2, 10, 0, "E <ij||ka> (i>j,ka)");
427
 
  two_energy += dpd_buf4_dot(&E, &G);
428
 
  dpd_buf4_close(&E);
429
 
 
430
 
  dpd_buf4_close(&G);
431
 
 
432
 
  dpd_file2_mat_close(&D1);
433
 
  dpd_file2_close(&D1);
434
 
  dpd_file2_mat_close(&D2);
435
 
  dpd_file2_close(&D2);
436
 
 
437
 
  dpd_file2_init(&D1, CC_OEI, 0, 0, 1, "DIA");
438
 
  dpd_file2_mat_init(&D1);
439
 
  dpd_file2_mat_rd(&D1);
440
 
  dpd_file2_init(&D2, CC_OEI, 0, 0, 1, "DAI");
441
 
  dpd_file2_mat_init(&D2);
442
 
  dpd_file2_mat_rd(&D2);
443
 
 
444
 
  dpd_buf4_init(&G, CC_GAMMA, 0, 0, 10, 0, 10, 0, "GiJkA");
445
 
  dpd_buf4_scm(&G,0.0);
446
 
  dpd_buf4_close(&G);
447
 
 
448
 
  dpd_buf4_init(&G, CC_GAMMA, 0, 0, 10, 0, 10, 0, "GiJkA");
449
 
  for(h=0; h < nirreps; h++) {
450
 
    dpd_buf4_mat_irrep_init(&G, h);
451
 
    dpd_buf4_mat_irrep_rd(&G, h);
452
 
 
453
 
    for(Gm=0; Gm < nirreps; Gm++) {
454
 
      Gi = Ga = h^Gm;
455
 
 
456
 
      for(i=0; i < occpi[Gi]; i++) {
457
 
        I = occ_off[Gi] + i;
458
 
        for(a=0; a < virtpi[Ga]; a++) {
459
 
          A = vir_off[Ga] + a;
460
 
          for(m=0; m < occpi[Gm]; m++) {
461
 
            M = occ_off[Gm] + m;
462
 
 
463
 
            MI = G.params->rowidx[M][I];
464
 
            MA = G.params->colidx[M][A];
465
 
 
466
 
            G.matrix[h][MI][MA] += 0.5 * (D1.matrix[Gi][i][a] +
467
 
                                          D2.matrix[Gi][i][a]);
468
 
          }
469
 
        }
470
 
      }
471
 
    }
472
 
 
473
 
    dpd_buf4_mat_irrep_wrt(&G, h);
474
 
    dpd_buf4_mat_irrep_close(&G, h);
475
 
  }
476
 
 
477
 
  dpd_buf4_init(&E, CC_EINTS, 0, 0, 10, 0, 10, 0, "E <ij|ka>");
478
 
  two_energy += 2 * dpd_buf4_dot(&E, &G);
479
 
  dpd_buf4_close(&E);
480
 
 
481
 
  dpd_buf4_close(&G);
482
 
 
483
 
  dpd_file2_mat_close(&D1);
484
 
  dpd_file2_close(&D1);
485
 
  dpd_file2_mat_close(&D2);
486
 
  dpd_file2_close(&D2);
487
 
 
488
 
  dpd_file2_init(&D1, CC_OEI, 0, 0, 1, "Dia");
489
 
  dpd_file2_mat_init(&D1);
490
 
  dpd_file2_mat_rd(&D1);
491
 
  dpd_file2_init(&D2, CC_OEI, 0, 0, 1, "Dai");
492
 
  dpd_file2_mat_init(&D2);
493
 
  dpd_file2_mat_rd(&D2);
494
 
 
495
 
  dpd_buf4_init(&G, CC_GAMMA, 0, 0, 10, 0, 10, 0, "GIjKa");
496
 
  dpd_buf4_scm(&G,0.0);
497
 
  dpd_buf4_close(&G);
498
 
 
499
 
  dpd_buf4_init(&G, CC_GAMMA, 0, 0, 10, 0, 10, 0, "GIjKa");
500
 
  for(h=0; h < nirreps; h++) {
501
 
    dpd_buf4_mat_irrep_init(&G, h);
502
 
    dpd_buf4_mat_irrep_rd(&G, h);
503
 
 
504
 
    for(Gm=0; Gm < nirreps; Gm++) {
505
 
      Gi = Ga = h^Gm;
506
 
 
507
 
      for(i=0; i < occpi[Gi]; i++) {
508
 
        I = occ_off[Gi] + i;
509
 
        for(a=0; a < virtpi[Ga]; a++) {
510
 
          A = vir_off[Ga] + a;
511
 
          for(m=0; m < occpi[Gm]; m++) {
512
 
            M = occ_off[Gm] + m;
513
 
 
514
 
            MI = G.params->rowidx[M][I];
515
 
            MA = G.params->colidx[M][A];
516
 
 
517
 
            G.matrix[h][MI][MA] += 0.5 * (D1.matrix[Gi][i][a] +
518
 
                                          D2.matrix[Gi][i][a]);
519
 
          }
520
 
        }
521
 
      }
522
 
    }
523
 
 
524
 
    dpd_buf4_mat_irrep_wrt(&G, h);
525
 
    dpd_buf4_mat_irrep_close(&G, h);
526
 
  }
527
 
 
528
 
  dpd_buf4_init(&E, CC_EINTS, 0, 0, 10, 0, 10, 0, "E <ij|ka>");
529
 
  two_energy += 2 * dpd_buf4_dot(&E, &G);
530
 
  dpd_buf4_close(&E);
531
 
 
532
 
  dpd_buf4_close(&G);
533
 
 
534
 
  total_two_energy += two_energy;
535
 
  fprintf(outfile, "\tIJKA energy                = %20.15f\n", two_energy);
536
 
  fflush(outfile);
537
 
 
538
 
  dpd_file2_mat_close(&D1);
539
 
  dpd_file2_close(&D1);
540
 
  dpd_file2_mat_close(&D2);
541
 
  dpd_file2_close(&D2);
542
 
 
543
 
  two_energy = 0.0;
544
 
  dpd_buf4_init(&DInts, CC_DINTS, 0, 2, 7, 2, 7, 0, "D <ij||ab> (i>j,a>b)");
545
 
  dpd_buf4_init(&G, CC_GAMMA, 0, 2, 7, 2, 7, 0, "GIJAB");
546
 
  two_energy += dpd_buf4_dot(&G, &DInts);
547
 
  dpd_buf4_close(&G);
548
 
  dpd_buf4_init(&G, CC_GAMMA, 0, 2, 7, 2, 7, 0, "Gijab");
549
 
  two_energy += dpd_buf4_dot(&G, &DInts); 
550
 
  dpd_buf4_close(&G);
551
 
  dpd_buf4_close(&DInts);
552
 
  dpd_buf4_init(&DInts, CC_DINTS, 0, 0, 5, 0, 5, 0, "D <ij|ab>");
553
 
  dpd_buf4_init(&G, CC_GAMMA, 0, 0, 5, 0, 5, 0, "GIjAb");
554
 
  two_energy += dpd_buf4_dot(&G, &DInts);
555
 
  dpd_buf4_close(&G);
556
 
  dpd_buf4_close(&DInts);
557
 
 
558
 
  two_energy *= 2;
559
 
  total_two_energy += two_energy;
560
 
  fprintf(outfile, "\tIJAB energy                = %20.15f\n", two_energy);
561
 
  fflush(outfile);
562
 
 
563
 
  dpd_file2_init(&D, CC_OEI, 0, 1, 1, "DAB");
564
 
  dpd_file2_mat_init(&D);
565
 
  dpd_file2_mat_rd(&D);
566
 
 
567
 
  dpd_buf4_init(&G, CC_GAMMA, 0, 10, 10, 10, 10, 0, "GIBJA");
568
 
  dpd_buf4_scm(&G, 0.0);
569
 
  dpd_buf4_close(&G);
570
 
 
571
 
  dpd_buf4_init(&G, CC_GAMMA, 0, 10, 10, 10, 10, 0, "GIBJA");
572
 
  for(h=0; h < nirreps; h++) {
573
 
    dpd_buf4_mat_irrep_init(&G, h);
574
 
    dpd_buf4_mat_irrep_rd(&G, h);
575
 
 
576
 
    for(Gm=0; Gm < nirreps; Gm++) {
577
 
      Ga = Gb = h^Gm;
578
 
 
579
 
      for(b=0; b < virtpi[Gb]; b++) {
580
 
        B = vir_off[Gb] + b;
581
 
        for(a=0; a < virtpi[Ga]; a++) {
582
 
          A = vir_off[Ga] + a;
583
 
          for(m=0; m < occpi[Gm]; m++) {
584
 
            M = occ_off[Gm] + m;
585
 
 
586
 
            MB = G.params->rowidx[M][B];
587
 
            MA = G.params->colidx[M][A];
588
 
 
589
 
            G.matrix[h][MB][MA] += D.matrix[Ga][a][b];
590
 
          }
591
 
        }
592
 
      }
593
 
    }
594
 
 
595
 
    dpd_buf4_mat_irrep_wrt(&G, h);
596
 
    dpd_buf4_mat_irrep_close(&G, h);
597
 
  }
598
 
 
599
 
  two_energy = 0.0;
600
 
  dpd_buf4_init(&C, CC_CINTS, 0, 10, 10, 10, 10, 0, "C <ia||jb>");
601
 
  two_energy += dpd_buf4_dot(&C, &G);
602
 
  dpd_buf4_close(&C);
603
 
 
604
 
  dpd_buf4_close(&G);
605
 
 
606
 
  dpd_file2_mat_close(&D);
607
 
  dpd_file2_close(&D);
608
 
 
609
 
 
610
 
  dpd_file2_init(&D, CC_OEI, 0, 1, 1, "Dab");
611
 
  dpd_file2_mat_init(&D);
612
 
  dpd_file2_mat_rd(&D);
613
 
 
614
 
  dpd_buf4_init(&G, CC_GAMMA, 0, 10, 10, 10, 10, 0, "Gibja");
615
 
  dpd_buf4_scm(&G, 0.0);
616
 
  dpd_buf4_close(&G);
617
 
 
618
 
  dpd_buf4_init(&G, CC_GAMMA, 0, 10, 10, 10, 10, 0, "Gibja");
619
 
  for(h=0; h < nirreps; h++) {
620
 
    dpd_buf4_mat_irrep_init(&G, h);
621
 
    dpd_buf4_mat_irrep_rd(&G, h);
622
 
 
623
 
    for(Gm=0; Gm < nirreps; Gm++) {
624
 
      Ga = Gb = h^Gm;
625
 
 
626
 
      for(b=0; b < virtpi[Gb]; b++) {
627
 
        B = vir_off[Gb] + b;
628
 
        for(a=0; a < virtpi[Ga]; a++) {
629
 
          A = vir_off[Ga] + a;
630
 
          for(m=0; m < occpi[Gm]; m++) {
631
 
            M = occ_off[Gm] + m;
632
 
 
633
 
            MB = G.params->rowidx[M][B];
634
 
            MA = G.params->colidx[M][A];
635
 
 
636
 
            G.matrix[h][MB][MA] += D.matrix[Ga][a][b];
637
 
          }
638
 
        }
639
 
      }
640
 
    }
641
 
 
642
 
    dpd_buf4_mat_irrep_wrt(&G, h);
643
 
    dpd_buf4_mat_irrep_close(&G, h);
644
 
  }
645
 
 
646
 
  dpd_buf4_init(&C, CC_CINTS, 0, 10, 10, 10, 10, 0, "C <ia||jb>");
647
 
  two_energy += dpd_buf4_dot(&C, &G);
648
 
  dpd_buf4_close(&C);
649
 
 
650
 
  dpd_buf4_close(&G);
651
 
 
652
 
  dpd_file2_mat_close(&D);
653
 
  dpd_file2_close(&D);
654
 
 
655
 
  dpd_file2_init(&D, CC_OEI, 0, 1, 1, "Dab");
656
 
  dpd_file2_mat_init(&D);
657
 
  dpd_file2_mat_rd(&D);
658
 
 
659
 
  dpd_buf4_init(&G, CC_GAMMA, 0, 10, 10, 10, 10, 0, "GIbJa");
660
 
  dpd_buf4_scm(&G, 0.0);
661
 
  dpd_buf4_close(&G);
662
 
 
663
 
  dpd_buf4_init(&G, CC_GAMMA, 0, 10, 10, 10, 10, 0, "GIbJa");
664
 
  for(h=0; h < nirreps; h++) {
665
 
    dpd_buf4_mat_irrep_init(&G, h);
666
 
    dpd_buf4_mat_irrep_rd(&G, h);
667
 
 
668
 
    for(Gm=0; Gm < nirreps; Gm++) {
669
 
      Ga = Gb = h^Gm;
670
 
 
671
 
      for(b=0; b < virtpi[Gb]; b++) {
672
 
        B = vir_off[Gb] + b;
673
 
        for(a=0; a < virtpi[Ga]; a++) {
674
 
          A = vir_off[Ga] + a;
675
 
          for(m=0; m < occpi[Gm]; m++) {
676
 
            M = occ_off[Gm] + m;
677
 
 
678
 
            MB = G.params->rowidx[M][B];
679
 
            MA = G.params->colidx[M][A];
680
 
 
681
 
            G.matrix[h][MB][MA] += D.matrix[Ga][a][b];
682
 
          }
683
 
        }
684
 
      }
685
 
    }
686
 
 
687
 
    dpd_buf4_mat_irrep_wrt(&G, h);
688
 
    dpd_buf4_mat_irrep_close(&G, h);
689
 
  }
690
 
 
691
 
  dpd_buf4_init(&C, CC_CINTS, 0, 10, 10, 10, 10, 0, "C <ia|jb>");
692
 
  two_energy += dpd_buf4_dot(&C, &G);
693
 
  dpd_buf4_close(&C);
694
 
 
695
 
  dpd_buf4_close(&G);
696
 
 
697
 
  dpd_file2_mat_close(&D);
698
 
  dpd_file2_close(&D);
699
 
 
700
 
  dpd_file2_init(&D, CC_OEI, 0, 1, 1, "DAB");
701
 
  dpd_file2_mat_init(&D);
702
 
  dpd_file2_mat_rd(&D);
703
 
 
704
 
  dpd_buf4_init(&G, CC_GAMMA, 0, 10, 10, 10, 10, 0, "GiBjA");
705
 
  dpd_buf4_scm(&G, 0.0);
706
 
  dpd_buf4_close(&G);
707
 
 
708
 
  dpd_buf4_init(&G, CC_GAMMA, 0, 10, 10, 10, 10, 0, "GiBjA");
709
 
  for(h=0; h < nirreps; h++) {
710
 
    dpd_buf4_mat_irrep_init(&G, h);
711
 
    dpd_buf4_mat_irrep_rd(&G, h);
712
 
 
713
 
    for(Gm=0; Gm < nirreps; Gm++) {
714
 
      Ga = Gb = h^Gm;
715
 
 
716
 
      for(b=0; b < virtpi[Gb]; b++) {
717
 
        B = vir_off[Gb] + b;
718
 
        for(a=0; a < virtpi[Ga]; a++) {
719
 
          A = vir_off[Ga] + a;
720
 
          for(m=0; m < occpi[Gm]; m++) {
721
 
            M = occ_off[Gm] + m;
722
 
 
723
 
            MB = G.params->rowidx[M][B];
724
 
            MA = G.params->colidx[M][A];
725
 
 
726
 
            G.matrix[h][MB][MA] += D.matrix[Ga][a][b];
727
 
          }
728
 
        }
729
 
      }
730
 
    }
731
 
    dpd_buf4_mat_irrep_wrt(&G, h);
732
 
    dpd_buf4_mat_irrep_close(&G, h);
733
 
  }
734
 
 
735
 
  dpd_buf4_init(&C, CC_CINTS, 0, 10, 10, 10, 10, 0, "C <ia|jb>");
736
 
  two_energy += dpd_buf4_dot(&C, &G); 
737
 
  dpd_buf4_close(&C);
738
 
 
739
 
  dpd_buf4_close(&G);
740
 
 
741
 
  dpd_file2_mat_close(&D);
742
 
  dpd_file2_close(&D);
743
 
 
744
 
  dpd_buf4_init(&DInts, CC_DINTS, 0, 10, 10, 10, 10, 0, "D <ij|ab> (ib,ja)");
745
 
  dpd_buf4_init(&G, CC_GAMMA, 0, 10, 10, 10, 10, 0, "GIbjA");
746
 
  two_energy -= dpd_buf4_dot(&G, &DInts);
747
 
  dpd_buf4_close(&G);
748
 
  dpd_buf4_init(&G, CC_GAMMA, 0, 10, 10, 10, 10, 0, "GiBJa");
749
 
  two_energy -= dpd_buf4_dot(&G, &DInts);
750
 
  dpd_buf4_close(&G);
751
 
  dpd_buf4_close(&DInts);
752
 
 
753
 
  total_two_energy += two_energy;
754
 
  fprintf(outfile, "\tIBJA energy                = %20.15f\n", two_energy);
755
 
  fflush(outfile);
756
 
 
757
 
  fprintf(outfile, "\tMP2 correlation energy    = %20.15f\n",
758
 
          one_energy + total_two_energy);
759
 
  fprintf(outfile, "\tTotal MP2 energy          = %20.15f\n",
760
 
          one_energy + total_two_energy + mo.Escf);
761
 
}
762
 
 
763
 
void uhf_sf_fold(void)
764
 
{
765
 
 
766
 
}