← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/ccsort/transform_two_rhf.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/ccsort/transform_two_rhf.c
 
1
/* transform_two_rhf(): Carry out a four-index transform of the
 
2
** SO-basis two-electron integrals using DPD constructs.
 
3
**
 
4
** NB: Two key DPD pairs for the SO/MO subspaces:
 
5
**  pair   spaces
 
6
**    0    All SO,SO
 
7
**    3    SO >= SO (+)
 
8
**    5    All SO,SO
 
9
**    8    MO >= MO (+)
 
10
**
 
11
** TDC, March 2006
 
12
*/
 
13
 
 
14
#include <stdio.h>
 
15
#include <stdlib.h>
 
16
#include <libciomr/libciomr.h>
 
17
#include <libiwl/iwl.h>
 
18
#include <libqt/qt.h>
 
19
#include <libchkpt/chkpt.h>
 
20
#include <psifiles.h>
 
21
#include <ccfiles.h>
 
22
#include <libdpd/dpd.h>
 
23
#include "MOInfo.h"
 
24
#include "Params.h"
 
25
#define EXTERN
 
26
#include "globals.h"
 
27
 
 
28
#define CC_TEI_HALFT0 91
 
29
#define CC_TEI_HALFT1 92
 
30
 
 
31
int file_build_presort(dpdfile4 *, int, double, int);
 
32
 
 
33
void transform_two_rhf(void)
 
