← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/cclambda/denom.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/cclambda/denom.c
1
 
#include <stdio.h>
2
 
#include <libdpd/dpd.h>
3
 
#define EXTERN
4
 
#include "globals.h"
5
 
 
6
 
void denom_rhf(struct L_Params);
7
 
void denom_rohf(struct L_Params);
8
 
void denom_uhf(struct L_Params);
9
 
 
10
 
void denom(struct L_Params L_params) {
11
 
  if(params.ref == 0) denom_rhf(L_params);
12
 
  else if(params.ref == 1) denom_rohf(L_params);
13
 
  else if(params.ref == 2) denom_uhf(L_params);
14
 
}
15
 
 
16
 
void denom_rhf(struct L_Params L_params)
17
 
{
18
 
  dpdfile2 FAE, FMI;
19
 
  dpdfile2 dIA;
20
 
  dpdfile4 dIjAb;
21
 
  dpdbuf4 d, bdIJAB, bdijab, bdIjAb;
22
 
  double tval;
23
 
  int nirreps,L_irr;
24
 
  int h, i, j, a, b, ij, ab;
25
 
  int I, J, A, B;
26
 
  int isym, jsym, asym, bsym;
27
 
  int *occpi, *virtpi;
28
 
  int *occ_off, *vir_off;
29
 
  int *openpi;
30
 
  double Fii, Fjj, Faa, Fbb;
31
 
 
32
 
  L_irr = L_params.irrep;
33
 
  nirreps = moinfo.nirreps;
34
 
  occpi = moinfo.occpi; virtpi = moinfo.virtpi;
35
 
  occ_off = moinfo.occ_off; vir_off = moinfo.vir_off;
36
 
 
37
 
  dpd_file2_init(&FMI, CC_OEI, 0, 0, 0, "FMI");
38
 
  dpd_file2_mat_init(&FMI);
39
 
  dpd_file2_mat_rd(&FMI);
40
 
 
41
 
  dpd_file2_init(&FAE, CC_OEI, 0, 1, 1, "FAE");
42
 
  dpd_file2_mat_init(&FAE);
43
 
  dpd_file2_mat_rd(&FAE);
44
 
 
45
 
  /* Alpha one-electron denominator */
46
 
  dpd_file2_init(&dIA, CC_DENOM, L_irr, 0, 1, "dIA");
47
 
  dpd_file2_mat_init(&dIA);
48
 
  for(h=0; h < nirreps; h++) { /* irreps of dIA and Fii */
49
 
    for(i=0; i < occpi[h]; i++) {
50
 
      Fii = FMI.matrix[h][i][i];
51
 
      for(a=0; a < virtpi[h^L_irr]; a++) {
52
 
        Faa = FAE.matrix[h^L_irr][a][a];
53
 
        dIA.matrix[h][i][a] = 1.0/(Fii - Faa + L_params.cceom_energy);
54
 
      }
55
 
    }
56
 
  }
57
 
  dpd_file2_mat_wrt(&dIA);
58
 
  dpd_file2_mat_close(&dIA);
59
 
  dpd_file2_close(&dIA);
60
 
 
61
 
  /* Alpha-beta two-electron denominator */
62
 
  dpd_file4_init(&dIjAb, CC_DENOM, L_irr, 0, 5, "dIjAb");
63
 
 
64
 
  for(h=0; h < nirreps; h++) {
65
 
    dpd_file4_mat_irrep_init(&dIjAb, h);
66
 
    /* Loop over the rows */
67
 
    for(ij=0; ij < dIjAb.params->rowtot[h]; ij++) {
68
 
          i = dIjAb.params->roworb[h][ij][0];
69
 
          j = dIjAb.params->roworb[h][ij][1];
70
 
          isym = dIjAb.params->psym[i];
71
 
          jsym = dIjAb.params->qsym[j];
72
 
 
73
 
          /* Convert to relative orbital index */
74
 
          I = i - occ_off[isym];
75
 
          J = j - occ_off[jsym];
76
 
          Fii = FMI.matrix[isym][I][I];
77
 
          Fjj = FMI.matrix[jsym][J][J];
78
 
 
79
 
          /* Loop over the columns */
80
 
          for(ab=0; ab < dIjAb.params->coltot[h^L_irr]; ab++) {
81
 
              a = dIjAb.params->colorb[h^L_irr][ab][0];
82
 
              b = dIjAb.params->colorb[h^L_irr][ab][1];
83
 
              asym = dIjAb.params->rsym[a];
84
 
              bsym = dIjAb.params->ssym[b];
85
 
 
86
 
              /* Convert to relative orbital index */
87
 
              A = a - vir_off[asym];
88
 
              B = b - vir_off[bsym];
89
 
 
90
 
              Faa = FAE.matrix[asym][A][A];
91
 
              Fbb = FAE.matrix[bsym][B][B];
92
 
 
93
 
              dIjAb.matrix[h][ij][ab] = 1.0/(Fii + Fjj - Faa - Fbb + L_params.cceom_energy);
94
 
            }
95
 
        }
96
 
    dpd_file4_mat_irrep_wrt(&dIjAb, h);
97
 
    dpd_file4_mat_irrep_close(&dIjAb, h);
98
 
  }
99
 
  dpd_file4_close(&dIjAb);
100
 
 
101
 
  dpd_file2_mat_close(&FMI);
102
 
  dpd_file2_mat_close(&FAE);
103
 
  dpd_file2_close(&FMI);
104
 
  dpd_file2_close(&FAE);
105
 
 
106
 
  return;
107
 
}
108
 
 
109
 
