← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/cints/Tools/compute_scf_opdm.c

  • 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/compute_scf_opdm.c
1
 
#include<stdio.h>
2
 
#include<string.h>
3
 
#include<libciomr/libciomr.h>
4
 
#include<libchkpt/chkpt.h>
5
 
#include<math.h>
6
 
#include<libint/libint.h>
7
 
 
8
 
#include"defines.h"
9
 
#define EXTERN
10
 
#include"global.h"
11
 
#include"prints.h"
12
 
#include"small_fns.h"
13
 
 
14
 
#define VIRT_SHELL -1000000       /* This should cause program to segfault if arrays are overrun */
15
 
 
16
 
void compute_scf_opdm()
17
 
{
18
 
  int i,j,k,l,mo,open_1st,openshell_count,closedmos;
19
 
  int sh_i, sh_j, ioffset, joffset, nao_i, nao_j;
20
 
  int *symblk_offset;                           /* The pointer to the first MO in the symmetry block */
21
 
  int *mo2moshell;                              /* array maping MO number to the mo shell */
22
 
  double *occ;                                  /* occupation vector - num_so long */
23
 
  double *shell_occ;                            /* occupation numbers for MO shells */
24
 
  double **mo_lagr;                                /* MO lagrangian */
25
 
  double tmp, max_elem;
26
 
  struct shell_pair* sp;
27
 
  
28
 
  /*--- Compute offsets of symmetry blocks ---*/
29
 
  symblk_offset = init_int_array(Symmetry.nirreps);
30
 
  symblk_offset[0] = 0;
31
 
  for(i=0;i<Symmetry.nirreps-1;i++)
32
 
    symblk_offset[i+1] = symblk_offset[i] + MOInfo.orbspi[i];
33
 
 
34
 
  /*----------------------------------------------------
35
 
    Compute occupation vector for spin-restricted cases
36
 
   ----------------------------------------------------*/
37
 
  if (UserOptions.reftype != uhf) {
38
 
    occ = init_array(MOInfo.num_mo);
39
 
    mo2moshell = init_int_array(MOInfo.num_mo);
40
 
    closedmos = (MOInfo.num_moshells != MOInfo.num_openmoshells);
41
 
    shell_occ = init_array(MOInfo.num_moshells);
42
 
    if (closedmos > 0)
43
 
      shell_occ[0] = 2.0;
44
 
    if (UserOptions.reftype != twocon) {        /* single-reference case */
45
 
      openshell_count = (closedmos > 0);
46
 
      for (i=0;i<Symmetry.nirreps;i++) {
47
 
        mo = symblk_offset[i];
48
 
        for (j=0;j<MOInfo.clsdpi[i];j++) {
49
 
          mo2moshell[mo] = 0;
50
 
          occ[mo] = 2.0;
51
 
          mo++;
52
 
        }
53
 
        if (MOInfo.openpi[i] > 0) {
54
 
          for (j=0;j<MOInfo.openpi[i];j++) {
55
 
            mo2moshell[mo] = openshell_count;
56
 
            occ[mo] = 1.0;
57
 
            mo++;
58
 
          }
59
 
          shell_occ[openshell_count] = 1.0;
60
 
          openshell_count++;
61
 
        }
62
 
        for (j=0;j<MOInfo.orbspi[i]-MOInfo.clsdpi[i]-MOInfo.openpi[i];j++) {
63
 
          mo2moshell[mo] = VIRT_SHELL;
64
 
          mo++;
65
 
        }
66
 
      }
67
 
    }
68
 
    else {              /* TCSCF for closed shells */
69
 
      openshell_count = (closedmos > 0);
70
 
      k = 0;
71
 
      for (i=0;i<Symmetry.nirreps;i++) {
72
 
        mo = symblk_offset[i];
73
 
        for (j=0;j<MOInfo.clsdpi[i];j++) {
74
 
          mo2moshell[mo] = 0;
75
 
          occ[mo] = 2.0;
76
 
          mo++;
77
 
        }
78
 
        if (MOInfo.openpi[i] > 0) {
79
 
          for (j=0;j<MOInfo.openpi[i];j++) {
80
 
            mo2moshell[mo] = openshell_count;
81
 
            occ[mo] = MOInfo.tcscf_occ[k];
82
 
            mo++;
83
 
          }
84
 
          shell_occ[openshell_count] = MOInfo.tcscf_occ[k];
85
 
          openshell_count++;
86
 
          k++;
87
 
        }
88
 
        for (j=0;j<MOInfo.orbspi[i]-MOInfo.clsdpi[i]-MOInfo.openpi[i];j++) {
89
 
          mo2moshell[mo] = VIRT_SHELL;
90
 
          mo++;
91
 
        }
92
 
      }
93
 
    }
94
 
  }  /*--- Done with occupations ---*/
95
 
 
96
 
 
97
 
  /*-----------------------------------
98
 
    Allocate and compute total density
99
 
   -----------------------------------*/
100
 
  Dens = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
101
 
  if (UserOptions.reftype != uhf) {
102
 
    for(i=0;i<BasisSet.num_ao;i++)
103
 
      for(j=0;j<=i;j++) {
104
 
        tmp = 0.0;
105
 
        for(k=0;k<Symmetry.nirreps;k++) {
106
 
          mo = symblk_offset[k];
107
 
          for(l=0; l < (MOInfo.clsdpi[k]+MOInfo.openpi[k]); l++) {
108
 
            tmp += occ[mo]*MOInfo.scf_evec[0][mo][i]*MOInfo.scf_evec[0][mo][j];
109
 
            mo++;
110
 
          }
111
 
        }
112
 
        Dens[j][i] = Dens[i][j] = tmp;
113
 
      }
114
 
  }
115
 
  else {
116
 
    /* UHF case : Alpha and Beta eigenvectors */
117
 
    for(i=0;i<BasisSet.num_ao;i++)
118
 
      for(j=0;j<=i;j++) {
119
 
        tmp = 0.0;
120
 
        for(k=0;k<Symmetry.nirreps;k++) {
121
 
          mo = symblk_offset[k];
122
 
          /*--- Doubly-occupied MOs : Alpha and Beta ---*/
123
 
          for(l=0; l < MOInfo.clsdpi[k]; l++) {
124
 
            tmp += MOInfo.scf_evec[0][mo][i]*
125
 
                   MOInfo.scf_evec[0][mo][j]
126
 
                +  MOInfo.scf_evec[1][mo][i]*
127
 
                   MOInfo.scf_evec[1][mo][j];
128
 
            mo++;
129
 
          }
130
 
          /*--- Singly-occupied MOs : Alpha only ---*/
131
 
          for(l=0; l < MOInfo.openpi[k]; l++) {
132
 
            tmp += MOInfo.scf_evec[0][mo][i]*
133
 
                   MOInfo.scf_evec[0][mo][j];
134
 
            mo++;
135
 
          }
136
 
        }
137
 
        Dens[j][i] = Dens[i][j] = tmp;
138
 
      }
139
 
  }
140
 
 
141
 
  /*---------------------------------
142
 
    Find maximal elements of Dens in
143
 
    each shell block
144
 
   ---------------------------------*/
145
 
  if (BasisSet.shell_pairs)
146
 
    for(sh_i=0;sh_i<BasisSet.num_shells;sh_i++) {
147
 
      ioffset = BasisSet.shells[sh_i].fao-1;
148
 
      nao_i = ioff[BasisSet.shells[sh_i].am];
149
 
      for(sh_j=0;sh_j<BasisSet.num_shells;sh_j++) {
150
 
        sp = &(BasisSet.shell_pairs[sh_i][sh_j]);
151
 
        joffset = BasisSet.shells[sh_j].fao-1;
152
 
        nao_j = ioff[BasisSet.shells[sh_j].am];
153
 
        max_elem = -1.0;
154
 
        for(i=0;i<nao_i;i++)
155
 
          for(j=0;j<nao_j;j++) {
156
 
            tmp = Dens[ioffset+i][joffset+j];
157
 
            tmp = fabs(tmp);
158
 
            if (max_elem < tmp)
159
 
              max_elem = tmp;
160
 
          }
161
 
        sp->Dmax = max_elem;
162
 
      }
163
 
    }
164
 
  
165
 
  /*--------------------------------
166
 
    Compute energy-weighted density
167
 
   --------------------------------*/
168
 
  Lagr = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
169
 
  if (UserOptions.reftype == rhf) {
170
 
    /*--- In closed-shell case compute energy-weighted density from SCF eigenvalues ---*/
171
 
    for(i=0;i<BasisSet.num_ao;i++)
172
 
      for(j=0;j<=i;j++) {
173
 
        tmp = 0.0;
174
 
        for(k=0;k<Symmetry.nirreps;k++) {
175
 
          mo = symblk_offset[k];
176
 
          for(l=0; l < (MOInfo.clsdpi[k]+MOInfo.openpi[k]); l++) {
177
 
            tmp += occ[mo]*MOInfo.scf_evals[0][mo]*
178
 
                   MOInfo.scf_evec[0][mo][i]*
179
 
                   MOInfo.scf_evec[0][mo][j];
180
 
            mo++;
181
 
          }
182
 
        }
183
 
        Lagr[j][i] = Lagr[i][j] = tmp;
184
 
      }
185
 
  }
186
 
  else if (UserOptions.reftype == uhf) {
187
 
    /*--- Use both alpha and beta eigenvalues in UHF case ---*/
188
 
    for(i=0;i<BasisSet.num_ao;i++)
189
 
      for(j=0;j<=i;j++) {
190
 
        tmp = 0.0;
191
 
        for(k=0;k<Symmetry.nirreps;k++) {
192
 
          mo = symblk_offset[k];
193
 
          /*--- Doubly-occupied MOs : Alpha and Beta ---*/
194
 
          for(l=0; l < MOInfo.clsdpi[k]; l++) {
195
 
            tmp += MOInfo.scf_evals[0][mo]*
196
 
                   MOInfo.scf_evec[0][mo][i]*
197
 
                   MOInfo.scf_evec[0][mo][j]
198
 
                +  MOInfo.scf_evals[1][mo]*
199
 
                   MOInfo.scf_evec[1][mo][i]*
200
 
                   MOInfo.scf_evec[1][mo][j];
201
 
            mo++;
202
 
          }
203
 
          /*--- Singly-occupied MOs : Alpha only ---*/
204
 
          for(l=0; l < MOInfo.openpi[k]; l++) {
205
 
            tmp += MOInfo.scf_evals[0][mo]*
206
 
                   MOInfo.scf_evec[0][mo][i]*
207
 
                   MOInfo.scf_evec[0][mo][j];
208
 
            mo++;
209
 
          }
210
 
        }
211
 
        Lagr[j][i] = Lagr[i][j] = tmp;
212
 
      }
213
 
  }
214
 
  else if (UserOptions.reftype == rohf || UserOptions.reftype == twocon) {
215
 
    /*--- In a restricted open-shell case just transform the lagrangian to AO basis ---*/
216
 
    mo_lagr = chkpt_rd_lagr();
217
 
    for(i=0;i<BasisSet.num_ao;i++)
218
 
      for(j=0;j<=i;j++) {
219
 
        tmp = 0.0;
220
 
        for(k=0;k<MOInfo.num_mo;k++)
221
 
          for(l=0;l<MOInfo.num_mo;l++)
222
 
            tmp += mo_lagr[k][l]*MOInfo.scf_evec[0][k][i]*MOInfo.scf_evec[0][l][j];
223
 
        Lagr[j][i] = Lagr[i][j] = tmp;
224
 
      }
225
 
  }
226
 
 
227
 
  /*----------------------------------------------
228
 
    Compute MO shell densities for ROHF and TCSCF
229
 
   ----------------------------------------------*/
230
 
  if (UserOptions.reftype == rohf || UserOptions.reftype == twocon) {
231
 
    ShDens = (double ***) malloc(sizeof(double **)*MOInfo.num_moshells);
232
 
    for(i=0;i<MOInfo.num_moshells;i++)
233
 
      ShDens[i] = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
234
 
    for(i=0;i<BasisSet.num_ao;i++)
235
 
      for(j=0;j<=i;j++) {
236
 
        for(k=0;k<Symmetry.nirreps;k++) {
237
 
          mo = symblk_offset[k];
238
 
          for(l=0; l < (MOInfo.clsdpi[k]+MOInfo.openpi[k]); l++) {
239
 
            ShDens[mo2moshell[mo]][j][i] = (ShDens[mo2moshell[mo]][i][j] +=
240
 
                                                MOInfo.scf_evec[0][mo][i]*MOInfo.scf_evec[0][mo][j]);
241
 
            mo++;
242
 
          }
243
 
        }
244
 
      }
245
 
    
246
 
    /* Check if the total density is a sum of shell densities - REMOVE AFTER DONE TESTING */
247
 
    for(i=0;i<BasisSet.num_ao;i++)
248
 
      for(j=0;j<BasisSet.num_ao;j++) {
249
 
        tmp = 0.0;
250
 
        for(k=0;k<MOInfo.num_moshells;k++)
251
 
          tmp += ShDens[k][i][j]*shell_occ[k];
252
 
        if (fabs(tmp - Dens[i][j]) > 1.0e-12)
253
 
          punt("Total density is not a sum of shell densities");
254
 
      }
255
 
 
256
 
    /*--- Do this dirty trick because some people decided to keep
257
 
          open-shells of diff irreps separate ---*/
258
 
    /*--- form a single open-shell if rohf and high-spin ---*/
259
 
/*    if (UserOptions.reftype == rohf) {
260
 
      i = (closedmos > 0) ? 1 : 0;
261
 
      for(j=i+1;j<MOInfo.num_moshells;j++) {
262
 
        for(k=0;k<BasisSet.num_ao;k++)
263
 
          for(l=0;l<BasisSet.num_ao;l++)
264
 
            ShDens[i][k][l] += ShDens[j][k][l];
265
 
        free_block(ShDens[j]);
266
 
      }
267
 
      BasisSet.num_moshells = 1 + i;
268
 
    }*/
269
 
 
270
 
    if (UserOptions.print_lvl >= PRINT_OPDM) {
271
 
      for(i=0;i<MOInfo.num_moshells;i++) {
272
 
        fprintf(outfile,"  -Density matrix for shell %d in AO basis :\n",i);
273
 
        print_mat(ShDens[i],BasisSet.num_ao,BasisSet.num_ao,outfile);
274
 
        fprintf(outfile,"\n\n");
275
 
      }
276
 
    }
277
 
  }
278
 
 
279
 
  print_opdm();
280
 
 
281
 
                /* Cleaning up */
282
 
 
283
 
  free(symblk_offset);
284
 
  if (UserOptions.reftype != uhf) {
285
 
    free(occ);
286
 
    free(mo2moshell);
287
 
    free(shell_occ);
288
 
  }
289
 
  if (UserOptions.reftype == rohf || UserOptions.reftype == twocon)
290
 
    free_block(mo_lagr);
291
 
  
292
 
  return;
293
 
}