← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/mcscf/scf_read_so_tei.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/mcscf/scf_read_so_tei.cc
 
1
#include <iostream>
 
2
#include <cmath>
 
3
 
 
4
#include <psifiles.h>
 
5
#include <libiwl/iwl.h>
 
6
#include <libciomr/libciomr.h>
 
7
#include <libmoinfo/libmoinfo.h>
 
8
 
 
9
#include "scf.h"
 
10
#include "memory_manager.h"
 
11
 
 
12
#define MAX(i,j) ((i>j) ? i : j)
 
13
#define MIN(i,j) ((i>j) ? j : i)
 
14
#define INDEX(i,j) ((i>j) ? (ioff[(i)]+(j)) : (ioff[(j)]+(i)))
 
15
 
 
16
extern FILE* outfile;
 
17
 
 
18
namespace psi{ namespace mcscf{
 
19
 
 
20
void SCF::read_so_tei()
 
21
{
 
22
  generate_pairs();
 
23
 
 
24
  total_symmetric_block_size = INDEX(pairpi[0]-1,pairpi[0]-1)+1;
 
25
 
 
26
  size_t free_memory = mem->get_FreeMemory();
 
27
  
 
28
  // Determine the number of matrix elements of the PK (and K) matrix to hold in core
 
29
  if(reference == rhf){
 
30
    nin_core        = min(free_memory / sizeof(double),total_symmetric_block_size);
 
31
  }else{
 
32
    nin_core        = min(free_memory / (2 * sizeof(double)),total_symmetric_block_size);
 
33
  }
 
34
  if(nin_core != total_symmetric_block_size)
 
35
    out_of_core = true;
 
36
 
 
37
  size_t total_symmetric_pairs = pairpi[0];
 
38
 
 
39
  // Determine how to split the two-electron operators
 
40
  nbatch            = 0;
 
41
  size_t pq_incore  = 0;
 
42
  size_t pqrs_index = 0;
 
43
 
 
44
  batch_pq_min[0]   = 0;
 
45
  batch_pq_max[0]   = 0;
 
46
  batch_index_min[0]= 0;
 
47
  batch_index_max[0]= 0;
 
48
 
 
49
  for(size_t pq = 0; pq < total_symmetric_pairs; ++pq){
 
50
    // Increment counters
 
51
    if(pq_incore + pq + 1 > nin_core){
 
52
      // The batch is full. Save info.
 
53
      batch_pq_max[nbatch]         = pq;
 
54
      batch_pq_min[nbatch + 1]     = pq;
 
55
      batch_index_max[nbatch]      = pqrs_index;
 
56
      batch_index_min[nbatch + 1]  = pqrs_index;
 
57
      pq_incore = 0;
 
58
      nbatch++;
 
59
    }
 
60
    pq_incore  += pq + 1;
 
61
    pqrs_index += pq + 1;
 
62
  }
 
63
  if(batch_pq_max[nbatch] != total_symmetric_pairs){
 
64
    batch_pq_max[nbatch]     = total_symmetric_pairs;
 
65
    batch_index_max[nbatch]  = total_symmetric_block_size;
 
66
    nbatch++;
 
67
  }
 
68
 
 
69
  for(int batch = 0; batch < nbatch; ++batch){
 
70
    batch_size[batch] = batch_index_max[ batch] - batch_index_min[batch];
 
71
    fprintf(outfile,"\n  batch %3d pq = [%8ld,%8ld] index = [%16ld,%16ld]",
 
72
                     batch,
 
73
                     batch_pq_min[batch],batch_pq_max[batch],
 
74
                     batch_index_min[batch],batch_index_max[batch]);
 
75
  }
 
76
  fflush(outfile);
 
77
  
 
78
  // Allocate the PK matrix
 
79
  allocate1(double,PK,nin_core);
 
80
  for(size_t i=0; i < nin_core; i++)
 
81
    PK[i]    =0.0;
 
82
  fprintf(outfile,"\n\n  Allocated the PK matrix (%d elements) ",nin_core);
 
83
  fflush(outfile);      
 
84
  
 
85
  if(reference != rhf){
 
86
    // Allocate the K matrix
 
87
    allocate1(double,K,nin_core);
 
88
    for(size_t i=0; i < nin_core; i++)
 
89
      K[i]    =0.0;
 
90
    fprintf(outfile,"\n  Allocated the  K matrix (%d elements) ",nin_core);
 
91
    fflush(outfile);
 
92
  }
 
93
  
 
94
  if(reference == rhf)
 
95
    read_so_tei_form_PK();
 
96
  else
 
97
    read_so_tei_form_PK_and_K();
 
98
}
 
99
 
 
100
void SCF::read_so_tei_form_PK()
 
101
{
 
102
  fprintf(outfile,"\n  Reading the two-electron integrals to form PK ... ");
 
103
  fflush(outfile);
 
104
 
 
105
  for(int batch = 0; batch < nbatch; ++batch){
 
106
        fprintf(outfile,"\n  batch %3d ... ",batch);
 
107
        fflush(outfile);
 
108
    // Compute the minimum and maximum indices
 
109
    size_t min_index   = batch_index_min[batch];
 
110
    size_t max_index   = batch_index_max[batch];
 
111
    size_t buffer_size = max_index - min_index;
 
112
 
 
113
    for(size_t pqrs = 0; pqrs < buffer_size; ++pqrs)
 
114
      PK[pqrs] = 0.0;
 
115
 
 
116
    double value;
 
117
    size_t p,q,r,s,four_index;
 
118
    int ilsti,nbuf,fi,index;
 
119
    struct iwlbuf ERIIN;
 
120
    iwl_buf_init(&ERIIN,PSIF_SO_TEI,0.0,1,1);
 
121
    do{
 
122
      ilsti = ERIIN.lastbuf;
 
123
      nbuf  = ERIIN.inbuf;
 
124
      fi=0;
 
125
      for(index=0;index<nbuf;index++){
 
126
        // Compute the [pq] index for this pqrs combination
 
127
        p = (short int)abs((float)ERIIN.labels[fi]);
 
128
        q = ERIIN.labels[fi+1];
 
129
        r = ERIIN.labels[fi+2];
 
130
        s = ERIIN.labels[fi+3];
 
131
        value = ERIIN.values[index];
 
132
  
 
133
        if(pair_sym[p][q] == 0)
 
134
        {
 
135
          four_index = INDEX(pair[p][q],pair[r][s]);
 
136
          if((four_index >= min_index) && (four_index < max_index)){
 
137
             PK[four_index - min_index]    += value;
 
138
          }
 
139
        }
 
140
        if(pair_sym[p][r] == 0)
 
141
        {
 
142
          four_index = INDEX(pair[p][r],pair[q][s]);
 
143
          if((four_index >= min_index) && (four_index < max_index)){
 
144
            if((p == r) || (q == s))
 
145
              PK[four_index - min_index]   -= 0.5 * value;
 
146
            else
 
147
              PK[four_index - min_index]   -= 0.25 * value;
 
148
          }
 
149
        }
 
150
        if(pair_sym[p][s] == 0){
 
151
          four_index = INDEX(pair[p][s],pair[q][r]);
 
152
          if((four_index >= min_index) && (four_index < max_index)){
 
153
            if((p != q) && (r != s)){
 
154
              if((p == s) || (q == r))
 
155
                PK[four_index - min_index] -= 0.5 * value;
 
156
              else
 
157
                PK[four_index - min_index] -= 0.25 * value;
 
158
            }
 
159
          }
 
160
        }
 
161
        fi+=4;
 
162
      }
 
163
      if(!ilsti)
 
164
        iwl_buf_fetch(&ERIIN);
 
165
    } while(!ilsti);
 
166
    iwl_buf_close(&ERIIN,1);
 
167
 
 
168
    // Half the diagonal elements held in core
 
169
    for(int pq = batch_pq_min[batch]; pq < batch_pq_max[batch]; ++pq){
 
170
      PK[INDEX(pq,pq) - min_index] *= 0.5;
 
171
    }
 
172
 
 
173
    // Write the PK matrix to disk
 
174
    write_Raffanetti("PK",PK,batch);
 
175
 
 
176
    fprintf(outfile,"done.");
 
177
    fflush(outfile);
 
178
  }
 
179
  fprintf(outfile,"\n");
 
180
  fflush(outfile);
 
181
}
 
182
 
 
183
void SCF::read_so_tei_form_PK_and_K()
 
184
{
 
185
  fprintf(outfile,"\n  Reading the two-electron integrals to form PK and K ... ");
 
186
  fflush(outfile);
 
187
 
 
188
  for(int batch = 0; batch < nbatch; ++batch){
 
189
        fprintf(outfile,"\n  batch %3d ... ",batch);
 
190
        fflush(outfile);
 
191
    // Compute the minimum and maximum indices
 
192
    size_t min_index   = batch_index_min[batch];
 
193
    size_t max_index   = batch_index_max[batch];
 
194
    size_t buffer_size = max_index - min_index;
 
195
 
 
196
    for(size_t pqrs = 0; pqrs < buffer_size; ++pqrs){
 
197
      PK[pqrs] = 0.0; K[pqrs] = 0.0;
 
198
    }
 
199
 
 
200
    double value;
 
201
    size_t p,q,r,s,four_index;
 
202
    int ilsti,nbuf,fi,index;
 
203
 
 
204
    struct iwlbuf ERIIN;
 
205
    iwl_buf_init(&ERIIN,PSIF_SO_TEI,0.0,1,1);
 
206
    do{
 
207
      ilsti = ERIIN.lastbuf;
 
208
      nbuf  = ERIIN.inbuf;
 
209
      fi=0;
 
210
      for(index=0;index<nbuf;index++){
 
211
        // Compute the [pq] index for this pqrs combination
 
212
        p = (short int)abs((float)ERIIN.labels[fi]);
 
213
        q = ERIIN.labels[fi+1];
 
214
        r = ERIIN.labels[fi+2];
 
215
        s = ERIIN.labels[fi+3];
 
216
        value = ERIIN.values[index];
 
217
  
 
218
        if(pair_sym[p][q] == 0)
 
219
        {
 
220
          four_index = INDEX(pair[p][q],pair[r][s]);
 
221
          if((four_index >= min_index) && (four_index < max_index)){
 
222
             PK[four_index - min_index]    += value;
 
223
          }
 
224
        }
 
225
        if(pair_sym[p][r] == 0)
 
226
        {
 
227
          four_index = INDEX(pair[p][r],pair[q][s]);
 
228
          if((four_index >= min_index) && (four_index < max_index)){
 
229
            if((p == r) || (q == s)){
 
230
              PK[four_index - min_index]   -= 0.5 * value;
 
231
               K[four_index - min_index]   -= 0.5 * value;
 
232
            }else{
 
233
              PK[four_index - min_index]   -= 0.25 * value;
 
234
               K[four_index - min_index]   -= 0.25 * value;
 
235
            }
 
236
          }
 
237
        }
 
238
        if(pair_sym[p][s] == 0){
 
239
          four_index = INDEX(pair[p][s],pair[q][r]);
 
240
          if((four_index >= min_index) && (four_index < max_index)){
 
241
            if((p != q) && (r != s)){
 
242
              if((p == s) || (q == r)){
 
243
                PK[four_index - min_index] -= 0.5 * value;
 
244
                 K[four_index - min_index] -= 0.5 * value;
 
245
              }else{
 
246
                PK[four_index - min_index] -= 0.25 * value;
 
247
                 K[four_index - min_index] -= 0.25 * value;
 
248
              }
 
249
            }
 
250
          }
 
251
        }
 
252
        fi+=4;
 
253
      }
 
254
      if(!ilsti)
 
255
        iwl_buf_fetch(&ERIIN);
 
256
    } while(!ilsti);
 
257
    iwl_buf_close(&ERIIN,1);
 
258
 
 
259
    // Half the diagonal elements held in core
 
260
    for(int pq = batch_pq_min[batch]; pq < batch_pq_max[batch]; ++pq){
 
261
      PK[INDEX(pq,pq) - min_index] *= 0.5;
 
262
       K[INDEX(pq,pq) - min_index] *= 0.5;
 
263
    }
 
264
 
 
265
    // Write the PK matrix to disk
 
266
    write_Raffanetti("PK",PK,batch);
 
267
    write_Raffanetti("K",K,batch);
 
268
    
 
269
    fprintf(outfile,"done.");
 
270
    fflush(outfile);
 
271
  }
 
272
  fprintf(outfile,"\n");
 
273
  fflush(outfile);
 
274
}
 
275
 
 
276
}} /* End Namespaces */