void denom_uhf(struct L_Params L_params)
110
 
{
111
 
  int nirreps, h, i, j, a, b, ij, ab, I, J, A, B, isym, jsym, asym, bsym, m, e;
112
 
  int *aoccpi, *boccpi, *avirtpi, *bvirtpi; 
113
 
  int *aocc_off, *bocc_off, *avir_off, *bvir_off, L_irr;
114
 
  dpdfile2 LFMIt, LFmit, LFaet, LFAEt;
115
 
  dpdfile2 FMI, Fmi, FAE, Fae;
116
 
  dpdfile2 dIA, dia;
117
 
  dpdfile4 dIJAB, dijab, dIjAb;
118
 
  double Fii, Fjj, Faa, Fbb;
119
 
 
120
 
  L_irr = L_params.irrep;
121
 
  nirreps = moinfo.nirreps;
122
 
  aoccpi = moinfo.aoccpi; 
123
 
  boccpi = moinfo.boccpi; 
124
 
  avirtpi = moinfo.avirtpi;
125
 
  bvirtpi = moinfo.bvirtpi;
126
 
  aocc_off = moinfo.aocc_off;
127
 
  bocc_off = moinfo.bocc_off;
128
 
  avir_off = moinfo.avir_off;
129
 
  bvir_off = moinfo.bvir_off;
130
 
 
131
 
  if((!strcmp(params.wfn,"CC2")) || (!strcmp(params.wfn,"EOM_CC2"))) {
132
 
 
133
 
    dpd_file2_init(&LFMIt, CC_OEI, 0, 0, 0, "fIJ");
134
 
    dpd_file2_mat_init(&LFMIt);
135
 
    dpd_file2_mat_rd(&LFMIt);
136
 
 
137
 
    dpd_file2_init(&LFmit, CC_OEI, 0, 2, 2, "fij");
138
 
    dpd_file2_mat_init(&LFmit);
139
 
    dpd_file2_mat_rd(&LFmit);
140
 
 
141
 
    dpd_file2_init(&LFaet, CC_OEI, 0, 3, 3, "fab");
142
 
    dpd_file2_mat_init(&LFaet);
143
 
    dpd_file2_mat_rd(&LFaet);
144
 
 
145
 
    dpd_file2_init(&LFAEt, CC_OEI, 0, 1, 1, "fAB");
146
 
    dpd_file2_mat_init(&LFAEt);
147
 
    dpd_file2_mat_rd(&LFAEt);
148
 
 
149
 
  }
150
 
  else {
151
 
    dpd_file2_init(&LFMIt, CC_OEI, 0, 0, 0, "FMI");
152
 
    dpd_file2_mat_init(&LFMIt);
153
 
    dpd_file2_mat_rd(&LFMIt);
154
 
 
155
 
    dpd_file2_init(&LFmit, CC_OEI, 0, 2, 2, "Fmi");
156
 
    dpd_file2_mat_init(&LFmit);
157
 
    dpd_file2_mat_rd(&LFmit);
158
 
 
159
 
    dpd_file2_init(&LFaet, CC_OEI, 0, 3, 3, "Fae");
160
 
    dpd_file2_mat_init(&LFaet);
161
 
    dpd_file2_mat_rd(&LFaet);
162
 
 
163
 
    dpd_file2_init(&LFAEt, CC_OEI, 0, 1, 1, "FAE");
164
 
    dpd_file2_mat_init(&LFAEt);
165
 
    dpd_file2_mat_rd(&LFAEt);
166
 
  }
167
 
 
168
 
  dpd_file2_init(&dIA, CC_DENOM, L_irr, 0, 1, "dIA");
169
 
  dpd_file2_mat_init(&dIA);
170
 
  for(h=0; h < nirreps; h++) {
171
 
    for(i=0; i < aoccpi[h]; i++) {
172
 
      Fii = LFMIt.matrix[h][i][i];
173
 
      for(a=0; a < avirtpi[h^L_irr]; a++) {
174
 
        Faa = LFAEt.matrix[h^L_irr][a][a];
175
 
        dIA.matrix[h][i][a] = 1.0/(Fii - Faa + L_params.cceom_energy);
176
 
      }
177
 
    }
178
 
  }
179
 
  dpd_file2_mat_wrt(&dIA);
180
 
  dpd_file2_mat_close(&dIA);
181
 
  dpd_file2_close(&dIA);
182
 
 
183
 
  dpd_file2_init(&dia, CC_DENOM, L_irr, 2, 3, "dia");
184
 
  dpd_file2_mat_init(&dia);
185
 
  for(h=0; h < nirreps; h++) {
186
 
    for(i=0; i < boccpi[h]; i++) {
187
 
      Fii = LFmit.matrix[h][i][i];
188
 
      for(a=0; a < bvirtpi[h^L_irr]; a++) {
189
 
        Faa = LFaet.matrix[h^L_irr][a][a];
190
 
        dia.matrix[h][i][a] = 1.0/(Fii - Faa + L_params.cceom_energy);
191
 
      }
192
 
    }
193
 
  }
194
 
  dpd_file2_mat_wrt(&dia);
195
 
  dpd_file2_mat_close(&dia);
196
 
  dpd_file2_close(&dia);
197
 
 
198
 
  dpd_file4_init(&dIJAB, CC_DENOM, L_irr, 1, 6, "dIJAB");
199
 
  for(h=0; h < nirreps; h++) {
200
 
    dpd_file4_mat_irrep_init(&dIJAB, h);
201
 
    for(ij=0; ij < dIJAB.params->rowtot[h]; ij++) {
202
 
      i = dIJAB.params->roworb[h][ij][0];
203
 
      j = dIJAB.params->roworb[h][ij][1];
204
 
      isym = dIJAB.params->psym[i];
205
 
      jsym = dIJAB.params->qsym[j];
206
 
      I = i - aocc_off[isym];
207
 
      J = j - aocc_off[jsym];
208
 
      Fii = LFMIt.matrix[isym][I][I];
209
 
      Fjj = LFMIt.matrix[jsym][J][J];
210
 
 
211
 
      for(ab=0; ab < dIJAB.params->coltot[h^L_irr]; ab++) {
212
 
        a = dIJAB.params->colorb[h^L_irr][ab][0];
213
 
        b = dIJAB.params->colorb[h^L_irr][ab][1];
214
 
        asym = dIJAB.params->rsym[a];
215
 
        bsym = dIJAB.params->ssym[b];
216
 
        A = a - avir_off[asym];
217
 
        B = b - avir_off[bsym];
218
 
        Faa = LFAEt.matrix[asym][A][A];
219
 
        Fbb = LFAEt.matrix[bsym][B][B];
220
 
 
221
 
        dIJAB.matrix[h][ij][ab] = 1.0/(Fii + Fjj - Faa - Fbb
222
 
                                       + L_params.cceom_energy);
223
 
      }
224
 
    }
225
 
    dpd_file4_mat_irrep_wrt(&dIJAB, h);
226
 
    dpd_file4_mat_irrep_close(&dIJAB, h);
227
 
  }
228
 
  dpd_file4_close(&dIJAB);
229
 
 
230
 
  dpd_file4_init(&dijab, CC_DENOM, L_irr, 11, 16, "dijab");
231
 
 
232
 
  for(h=0; h < nirreps; h++) {
233
 
    dpd_file4_mat_irrep_init(&dijab, h);
234
 
    for(ij=0; ij < dijab.params->rowtot[h]; ij++) {
235
 
      i = dijab.params->roworb[h][ij][0];
236
 
      j = dijab.params->roworb[h][ij][1];
237
 
      isym = dijab.params->psym[i];
238
 
      jsym = dijab.params->qsym[j];
239
 
      I = i - bocc_off[isym];
240
 
      J = j - bocc_off[jsym];
241
 
      Fii = LFmit.matrix[isym][I][I];
242
 
      Fjj = LFmit.matrix[jsym][J][J];
243
 
 
244
 
      for(ab=0; ab < dijab.params->coltot[h^L_irr]; ab++) {
245
 
        a = dijab.params->colorb[h^L_irr][ab][0];
246
 
        b = dijab.params->colorb[h^L_irr][ab][1];
247
 
        asym = dijab.params->rsym[a];
248
 
        bsym = dijab.params->ssym[b];
249
 
        A = a - bvir_off[asym];
250
 
        B = b - bvir_off[bsym];
251
 
        Faa = LFaet.matrix[asym][A][A];
252
 
        Fbb = LFaet.matrix[bsym][B][B];
253
 
 
254
 
        dijab.matrix[h][ij][ab] = 1.0/(Fii + Fjj - Faa - Fbb
255
 
                                       + L_params.cceom_energy);
256
 
      }
257
 
    }
258
 
    dpd_file4_mat_irrep_wrt(&dijab, h);
259
 
    dpd_file4_mat_irrep_close(&dijab, h);
260
 
  }
261
 
  dpd_file4_close(&dijab);
262
 
 
263
 
  dpd_file4_init(&dIjAb, CC_DENOM, L_irr, 22, 28, "dIjAb");
264
 
 
265
 
  for(h=0; h < nirreps; h++) {
266
 
    dpd_file4_mat_irrep_init(&dIjAb, h);
267
 
    for(ij=0; ij < dIjAb.params->rowtot[h]; ij++) {
268
 
      i = dIjAb.params->roworb[h][ij][0];
269
 
      j = dIjAb.params->roworb[h][ij][1];
270
 
      isym = dIjAb.params->psym[i];
271
 
      jsym = dIjAb.params->qsym[j];
272
 
      I = i - aocc_off[isym];
273
 
      J = j - bocc_off[jsym];
274
 
      Fii = LFMIt.matrix[isym][I][I];
275
 
      Fjj = LFmit.matrix[jsym][J][J];
276
 
 
277
 
      for(ab=0; ab < dIjAb.params->coltot[h^L_irr]; ab++) {
278
 
        a = dIjAb.params->colorb[h^L_irr][ab][0];
279
 
        b = dIjAb.params->colorb[h^L_irr][ab][1];
280
 
        asym = dIjAb.params->rsym[a];
281
 
        bsym = dIjAb.params->ssym[b];
282
 
        A = a - avir_off[asym];
283
 
        B = b - bvir_off[bsym];
284
 
        Faa = LFAEt.matrix[asym][A][A];
285
 
        Fbb = LFaet.matrix[bsym][B][B];
286
 
 
287
 
        dIjAb.matrix[h][ij][ab] = 1.0/(Fii + Fjj - Faa - Fbb
288
 
                                       + L_params.cceom_energy);
289
 
      }
290
 
    }
291
 
    dpd_file4_mat_irrep_wrt(&dIjAb, h);
292
 
    dpd_file4_mat_irrep_close(&dIjAb, h);
293
 
  }
294
 
  dpd_file4_close(&dIjAb);
295
 
 
296
 
  dpd_file2_mat_close(&LFMIt);
297
 
  dpd_file2_mat_close(&LFmit);
298
 
  dpd_file2_mat_close(&LFAEt);
299
 
  dpd_file2_mat_close(&LFaet);
300
 
  dpd_file2_close(&LFMIt);
301
 
  dpd_file2_close(&LFmit);
302
 
  dpd_file2_close(&LFAEt);
303
 
  dpd_file2_close(&LFaet);
304
 
 
305
 
  /*     if((!strcmp(params.wfn,"CC2")) || (!strcmp(params.wfn,"EOM_CC2"))) { */
306
 
  /*       dpd_file2_init(&FMI, CC_OEI, 0, 0, 0, "FMI"); */
307
 
  /*       dpd_file2_init(&Fmi, CC_OEI, 0, 2, 2, "Fmi"); */
308
 
 
309
 
  /*       dpd_file2_mat_init(&FMI); */
310
 
  /*       dpd_file2_mat_rd(&FMI); */
311
 
  /*       dpd_file2_mat_init(&Fmi); */
312
 
  /*       dpd_file2_mat_rd(&Fmi); */
313
 
 
314
 
  /*       for(h=0; h < moinfo.nirreps; h++) { */
315
 
  /*    for(m=0; m < FMI.params->rowtot[h]; m++)  */
316
 
  /*      FMI.matrix[h][m][m] = 0; */
317
 
  /*    for(m=0; m < Fmi.params->rowtot[h]; m++)  */
318
 
  /*      Fmi.matrix[h][m][m] = 0; */
319
 
  /*       } */
320
 
 
321
 
  /*       dpd_file2_mat_wrt(&FMI); */
322
 
  /*       dpd_file2_mat_close(&FMI); */
323
 
  /*       dpd_file2_mat_wrt(&Fmi); */
324
 
  /*       dpd_file2_mat_close(&Fmi); */
325
 
 
326
 
  /*       dpd_file2_close(&FMI); */
327
 
  /*       dpd_file2_close(&Fmi); */
328
 
 
329
 
  /*       dpd_file2_init(&FAE, CC_OEI, 0, 1, 1, "FAE"); */
330
 
  /*       dpd_file2_init(&Fae, CC_OEI, 0, 3, 3, "Fae"); */
331
 
 
332
 
  /*       dpd_file2_mat_init(&FAE); */
333
 
  /*       dpd_file2_mat_rd(&FAE); */
334
 
  /*       dpd_file2_mat_init(&Fae); */
335
 
  /*       dpd_file2_mat_rd(&Fae); */
336
 
  
337
 
  /*       for(h=0; h < moinfo.nirreps; h++) { */
338
 
  /*    for(e=0; e < FAE.params->coltot[h]; e++)  */
339
 
  /*      FAE.matrix[h][e][e] = 0; */
340
 
  /*    for(e=0; e < Fae.params->coltot[h]; e++) */
341
 
  /*      Fae.matrix[h][e][e] = 0; */
342
 
  /*       } */
343
 
 
344
 
  /*       dpd_file2_mat_wrt(&FAE); */
345
 
  /*       dpd_file2_mat_close(&FAE); */
346
 
  /*       dpd_file2_mat_wrt(&Fae); */
347
 
  /*       dpd_file2_mat_close(&Fae); */
348
 
 
349
 
  /*       dpd_file2_close(&FAE); */
350
 
  /*       dpd_file2_close(&Fae); */
351
 
  /*     } */
352
 
 
353
 
  return;
354
 
}
355
 
 
356
 