34
{
 
35
  dpdfile4 I;
 
36
  dpdbuf4 J, K;
 
37
  int h, pq, p, q;
 
38
  double **TMP;
 
39
  int *sopi, *mopi;
 
40
  int nso, nmo, nrows, ncols, nlinks, J_rs, K_rs, Gs, Gr;
 
41
  double ***C;
 
42
  int *pitz2dpd_occ, *pitz2dpd_vir;
 
43
 
 
44
  nso = moinfo.nso;
 
45
  nmo = moinfo.nmo;
 
46
  sopi = moinfo.sopi;
 
47
  mopi = moinfo.mopi;
 
48
  C = (double ***) malloc(moinfo.nirreps * sizeof(double **));
 
49
  chkpt_init(PSIO_OPEN_OLD);
 
50
  for(h=0; h < moinfo.nirreps; h++)
 
51
    C[h] = chkpt_rd_scf_irrep(h);
 
52
  chkpt_close();
 
53
 
 
54
  /* We'll need to map from Pitzer directly to CC/DPD ordering, so let's build the lookups */
 
55
  pitz2dpd_occ = init_int_array(moinfo.nactive);
 
56
  pitz2dpd_vir = init_int_array(moinfo.nactive);
 
57
  for(i=0; i < moinfo.nactive; i++) {
 
58
    pitz2dpd_occ[i] = moinfo.cc_occ[moinfo.pitz2qt[i]];
 
59
    pitz2dpd_vir[i] = moinfo.cc_vir[moinfo.pitz2qt[i]];
 
60
  }
 
61
 
 
62
  psio_open(PSIF_SO_PRESORT, 0); /* presorted integrals */
 
63
 
 
64
  dpd_file4_init(&I, PSIF_SO_PRESORT, 0, 3, 3, "SO Ints (pq,rs)");
 
65
  file_build_presort(&I, PSIF_SO_TEI, params.tolerance, 1);  /* still need to add frozen-core operator in here */
 
66
  dpd_file4_close(&I);
 
67
 
 
68
  TMP = block_matrix(nso,nso);
 
69
 
 
70
  psio_open(CC_TEI_HALFT0, 0); /* half-transformed integrals */
 
71
 
 
72
  dpd_buf4_init(&J, PSIF_SO_PRESORT, 0, 3, 0, 3, 3, 0, "SO Ints (pq,rs)");
 
73
  dpd_buf4_init(&K, CC_TEI_HALFT0, 0, 3, 5, 3, 8, 0, "Half-Transformed Ints (pq,ij)");
 
74
  for(h=0; h < moinfo.nirreps; h++) {
 
75
    dpd_buf4_mat_irrep_row_init(&J, h);
 
76
    dpd_buf4_mat_irrep_row_init(&K, h);
 
77
 
 
78
    for(pq=0; pq < J.params->rowtot[h]; pq++) {
 
79
      dpd_buf4_mat_irrep_row_rd(&J, h, pq);
 
80
 
 
81
      for(Gr=0; Gr < moinfo.nirreps; Gr++) {
 
82
        Gs = h^Gr;
 
83
 
 
84
        nrows = sopi[Gr];
 
85
        ncols = mopi[Gs];
 
86
        nlinks = sopi[Gs];
 
87
        J_rs = J.col_offset[h][Gr];
 
88
        if(nrows && ncols && nlinks)
 
89
          C_DGEMM('n','n',nrows,ncols,nlinks,1.0,&J.matrix[h][0][J_rs],nlinks,C[Gs][0],
 
90
                  nlinks,0.0,TMP[0],nso);
 
91
 
 
92
        nrows = mopi[Gr];
 
93
        ncols = mopi[Gs];
 
94
        nlinks = sopi[Gr];
 
95
        K_rs = K.col_offset[h][Gr];
 
96
        if(nrows && ncols && nlinks)
 
97
          C_DGEMM('t','n',nrows,ncols,nlinks,1.0,C[Gr][0],nlinks,TMP[0],nso,
 
98
                  0.0,&K.matrix[h][0][K_rs],ncols);
 
99
      }
 
100
 
 
101
      dpd_buf4_mat_irrep_row_wrt(&K, h, pq);
 
102
    }
 
103
 
 
104
    dpd_buf4_mat_irrep_row_close(&J, h);
 
105
    dpd_buf4_mat_irrep_row_close(&K, h);
 
106
  }
 
107
  dpd_buf4_close(&K);
 
108
  dpd_buf4_close(&J);
 
109
 
 
110
  psio_close(PSIF_SO_PRESORT, 0); /* delete the presorted integrals */
 
111
  psio_open(CC_TEI_HALFT1, 0); /* transposed half-transformed integrals */
 
112
 
 
113
  dpd_buf4_init(&K, CC_TEI_HALFT0, 0, 3, 5, 3, 8, 0, "Half-Transformed Ints (pq,ij)");
 
114
  dpd_buf4_sort(&K, CC_TEI_HALFT1, rspq, 8, 3, "Half-Transformed Ints (ij,pq)");
 
115
  dpd_buf4_close(&K);
 
116
 
 
117
  psio_close(CC_TEI_HALFT0, 0); /* delete the original half-transformed integrals */
 
118
 
 
119
  dpd_buf4_init(&J, CC_TEI_HALFT1, 0, 8, 0, 8, 3, 0, "Half-Transformed Ints (ij,pq)");
 
120
  dpd_buf4_init(&K, CC_MISC, 0, 8, 5, 8, 8, 0, "MO Ints (ij,kl)"); /* place-holder buffer */
 
121
  for(h=0; h < moinfo.nirreps; h++) {
 
122
    dpd_buf4_mat_irrep_row_init(&J, h);
 
123
    buffer = init_array(K.params->coltot[h]);
 
124
 
 
125
    for(pq=0; pq < J.params->rowtot[h]; pq++) {
 
126
      if(occ[] && occ[])
 
127
        p = J.params->roworb[h][pq][0];
 
128
      q = J.params->roworb[h][pq][1];
 
129
 
 
130
      dpd_buf4_mat_irrep_row_rd(&J, h, pq);
 
131
 
 
132
      for(Gr=0; Gr < moinfo.nirreps; Gr++) {
 
133
        Gs = h^Gr;
 
134
 
 
135
        nrows = sopi[Gr];
 
136
        ncols = mopi[Gs];
 
137
        nlinks = sopi[Gs];
 
138
        J_rs = J.col_offset[h][Gr];
 
139
        if(nrows && ncols && nlinks)
 
140
          C_DGEMM('n','n',nrows,ncols,nlinks,1.0,&J.matrix[h][0][J_rs],nlinks,C[Gs][0],
 
141
                  nlinks,0.0,TMP[0],nso);
 
142
 
 
143
        nrows = mopi[Gr];
 
144
        ncols = mopi[Gs];
 
145
        nlinks = sopi[Gr];
 
146
        K_rs = K.col_offset[h][Gr];
 
147
        if(nrows && ncols && nlinks)
 
148
          C_DGEMM('t','n',nrows,ncols,nlinks,1.0,C[Gr][0],nlinks,TMP[0],nso,
 
149
                  0.0,&buffer[K_rs],ncols);
 
150
 
 
151
        /* dump the result straight into the distribute function */
 
152
      }
 
153
      for(rs=0; rs < J.params->coltot[h]; rs++) {
 
154
        r = J.params->colworb[h][rs][0];
 
155
        s = J.params->colworb[h][rs][1];
 
156
 
 
157
 
 
158
      }
 
159
    }
 
160
 
 
161
    dpd_buf4_mat_irrep_row_close(&J, h);
 
162
    dpd_buf4_mat_irrep_row_close(&K, h);
 
163
  }
 
164
  dpd_buf4_print(&K, outfile, 1);
 
165
  dpd_buf4_close(&K);
 
166
  dpd_buf4_close(&J);
 
167
 
 
168
  psio_close(CC_TEI_HALFT1, 0); /* delete the transposed half-transformed integrals */
 
169
 
 
170
  free_block(TMP);
 
171
  for(h=0; h < moinfo.nirreps; h++)
 
172
    free_block(C[h]);
 
173
  free(C);
 
174
 
 
175
  free(pitz2dpd_occ);
 
176
  free(pitz2dpd_vir);
 
177
}