← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/cints/Tools/read_scf_opdm.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/cints/Tools/read_scf_opdm.cc
 
1
/*! \file read_scf_opdm.cc
 
2
    \ingroup CINTS
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
#include <cstdio>
 
6
#include <cmath>
 
7
#include <cstdlib>
 
8
#include <libciomr/libciomr.h>
 
9
#include <libpsio/psio.h>
 
10
#include <libint/libint.h>
 
11
 
 
12
#include "defines.h"
 
13
#define EXTERN
 
14
#include "global.h"
 
15
 
 
16
#define OFFDIAG_DENS_FACTOR 0.5
 
17
 
 
18
namespace psi { namespace CINTS {
 
19
 
 
20
void read_scf_opdm()
 
21
{ 
 
22
  int i, j, ij, nao_i, nao_j, sh_i, sh_j;
 
23
  int ioffset, joffset;
 
24
  int nstri = ioff[Symmetry.num_so];
 
25
  double *dens, *denso;
 
26
  double **sq_dens, **sq_denso;
 
27
  double **tmp_mat;
 
28
  double pfac;
 
29
  double max_elem, temp;
 
30
  struct shell_pair* sp;
 
31
 
 
32
  psio_open(IOUnits.itapDSCF, PSIO_OPEN_OLD);
 
33
  psio_read_entry(IOUnits.itapDSCF, "Integrals cutoff", (char *) &(UserOptions.cutoff), sizeof(double));
 
34
  dens = init_array(nstri);
 
35
  if (UserOptions.reftype == rohf || UserOptions.reftype == uhf)
 
36
    denso = init_array(nstri);
 
37
 
 
38
  if (UserOptions.reftype == uhf) {
 
39
      denso = init_array(nstri);
 
40
      psio_read_entry(IOUnits.itapDSCF, "Difference Alpha Density", (char *) dens, sizeof(double)*nstri);
 
41
      psio_read_entry(IOUnits.itapDSCF, "Difference Beta Density", (char *) denso, sizeof(double)*nstri);
 
42
      /*--- Form total and open-shell (spin) density matrices first ---*/
 
43
      for(i=0;i<nstri;i++) {
 
44
          temp = dens[i] + denso[i];
 
45
        denso[i] = dens[i] - denso[i];
 
46
        dens[i] = temp;
 
47
      }
 
48
  }
 
49
  else {
 
50
      psio_read_entry(IOUnits.itapDSCF, "Difference Density", (char *) dens, sizeof(double)*nstri);
 
51
      if (UserOptions.reftype == rohf) {
 
52
          denso = init_array(nstri);
 
53
          psio_read_entry(IOUnits.itapDSCF, "Difference Open-Shell Density", (char *) denso, sizeof(double)*nstri);
 
54
      }
 
55
  }
 
56
  psio_close(IOUnits.itapDSCF, 1);
 
57
 
 
58
  
 
59
  /*--------------------------
 
60
    Remove after done testing
 
61
   --------------------------*/
 
62
/*  fprintf(outfile,"  Total density matrix in SO basis :\n");
 
63
  print_array(dens,Symmetry.num_so,outfile);
 
64
  fprintf(outfile,"\n\n");
 
65
  if (UserOptions.reftype == rohf) {
 
66
    fprintf(outfile,"  Total open-shell density matrix in SO basis :\n");
 
67
    print_array(denso,Symmetry.num_so,outfile);
 
68
    fprintf(outfile,"\n\n");
 
69
  }*/
 
70
 
 
71
  
 
72
  /*------------------------
 
73
    convert to square forms
 
74
   ------------------------*/
 
75
  sq_dens = init_matrix(BasisSet.num_ao,BasisSet.num_ao);
 
76
  ij = 0;
 
77
  for(i=0;i<Symmetry.num_so;i++) {
 
78
    for(j=0;j<i;j++) {
 
79
      sq_dens[i][j] = sq_dens[j][i] =  OFFDIAG_DENS_FACTOR*dens[ij];
 
80
      ij++;
 
81
    }
 
82
    sq_dens[i][i] = dens[ij];
 
83
    ij++;
 
84
  }
 
85
  free(dens);
 
86
  if (UserOptions.reftype == rohf || UserOptions.reftype == uhf) {
 
87
    sq_denso = init_matrix(BasisSet.num_ao,BasisSet.num_ao);
 
88
    ij = 0;
 
89
    for(i=0;i<Symmetry.num_so;i++) {
 
90
      for(j=0;j<i;j++) {
 
91
        sq_denso[i][j] = sq_denso[j][i] = OFFDIAG_DENS_FACTOR*denso[ij];
 
92
        ij++;
 
93
      }
 
94
      sq_denso[i][i] = denso[ij];
 
95
      ij++;
 
96
    }
 
97
    free(denso);
 
98
  }
 
