← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/transqt2/transtwo_rhf.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/transqt2/transtwo_rhf.cc
 
1
/*! \file
 
2
    \ingroup TRANSQT2
 
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 <libpsio/psio.h>
 
10
#include <libqt/qt.h>
 
11
#include <libiwl/iwl.h>
 
12
#include <libchkpt/chkpt.h>
 
13
#include <libdpd/dpd.h>
 
14
#define EXTERN
 
15
#include "globals.h"
 
16
 
 
17
namespace psi {
 
18
  namespace transqt2 {
 
19
 
 
20
void transtwo_rhf(void)
 
21
{
 
22
  int nirreps, nso, nmo;
 
23
  double ***C, **scratch, **TMP;
 
24
  int h, p, q, r, s, pq, rs, Gs, Gr, PQ, RS;
 
25
  int nrows, ncols, nlinks;
 
26
  unsigned long int memfree, rows_per_bucket, rows_left, this_bucket_rows; 
 
27
  int nbuckets, n;
 
28
  dpdbuf4 J, K;
 
29
  struct iwlbuf MBuff;
 
30
  int stat;
 
31
  int *C_offset;
 
32
 
 
33
  nirreps = moinfo.nirreps;
 
34
  nso = moinfo.nso;
 
35
  nmo = moinfo.nmo;
 
36
  C = moinfo.C;
 
37
  C_offset = moinfo.C_offset;
 
38
 
 
39
  TMP = block_matrix(nso,nso);
 
40
 
 
41
  if(params.print_lvl) {
 
42
    fprintf(outfile, "\tStarting first half-transformation.\n");
 
43
    fflush(outfile);
 
44
  }
 
45
 
 
46
  psio_open(PSIF_SO_PRESORT, PSIO_OPEN_OLD);
 
47
  psio_open(PSIF_HALFT0, PSIO_OPEN_NEW);
 
48
 
 
49
  timer_on("RHF:1sthalf");
 
50
  dpd_buf4_init(&J, PSIF_SO_PRESORT, 0, 3, 0, 3, 3, 0, "SO Ints (pq,rs)");
 
51
  dpd_buf4_init(&K, PSIF_HALFT0, 0, 3, 5, 3, 8, 0, "Half-Transformed Ints (pq,ij)");
 
52
  for(h=0; h < nirreps; h++) {
 
53
 
 
54
    memfree = (unsigned long int) (dpd_memfree() - J.params->coltot[h] - K.params->coltot[h]);
 
55
    if(J.params->coltot[h]) {
 
56
      rows_per_bucket = memfree/(2 * J.params->coltot[h]);
 
57
      if(rows_per_bucket > J.params->rowtot[h]) rows_per_bucket = (unsigned long int) J.params->rowtot[h];
 
58
      nbuckets = (int) ceil(((double) J.params->rowtot[h])/((double) rows_per_bucket));
 
59
      rows_left = (unsigned long int) (J.params->rowtot[h] % rows_per_bucket);
 
60
    }
 
61
    else {
 
62
      nbuckets = 0;
 
63
      rows_per_bucket = 0;
 
64
      rows_left = 0;
 
65
    }
 
66
 
 
67
    if(params.print_lvl > 1) {
 
68
      fprintf(outfile, "\th = %d; memfree         = %lu\n", h, memfree);
 
69
      fprintf(outfile, "\th = %d; rows_per_bucket = %lu\n", h, rows_per_bucket);
 
70
      fprintf(outfile, "\th = %d; rows_left       = %lu\n", h, rows_left);
 
71
      fprintf(outfile, "\th = %d; nbuckets        = %d\n", h, nbuckets);
 
72
      fflush(outfile);
 
73
    }
 
74
 
 
75
    dpd_buf4_mat_irrep_init_block(&J, h, rows_per_bucket);
 
76
    dpd_buf4_mat_irrep_init_block(&K, h, rows_per_bucket);
 
77
 
 
78
    for(n=0; n < nbuckets; n++) {
 
79
      if(nbuckets == 1) this_bucket_rows = rows_per_bucket;
 
80
      else this_bucket_rows = (n < nbuckets-1) ? rows_per_bucket : rows_left;
 
81
      dpd_buf4_mat_irrep_rd_block(&J, h, n*rows_per_bucket, this_bucket_rows);
 
82
      for(pq=0; pq < this_bucket_rows; pq++) {
 
83
        for(Gr=0; Gr < nirreps; Gr++) {
 
84
          Gs = h^Gr;
 
85
          nrows = moinfo.sopi[Gr];
 
86
          ncols = moinfo.actpi[Gs];
 
87
          nlinks = moinfo.sopi[Gs];
 
88
          rs = J.col_offset[h][Gr];
 
89
          if(nrows && ncols && nlinks)
 
90
            C_DGEMM('n','n',nrows,ncols,nlinks,1.0,&J.matrix[h][pq][rs],nlinks,
 
91
                    &(C[Gs][0][C_offset[Gs]]),moinfo.mopi[Gs],0.0,TMP[0],nso);
 
92
 
 
93
          nrows = moinfo.actpi[Gr];
 
94
          ncols = moinfo.actpi[Gs];
 
95
          nlinks = moinfo.sopi[Gr];
 
96
          rs = K.col_offset[h][Gr];
 
97
          if(nrows && ncols && nlinks)
 
98
            C_DGEMM('t','n',nrows,ncols,nlinks,1.0,&(C[Gr][0][C_offset[Gr]]),moinfo.mopi[Gr],TMP[0],nso,
 
99
                    0.0,&K.matrix[h][pq][rs],ncols);
 
100
        } /* Gr */
 
