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

« back to all changes in this revision

Viewing changes to src/bin/cscf/dmatuhf.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:
 
1
/*! \file
 
2
    \ingroup CSCF
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
/* $Log$
 
6
 * Revision 1.6  2004/05/03 04:32:40  crawdad
 
7
 * Major mods based on merge with stable psi-3-2-1 release.  Note that this
 
8
 * version has not been fully tested and some scf-optn test cases do not run
 
9
 * correctly beccause of changes in mid-March 2004 to optking.
 
10
 * -TDC
 
11
 *
 
12
/* Revision 1.5.8.1  2004/04/10 19:41:32  crawdad
 
13
/* Fixed the DIIS code for UHF cases.  The new version uses the Pulay scheme of
 
14
/* building the error vector in the AO basis as FDS-SDF, followed by xformation
 
15
/* into the orthogonal AO basis.   This code converges faster for test cases
 
16
/* like cc8, but fails for linearly dependent basis sets for unknown reasons.
 
17
/* -TDC
 
18
/*
 
19
/* Revision 1.5  2002/04/03 02:06:01  janssen
 
20
/* Finish changes to use new include paths for libraries.
 
21
/*
 
22
/* Revision 1.4  2000/12/05 19:40:03  sbrown
 
23
/* Added Unrestricted Kohn-Sham DFT.
 
24
/*
 
25
/* Revision 1.3  2000/06/22 22:15:00  evaleev
 
26
/* Modifications for KS DFT. Reading in XC Fock matrices and XC energy in formg_direct need to be uncommented (at present those are not produced by CINTS yet).
 
27
/*
 
28
/* Revision 1.2  2000/06/02 13:32:15  kenny
 
29
/*
 
30
/*
 
31
/* Added dynamic integral accuracy cutoffs for direct scf.  Added a few global
 
32
/* variables.  Added keyword 'dyn_acc'; true--use dynamic cutoffs.  Use of
 
33
/* 'dconv' and 'delta' to keep track of density convergence somewhat awkward,
 
34
/* but avoids problems when accuracy is switched and we have to wipe out density
 
35
/* matrices.  Also added error message and exit if direct rohf singlet is
 
36
/* attempted since it doesn't work.
 
37
/* --Joe Kenny
 
38
/*
 
39
/* Revision 1.1.1.1  2000/02/04 22:52:33  evaleev
 
40
/* Started PSI 3 repository
 
41
/*
 
42
/* Revision 1.2  1999/11/17 19:40:46  evaleev
 
43
/* Made all the adjustments necessary to have direct UHF working. Still doesn't work though..
 
44
/*
 
45
/* Revision 1.1  1999/11/02 23:55:56  sbrown
 
46
/* Shawn Brown - (11/2/99) Modified to the code in a few major ways.
 
47
/*
 
48
/* 1.  Added the capability to do UHF.  All of the features available with the
 
49
/* other refrences have been added for UHF.
 
50
/*
 
51
/* 2.  For UHF, I had to alter the structure of file30. (See cleanup.c for a
 
52
/* map)  This entailed adding a pointer array right after the header in the SCF
 
53
/* section of file30 that pointed to all of the data for the SCF caclulation.
 
54
/* Functions were added to libfile30 to account for this and they are
 
55
/* incorporated in this code.
 
56
/*
 
57
/* 3.  Updated and fixed all of the problems associated with my previous
 
58
/* guessing code.  The code no longer uses OPENTYPE to specify the type of
 
59
/* occupation.  The keword REFERENCE and MULTP can now be used to indicate any
 
60
/* type of calculation.  (e.g. ROHF with MULTP of 1 is an open shell singlet
 
61
/* ROHF calculation)  This code was moved to occ_fun.c.  The code can also
 
62
/* guess at any multplicity in a highspin case, provided enough electrons.
 
63
/*
 
64
/* Revision 1.1.1.1  1999/04/12 16:59:25  evaleev
 
65
/* Added a version of CSCF that can work with CINTS.
 
66
/* -Ed
 
67
/*
 
68
 * Revision 1.1  1991/06/15  20:22:21  seidl
 
69
 * Initial revision
 
70
 * */
 
71
 
 
72
#define EXTERN
 
73
#include <libpsio/psio.h>
 
74
#include "includes.h"
 
75
#include "common.h"
 