99
 
 
100
  /*----------------------
 
101
    transform to AO basis
 
102
   ----------------------*/
 
103
  tmp_mat = init_matrix(Symmetry.num_so,BasisSet.num_ao);
 
104
  mmult(sq_dens,0,Symmetry.usotao,0,tmp_mat,0,Symmetry.num_so,Symmetry.num_so,BasisSet.num_ao,0);
 
105
  mmult(Symmetry.usotao,1,tmp_mat,0,sq_dens,0,BasisSet.num_ao,Symmetry.num_so,BasisSet.num_ao,0);
 
106
  if (UserOptions.reftype == rohf || UserOptions.reftype == uhf) {
 
107
    mmult(sq_denso,0,Symmetry.usotao,0,tmp_mat,0,Symmetry.num_so,Symmetry.num_so,BasisSet.num_ao,0);
 
108
    mmult(Symmetry.usotao,1,tmp_mat,0,sq_denso,0,BasisSet.num_ao,Symmetry.num_so,BasisSet.num_ao,0);
 
109
  }
 
110
  free_matrix(tmp_mat,Symmetry.num_so);
 
111
 
 
112
 
 
113
  /*--------------------------
 
114
    Remove after done testing
 
115
   --------------------------*/
 
116
  /*fprintf(outfile,"  Total density matrix in AO basis :\n");
 
117
  print_mat(sq_dens,BasisSet.num_ao,BasisSet.num_ao,outfile);
 
118
  fprintf(outfile,"\n\n");
 
119
  if (UserOptions.reftype == rohf || UserOptions.reftype == uhf) {
 
120
    fprintf(outfile,"  Total open-shell density matrix in AO basis :\n");
 
121
    print_mat(sq_denso,BasisSet.num_ao,BasisSet.num_ao,outfile);
 
122
    fprintf(outfile,"\n\n");
 
123
    }*/
 
124
  
 
125
  /*--------------------------------
 
126
    transform to shell-blocked form
 
127
   --------------------------------*/
 