101
      } /* pq */
 
102
      dpd_buf4_mat_irrep_wrt_block(&K, h, n*rows_per_bucket, this_bucket_rows);
 
103
    }
 
104
 
 
105
    dpd_buf4_mat_irrep_close_block(&J, h, rows_per_bucket);
 
106
    dpd_buf4_mat_irrep_close_block(&K, h, rows_per_bucket);
 
107
  }
 
108
  dpd_buf4_close(&K);
 
109
  dpd_buf4_close(&J);
 
110
 
 
111
  psio_close(PSIF_SO_PRESORT, 0);
 
112
 
 
113
  timer_off("RHF:1sthalf");
 
114
  timer_on("RHF:midsort");
 
115
  if(params.print_lvl) {
 
116
    fprintf(outfile, "\tSorting half-transformed integrals.\n");
 
117
    fflush(outfile);
 
118
  }
 
119
 
 
120
  psio_open(PSIF_HALFT1, PSIO_OPEN_NEW);
 
121
 
 
122
  dpd_buf4_init(&K, PSIF_HALFT0, 0, 3, 8, 3, 8, 0, "Half-Transformed Ints (pq,ij)");
 
123
  dpd_buf4_sort(&K, PSIF_HALFT1, rspq, 8, 3, "Half-Transformed Ints (ij,pq)");
 
124
  dpd_buf4_close(&K);
 
125
 
 
126
  psio_close(PSIF_HALFT0, 0);
 
127
  timer_off("RHF:midsort");
 
128
  timer_on("RHF:2ndhalf");
 
129
 
 
130
  if(params.print_lvl) {
 
131
    fprintf(outfile, "\tStarting second half-transformation.\n");
 
132
    fflush(outfile);
 
133
  }
 
134
  iwl_buf_init(&MBuff, PSIF_MO_TEI, params.tolerance, 0, 0);
 
135
 
 
136
  dpd_buf4_init(&J, PSIF_HALFT1, 0, 8, 0, 8, 3, 0, "Half-Transformed Ints (ij,pq)");
 
137
  dpd_buf4_init(&K, CC_MISC, 0, 8, 5, 8, 8, 0, "MO Ints (ij,kl)");
 