76
 
 
77
namespace psi { namespace cscf {
 
78
 
 
79
void dmatuhf()
 
80
{
 
81
  int i,j,k,l,ij,n,m,jj,kk;
 
82
  int max, off, ntri;
 
83
  int ndocc;
 
84
  int nn;
 
85
  double ptemp,ctmp;
 
86
  struct symm *s;
 
87
  struct spin *sp;
 
88
  double *dmat;
 
89
  double **cmat;
 
90
   
 
91
  for(m = 0; m < 2; m++){
 
92
    for (l=0; l < num_ir ; l++) {
 
93
      sp = &spin_info[m];
 
94
      if (n=scf_info[l].num_so) {
 
95
        ndocc = sp->scf_spin[l].noccup;
 
96
        for (i=ij=0; i < n ; i++ ) {
 
97
          /*--------------------------------------
 
98
 
 
99
          OFF DIAGONAL ELEMENTS
 
100
 
 
101
          -------------------------------------*/
 
102
 
 
103
          for (j=0; j <i; j++,ij++) {
 
104
            ptemp=0.0;
 
105
            for (k=0; k < ndocc ; k++)
 
106
              ptemp += 2.0*sp->scf_spin[l].cmat[i][k]
 
107
                *sp->scf_spin[l].cmat[j][k];
 
108
            sp->scf_spin[l].dpmat[ij] = ptemp - sp->scf_spin[l].pmat[ij];
 
109
            sp->scf_spin[l].pmat[ij] = ptemp; 
 
110
          }
 
111
          /*----------------------------------
 
112
 
 
113
          DIAGONAL ELEMENTS
 
114
 
 
115
          --------------------------------*/
 
116
          ptemp = 0.0;
 
117
          for (k=0; k < ndocc ; k++) {
 
118
            ctmp=sp->scf_spin[l].cmat[i][k];
 
119
            ptemp += ctmp*ctmp;
 
120
          }
 
121
          sp->scf_spin[l].dpmat[ij] = ptemp - sp->scf_spin[l].pmat[ij];
 
122
          sp->scf_spin[l].pmat[ij] = ptemp;
 
123
          ij++;
 
124
        } 
 
125
        if(print & 4) {
 
126
          fprintf(outfile,
 
127
                  "\nSpin case %d density matrix for irrep %s",m,scf_info[l].irrep_label);
 
128
          print_array(spin_info[m].scf_spin[l].pmat,
 
129
                      scf_info[l].num_so,outfile);
 
130
        }
 
131
      }
 
132
    }
 
133
  }
 
134
   
 
135
  for(l=0;l < num_ir; l++){
 
136
    if(nn=scf_info[l].num_so) {
 
137
      for(ij=0;ij<ioff[nn];ij++) {
 
138
        ptemp = spin_info[0].scf_spin[l].pmat[ij] +
 
139
          spin_info[1].scf_spin[l].pmat[ij];
 
140
        scf_info[l].dpmat[ij] = ptemp - scf_info[l].pmat[ij];
 
141
        scf_info[l].pmat[ij] = ptemp;
 
142
      }
 
143
    }
 
144
  }
 
145
 
 
146
  /*-----------------------
 
147
    Handle direct SCF here
 
148
    -----------------------*/
 
149
  if(direct_scf) {
 
150
    /*decide what accuracy to request for direct_scf*/
 
151
    if (dyn_acc) {
 
152
      if((iter<30)&&(tight_ints==0)&&(delta>1.0E-5)) {
 
153
        eri_cutoff=1.0E-6;
 
154
      }
 
155
      if((tight_ints==0)&&(delta<=1.0E-5)){
 
156
        fprintf(outfile,"  Switching to full integral accuracy\n");
 
157
        acc_switch=1;
 
158
        tight_ints=1;
 
159
        eri_cutoff=1.0E-14;
 
160
      }
 
161
    }
 
162
 
 
163
    psio_open(itapDSCF, PSIO_OPEN_NEW);
 
164
    psio_write_entry(itapDSCF, "Integrals cutoff", (char *) &eri_cutoff, sizeof(double));
 
165
 
 
166
    ntri = nbasis*(nbasis+1)/2;
 
167
    dmat = init_array(ntri);
 
168
 
 
169
    /*--- Get full dpmata ---*/
 
170
    for(i=0;i<num_ir;i++) {
 
171
      max = scf_info[i].num_so;
 
172
      off = scf_info[i].ideg;
 
173
      for(j=0;j<max;j++) {
 
174
        jj = j + off;
 
175
        for(k=0;k<=j;k++) {
 
176
          kk = k + off;
 
177
          if(acc_switch) {
 
178
            dmat[ioff[jj]+kk] = spin_info[0].scf_spin[i].pmat[ioff[j]+k];
 
179
            spin_info[0].scf_spin[i].dpmat[ioff[j]+k] = 0.0;
 
180
          }
 
181
          else {
 
182
            dmat[ioff[jj]+kk] = spin_info[0].scf_spin[i].dpmat[ioff[j]+k];
 
183
          }
 
184
        }
 
185
      }
 
186
    }
 
187
    psio_write_entry(itapDSCF, "Difference Alpha Density", (char *) dmat, sizeof(double)*ntri);
 
188
 
 
189
    /*--- Get full dpmatb ---*/
 
190
    for(i=0;i<num_ir;i++) {
 
191
      max = scf_info[i].num_so;
 
192
      off = scf_info[i].ideg;
 
193
      for(j=0;j<max;j++) {
 
194
        jj = j + off;
 
195
        for(k=0;k<=j;k++) {
 
196
          kk = k + off;
 
197
          if(acc_switch) {
 
198
            dmat[ioff[jj]+kk] = spin_info[1].scf_spin[i].pmat[ioff[j]+k];
 
199
            spin_info[1].scf_spin[i].dpmat[ioff[j]+k] = 0.0;
 
200
          }
 
201
          else {
 
202
            dmat[ioff[jj]+kk] = spin_info[1].scf_spin[i].dpmat[ioff[j]+k];
 
203
          }
 
204
        }
 
205
      }
 
206
    }
 
207
    psio_write_entry(itapDSCF, "Difference Beta Density", (char *) dmat, sizeof(double)*ntri);
 
208
     
 
209
    if (ksdft) {
 
210
      /* ---- Get Occupied Eigenvector matrix for DFT ----*/
 
211
         
 
212
      /*------ Get the Alpha Occupied Eigenvector --------*/
 
213
      ntri = nbfso*a_elec;
 
214
      cmat = block_matrix(nbfso,a_elec);
 
215
      for(i=l=0;i<num_ir;i++){
 
216
        max = spin_info[0].scf_spin[i].nclosed;
 
217
        off = scf_info[i].ideg;
 
218
        for(j=0;j<max;j++){
 
219
          for(k=0;k<scf_info[i].num_mo;k++) {
 
220
            kk = k + off;
 
221
            cmat[kk][l] = spin_info[0].scf_spin[i].cmat[k][j];
 
222
          }
 
223
          l++;
 
224
        }
 
225
      }
 
226
      /*fprintf(outfile,"\na_elec = %d",a_elec);
 
227
        fprintf(outfile,"\nAlpha Occupied Eigenvector from CSCF");
 
228
        print_mat(cmat,nbfso,a_elec,outfile);*/
 
229
         
 
230
      psio_write_entry(itapDSCF, "Number of MOs", 
 
231
                       (char *) &(nmo),sizeof(int)); 
 
232
        
 
233
      psio_write_entry(itapDSCF, "Number of Alpha DOCC",
 
234
                       (char *) &(a_elec),sizeof(int));
 
235
       
 
236
      psio_write_entry(itapDSCF, "Alpha Occupied SCF Eigenvector", 
 
237
                       (char *) &(cmat[0][0]),sizeof(double)*ntri);
 
238
      /*fprintf(outfile,"\na_elec = %d",a_elec);*/
 
239
      free_block(cmat);
 
240
      /*------ Get the Alpha Occupied Eigenvector --------*/
 
241
      ntri = nbfso*b_elec;
 
242
      cmat = block_matrix(nbfso,b_elec);
 
243
      for(i=l=0;i<num_ir;i++){
 
244
        max = spin_info[1].scf_spin[i].nclosed;
 
245
        off = scf_info[i].ideg;
 
246
        for(j=0;j<max;j++){
 
247
          for(k=0;k<scf_info[i].num_mo;k++) {
 
248
            kk = k + off;
 
249
            cmat[kk][l] = spin_info[1].scf_spin[i].cmat[k][j];
 
250
          }
 
251
          l++;
 
252
        }
 
253
      }
 
254
      /*fprintf(outfile,"\nBeta Occupied Eigenvector from CSCF");
 
255
        print_mat(cmat,nbfso,b_elec,outfile);*/
 
256
 
 
257
      psio_write_entry(itapDSCF, "Number of Beta DOCC",
 
258
                       (char *) &(b_elec),sizeof(int));
 
259
      psio_write_entry(itapDSCF, "Beta Occupied SCF Eigenvector", 
 
260
                       (char *) &(cmat[0][0]),sizeof(double)*ntri);
 
261
      free_block(cmat);
 
262
         
 
263
      /*--- Get full dpmata ---*/
 
264
      for(i=0;i<num_ir;i++) {
 
265
        max = scf_info[i].num_so;
 
266
        off = scf_info[i].ideg;
 
267
        for(j=0;j<max;j++) {
 
268
          jj = j + off;
 
269
          for(k=0;k<=j;k++) {
 
270
            kk = k + off;
 
271
            dmat[ioff[jj]+kk] = spin_info[0].scf_spin[i].pmat[ioff[j]+k];
 
272
          }
 
273
        }
 
274
      }
 
275
      psio_write_entry(itapDSCF, "Total Alpha Density", (char *) dmat, sizeof(double)*ntri);
 
276
 
 
277
      /*--- Get full dpmatb ---*/
 
278
      for(i=0;i<num_ir;i++) {
 
279
        max = scf_info[i].num_so;
 
280
        off = scf_info[i].ideg;
 
281
        for(j=0;j<max;j++) {
 
282
          jj = j + off;
 
283
          for(k=0;k<=j;k++) {
 
284
            kk = k + off;
 
285
            dmat[ioff[jj]+kk] = spin_info[1].scf_spin[i].pmat[ioff[j]+k];
 
286
          }
 
287
        }
 
288
      }
 
289
      psio_write_entry(itapDSCF, "Total Beta Density", (char *) dmat, sizeof(double)*ntri);
 
290
    }
 
291
    free(dmat);
 
292
    psio_close(itapDSCF, 1);
 
293
  }
 
294
 
 
295
  return;
 
296
}
 
297
 
 
298
}} // namespace psi::cscf