← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/ccenergy/fock_build.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/ccenergy/fock_build.cc
 
1
/*! \file
 
2
    \ingroup CCENERGY
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
#include <cstdio>
 
6
#include <cstdlib>
 
7
#include <cmath>
 
8
#include <libciomr/libciomr.h>
 
9
#include <libiwl/iwl.h>
 
10
#include <libqt/qt.h>
 
11
#include <psifiles.h>
 
12
#include "MOInfo.h"
 
13
#define EXTERN
 
14
#include "globals.h"
 
15
 
 
16
namespace psi { namespace ccenergy {
 
17
 
 
18
#define INDEX(i,j) ((i>j) ? (ioff[(i)]+(j)) : (ioff[(j)]+(i)))
 
19
 
 
20
void rhf_fock_build(double **fock, double  **D)
 
21
{
 
22
  int i, j, ij;
 
23
  int nso, ntri, stat;
 
24
  double *scratch;
 
25
  int lastbuf, idx, p, q, r, s, pq, rs;
 
26
  double value;
 
27
  Value *valptr;
 
28
  Label *lblptr;
 
29
  struct iwlbuf InBuf;
 
30
 
 
31
  nso = moinfo.nso;
 
32
  ntri = nso * (nso+1)/2;
 
33
 
 
34
  /* one-electron contributions */
 
35
  scratch = init_array(ntri);
 
36
  stat = iwl_rdone(PSIF_OEI, PSIF_SO_T, scratch, ntri, 0, 0, outfile);
 
37
  for(i=0, ij=0; i < nso; i++)
 
38
    for(j=0; j <= i; j++, ij++)
 
39
      fock[i][j] = fock[j][i] = scratch[ij];
 
40
  stat = iwl_rdone(PSIF_OEI, PSIF_SO_V, scratch, ntri, 0, 0, outfile);
 
41
  for(i=0, ij=0; i < nso; i++)
 
42
    for(j=0; j <= i; j++, ij++) {
 
43
      fock[i][j] += scratch[ij];
 
44
      if(i!=j) fock[j][i] += scratch[ij];
 
45
    }
 
46
  free(scratch);
 
47
 
 
48
  iwl_buf_init(&InBuf, PSIF_SO_TEI, 0.0, 1, 1);
 