138
  for(h=0; h < nirreps; h++) {
 
139
 
 
140
    memfree = (unsigned long int) (dpd_memfree() - J.params->coltot[h] - K.params->coltot[h]);
 
141
    if(J.params->coltot[h]) {
 
142
      rows_per_bucket = memfree/(2 * J.params->coltot[h]);
 
143
      if(rows_per_bucket > J.params->rowtot[h]) rows_per_bucket = (unsigned long int) J.params->rowtot[h];
 
144
      nbuckets = (int) ceil(((double) J.params->rowtot[h])/((double) rows_per_bucket));
 
145
      rows_left = (unsigned long int) (J.params->rowtot[h] % rows_per_bucket);
 
146
    }
 
147
    else {
 
148
      nbuckets = 0;
 
149
      rows_per_bucket = 0;
 
150
      rows_left = 0;
 
151
    }
 
152
 
 
153
    if(params.print_lvl > 1) {
 
154
      fprintf(outfile, "\th = %d; memfree         = %lu\n", h, memfree);
 
155
      fprintf(outfile, "\th = %d; rows_per_bucket = %lu\n", h, rows_per_bucket);
 
156
      fprintf(outfile, "\th = %d; rows_left       = %lu\n", h, rows_left);
 
157
      fprintf(outfile, "\th = %d; nbuckets        = %d\n", h, nbuckets);
 
158
      fflush(outfile);
 
159
    }
 
160
 
 
161
    dpd_buf4_mat_irrep_init_block(&J, h, rows_per_bucket);
 
162
    dpd_buf4_mat_irrep_init_block(&K, h, rows_per_bucket);
 
163
 
 
164
    for(n=0; n < nbuckets; n++) {
 
165
      if(nbuckets == 1) this_bucket_rows = rows_per_bucket;
 
166
      else this_bucket_rows = (n < nbuckets-1) ? rows_per_bucket : rows_left;
 
167
      dpd_buf4_mat_irrep_rd_block(&J, h, n*rows_per_bucket, this_bucket_rows);
 
168
      for(pq=0; pq < this_bucket_rows; pq++) {
 
169
        for(Gr=0; Gr < nirreps; Gr++) {
 
170
          Gs = h^Gr;
 
171
          nrows = moinfo.sopi[Gr];
 
172
          ncols = moinfo.actpi[Gs];
 
173
          nlinks = moinfo.sopi[Gs];
 
174
          rs = J.col_offset[h][Gr];
 
175
          if(nrows && ncols && nlinks)
 
176
            C_DGEMM('n','n',nrows,ncols,nlinks,1.0,&J.matrix[h][pq][rs],nlinks,
 
177
                    &(C[Gs][0][C_offset[Gs]]),moinfo.mopi[Gs],0.0,TMP[0],nso);
 
178
 
 
179
          nrows = moinfo.actpi[Gr];
 
180
          ncols = moinfo.actpi[Gs];
 
181
          nlinks = moinfo.sopi[Gr];
 
182
          rs = K.col_offset[h][Gr];
 
183
          if(nrows && ncols && nlinks)
 
184
            C_DGEMM('t','n',nrows,ncols,nlinks,1.0,&(C[Gr][0][C_offset[Gr]]),moinfo.mopi[Gr],TMP[0],nso,
 
185
                    0.0,&K.matrix[h][pq][rs],ncols);
 
186
        } /* Gr */
 
187
 
 
188
        p = moinfo.pitz2corr_two[K.params->roworb[h][pq+n*rows_per_bucket][0]];
 
189
        q = moinfo.pitz2corr_two[K.params->roworb[h][pq+n*rows_per_bucket][1]];
 
190
        PQ = INDEX(p,q);
 
191
        for(rs=0; rs < K.params->coltot[h]; rs++) {
 
192
          r = moinfo.pitz2corr_two[K.params->colorb[h][rs][0]];
 
193
          s = moinfo.pitz2corr_two[K.params->colorb[h][rs][1]];
 
194
          RS = INDEX(r,s);
 
195
          if(r >= s && RS <= PQ)
 
196
            iwl_buf_wrt_val(&MBuff, p, q, r, s, K.matrix[h][pq][rs], params.print_tei, outfile, 0);
 
197
        } /* rs */
 
198
      } /* pq */
 
199
    }
 
200
    dpd_buf4_mat_irrep_close_block(&J, h, rows_per_bucket);
 
201
    dpd_buf4_mat_irrep_close_block(&K, h, rows_per_bucket);
 
202
  }
 
203
  dpd_buf4_close(&K);
 
204
  dpd_buf4_close(&J);
 
205
 
 
206
  iwl_buf_flush(&MBuff, 1);
 
207
  iwl_buf_close(&MBuff, 1);
 
208
 
 
209
  timer_off("RHF:2ndhalf");
 
210
 
 
211
  free_block(TMP);
 
212
 
 
213
  psio_close(PSIF_HALFT1, 0);
 
214
 
 
215
  if(params.print_lvl) {
 
216
    fprintf(outfile, "\tTwo-electron integral transformation complete.\n");
 
217
    fflush(outfile);
 
218
  }
 
219
}
 
220
 
 
221
  } // namespace transqt2
 
222
} // namespace psi