← Back to branch summary

~ubuntu-branches/ubuntu/vivid/psicode/vivid

« back to all changes in this revision

Viewing changes to src/bin/transqt2/transtwo_uhf.c

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck
  • Date: 2008-06-07 16:49:57 UTC
  • mfrom: (2.1.2 hardy)
  • Revision ID: james.westby@ubuntu.com-20080607164957-8pifvb133yjlkagn
Tags: 3.3.0-3
http://bugs.debian.org/464867
* debian/rules (DEB_MAKE_CHECK_TARGET): Do not abort test suite on
  failures.
* debian/rules (DEB_CONFIGURE_EXTRA_FLAGS): Set ${bindir} to /usr/lib/psi.
* debian/rules (install/psi3): Move psi3 file to /usr/bin.
* debian/patches/07_464867_move_executables.dpatch: New patch, add
  /usr/lib/psi to the $PATH, so that the moved executables are found.
  (closes: #464867)
* debian/patches/00list: Adjusted.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/bin/transqt2/transtwo_uhf.c
 
1
#include <stdio.h>
 
2
#include <stdlib.h>
 
3
#include <math.h>
 
4
#include <libciomr/libciomr.h>
 
5
#include <libpsio/psio.h>
 
6
#include <libqt/qt.h>
 
7
#include <libiwl/iwl.h>
 
8
#include <libchkpt/chkpt.h>
 
9
#include <libdpd/dpd.h>
 
10
#define EXTERN
 
11
#include "globals.h"
 
12
 
 
13
void transtwo_uhf(void)
 
14
{
 
15
  int nirreps, nso, nmo;
 
16
  double ***C, ***C_a, ***C_b, **scratch, **TMP;
 
17
  int h, p, q, r, s, pq, rs, Gs, Gr, PQ, RS;
 
18
  int nrows, ncols, nlinks;
 
19
  unsigned long int memfree, rows_per_bucket, rows_left; 
 
20
  int nbuckets, n;
 
21
  dpdbuf4 J, K;
 
22
  struct iwlbuf MBuff;
 
23
  int stat;
 
24
 
 
25
  nirreps = moinfo.nirreps;
 
26
  nso = moinfo.nso;
 
27
  nmo = moinfo.nmo;
 
28
 
 
29
  /* grab MOs and remove frozen core/virt */
 
30
  C_a = (double ***) malloc(nirreps * sizeof(double **));
 
31
  C_b = (double ***) malloc(nirreps * sizeof(double **));
 
32
  chkpt_init(PSIO_OPEN_OLD);
 
33
  for(h=0; h < nirreps; h++) {
 
34
    scratch = chkpt_rd_alpha_scf_irrep(h);
 
35
    C_a[h] = block_matrix(moinfo.sopi[h],moinfo.mopi[h]);
 
36
    for(q=0; q < moinfo.mopi[h]; q++)
 
37
      for(p=0; p < moinfo.sopi[h]; p++)
 
38
        C_a[h][p][q] = scratch[p][q+moinfo.frdocc[h]];
 
39
    if(params.print_lvl > 2) {
 
40
      fprintf(outfile, "\n\tAlpha MOs for irrep %d:\n",h);
 
41
      mat_print(C_a[h], moinfo.sopi[h], moinfo.mopi[h], outfile);
 
42
    }
 
43
    free_block(scratch);
 
44
 
 
45
    scratch = chkpt_rd_beta_scf_irrep(h);
 
46
    C_b[h] = block_matrix(moinfo.sopi[h],moinfo.mopi[h]);
 
47
    for(q=0; q < moinfo.mopi[h]; q++)
 
48
      for(p=0; p < moinfo.sopi[h]; p++)
 
49
        C_b[h][p][q] = scratch[p][q+moinfo.frdocc[h]];
 
50
    if(params.print_lvl > 2) {
 
51
      fprintf(outfile, "\n\tBeta MOs for irrep %d:\n",h);
 
52
      mat_print(C_b[h], moinfo.sopi[h], moinfo.mopi[h], outfile);
 
53
    }
 
54
    free_block(scratch);
 
55
  }
 
56
  chkpt_close();
 
57
 
 
58
  TMP = block_matrix(nso,nso);
 
59
 
 
60
  /*** AA/AB two-electron integral transformation ***/
 
61
 
 
62
  if(params.print_lvl) {
 
63
    fprintf(outfile, "\tStarting AA/AB first half-transformation.\n");
 
64
    fflush(outfile);
 
65
  }
 
66
 
 
67
  psio_open(PSIF_SO_PRESORT, PSIO_OPEN_OLD);
 
68
  psio_open(PSIF_HALFT0, PSIO_OPEN_NEW);
 
69
 
 
70
  C = C_a; /* Use alpha MOs for this half-transformation */
 
71
 
 
72
  dpd_buf4_init(&J, PSIF_SO_PRESORT, 0, 3, 0, 3, 3, 0, "SO Ints (pq,rs)");
 
73
  dpd_buf4_init(&K, PSIF_HALFT0, 0, 3, 5, 3, 8, 0, "Half-Transformed Ints (pq,ij)");
 
74
  for(h=0; h < nirreps; h++) {
 
75
 
 
76
    memfree = (unsigned long int) (dpd_memfree() - J.params->coltot[h] - K.params->coltot[h]);
 
77
    rows_per_bucket = memfree/(2 * J.params->coltot[h]);
 
78
    if(rows_per_bucket > J.params->rowtot[h]) rows_per_bucket = (unsigned long int) J.params->rowtot[h];
 
79
    nbuckets = (int) ceil(((double) J.params->rowtot[h])/((double) rows_per_bucket));
 
80
    rows_left = (unsigned long int) (J.params->rowtot[h] % rows_per_bucket);
 
81
    if(params.print_lvl > 1) {
 
82
      fprintf(outfile, "\th = %d; memfree         = %lu\n", h, memfree);
 
83
      fprintf(outfile, "\th = %d; rows_per_bucket = %lu\n", h, rows_per_bucket);
 
84
      fprintf(outfile, "\th = %d; rows_left       = %lu\n", h, rows_left);
 
85
      fprintf(outfile, "\th = %d; nbuckets        = %d\n", h, nbuckets);
 
86
      fflush(outfile);
 
87
    }
 
88
 
 
89
    dpd_buf4_mat_irrep_init_block(&J, h, rows_per_bucket);
 
90
    dpd_buf4_mat_irrep_init_block(&K, h, rows_per_bucket);
 
91
 
 
92
    for(n=0; n < (rows_left ? nbuckets-1 : nbuckets); n++) {
 
93
      dpd_buf4_mat_irrep_rd_block(&J, h, n*rows_per_bucket, rows_per_bucket);
 
94
      for(pq=0; pq < rows_per_bucket; pq++) {
 
95
        for(Gr=0; Gr < nirreps; Gr++) {
 
96
          Gs = h^Gr;
 
97
          nrows = moinfo.sopi[Gr];
 
98
          ncols = moinfo.mopi[Gs];
 
99
          nlinks = moinfo.sopi[Gs];
 
100
          rs = J.col_offset[h][Gr];
 
101
          if(nrows && ncols && nlinks)
 
102
            C_DGEMM('n','n',nrows,ncols,nlinks,1.0,&J.matrix[h][pq][rs],nlinks,
 
103
                    C[Gs][0],ncols,0.0,TMP[0],nso);
 
104
 
 
105
          nrows = moinfo.mopi[Gr];
 
106
          ncols = moinfo.mopi[Gs];
 
107
          nlinks = moinfo.sopi[Gr];
 
108
          rs = K.col_offset[h][Gr];
 
109
          if(nrows && ncols && nlinks)
 
110
            C_DGEMM('t','n',nrows,ncols,nlinks,1.0,C[Gr][0],nrows,TMP[0],nso,
 
111
                    0.0,&K.matrix[h][pq][rs],ncols);
 
112
        } /* Gr */
 
113
      } /* pq */
 
114
      dpd_buf4_mat_irrep_wrt_block(&K, h, n*rows_per_bucket, rows_per_bucket);
 
115
    }
 
116
    if(rows_left) {
 
117
      dpd_buf4_mat_irrep_rd_block(&J, h, n*rows_per_bucket, rows_left);
 
118
      for(pq=0; pq < rows_left; pq++) {
 
119
        for(Gr=0; Gr < nirreps; Gr++) {
 
120
          Gs = h^Gr;
 
121
 
 
122
          nrows = moinfo.sopi[Gr];
 
123
          ncols = moinfo.mopi[Gs];
 
124
          nlinks = moinfo.sopi[Gs];
 
125
          rs = J.col_offset[h][Gr];
 
126
          if(nrows && ncols && nlinks)
 
127
            C_DGEMM('n','n',nrows,ncols,nlinks,1.0,&J.matrix[h][pq][rs],nlinks,
 
128
                    C[Gs][0],ncols,0.0,TMP[0],nso);
 
129
 
 
130
          nrows = moinfo.mopi[Gr];
 
131
          ncols = moinfo.mopi[Gs];
 
132
          nlinks = moinfo.sopi[Gr];
 
133
          rs = K.col_offset[h][Gr];
 
134
          if(nrows && ncols && nlinks)
 
135
            C_DGEMM('t','n',nrows,ncols,nlinks,1.0,C[Gr][0],nrows,TMP[0],nso,
 
136
                    0.0,&K.matrix[h][pq][rs],ncols);
 
137
        } /* Gr */
 
138
      } /* pq */
 
139
 
 
140
      dpd_buf4_mat_irrep_wrt_block(&K, h, n*rows_per_bucket, rows_left);
 
141
    }
 
142
 
 
143
    dpd_buf4_mat_irrep_close_block(&J, h, rows_per_bucket);
 
144
    dpd_buf4_mat_irrep_close_block(&K, h, rows_per_bucket);
 
145
  }
 
146
  dpd_buf4_close(&K);
 
147
  dpd_buf4_close(&J);
 
148
 
 
149
  psio_close(PSIF_SO_PRESORT, 1); /* must keep the presort file for the upcoming BB transformation */
 
150
 
 
151
  if(params.print_lvl) {
 
152
    fprintf(outfile, "\tSorting AA/AB half-transformed integrals.\n");
 
153
    fflush(outfile);
 
154
  }
 
155
 
 
156
  psio_open(PSIF_HALFT1, PSIO_OPEN_NEW);
 
157
 
 
158
  dpd_buf4_init(&K, PSIF_HALFT0, 0, 3, 8, 3, 8, 0, "Half-Transformed Ints (pq,ij)");
 
159
  dpd_buf4_sort(&K, PSIF_HALFT1, rspq, 8, 3, "Half-Transformed Ints (ij,pq)");
 
160
  dpd_buf4_close(&K);
 
161
 
 
162
  psio_close(PSIF_HALFT0, 0);
 
163
 
 
164
  if(params.print_lvl) {
 
165
    fprintf(outfile, "\tStarting AA second half-transformation.\n");
 
166
    fflush(outfile);
 
167
  }
 
168
  iwl_buf_init(&MBuff, PSIF_MO_AA_TEI, params.tolerance, 0, 0);
 
169
 
 
170
  C = C_a; /* Usa alpha MOs for this half-transformation */
 
171
 
 
172
  dpd_buf4_init(&J, PSIF_HALFT1, 0, 8, 0, 8, 3, 0, "Half-Transformed Ints (ij,pq)");
 
173
  dpd_buf4_init(&K, CC_MISC, 0, 8, 5, 8, 8, 0, "MO Ints (ij,kl)");
 
174
  for(h=0; h < nirreps; h++) {
 
175
 
 
176
    memfree = (unsigned long int) (dpd_memfree() - J.params->coltot[h] - K.params->coltot[h]);
 
177
    rows_per_bucket = memfree/(2 * J.params->coltot[h]);
 
178
    if(rows_per_bucket > J.params->rowtot[h]) rows_per_bucket = (unsigned long int) J.params->rowtot[h];
 
179
    nbuckets = (int) ceil(((double) J.params->rowtot[h])/((double) rows_per_bucket));
 
180
    rows_left = (unsigned long int) (J.params->rowtot[h] % rows_per_bucket);
 
181
 
 
182
    if(params.print_lvl > 1) {
 
183
      fprintf(outfile, "\th = %d; memfree         = %lu\n", h, memfree);
 
184
      fprintf(outfile, "\th = %d; rows_per_bucket = %lu\n", h, rows_per_bucket);
 
185
      fprintf(outfile, "\th = %d; rows_left       = %lu\n", h, rows_left);
 
186
      fprintf(outfile, "\th = %d; nbuckets        = %d\n", h, nbuckets);
 
187
      fflush(outfile);
 
188
    }
 
189
 
 
190
    dpd_buf4_mat_irrep_init_block(&J, h, rows_per_bucket);
 
191
    dpd_buf4_mat_irrep_init_block(&K, h, rows_per_bucket);
 
192
 
 
193
    for(n=0; n < (rows_left ? nbuckets-1 : nbuckets); n++) {
 
194
      dpd_buf4_mat_irrep_rd_block(&J, h, n*rows_per_bucket, rows_per_bucket);
 
195
      for(pq=0; pq < rows_per_bucket; pq++) {
 
196
        for(Gr=0; Gr < nirreps; Gr++) {
 
197
          Gs = h^Gr;
 
198
          nrows = moinfo.sopi[Gr];
 
199
          ncols = moinfo.mopi[Gs];
 
200
          nlinks = moinfo.sopi[Gs];
 
201
          rs = J.col_offset[h][Gr];
 
202
          if(nrows && ncols && nlinks)
 
203
            C_DGEMM('n','n',nrows,ncols,nlinks,1.0,&J.matrix[h][pq][rs],nlinks,
 
204
                    C[Gs][0],ncols,0.0,TMP[0],nso);
 
205
 
 
206
          nrows = moinfo.mopi[Gr];
 
207
          ncols = moinfo.mopi[Gs];
 
208
          nlinks = moinfo.sopi[Gr];
 
209
          rs = K.col_offset[h][Gr];
 
210
          if(nrows && ncols && nlinks)
 
211
            C_DGEMM('t','n',nrows,ncols,nlinks,1.0,C[Gr][0],nrows,TMP[0],nso,
 
212
                    0.0,&K.matrix[h][pq][rs],ncols);
 
213
        } /* Gr */
 
214
 
 
215
        p = moinfo.act2qt_A[K.params->roworb[h][pq+n*rows_per_bucket][0]];
 
216
        q = moinfo.act2qt_A[K.params->roworb[h][pq+n*rows_per_bucket][1]];
 
217
        PQ = INDEX(p,q);
 
218
        for(rs=0; rs < K.params->coltot[h]; rs++) {
 
219
          r = moinfo.act2qt_A[K.params->colorb[h][rs][0]];
 
220
          s = moinfo.act2qt_A[K.params->colorb[h][rs][1]];
 
221
          RS = INDEX(r,s);
 
222
          if(r >= s && RS <= PQ)
 
223
            iwl_buf_wrt_val(&MBuff, p, q, r, s, K.matrix[h][pq][rs], params.print_tei, outfile, 0);
 
224
        } /* rs */
 
225
      } /* pq */
 
226
    }
 
227
    if(rows_left) {
 
228
      dpd_buf4_mat_irrep_rd_block(&J, h, n*rows_per_bucket, rows_left);
 
229
      for(pq=0; pq < rows_left; pq++) {
 
230
        for(Gr=0; Gr < nirreps; Gr++) {
 
231
          Gs = h^Gr;
 
232
          nrows = moinfo.sopi[Gr];
 
233
          ncols = moinfo.mopi[Gs];
 
234
          nlinks = moinfo.sopi[Gs];
 
235
          rs = J.col_offset[h][Gr];
 
236
          if(nrows && ncols && nlinks)
 
237
            C_DGEMM('n','n',nrows,ncols,nlinks,1.0,&J.matrix[h][pq][rs],nlinks,
 
238
                    C[Gs][0],ncols,0.0,TMP[0],nso);
 
239
 
 
240
          nrows = moinfo.mopi[Gr];
 
241
          ncols = moinfo.mopi[Gs];
 
242
          nlinks = moinfo.sopi[Gr];
 
243
          rs = K.col_offset[h][Gr];
 
244
          if(nrows && ncols && nlinks)
 
245
            C_DGEMM('t','n',nrows,ncols,nlinks,1.0,C[Gr][0],nrows,TMP[0],nso,
 
246
                    0.0,&K.matrix[h][pq][rs],ncols);
 
247
        } /* Gr */
 
248
 
 
249
        p = moinfo.act2qt_A[K.params->roworb[h][pq+n*rows_per_bucket][0]];
 
250
        q = moinfo.act2qt_A[K.params->roworb[h][pq+n*rows_per_bucket][1]];
 
251
        PQ = INDEX(p,q);
 
252
        for(rs=0; rs < K.params->coltot[h]; rs++) {
 
253
          r = moinfo.act2qt_A[K.params->colorb[h][rs][0]];
 
254
          s = moinfo.act2qt_A[K.params->colorb[h][rs][1]];
 
255
          RS = INDEX(r,s);
 
256
          if(r >= s && RS <= PQ)
 
257
            iwl_buf_wrt_val(&MBuff, p, q, r, s, K.matrix[h][pq][rs], params.print_tei, outfile, 0);
 
258
        } /* rs */
 
259
      } /* pq */
 
260
    }
 
261
    dpd_buf4_mat_irrep_close_block(&J, h, rows_per_bucket);
 
262
    dpd_buf4_mat_irrep_close_block(&K, h, rows_per_bucket);
 
263
  }
 
264
  dpd_buf4_close(&K);
 
265
  dpd_buf4_close(&J);
 
266
 
 
267
  iwl_buf_flush(&MBuff, 1);
 
268
  iwl_buf_close(&MBuff, 1);
 
269
 
 
270
  if(params.print_lvl) {
 
271
    fprintf(outfile, "\tStarting AB second half-transformation.\n");
 
272
    fflush(outfile);
 
273
  }
 
274
  iwl_buf_init(&MBuff, PSIF_MO_AB_TEI, params.tolerance, 0, 0);
 
275
 
 
276
  C = C_b; /* Usa beta MOs for this half-transformation */
 
277
 
 
278
  dpd_buf4_init(&J, PSIF_HALFT1, 0, 8, 0, 8, 3, 0, "Half-Transformed Ints (ij,pq)");
 
279
  dpd_buf4_init(&K, CC_MISC, 0, 8, 5, 8, 8, 0, "MO Ints (ij,kl)");
 
280
  for(h=0; h < nirreps; h++) {
 
281
 
 
282
    memfree = (unsigned long int) (dpd_memfree() - J.params->coltot[h] - K.params->coltot[h]);
 
283
    rows_per_bucket = memfree/(2 * J.params->coltot[h]);
 
284
    if(rows_per_bucket > J.params->rowtot[h]) rows_per_bucket = (unsigned long int) J.params->rowtot[h];
 
285
    nbuckets = (int) ceil(((double) J.params->rowtot[h])/((double) rows_per_bucket));
 
286
    rows_left = (unsigned long int) (J.params->rowtot[h] % rows_per_bucket);
 
287
 
 
288
    if(params.print_lvl > 1) {
 
289
      fprintf(outfile, "\th = %d; memfree         = %lu\n", h, memfree);
 
290
      fprintf(outfile, "\th = %d; rows_per_bucket = %lu\n", h, rows_per_bucket);
 
291
      fprintf(outfile, "\th = %d; rows_left       = %lu\n", h, rows_left);
 
292
      fprintf(outfile, "\th = %d; nbuckets        = %d\n", h, nbuckets);
 
293
      fflush(outfile);
 
294
    }
 
295
 
 
296
    dpd_buf4_mat_irrep_init_block(&J, h, rows_per_bucket);
 
297
    dpd_buf4_mat_irrep_init_block(&K, h, rows_per_bucket);
 
298
 
 
299
    for(n=0; n < (rows_left ? nbuckets-1 : nbuckets); n++) {
 
300
      dpd_buf4_mat_irrep_rd_block(&J, h, n*rows_per_bucket, rows_per_bucket);
 
301
      for(pq=0; pq < rows_per_bucket; pq++) {
 
302
        for(Gr=0; Gr < nirreps; Gr++) {
 
303
          Gs = h^Gr;
 
304
          nrows = moinfo.sopi[Gr];
 
305
          ncols = moinfo.mopi[Gs];
 
306
          nlinks = moinfo.sopi[Gs];
 
307
          rs = J.col_offset[h][Gr];
 
308
          if(nrows && ncols && nlinks)
 
309
            C_DGEMM('n','n',nrows,ncols,nlinks,1.0,&J.matrix[h][pq][rs],nlinks,
 
310
                    C[Gs][0],ncols,0.0,TMP[0],nso);
 
311
 
 
312
          nrows = moinfo.mopi[Gr];
 
313
          ncols = moinfo.mopi[Gs];
 
314
          nlinks = moinfo.sopi[Gr];
 
315
          rs = K.col_offset[h][Gr];
 
316
          if(nrows && ncols && nlinks)
 
317
            C_DGEMM('t','n',nrows,ncols,nlinks,1.0,C[Gr][0],nrows,TMP[0],nso,
 
318
                    0.0,&K.matrix[h][pq][rs],ncols);
 
319
        } /* Gr */
 
320
 
 
321
        p = moinfo.act2qt_A[K.params->roworb[h][pq+n*rows_per_bucket][0]];
 
322
        q = moinfo.act2qt_A[K.params->roworb[h][pq+n*rows_per_bucket][1]];
 
323
        PQ = INDEX(p,q);
 
324
        for(rs=0; rs < K.params->coltot[h]; rs++) {
 
325
          r = moinfo.act2qt_B[K.params->colorb[h][rs][0]];
 
326
          s = moinfo.act2qt_B[K.params->colorb[h][rs][1]];
 
327
          RS = INDEX(r,s);
 
328
          if(r >= s)
 
329
            iwl_buf_wrt_val(&MBuff, p, q, r, s, K.matrix[h][pq][rs], params.print_tei, outfile, 0);
 
330
        } /* rs */
 
331
      } /* pq */
 
332
    }
 
333
    if(rows_left) {
 
334
      dpd_buf4_mat_irrep_rd_block(&J, h, n*rows_per_bucket, rows_left);
 
335
      for(pq=0; pq < rows_left; pq++) {
 
336
        for(Gr=0; Gr < nirreps; Gr++) {
 
337
          Gs = h^Gr;
 
338
          nrows = moinfo.sopi[Gr];
 
339
          ncols = moinfo.mopi[Gs];
 
340
          nlinks = moinfo.sopi[Gs];
 
341
          rs = J.col_offset[h][Gr];
 
342
          if(nrows && ncols && nlinks)
 
343
            C_DGEMM('n','n',nrows,ncols,nlinks,1.0,&J.matrix[h][pq][rs],nlinks,
 
344
                    C[Gs][0],ncols,0.0,TMP[0],nso);
 
345
 
 
346
          nrows = moinfo.mopi[Gr];
 
347
          ncols = moinfo.mopi[Gs];
 
348
          nlinks = moinfo.sopi[Gr];
 
349
          rs = K.col_offset[h][Gr];
 
350
          if(nrows && ncols && nlinks)
 
351
            C_DGEMM('t','n',nrows,ncols,nlinks,1.0,C[Gr][0],nrows,TMP[0],nso,
 
352
                    0.0,&K.matrix[h][pq][rs],ncols);
 
353
        } /* Gr */
 
354
 
 
355
        p = moinfo.act2qt_A[K.params->roworb[h][pq+n*rows_per_bucket][0]];
 
356
        q = moinfo.act2qt_A[K.params->roworb[h][pq+n*rows_per_bucket][1]];
 
357
        PQ = INDEX(p,q);
 
358
        for(rs=0; rs < K.params->coltot[h]; rs++) {
 
359
          r = moinfo.act2qt_B[K.params->colorb[h][rs][0]];
 
360
          s = moinfo.act2qt_B[K.params->colorb[h][rs][1]];
 
361
          RS = INDEX(r,s);
 
362
          if(r >= s)
 
363
            iwl_buf_wrt_val(&MBuff, p, q, r, s, K.matrix[h][pq][rs], params.print_tei, outfile, 0);
 
364
        } /* rs */
 
365
      } /* pq */
 
366
    }
 
367
    dpd_buf4_mat_irrep_close_block(&J, h, rows_per_bucket);
 
368
    dpd_buf4_mat_irrep_close_block(&K, h, rows_per_bucket);
 
369
  }
 
370
  dpd_buf4_close(&K);
 
371
  dpd_buf4_close(&J);
 
372
 
 
373
  iwl_buf_flush(&MBuff, 1);
 
374
  iwl_buf_close(&MBuff, 1);
 
375
 
 
376
  psio_close(PSIF_HALFT1, 0);
 
377
 
 
378
  /*** AA/AB two-electron integral transformation complete ***/
 
379
 
 
380
  /*** BB two-electron integral transformation ***/
 
381
 
 
382
  if(params.print_lvl) {
 
383
    fprintf(outfile, "\tStarting BB first half-transformation.\n");
 
384
    fflush(outfile);
 
385
  }
 
386
 
 
387
  psio_open(PSIF_SO_PRESORT, PSIO_OPEN_OLD);
 
388
  psio_open(PSIF_HALFT0, PSIO_OPEN_NEW);
 
389
 
 
390
  C = C_b; /* Use beta MOs for this half-transformation */
 
391
 
 
392
  dpd_buf4_init(&J, PSIF_SO_PRESORT, 0, 3, 0, 3, 3, 0, "SO Ints (pq,rs)");
 
393
  dpd_buf4_init(&K, PSIF_HALFT0, 0, 3, 5, 3, 8, 0, "Half-Transformed Ints (pq,ij)");
 
394
  for(h=0; h < nirreps; h++) {
 
395
 
 
396
    memfree = (unsigned long int) (dpd_memfree() - J.params->coltot[h] - K.params->coltot[h]);
 
397
    rows_per_bucket = memfree/(2 * J.params->coltot[h]);
 
398
    if(rows_per_bucket > J.params->rowtot[h]) rows_per_bucket = (unsigned long int) J.params->rowtot[h];
 
399
    nbuckets = (int) ceil(((double) J.params->rowtot[h])/((double) rows_per_bucket));
 
400
    rows_left = (unsigned long int) (J.params->rowtot[h] % rows_per_bucket);
 
401
    if(params.print_lvl > 1) {
 
402
      fprintf(outfile, "\th = %d; memfree         = %lu\n", h, memfree);
 
403
      fprintf(outfile, "\th = %d; rows_per_bucket = %lu\n", h, rows_per_bucket);
 
404
      fprintf(outfile, "\th = %d; rows_left       = %lu\n", h, rows_left);
 
405
      fprintf(outfile, "\th = %d; nbuckets        = %d\n", h, nbuckets);
 
406
      fflush(outfile);
 
407
    }
 
408
 
 
409
    dpd_buf4_mat_irrep_init_block(&J, h, rows_per_bucket);
 
410
    dpd_buf4_mat_irrep_init_block(&K, h, rows_per_bucket);
 
411
 
 
412
    for(n=0; n < (rows_left ? nbuckets-1 : nbuckets); n++) {
 
413
      dpd_buf4_mat_irrep_rd_block(&J, h, n*rows_per_bucket, rows_per_bucket);
 
414
      for(pq=0; pq < rows_per_bucket; pq++) {
 
415
        for(Gr=0; Gr < nirreps; Gr++) {
 
416
          Gs = h^Gr;
 
417
          nrows = moinfo.sopi[Gr];
 
418
          ncols = moinfo.mopi[Gs];
 
419
          nlinks = moinfo.sopi[Gs];
 
420
          rs = J.col_offset[h][Gr];
 
421
          if(nrows && ncols && nlinks)
 
422
            C_DGEMM('n','n',nrows,ncols,nlinks,1.0,&J.matrix[h][pq][rs],nlinks,
 
423
                    C[Gs][0],ncols,0.0,TMP[0],nso);
 
424
 
 
425
          nrows = moinfo.mopi[Gr];
 
426
          ncols = moinfo.mopi[Gs];
 
427
          nlinks = moinfo.sopi[Gr];
 
428
          rs = K.col_offset[h][Gr];
 
429
          if(nrows && ncols && nlinks)
 
430
            C_DGEMM('t','n',nrows,ncols,nlinks,1.0,C[Gr][0],nrows,TMP[0],nso,
 
431
                    0.0,&K.matrix[h][pq][rs],ncols);
 
432
        } /* Gr */
 
433
      } /* pq */
 
434
      dpd_buf4_mat_irrep_wrt_block(&K, h, n*rows_per_bucket, rows_per_bucket);
 
435
    }
 
436
    if(rows_left) {
 
437
      dpd_buf4_mat_irrep_rd_block(&J, h, n*rows_per_bucket, rows_left);
 
438
      for(pq=0; pq < rows_left; pq++) {
 
439
        for(Gr=0; Gr < nirreps; Gr++) {
 
440
          Gs = h^Gr;
 
441
 
 
442
          nrows = moinfo.sopi[Gr];
 
443
          ncols = moinfo.mopi[Gs];
 
444
          nlinks = moinfo.sopi[Gs];
 
445
          rs = J.col_offset[h][Gr];
 
446
          if(nrows && ncols && nlinks)
 
447
            C_DGEMM('n','n',nrows,ncols,nlinks,1.0,&J.matrix[h][pq][rs],nlinks,
 
448
                    C[Gs][0],ncols,0.0,TMP[0],nso);
 
449
 
 
450
          nrows = moinfo.mopi[Gr];
 
451
          ncols = moinfo.mopi[Gs];
 
452
          nlinks = moinfo.sopi[Gr];
 
453
          rs = K.col_offset[h][Gr];
 
454
          if(nrows && ncols && nlinks)
 
455
            C_DGEMM('t','n',nrows,ncols,nlinks,1.0,C[Gr][0],nrows,TMP[0],nso,
 
456
                    0.0,&K.matrix[h][pq][rs],ncols);
 
457
        } /* Gr */
 
458
      } /* pq */
 
459
 
 
460
      dpd_buf4_mat_irrep_wrt_block(&K, h, n*rows_per_bucket, rows_left);
 
461
    }
 
462
 
 
463
    dpd_buf4_mat_irrep_close_block(&J, h, rows_per_bucket);
 
464
    dpd_buf4_mat_irrep_close_block(&K, h, rows_per_bucket);
 
465
  }
 
466
  dpd_buf4_close(&K);
 
467
  dpd_buf4_close(&J);
 
468
 
 
469
  psio_close(PSIF_SO_PRESORT, 0); 
 
470
 
 
471
  if(params.print_lvl) {
 
472
    fprintf(outfile, "\tSorting BB half-transformed integrals.\n");
 
473
    fflush(outfile);
 
474
  }
 
475
 
 
476
  psio_open(PSIF_HALFT1, PSIO_OPEN_NEW);
 
477
 
 
478
  dpd_buf4_init(&K, PSIF_HALFT0, 0, 3, 8, 3, 8, 0, "Half-Transformed Ints (pq,ij)");
 
479
  dpd_buf4_sort(&K, PSIF_HALFT1, rspq, 8, 3, "Half-Transformed Ints (ij,pq)");
 
480
  dpd_buf4_close(&K);
 
481
 
 
482
  psio_close(PSIF_HALFT0, 0);
 
483
 
 
484
  if(params.print_lvl) {
 
485
    fprintf(outfile, "\tStarting BB second half-transformation.\n");
 
486
    fflush(outfile);
 
487
  }
 
488
  iwl_buf_init(&MBuff, PSIF_MO_BB_TEI, params.tolerance, 0, 0);
 
489
 
 
490
  C = C_b; /* Usa beta MOs for this half-transformation */
 
491
 
 
492
  dpd_buf4_init(&J, PSIF_HALFT1, 0, 8, 0, 8, 3, 0, "Half-Transformed Ints (ij,pq)");
 
493
  dpd_buf4_init(&K, CC_MISC, 0, 8, 5, 8, 8, 0, "MO Ints (ij,kl)");
 
494
  for(h=0; h < nirreps; h++) {
 
495
 
 
496
    memfree = (unsigned long int) (dpd_memfree() - J.params->coltot[h] - K.params->coltot[h]);
 
497
    rows_per_bucket = memfree/(2 * J.params->coltot[h]);
 
498
    if(rows_per_bucket > J.params->rowtot[h]) rows_per_bucket = (unsigned long int) J.params->rowtot[h];
 
499
    nbuckets = (int) ceil(((double) J.params->rowtot[h])/((double) rows_per_bucket));
 
500
    rows_left = (unsigned long int) (J.params->rowtot[h] % rows_per_bucket);
 
501
 
 
502
    if(params.print_lvl > 1) {
 
503
      fprintf(outfile, "\th = %d; memfree         = %lu\n", h, memfree);
 
504
      fprintf(outfile, "\th = %d; rows_per_bucket = %lu\n", h, rows_per_bucket);
 
505
      fprintf(outfile, "\th = %d; rows_left       = %lu\n", h, rows_left);
 
506
      fprintf(outfile, "\th = %d; nbuckets        = %d\n", h, nbuckets);
 
507
      fflush(outfile);
 
508
    }
 
509
 
 
510
    dpd_buf4_mat_irrep_init_block(&J, h, rows_per_bucket);
 
511
    dpd_buf4_mat_irrep_init_block(&K, h, rows_per_bucket);
 
512
 
 
513
    for(n=0; n < (rows_left ? nbuckets-1 : nbuckets); n++) {
 
514
      dpd_buf4_mat_irrep_rd_block(&J, h, n*rows_per_bucket, rows_per_bucket);
 
515
      for(pq=0; pq < rows_per_bucket; pq++) {
 
516
        for(Gr=0; Gr < nirreps; Gr++) {
 
517
          Gs = h^Gr;
 
518
          nrows = moinfo.sopi[Gr];
 
519
          ncols = moinfo.mopi[Gs];
 
520
          nlinks = moinfo.sopi[Gs];
 
521
          rs = J.col_offset[h][Gr];
 
522
          if(nrows && ncols && nlinks)
 
523
            C_DGEMM('n','n',nrows,ncols,nlinks,1.0,&J.matrix[h][pq][rs],nlinks,
 
524
                    C[Gs][0],ncols,0.0,TMP[0],nso);
 
525
 
 
526
          nrows = moinfo.mopi[Gr];
 
527
          ncols = moinfo.mopi[Gs];
 
528
          nlinks = moinfo.sopi[Gr];
 
529
          rs = K.col_offset[h][Gr];
 
530
          if(nrows && ncols && nlinks)
 
531
            C_DGEMM('t','n',nrows,ncols,nlinks,1.0,C[Gr][0],nrows,TMP[0],nso,
 
532
                    0.0,&K.matrix[h][pq][rs],ncols);
 
533
        } /* Gr */
 
534
 
 
535
        p = moinfo.act2qt_B[K.params->roworb[h][pq+n*rows_per_bucket][0]];
 
536
        q = moinfo.act2qt_B[K.params->roworb[h][pq+n*rows_per_bucket][1]];
 
537
        PQ = INDEX(p,q);
 
538
        for(rs=0; rs < K.params->coltot[h]; rs++) {
 
539
          r = moinfo.act2qt_B[K.params->colorb[h][rs][0]];
 
540
          s = moinfo.act2qt_B[K.params->colorb[h][rs][1]];
 
541
          RS = INDEX(r,s);
 
542
          if(r >= s && RS <= PQ)
 
543
            iwl_buf_wrt_val(&MBuff, p, q, r, s, K.matrix[h][pq][rs], params.print_tei, outfile, 0);
 
544
        } /* rs */
 
545
      } /* pq */
 
546
    }
 
547
    if(rows_left) {
 
548
      dpd_buf4_mat_irrep_rd_block(&J, h, n*rows_per_bucket, rows_left);
 
549
      for(pq=0; pq < rows_left; pq++) {
 
550
        for(Gr=0; Gr < nirreps; Gr++) {
 
551
          Gs = h^Gr;
 
552
          nrows = moinfo.sopi[Gr];
 
553
          ncols = moinfo.mopi[Gs];
 
554
          nlinks = moinfo.sopi[Gs];
 
555
          rs = J.col_offset[h][Gr];
 
556
          if(nrows && ncols && nlinks)
 
557
            C_DGEMM('n','n',nrows,ncols,nlinks,1.0,&J.matrix[h][pq][rs],nlinks,
 
558
                    C[Gs][0],ncols,0.0,TMP[0],nso);
 
559
 
 
560
          nrows = moinfo.mopi[Gr];
 
561
          ncols = moinfo.mopi[Gs];
 
562
          nlinks = moinfo.sopi[Gr];
 
563
          rs = K.col_offset[h][Gr];
 
564
          if(nrows && ncols && nlinks)
 
565
            C_DGEMM('t','n',nrows,ncols,nlinks,1.0,C[Gr][0],nrows,TMP[0],nso,
 
566
                    0.0,&K.matrix[h][pq][rs],ncols);
 
567
        } /* Gr */
 
568
 
 
569
        p = moinfo.act2qt_A[K.params->roworb[h][pq+n*rows_per_bucket][0]];
 
570
        q = moinfo.act2qt_A[K.params->roworb[h][pq+n*rows_per_bucket][1]];
 
571
        PQ = INDEX(p,q);
 
572
        for(rs=0; rs < K.params->coltot[h]; rs++) {
 
573
          r = moinfo.act2qt_A[K.params->colorb[h][rs][0]];
 
574
          s = moinfo.act2qt_A[K.params->colorb[h][rs][1]];
 
575
          RS = INDEX(r,s);
 
576
          if(r >= s && RS <= PQ)
 
577
            iwl_buf_wrt_val(&MBuff, p, q, r, s, K.matrix[h][pq][rs], params.print_tei, outfile, 0);
 
578
        } /* rs */
 
579
      } /* pq */
 
580
    }
 
581
    dpd_buf4_mat_irrep_close_block(&J, h, rows_per_bucket);
 
582
    dpd_buf4_mat_irrep_close_block(&K, h, rows_per_bucket);
 
583
  }
 
584
  dpd_buf4_close(&K);
 
585
  dpd_buf4_close(&J);
 
586
 
 
587
  iwl_buf_flush(&MBuff, 1);
 
588
  iwl_buf_close(&MBuff, 1);
 
589
 
 
590
  /*** BB two-electron integral transformation complete ***/
 
591
 
 
592
  free_block(TMP);
 
593
 
 
594
  for(h=0; h < nirreps; h++) {
 
595
    free_block(C_a[h]);
 
596
    free_block(C_b[h]);
 
597
  }
 
598
  free(C_a);
 
599
  free(C_b);
 
600
 
 
601
  if(params.print_lvl) {
 
602
    fprintf(outfile, "\tTwo-electron integral transformation complete.\n");
 
603
    fflush(outfile);
 
604
  }
 
605
}