← Back to branch summary

~ubuntu-branches/ubuntu/precise/psicode/precise

« back to all changes in this revision

Viewing changes to src/bin/ccsort/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/ccsort/denom.c
1
 
#include <stdio.h>
2
 
#include <stdlib.h>
3
 
#include <libdpd/dpd.h>
4
 
#include "MOInfo.h"
5
 
#include "Params.h"
6
 
#define EXTERN
7
 
#include "globals.h"
8
 
 
9
 
void denom_rhf(void);
10
 
void denom_uhf(void);
11
 
 
12
 
void denom(void)
13
 
{
14
 
  if(params.ref == 2) denom_uhf();
15
 
  else denom_rhf();
16
 
}
17
 
 
18
 
void denom_uhf(void)
19
 
{
20
 
  int h, nirreps;
21
 
  int i, j, a, b, ij, ab, I, J, A, B, isym, jsym, asym, bsym;
22
 
  int *aoccpi, *avirtpi, *boccpi, *bvirtpi;
23
 
  int *aocc_off, *bocc_off, *avir_off, *bvir_off;
24
 
  double fii, fjj, faa, fbb;
25
 
  dpdfile2 dIA, fIJ, fij, fAB, fab;
26
 
  dpdfile4 dIJAB;
27
 
 
28
 
  nirreps = moinfo.nirreps;
29
 
  aoccpi = moinfo.aoccpi;
30
 
  avirtpi = moinfo.avirtpi;
31
 
  boccpi = moinfo.boccpi;
32
 
  bvirtpi = moinfo.bvirtpi;
33
 
  aocc_off = moinfo.aocc_off;
34
 
  bocc_off = moinfo.bocc_off;
35
 
  avir_off = moinfo.avir_off;
36
 
  bvir_off = moinfo.bvir_off;
37
 
 
38
 
  dpd_file2_init(&fIJ, CC_OEI, 0, 0, 0, "fIJ");
39
 
  dpd_file2_mat_init(&fIJ);
40
 
  dpd_file2_mat_rd(&fIJ);
41
 
  
42
 
  dpd_file2_init(&fij, CC_OEI, 0, 2, 2, "fij");
43
 
  dpd_file2_mat_init(&fij);
44
 
  dpd_file2_mat_rd(&fij);
45
 
  
46
 
  dpd_file2_init(&fAB, CC_OEI, 0, 1, 1, "fAB");
47
 
  dpd_file2_mat_init(&fAB);
48
 
  dpd_file2_mat_rd(&fAB);
49
 
  
50
 
  dpd_file2_init(&fab, CC_OEI, 0, 3, 3, "fab");
51
 
  dpd_file2_mat_init(&fab);
52
 
  dpd_file2_mat_rd(&fab);
53
 
 
54
 
  dpd_file2_init(&dIA, CC_OEI, 0, 0, 1, "dIA");
55
 
  dpd_file2_mat_init(&dIA);
56
 
  for(h=0; h < nirreps; h++) {
57
 
      for(i=0; i < aoccpi[h]; i++) {
58
 
          fii = fIJ.matrix[h][i][i];
59
 
          for(a=0; a < avirtpi[h]; a++) {
60
 
              faa = fAB.matrix[h][a][a];
61
 
              dIA.matrix[h][i][a] = 1.0/(fii - faa);
62
 
            }
63
 
        }
64
 
    }
65
 
  dpd_file2_mat_wrt(&dIA);
66
 
  dpd_file2_mat_close(&dIA);
67
 
  dpd_file2_close(&dIA);
68
 
 
69
 
  dpd_file2_init(&dIA, CC_OEI, 0, 2, 3, "dia");
70
 
  dpd_file2_mat_init(&dIA);
71
 
  for(h=0; h < nirreps; h++) {
72
 
      for(i=0; i < boccpi[h]; i++) {
73
 
          fii = fij.matrix[h][i][i];
74
 
          for(a=0; a < bvirtpi[h]; a++) {
75
 
              faa = fab.matrix[h][a][a];
76
 
              dIA.matrix[h][i][a] = 1.0/(fii - faa);
77
 
            }
78
 
        }
79
 
    }
80
 
  dpd_file2_mat_wrt(&dIA);
81
 
  dpd_file2_mat_close(&dIA);
82
 
  dpd_file2_close(&dIA);
83
 
 
84
 
  dpd_file4_init(&dIJAB, CC_DENOM, 0, 1, 6, "dIJAB");
85
 
  for(h=0; h < nirreps; h++) {
86
 
      dpd_file4_mat_irrep_init(&dIJAB, h);
87
 
      for(ij=0; ij < dIJAB.params->rowtot[h]; ij++) {
88
 
          i = dIJAB.params->roworb[h][ij][0];
89
 
          j = dIJAB.params->roworb[h][ij][1];
90
 
          isym = dIJAB.params->psym[i];
91
 
          jsym = dIJAB.params->qsym[j];
92
 
          I = i - aocc_off[isym];
93
 
          J = j - aocc_off[jsym];
94
 
          fii = fIJ.matrix[isym][I][I];
95
 
          fjj = fIJ.matrix[jsym][J][J];
96
 
 
97
 
          for(ab=0; ab < dIJAB.params->coltot[h]; ab++) {
98
 
              a = dIJAB.params->colorb[h][ab][0];
99
 
              b = dIJAB.params->colorb[h][ab][1];
100
 
              asym = dIJAB.params->rsym[a];
101
 
              bsym = dIJAB.params->ssym[b];
102
 
              A = a - avir_off[asym];
103
 
              B = b - avir_off[bsym];
104
 
              faa = fAB.matrix[asym][A][A];
105
 
              fbb = fAB.matrix[bsym][B][B];
106
 
 
107
 
              dIJAB.matrix[h][ij][ab] = 1.0/(fii + fjj - faa - fbb);
108
 
            }
109
 
        }
110
 
      dpd_file4_mat_irrep_wrt(&dIJAB, h);
111
 
      dpd_file4_mat_irrep_close(&dIJAB, h);
112
 
    }
113
 
  dpd_file4_close(&dIJAB);
114
 
 
115
 
  dpd_file4_init(&dIJAB, CC_DENOM, 0, 11, 16, "dijab");
116
 
  for(h=0; h < nirreps; h++) {
117
 
      dpd_file4_mat_irrep_init(&dIJAB, h);
118
 
      for(ij=0; ij < dIJAB.params->rowtot[h]; ij++) {
119
 
          i = dIJAB.params->roworb[h][ij][0];
120
 
          j = dIJAB.params->roworb[h][ij][1];
121
 
          isym = dIJAB.params->psym[i];
122
 
          jsym = dIJAB.params->qsym[j];
123
 
          I = i - bocc_off[isym];
124
 
          J = j - bocc_off[jsym];
125
 
          fii = fij.matrix[isym][I][I];
126
 
          fjj = fij.matrix[jsym][J][J];
127
 
 
128
 
          for(ab=0; ab < dIJAB.params->coltot[h]; ab++) {
129
 
              a = dIJAB.params->colorb[h][ab][0];
130
 
              b = dIJAB.params->colorb[h][ab][1];
131
 
              asym = dIJAB.params->rsym[a];
132
 
              bsym = dIJAB.params->ssym[b];
133
 
              A = a - bvir_off[asym];
134
 
              B = b - bvir_off[bsym];
135
 
              faa = fab.matrix[asym][A][A];
136
 
              fbb = fab.matrix[bsym][B][B];
137
 
 
138
 
              dIJAB.matrix[h][ij][ab] = 1.0/(fii + fjj - faa - fbb);
139
 
            }
140
 
        }
141
 
      dpd_file4_mat_irrep_wrt(&dIJAB, h);
142
 
      dpd_file4_mat_irrep_close(&dIJAB, h);
143
 
    }
144
 
  dpd_file4_close(&dIJAB);
145
 
 
146
 
  dpd_file4_init(&dIJAB, CC_DENOM, 0, 22, 28, "dIjAb");
147
 
  for(h=0; h < nirreps; h++) {
148
 
      dpd_file4_mat_irrep_init(&dIJAB, h);
149
 
      for(ij=0; ij < dIJAB.params->rowtot[h]; ij++) {
150
 
          i = dIJAB.params->roworb[h][ij][0];
151
 
          j = dIJAB.params->roworb[h][ij][1];
152
 
          isym = dIJAB.params->psym[i];
153
 
          jsym = dIJAB.params->qsym[j];
154
 
          I = i - aocc_off[isym];
155
 
          J = j - bocc_off[jsym];
156
 
          fii = fIJ.matrix[isym][I][I];
157
 
          fjj = fij.matrix[jsym][J][J];
158
 
 
159
 
          for(ab=0; ab < dIJAB.params->coltot[h]; ab++) {
160
 
              a = dIJAB.params->colorb[h][ab][0];
161
 
              b = dIJAB.params->colorb[h][ab][1];
162
 
              asym = dIJAB.params->rsym[a];
163
 
              bsym = dIJAB.params->ssym[b];
164
 
              A = a - avir_off[asym];
165
 
              B = b - bvir_off[bsym];
166
 
              faa = fAB.matrix[asym][A][A];
167
 
              fbb = fab.matrix[bsym][B][B];
168
 
 
169
 
              dIJAB.matrix[h][ij][ab] = 1.0/(fii + fjj - faa - fbb);
170
 
            }
171
 
        }
172
 
      dpd_file4_mat_irrep_wrt(&dIJAB, h);
173
 
      dpd_file4_mat_irrep_close(&dIJAB, h);
174
 
    }
175
 
  dpd_file4_close(&dIJAB);
176
 
 
177
 
  dpd_file2_mat_close(&fIJ);
178
 
  dpd_file2_mat_close(&fij);
179
 
  dpd_file2_mat_close(&fAB);
180
 
  dpd_file2_mat_close(&fab);
181
 
  dpd_file2_close(&fIJ);
182
 
  dpd_file2_close(&fij);
183
 
  dpd_file2_close(&fAB);
184
 
  dpd_file2_close(&fab);
185
 
}
186
 
 
187
 