128
  switch (UserOptions.reftype) {
 
129
    case rohf:
 
130
      for(sh_i=0;sh_i<BasisSet.num_shells;sh_i++) {
 
131
        ioffset = BasisSet.shells[sh_i].fao-1;
 
132
        nao_i = ioff[BasisSet.shells[sh_i].am];
 
133
        for(sh_j=0;sh_j<BasisSet.num_shells;sh_j++) {
 
134
          sp = &(BasisSet.shell_pairs[sh_i][sh_j]);
 
135
          joffset = BasisSet.shells[sh_j].fao-1;
 
136
          nao_j = ioff[BasisSet.shells[sh_j].am];
 
137
          if (sp->dmato == NULL)
 
138
            sp->dmato = block_matrix(nao_i,nao_j);
 
139
 
 
140
          max_elem = sp->Dmax;
 
141
          for(i=0;i<nao_i;i++)
 
142
            for(j=0;j<nao_j;j++) {
 
143
              temp = sp->dmato[i][j] = sq_denso[ioffset+i][joffset+j];
 
144
              temp = fabs(temp);
 
145
              if (max_elem < temp)
 
146
                max_elem = temp;
 
147
            }
 
148
          sp->Dmax = max_elem;
 
149
        }
 
150
      }
 
151
      free_matrix(sq_denso,BasisSet.num_ao);
 
152
 
 
153
      for(sh_i=0;sh_i<BasisSet.num_shells;sh_i++) {
 
154
        ioffset = BasisSet.shells[sh_i].fao-1;
 
155
        nao_i = ioff[BasisSet.shells[sh_i].am];
 
156
        for(sh_j=0;sh_j<BasisSet.num_shells;sh_j++) {
 
157
          sp = &(BasisSet.shell_pairs[sh_i][sh_j]);
 
158
          joffset = BasisSet.shells[sh_j].fao-1;
 
159
          nao_j = ioff[BasisSet.shells[sh_j].am];
 
160
          if (sp->dmat == NULL)
 
161
            sp->dmat = block_matrix(nao_i,nao_j);
 
162
          max_elem = -1.0;
 
163
          for(i=0;i<nao_i;i++)
 
164
            for(j=0;j<nao_j;j++) {
 
165
              temp = sp->dmat[i][j] = sq_dens[ioffset+i][joffset+j];
 
166
              temp = fabs(temp);
 
167
              if (max_elem < temp)
 
168
                max_elem = temp;
 
169
            }
 
170
          sp->Dmax = MAX(sp->Dmax,max_elem);
 
171
        }
 
172
      }
 
173
      free_matrix(sq_dens,BasisSet.num_ao);
 
174
      break;
 
175
 
 
176
    case uhf:
 
177
      for(sh_i=0;sh_i<BasisSet.num_shells;sh_i++) {
 
178
        ioffset = BasisSet.shells[sh_i].fao-1;
 
179
        nao_i = ioff[BasisSet.shells[sh_i].am];
 
180
        for(sh_j=0;sh_j<BasisSet.num_shells;sh_j++) {
 
181
          sp = &(BasisSet.shell_pairs[sh_i][sh_j]);
 
182
          joffset = BasisSet.shells[sh_j].fao-1;
 
183
          nao_j = ioff[BasisSet.shells[sh_j].am];
 
184
          /*--- Check if can reuse space from dmata ---*/
 
185
          if (sp->dmato == NULL)
 
186
              if (sp->dmata == NULL)
 
187
                  sp->dmato = block_matrix(nao_i,nao_j);
 
188
              else {
 
189
                  sp->dmato = sp->dmata;
 
190
                  sp->dmata = NULL;
 
191
              }
 
192
          
 
193
          max_elem = sp->Dmax;
 
194
          for(i=0;i<nao_i;i++)
 
195
              for(j=0;j<nao_j;j++) {
 
196
                  temp = sp->dmato[i][j] = sq_denso[ioffset+i][joffset+j];
 
197
                  temp = fabs(temp);
 
198
                  if (max_elem < temp)
 
199
                      max_elem = temp;
 
200
              }
 
201
          sp->Dmax = max_elem;
 
202
        }
 
203
      }
 
204
      free_matrix(sq_denso,BasisSet.num_ao);
 
205
      
 
206
      for(sh_i=0;sh_i<BasisSet.num_shells;sh_i++) {
 
207
          ioffset = BasisSet.shells[sh_i].fao-1;
 
208
          nao_i = ioff[BasisSet.shells[sh_i].am];
 
209
          for(sh_j=0;sh_j<BasisSet.num_shells;sh_j++) {
 
210
              sp = &(BasisSet.shell_pairs[sh_i][sh_j]);
 
211
              joffset = BasisSet.shells[sh_j].fao-1;
 
212
              nao_j = ioff[BasisSet.shells[sh_j].am];
 
213
              /*--- Check if can reuse space from dmatb ---*/
 
214
              if (sp->dmat == NULL)
 
215
                  if (sp->dmatb == NULL)
 
216
                      sp->dmat = block_matrix(nao_i,nao_j);
 
217
                  else {
 
218
                      sp->dmat = sp->dmatb;
 
219
                      sp->dmatb = NULL;
 
220
                  }
 
221
              
 
222
              max_elem = sp->Dmax;
 
223
              for(i=0;i<nao_i;i++)
 
224
                  for(j=0;j<nao_j;j++) {
 
225
                      temp = sp->dmat[i][j] = sq_dens[ioffset+i][joffset+j];
 
226
                      temp = fabs(temp);
 
227
                      if (max_elem < temp)
 
228
                          max_elem = temp;
 
229
                  }
 
230
              sp->Dmax = MAX(sp->Dmax,max_elem);      
 
231
          }
 
232
      }
 
233
  
 
234
      free_matrix(sq_dens,BasisSet.num_ao);
 
235
      break;
 
236
 
 
237
    case rhf:
 
238
      for(sh_i=0;sh_i<BasisSet.num_shells;sh_i++) {
 
239
        ioffset = BasisSet.shells[sh_i].fao-1;
 
240
        nao_i = ioff[BasisSet.shells[sh_i].am];
 
241
        for(sh_j=0;sh_j<BasisSet.num_shells;sh_j++) {
 
242
          sp = &(BasisSet.shell_pairs[sh_i][sh_j]);
 
243
          joffset = BasisSet.shells[sh_j].fao-1;
 
244
          nao_j = ioff[BasisSet.shells[sh_j].am];
 
245
          if (sp->dmat == NULL)
 
246
            sp->dmat = block_matrix(nao_i,nao_j);
 
247
          max_elem = -1.0;
 
248
          for(i=0;i<nao_i;i++)
 
249
            for(j=0;j<nao_j;j++) {
 
250
              temp = sp->dmat[i][j] = sq_dens[ioffset+i][joffset+j];
 
251
              temp = fabs(temp);
 
252
              if (max_elem < temp)
 
253
                max_elem = temp;
 
254
            }
 
255
          sp->Dmax = max_elem;
 
256
        }
 
257
      }
 
258
      free_matrix(sq_dens,BasisSet.num_ao);
 
259
      break;
 
260
  }
 
261
 
 
262
  return;
 
263
}
 
264
 
 
265
};};