void denom_rohf(struct L_Params L_params)
357
 
{
358
 
  dpdfile2 LFAEt, LFaet, LFMIt, LFmit;
359
 
  dpdfile2 dIA, dia;
360
 
  dpdfile4 dIJAB, dijab, dIjAb;
361
 
  dpdbuf4 d, bdIJAB, bdijab, bdIjAb;
362
 
  double tval;
363
 
  int nirreps,L_irr;
364
 
  int h, i, j, a, b, ij, ab;
365
 
  int I, J, A, B;
366
 
  int isym, jsym, asym, bsym;
367
 
  int *occpi, *virtpi;
368
 
  int *occ_off, *vir_off;
369
 
  int *openpi;
370
 
  double Fii, Fjj, Faa, Fbb;
371
 
 
372
 
  L_irr = L_params.irrep;
373
 
  nirreps = moinfo.nirreps;
374
 
  occpi = moinfo.occpi; virtpi = moinfo.virtpi;
375
 
  openpi = moinfo.openpi;
376
 
  occ_off = moinfo.occ_off; vir_off = moinfo.vir_off;
377
 
 
378
 
  dpd_file2_init(&LFMIt, CC_OEI, 0, 0, 0, "FMI");
379
 
  dpd_file2_mat_init(&LFMIt);
380
 
  dpd_file2_mat_rd(&LFMIt);
381
 
 
382
 
  dpd_file2_init(&LFmit, CC_OEI, 0, 0, 0, "Fmi");
383
 
  dpd_file2_mat_init(&LFmit);
384
 
  dpd_file2_mat_rd(&LFmit);
385
 
 
386
 
  dpd_file2_init(&LFaet, CC_OEI, 0, 1, 1, "Fae");
387
 
  dpd_file2_mat_init(&LFaet);
388
 
  dpd_file2_mat_rd(&LFaet);
389
 
 
390
 
  dpd_file2_init(&LFAEt, CC_OEI, 0, 1, 1, "FAE");
391
 
  dpd_file2_mat_init(&LFAEt);
392
 
  dpd_file2_mat_rd(&LFAEt);
393
 
 
394
 
  /* Alpha one-electron denominator */
395
 
  dpd_file2_init(&dIA, CC_DENOM, L_irr, 0, 1, "dIA");
396
 
  dpd_file2_mat_init(&dIA);
397
 
  for(h=0; h < nirreps; h++) { /* irreps of dIA and Fii */
398
 
    for(i=0; i < occpi[h]; i++) {
399
 
      Fii = LFMIt.matrix[h][i][i];
400
 
      for(a=0; a < (virtpi[h^L_irr] - openpi[h^L_irr]); a++) {
401
 
        Faa = LFAEt.matrix[h^L_irr][a][a];
402
 
        dIA.matrix[h][i][a] = 1.0/(Fii - Faa + L_params.cceom_energy);
403
 
      }
404
 
    }
405
 
  }
406
 
  dpd_file2_mat_wrt(&dIA);
407
 
  dpd_file2_mat_close(&dIA);
408
 
  dpd_file2_close(&dIA);
409
 
 
410
 
  /* Beta one-electron denominator */
411
 
  dpd_file2_init(&dia, CC_DENOM, L_irr, 0, 1, "dia");
412
 
  dpd_file2_mat_init(&dia);
413
 
  for(h=0; h < nirreps; h++) {
414
 
    for(i=0; i < (occpi[h] - openpi[h]); i++) {
415
 
      Fii = LFmit.matrix[h][i][i];
416
 
      for(a=0; a < virtpi[h^L_irr]; a++) {
417
 
        Faa = LFaet.matrix[h^L_irr][a][a];
418
 
        dia.matrix[h][i][a] = 1.0/(Fii - Faa + L_params.cceom_energy);
419
 
      }
420
 
    }
421
 
  }
422
 
  dpd_file2_mat_wrt(&dia);
423
 
  dpd_file2_mat_close(&dia);
424
 
  dpd_file2_close(&dia);
425
 
 
426
 
  /* Alpha-alpha two-electron denominator */
427
 
  dpd_file4_init(&dIJAB, CC_DENOM, L_irr, 1, 6, "dIJAB");
428
 
 
429
 
  for(h=0; h < nirreps; h++) {
430
 
    dpd_file4_mat_irrep_init(&dIJAB, h);
431
 
      /* Loop over the rows */
432
 
      for(ij=0; ij < dIJAB.params->rowtot[h]; ij++) {
433
 
          i = dIJAB.params->roworb[h][ij][0];
434
 
          j = dIJAB.params->roworb[h][ij][1];
435
 
          isym = dIJAB.params->psym[i];
436
 
          jsym = dIJAB.params->qsym[j];
437
 
 
438
 
          /* Convert to relative orbital index */
439
 
          I = i - occ_off[isym];
440
 
          J = j - occ_off[jsym];
441
 
 
442
 
          Fii = LFMIt.matrix[isym][I][I];
443
 
          Fjj = LFMIt.matrix[jsym][J][J];
444
 
 
445
 
          /* Loop over the columns */
446
 
          for(ab=0; ab < dIJAB.params->coltot[h^L_irr]; ab++) {
447
 
              a = dIJAB.params->colorb[h^L_irr][ab][0];
448
 
              b = dIJAB.params->colorb[h^L_irr][ab][1];
449
 
              asym = dIJAB.params->rsym[a];
450
 
              bsym = dIJAB.params->ssym[b];
451
 
 
452
 
              /* Convert to relative orbital index */
453
 
              A = a - vir_off[asym];
454
 
              B = b - vir_off[bsym];
455
 
 
456
 
              Faa = LFAEt.matrix[asym][A][A];
457
 
              Fbb = LFAEt.matrix[bsym][B][B];
458
 
 
459
 
              dIJAB.matrix[h][ij][ab] =
460
 
                ((A >= (virtpi[asym] - openpi[asym])) ||
461
 
                 (B >= (virtpi[bsym] - openpi[bsym])) ?
462
 
                 0.0 : 1.0/(Fii + Fjj - Faa - Fbb
463
 
                   + L_params.cceom_energy));
464
 
          }
465
 
      }
466
 
    dpd_file4_mat_irrep_wrt(&dIJAB, h);
467
 
    dpd_file4_mat_irrep_close(&dIJAB, h);
468
 
  }
469
 
  dpd_file4_close(&dIJAB);
470
 
 
471
 
  /* Beta-beta two-electron denominator */
472
 
  dpd_file4_init(&dijab, CC_DENOM, L_irr, 1, 6, "dijab");
473
 
 
474
 
  for(h=0; h < nirreps; h++) {
475
 
    dpd_file4_mat_irrep_init(&dijab, h);
476
 
    /* Loop over the rows */
477
 
    for(ij=0; ij < dijab.params->rowtot[h]; ij++) {
478
 
          i = dijab.params->roworb[h][ij][0];
479
 
          j = dijab.params->roworb[h][ij][1];
480
 
          isym = dijab.params->psym[i];
481
 
          jsym = dijab.params->qsym[j];
482
 
 
483
 
          /* Convert to relative orbital index */
484
 
          I = i - occ_off[isym];
485
 
          J = j - occ_off[jsym];
486
 
 
487
 
          Fii = LFmit.matrix[isym][I][I];
488
 
          Fjj = LFmit.matrix[jsym][J][J];
489
 
 
490
 
          /* Loop over the columns */
491
 
          for(ab=0; ab < dijab.params->coltot[h^L_irr]; ab++) {
492
 
              a = dijab.params->colorb[h^L_irr][ab][0];
493
 
              b = dijab.params->colorb[h^L_irr][ab][1];
494
 
              asym = dijab.params->rsym[a];
495
 
              bsym = dijab.params->ssym[b];
496
 
 
497
 
              /* Convert to relative orbital index */
498
 
              A = a - vir_off[asym];
499
 
              B = b - vir_off[bsym];
500
 
 
501
 
              Faa = LFaet.matrix[asym][A][A];
502
 
              Fbb = LFaet.matrix[bsym][B][B];
503
 
 
504
 
              dijab.matrix[h][ij][ab] =
505
 
                ((I >= (occpi[isym] - openpi[isym])) ||
506
 
                 (J >= (occpi[jsym] - openpi[jsym])) ?
507
 
                 0.0 : 1.0/(Fii + Fjj - Faa - Fbb
508
 
                   + L_params.cceom_energy));
509
 
      }
510
 
    }
511
 
    dpd_file4_mat_irrep_wrt(&dijab, h);
512
 
    dpd_file4_mat_irrep_close(&dijab, h);
513
 
  }
514
 
  dpd_file4_close(&dijab);
515
 
 
516
 
 
517
 
  /* Alpha-beta two-electron denominator */
518
 
  dpd_file4_init(&dIjAb, CC_DENOM, L_irr, 0, 5, "dIjAb");
519
 
 
520
 
  for(h=0; h < nirreps; h++) {
521
 
    dpd_file4_mat_irrep_init(&dIjAb, h);
522
 
    /* Loop over the rows */
523
 
    for(ij=0; ij < dIjAb.params->rowtot[h]; ij++) {
524
 
          i = dIjAb.params->roworb[h][ij][0];
525
 
          j = dIjAb.params->roworb[h][ij][1];
526
 
          isym = dIjAb.params->psym[i];
527
 
          jsym = dIjAb.params->qsym[j];
528
 
 
529
 
          /* Convert to relative orbital index */
530
 
          I = i - occ_off[isym];
531
 
          J = j - occ_off[jsym];
532
 
          Fii = LFMIt.matrix[isym][I][I];
533
 
          Fjj = LFmit.matrix[jsym][J][J];
534
 
 
535
 
          /* Loop over the columns */
536
 
          for(ab=0; ab < dIjAb.params->coltot[h^L_irr]; ab++) {
537
 
              a = dIjAb.params->colorb[h^L_irr][ab][0];
538
 
              b = dIjAb.params->colorb[h^L_irr][ab][1];
539
 
              asym = dIjAb.params->rsym[a];
540
 
              bsym = dIjAb.params->ssym[b];
541
 
 
542
 
              /* Convert to relative orbital index */
543
 
              A = a - vir_off[asym];
544
 
              B = b - vir_off[bsym];
545
 
 
546
 
              Faa = LFAEt.matrix[asym][A][A];
547
 
              Fbb = LFaet.matrix[bsym][B][B];
548
 
 
549
 
              dIjAb.matrix[h][ij][ab] =
550
 
                ((A >= (virtpi[asym] - openpi[asym])) ||
551
 
                 (J >= (occpi[jsym] - openpi[jsym])) ?
552
 
                 0.0 : 1.0/(Fii + Fjj - Faa - Fbb
553
 
                   + L_params.cceom_energy));
554
 
            }
555
 
        }
556
 
    dpd_file4_mat_irrep_wrt(&dIjAb, h);
557
 
    dpd_file4_mat_irrep_close(&dIjAb, h);
558
 
  }
559
 
  dpd_file4_close(&dIjAb);
560
 
 
561
 
  dpd_file2_mat_close(&LFMIt);
562
 
  dpd_file2_mat_close(&LFmit);
563
 
  dpd_file2_mat_close(&LFAEt);
564
 
  dpd_file2_mat_close(&LFaet);
565
 
  dpd_file2_close(&LFMIt);
566
 
  dpd_file2_close(&LFmit);
567
 
  dpd_file2_close(&LFAEt);
568
 
  dpd_file2_close(&LFaet);
569
 
 
570
 
  return;
571
 
}
572
 
 
573
 
 
574