49
  do {
 
50
 
 
51
    lastbuf = InBuf.lastbuf;
 
52
    lblptr = InBuf.labels;
 
53
    valptr = InBuf.values;
 
54
 
 
55
    for(idx=4*InBuf.idx; InBuf.idx < InBuf.inbuf; InBuf.idx++) {
 
56
      p = abs((int) lblptr[idx++]);
 
57
      q = (int) lblptr[idx++];
 
58
      r = (int) lblptr[idx++];
 
59
      s = (int) lblptr[idx++];
 
60
      value = (double) valptr[InBuf.idx];
 
61
 
 
62
      pq = INDEX(p,q);
 
63
      rs = INDEX(r,s);
 
64
 
 
65
      /*
 
66
      fprintf(outfile, "%d %d %d %d [%d] [%d] %20.15f\n", p, q, r, s, pq, rs, value);
 
67
      */
 
68
 
 
69
      /* (pq|rs) */
 
70
      fock[p][q] += 2.0 * D[r][s] * value;
 
71
      fock[p][r] -= D[q][s] * value;
 
72
 
 
73
      if(p!=q && r!=s && pq != rs) {
 
74
 
 
75
        /* (pq|sr) */
 
76
        fock[p][q] += 2.0 * D[s][r] * value;
 
77
        fock[p][s] -= D[q][r] * value;
 
78
 
 
79
        /* (qp|rs) */
 
80
        fock[q][p] += 2.0 * D[r][s] * value;
 
81
        fock[q][r] -= D[p][s] * value;
 
82
 
 
83
        /* (qp|sr) */
 
84
        fock[q][p] += 2.0 * D[s][r] * value;
 
85
        fock[q][s] -= D[p][r] * value;
 
86
 
 
87
        /* (rs|pq) */
 
88
        fock[r][s] += 2.0 * D[p][q] * value;
 
89
        fock[r][p] -= D[s][q] * value;
 
90
 
 
91
        /* (rs|qp) */
 
92
        fock[r][s] += 2.0 * D[q][p] * value;
 
93
        fock[r][q] -= D[s][p] * value;
 
94
 
 
95
        /* (sr|pq) */
 
96
        fock[s][r] += 2.0 * D[p][q] * value;
 
97
        fock[s][p] -= D[r][q] * value;
 
98
 
 
99
        /* (sr|qp) */
 
100
        fock[s][r] += 2.0 * D[q][p] * value;
 
101
        fock[s][q] -= D[r][p] * value;
 
102
      }
 
103
      else if(p!=q && r!=s && pq==rs) {
 
104
 
 
105
        /* (pq|sr) */
 
106
        fock[p][q] += 2.0 * D[s][r] * value;
 
107
        fock[p][s] -= D[q][r] * value;
 
108
 
 
109
        /* (qp|rs) */
 
110
        fock[q][p] += 2.0 * D[r][s] * value;
 
111
        fock[q][r] -= D[p][s] * value;
 
112
 
 
113
        /* (qp|sr) */
 
114
        fock[q][p] += 2.0 * D[s][r] * value;
 
115
        fock[q][s] -= D[p][r] * value;
 
116
 
 
117
      }
 
118
      else if(p!=q && r==s) {
 
119
 
 
120
        /* (qp|rs) */
 
121
        fock[q][p] += 2.0 * D[r][s] * value;
 
122
        fock[q][r] -= D[p][s] * value;
 
123
 
 
124
        /* (rs|pq) */
 
125
        fock[r][s] += 2.0 * D[p][q] * value;
 
126
        fock[r][p] -= D[s][q] * value;
 
127
 
 
128
        /* (rs|qp) */
 
129
        fock[r][s] += 2.0 * D[q][p] * value;
 
130
        fock[r][q] -= D[s][p] * value;
 
131
 
 
132
      }
 
133
      else if(p==q && r!=s) {
 
134
 
 
135
        /* (pq|sr) */
 
136
        fock[p][q] += 2.0 * D[s][r] * value;
 
137
        fock[p][s] -= D[q][r] * value;
 
138
 
 
139
        /* (rs|pq) */
 
140
        fock[r][s] += 2.0 * D[p][q] * value;
 
141
        fock[r][p] -= D[s][q] * value;
 
142
 
 
143
        /* (sr|pq) */
 
144
        fock[s][r] += 2.0 * D[p][q] * value;
 
145
        fock[s][p] -= D[r][q] * value;
 
146
 
 
147
      }
 
148
      else if(p==q && r==s && pq!=rs) {
 
149
 
 
150
        /* (rs|pq) */
 
151
        fock[r][s] += 2.0 * D[p][q] * value;
 
152
        fock[r][p] -= D[s][q] * value;
 
153
 
 
154
      }
 
155
    }
 
156
 
 
157
    if(!lastbuf) iwl_buf_fetch(&InBuf);
 
158
 
 
159
  } while (!lastbuf);
 
160
  iwl_buf_close(&InBuf, 1);
 
161
}
 
162
 
 
163
void uhf_fock_build(double **fock_a, double **fock_b, double **D_a, double **D_b)
 