void denom_rhf(void)
188
 
{
189
 
  int nirreps;
190
 
  int h, i, j, a, b, ij, ab;
191
 
  int I, J, A, B;
192
 
  int isym, jsym, asym, bsym;
193
 
  int *occpi, *virtpi;
194
 
  int *occ_off, *vir_off;
195
 
  int *openpi;
196
 
  double fii, fjj, faa, fbb;
197
 
  dpdfile2 fIJ, fij, fAB, fab;
198
 
  dpdfile2 dIA, dia;
199
 
  dpdfile4 dIJAB, dijab, dIjAb;
200
 
 
201
 
  nirreps = moinfo.nirreps;
202
 
  occpi = moinfo.occpi; virtpi = moinfo.virtpi;
203
 
  openpi = moinfo.openpi;
204
 
  occ_off = moinfo.occ_off; vir_off = moinfo.vir_off;
205
 
 
206
 
  /* Grab Fock matrices from disk */
207
 
  dpd_file2_init(&fIJ, CC_OEI, 0, 0, 0, "fIJ");
208
 
  dpd_file2_mat_init(&fIJ);
209
 
  dpd_file2_mat_rd(&fIJ);
210
 
  
211
 
  dpd_file2_init(&fij, CC_OEI, 0, 0, 0, "fij");
212
 
  dpd_file2_mat_init(&fij);
213
 
  dpd_file2_mat_rd(&fij);
214
 
  
215
 
  dpd_file2_init(&fAB, CC_OEI, 0, 1, 1, "fAB");
216
 
  dpd_file2_mat_init(&fAB);
217
 
  dpd_file2_mat_rd(&fAB);
218
 
  
219
 
  dpd_file2_init(&fab, CC_OEI, 0, 1, 1, "fab");
220
 
  dpd_file2_mat_init(&fab);
221
 
  dpd_file2_mat_rd(&fab);
222
 
 
223
 
  /* Alpha one-electron denominator */
224
 
  dpd_file2_init(&dIA, CC_OEI, 0, 0, 1, "dIA");
225
 
  dpd_file2_mat_init(&dIA);
226
 
 
227
 
  for(h=0; h < nirreps; h++) {
228
 
      
229
 
      for(i=0; i < occpi[h]; i++) {
230
 
          fii = fIJ.matrix[h][i][i];
231
 
 
232
 
          for(a=0; a < (virtpi[h] - openpi[h]); a++) {
233
 
              faa = fAB.matrix[h][a][a];
234
 
 
235
 
              dIA.matrix[h][i][a] = 1.0/(fii - faa);
236
 
            }
237
 
        }
238
 
      
239
 
    }
240
 
 
241
 
  dpd_file2_mat_wrt(&dIA);
242
 
  dpd_file2_mat_close(&dIA);
243
 
  dpd_file2_close(&dIA);
244
 
 
245
 
  /* Beta one-electron denominator */
246
 
  dpd_file2_init(&dia, CC_OEI, 0, 0, 1, "dia");
247
 
  dpd_file2_mat_init(&dia);
248
 
 
249
 
  for(h=0; h < nirreps; h++) {
250
 
      
251
 
      for(i=0; i < (occpi[h] - openpi[h]); i++) {
252
 
          fii = fij.matrix[h][i][i];
253
 
 
254
 
          for(a=0; a < virtpi[h]; a++) {
255
 
              faa = fab.matrix[h][a][a];
256
 
 
257
 
              dia.matrix[h][i][a] = 1.0/(fii - faa);
258
 
            }
259
 
        }
260
 
      
261
 
    }
262
 
  
263
 
  dpd_file2_mat_wrt(&dia);
264
 
  dpd_file2_mat_close(&dia);
265
 
  dpd_file2_close(&dia);
266
 
 
267
 
  /* Alpha-alpha two-electron denominator */
268
 
  dpd_file4_init(&dIJAB, CC_DENOM, 0, 1, 6, "dIJAB");
269
 
 
270
 
  for(h=0; h < nirreps; h++) {
271
 
 
272
 
      dpd_file4_mat_irrep_init(&dIJAB, h);
273
 
 
274
 
      /* Loop over the rows */
275
 
      for(ij=0; ij < dIJAB.params->rowtot[h]; ij++) {
276
 
          i = dIJAB.params->roworb[h][ij][0];
277
 
          j = dIJAB.params->roworb[h][ij][1];
278
 
          isym = dIJAB.params->psym[i];
279
 
          jsym = dIJAB.params->qsym[j];
280
 
 
281
 
          /* Convert to relative orbital index */
282
 
          I = i - occ_off[isym];
283
 
          J = j - occ_off[jsym];
284
 
 
285
 
          fii = fIJ.matrix[isym][I][I];
286
 
          fjj = fIJ.matrix[jsym][J][J];
287
 
          
288
 
          /* Loop over the columns */
289
 
          for(ab=0; ab < dIJAB.params->coltot[h]; ab++) {
290
 
              a = dIJAB.params->colorb[h][ab][0];
291
 
              b = dIJAB.params->colorb[h][ab][1];
292
 
              asym = dIJAB.params->rsym[a];
293
 
              bsym = dIJAB.params->ssym[b];
294
 
 
295
 
              /* Convert to relative orbital index */
296
 
              A = a - vir_off[asym];
297
 
              B = b - vir_off[bsym];
298
 
 
299
 
              faa = fAB.matrix[asym][A][A];
300
 
              fbb = fAB.matrix[bsym][B][B];
301
 
 
302
 
              dIJAB.matrix[h][ij][ab] =
303
 
                ((A >= (virtpi[asym] - openpi[asym])) ||
304
 
                 (B >= (virtpi[bsym] - openpi[bsym])) ?
305
 
                 0.0 : 1.0/(fii + fjj - faa - fbb));
306
 
            }
307
 
        }
308
 
 
309
 
      dpd_file4_mat_irrep_wrt(&dIJAB, h);
310
 
      dpd_file4_mat_irrep_close(&dIJAB, h);
311
 
 
312
 
    }
313
 
 
314
 
  dpd_file4_close(&dIJAB);
315
 
 
316
 
  /* Beta-beta two-electron denominator */
317
 
  dpd_file4_init(&dijab, CC_DENOM, 0, 1, 6, "dijab");
318
 
 
319
 
  for(h=0; h < nirreps; h++) {
320
 
 
321
 
      dpd_file4_mat_irrep_init(&dijab, h);
322
 
 
323
 
      /* Loop over the rows */
324
 
      for(ij=0; ij < dijab.params->rowtot[h]; ij++) {
325
 
          i = dijab.params->roworb[h][ij][0];
326
 
          j = dijab.params->roworb[h][ij][1];
327
 
          isym = dijab.params->psym[i];
328
 
          jsym = dijab.params->qsym[j];
329
 
 
330
 
          /* Convert to relative orbital index */
331
 
          I = i - occ_off[isym];
332
 
          J = j - occ_off[jsym];
333
 
 
334
 
          fii = fij.matrix[isym][I][I];
335
 
          fjj = fij.matrix[jsym][J][J];
336
 
          
337
 
          /* Loop over the columns */
338
 
          for(ab=0; ab < dijab.params->coltot[h]; ab++) {
339
 
              a = dijab.params->colorb[h][ab][0];
340
 
              b = dijab.params->colorb[h][ab][1];
341
 
              asym = dijab.params->rsym[a];
342
 
              bsym = dijab.params->ssym[b];
343
 
 
344
 
              /* Convert to relative orbital index */
345
 
              A = a - vir_off[asym];
346
 
              B = b - vir_off[bsym];
347
 
 
348
 
              faa = fab.matrix[asym][A][A];
349
 
              fbb = fab.matrix[bsym][B][B];
350
 
 
351
 
              dijab.matrix[h][ij][ab] =
352
 
                ((I >= (occpi[isym] - openpi[isym])) ||
353
 
                 (J >= (occpi[jsym] - openpi[jsym])) ?
354
 
                 0.0 : 1.0/(fii + fjj - faa - fbb));
355
 
            }
356
 
        }
357
 
 
358
 
      dpd_file4_mat_irrep_wrt(&dijab, h);
359
 
      dpd_file4_mat_irrep_close(&dijab, h);
360
 
 
361
 
    }
362
 
 
363
 
  dpd_file4_close(&dijab);
364
 
 
365
 
  /* Alpha-beta two-electron denominator */
366
 
  dpd_file4_init(&dIjAb, CC_DENOM, 0, 0, 5, "dIjAb");
367
 
 
368
 
  for(h=0; h < nirreps; h++) {
369
 
 
370
 
      dpd_file4_mat_irrep_init(&dIjAb, h);
371
 
 
372
 
      /* Loop over the rows */
373
 
      for(ij=0; ij < dIjAb.params->rowtot[h]; ij++) {
374
 
          i = dIjAb.params->roworb[h][ij][0];
375
 
          j = dIjAb.params->roworb[h][ij][1];
376
 
          isym = dIjAb.params->psym[i];
377
 
          jsym = dIjAb.params->qsym[j];
378
 
 
379
 
          /* Convert to relative orbital index */
380
 
          I = i - occ_off[isym];
381
 
          J = j - occ_off[jsym];
382
 
 
383
 
          fii = fIJ.matrix[isym][I][I];
384
 
          fjj = fij.matrix[jsym][J][J];
385
 
          
386
 
          /* Loop over the columns */
387
 
          for(ab=0; ab < dIjAb.params->coltot[h]; ab++) {
388
 
              a = dIjAb.params->colorb[h][ab][0];
389
 
              b = dIjAb.params->colorb[h][ab][1];
390
 
              asym = dIjAb.params->rsym[a];
391
 
              bsym = dIjAb.params->ssym[b];
392
 
 
393
 
              /* Convert to relative orbital index */
394
 
              A = a - vir_off[asym];
395
 
              B = b - vir_off[bsym];
396
 
 
397
 
              faa = fAB.matrix[asym][A][A];
398
 
              fbb = fab.matrix[bsym][B][B];
399
 
 
400
 
              dIjAb.matrix[h][ij][ab] =
401
 
                ((A >= (virtpi[asym] - openpi[asym])) ||
402
 
                 (J >= (occpi[jsym] - openpi[jsym])) ?
403
 
                 0.0 : 1.0/(fii + fjj - faa - fbb));
404
 
            }
405
 
        }
406
 
 
407
 
      dpd_file4_mat_irrep_wrt(&dIjAb, h);
408
 
      dpd_file4_mat_irrep_close(&dIjAb, h);
409
 
 
410
 
    }
411
 
 
412
 
  dpd_file4_close(&dIjAb);
413
 
 
414
 
  dpd_file2_mat_close(&fIJ);
415
 
  dpd_file2_mat_close(&fij);
416
 
  dpd_file2_mat_close(&fAB);
417
 
  dpd_file2_mat_close(&fab);
418
 
  dpd_file2_close(&fIJ);
419
 
  dpd_file2_close(&fij);
420
 
  dpd_file2_close(&fAB);
421
 
  dpd_file2_close(&fab);
422
 
 
423
 
}