164
{
 
165
  int i, j, ij;
 
166
  int nso, ntri, ntei, stat;
 
167
  double *scratch;
 
168
  int lastbuf, idx, p, q, r, s, pq, rs;
 
169
  double value;
 
170
  Value *valptr;
 
171
  Label *lblptr;
 
172
  struct iwlbuf InBuf;
 
173
  double **Dt;
 
174
 
 
175
  nso = moinfo.nso;
 
176
  ntri = nso * (nso+1)/2;
 
177
  ntei = ntri * (ntri+1)/2;
 
178
 
 
179
  Dt = block_matrix(nso, nso);
 
180
  for(p=0; p < nso; p++)
 
181
    for(q=0; q < nso; q++)
 
182
      Dt[p][q] = D_a[p][q] + D_b[p][q];
 
183
 
 
184
  /* one-electron contributions */
 
185
  scratch = init_array(ntri);
 
186
  stat = iwl_rdone(PSIF_OEI, PSIF_SO_T, scratch, ntri, 0, 0, outfile);
 
187
  for(i=0, ij=0; i < nso; i++)
 
188
    for(j=0; j <= i; j++, ij++) {
 
189
      fock_a[i][j] = fock_a[j][i] = scratch[ij];
 
190
      fock_b[i][j] = fock_b[j][i] = scratch[ij];
 
191
    }
 
192
  stat = iwl_rdone(PSIF_OEI, PSIF_SO_V, scratch, ntri, 0, 0, outfile);
 
193
  for(i=0, ij=0; i < nso; i++)
 
194
    for(j=0; j <= i; j++, ij++) {
 
195
      fock_a[i][j] += scratch[ij];
 
196
      if(i!=j) fock_a[j][i] += scratch[ij];
 
197
      fock_b[i][j] += scratch[ij];
 
198
      if(i!=j) fock_b[j][i] += scratch[ij];
 
199
    }
 
200
  free(scratch);
 
201
 
 
202
  iwl_buf_init(&InBuf, PSIF_SO_TEI, 0.0, 1, 1);
 
203
  do {
 
204
 
 
205
    lastbuf = InBuf.lastbuf;
 
206
    lblptr = InBuf.labels;
 
207
    valptr = InBuf.values;
 
208
 
 
209
    for(idx=4*InBuf.idx; InBuf.idx < InBuf.inbuf; InBuf.idx++) {
 
210
      p = abs((int) lblptr[idx++]);
 
211
      q = (int) lblptr[idx++];
 
212
      r = (int) lblptr[idx++];
 
213
      s = (int) lblptr[idx++];
 
214
      value = (double) valptr[InBuf.idx];
 
215
 
 
216
      pq = INDEX(p,q);
 
217
      rs = INDEX(r,s);
 
218
 
 
219
      /*
 
220
      fprintf(outfile, "%d %d %d %d [%d] [%d] %20.15f\n", p, q, r, s, pq, rs, value);
 
221
      */
 
222
 
 
223
      /* (pq|rs) */
 
224
      fock_a[p][q] += Dt[r][s] * value;
 
225
      fock_a[p][r] -= D_a[q][s] * value;
 
226
      fock_b[p][q] += Dt[r][s] * value;
 
227
      fock_b[p][r] -= D_b[q][s] * value;
 
228
 
 
229
      if(p!=q && r!=s && pq != rs) {
 
230
 
 
231
        /* (pq|sr) */
 
232
        fock_a[p][q] += Dt[s][r] * value;
 
233
        fock_a[p][s] -= D_a[q][r] * value;
 
234
        fock_b[p][q] += Dt[s][r] * value;
 
235
        fock_b[p][s] -= D_b[q][r] * value;
 
236
 
 
237
        /* (qp|rs) */
 
238
        fock_a[q][p] += Dt[r][s] * value;
 
239
        fock_a[q][r] -= D_a[p][s] * value;
 
240
        fock_b[q][p] += Dt[r][s] * value;
 
241
        fock_b[q][r] -= D_b[p][s] * value;
 
242
 
 
243
        /* (qp|sr) */
 
244
        fock_a[q][p] += Dt[s][r] * value;
 
245
        fock_a[q][s] -= D_a[p][r] * value;
 
246
        fock_b[q][p] += Dt[s][r] * value;
 
247
        fock_b[q][s] -= D_b[p][r] * value;
 
248
 
 
249
        /* (rs|pq) */
 
250
        fock_a[r][s] += Dt[p][q] * value;
 
251
        fock_a[r][p] -= D_a[s][q] * value;
 
252
        fock_b[r][s] += Dt[p][q] * value;
 
253
        fock_b[r][p] -= D_b[s][q] * value;
 
254
 
 
255
        /* (rs|qp) */
 
256
        fock_a[r][s] += Dt[q][p] * value;
 
257
        fock_a[r][q] -= D_a[s][p] * value;
 
258
        fock_b[r][s] += Dt[q][p] * value;
 
259
        fock_b[r][q] -= D_b[s][p] * value;
 
260
 
 
261
        /* (sr|pq) */
 
262
        fock_a[s][r] += Dt[p][q] * value;
 
263
        fock_a[s][p] -= D_a[r][q] * value;
 
264
        fock_b[s][r] += Dt[p][q] * value;
 
265
        fock_b[s][p] -= D_b[r][q] * value;
 
266
 
 
267
        /* (sr|qp) */
 
268
        fock_a[s][r] += Dt[q][p] * value;
 
269
        fock_a[s][q] -= D_a[r][p] * value;
 
270
        fock_b[s][r] += Dt[q][p] * value;
 
271
        fock_b[s][q] -= D_b[r][p] * value;
 
272
      }
 
273
      else if(p!=q && r!=s && pq==rs) {
 
274
 
 
275
        /* (pq|sr) */
 
276
        fock_a[p][q] += Dt[s][r] * value;
 
277
        fock_a[p][s] -= D_a[q][r] * value;
 
278
        fock_b[p][q] += Dt[s][r] * value;
 
279
        fock_b[p][s] -= D_b[q][r] * value;
 
280
 
 
281
        /* (qp|rs) */
 
282
        fock_a[q][p] += Dt[r][s] * value;
 
283
        fock_a[q][r] -= D_a[p][s] * value;
 
284
        fock_b[q][p] += Dt[r][s] * value;
 
285
        fock_b[q][r] -= D_b[p][s] * value;
 
286
 
 
287
        /* (qp|sr) */
 
288
        fock_a[q][p] += Dt[s][r] * value;
 
289
        fock_a[q][s] -= D_a[p][r] * value;
 
290
        fock_b[q][p] += Dt[s][r] * value;
 
291
        fock_b[q][s] -= D_b[p][r] * value;
 
292
 
 
293
      }
 
294
      else if(p!=q && r==s) {
 
295
 
 
296
        /* (qp|rs) */
 
297
        fock_a[q][p] += Dt[r][s] * value;
 
298
        fock_a[q][r] -= D_a[p][s] * value;
 
299
        fock_b[q][p] += Dt[r][s] * value;
 
300
        fock_b[q][r] -= D_b[p][s] * value;
 
301
 
 
302
        /* (rs|pq) */
 
303
        fock_a[r][s] += Dt[p][q] * value;
 
304
        fock_a[r][p] -= D_a[s][q] * value;
 
305
        fock_b[r][s] += Dt[p][q] * value;
 
306
        fock_b[r][p] -= D_b[s][q] * value;
 
307
 
 
308
        /* (rs|qp) */
 
309
        fock_a[r][s] += Dt[q][p] * value;
 
310
        fock_a[r][q] -= D_a[s][p] * value;
 
311
        fock_b[r][s] += Dt[q][p] * value;
 
312
        fock_b[r][q] -= D_b[s][p] * value;
 
313
 
 
314
      }
 
315
      else if(p==q && r!=s) {
 
316
 
 
317
        /* (pq|sr) */
 
318
        fock_a[p][q] += Dt[s][r] * value;
 
319
        fock_a[p][s] -= D_a[q][r] * value;
 
320
        fock_b[p][q] += Dt[s][r] * value;
 
321
        fock_b[p][s] -= D_b[q][r] * value;
 
322
 
 
323
        /* (rs|pq) */
 
324
        fock_a[r][s] += Dt[p][q] * value;
 
325
        fock_a[r][p] -= D_a[s][q] * value;
 
326
        fock_b[r][s] += Dt[p][q] * value;
 
327
        fock_b[r][p] -= D_b[s][q] * value;
 
328
 
 
329
        /* (sr|pq) */
 
330
        fock_a[s][r] += Dt[p][q] * value;
 
331
        fock_a[s][p] -= D_a[r][q] * value;
 
332
        fock_b[s][r] += Dt[p][q] * value;
 
333
        fock_b[s][p] -= D_b[r][q] * value;
 
334
 
 
335
      }
 
336
      else if(p==q && r==s && pq!=rs) {
 
337
 
 
338
        /* (rs|pq) */
 
339
        fock_a[r][s] += Dt[p][q] * value;
 
340
        fock_a[r][p] -= D_a[s][q] * value;
 
341
        fock_b[r][s] += Dt[p][q] * value;
 
342
        fock_b[r][p] -= D_b[s][q] * value;
 
343
 
 
344
      }
 
345
    }
 
346
 
 
347
    if(!lastbuf) iwl_buf_fetch(&InBuf);
 
348
 
 
349
  } while (!lastbuf);
 
350
  iwl_buf_close(&InBuf, 1);
 
351
 
 
352
  free_block(Dt);
 
353
}
 
354
}} // namespace psi::ccenergy