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

« back to all changes in this revision

Viewing changes to src/bin/cscf/uhf_iter.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.4.3  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.4.2  2004/04/09 00:17:37  evaleev
 
20
/* Corrected dimensions of the matrix.
 
21
/*
 
22
/* Revision 1.5.4.1  2004/04/07 03:23:32  crawdad
 
23
/* Working to fix UHF-based DIIS.
 
24
/* -TDC
 
25
/*
 
26
/* Revision 1.5  2002/12/06 15:50:32  crawdad
 
27
/* Changed all exit values to PSI_RETURN_SUCCESS or PSI_RETURN_FAILURE as
 
28
/* necessary.  This is new for the PSI3 execution driver.
 
29
/* -TDC
 
30
/*
 
31
/* Revision 1.4  2000/12/05 19:40:04  sbrown
 
32
/* Added Unrestricted Kohn-Sham DFT.
 
33
/*
 
34
/* Revision 1.3  2000/10/13 19:51:22  evaleev
 
35
/* Cleaned up a lot of stuff in order to get CSCF working with the new "Mo-projection-capable" INPUT.
 
36
/*
 
37
/* Revision 1.2  2000/06/22 22:15:02  evaleev
 
38
/* 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).
 
39
/*
 
40
/* Revision 1.1.1.1  2000/02/04 22:52:34  evaleev
 
41
/* Started PSI 3 repository
 
42
/*
 
43
/* Revision 1.3  1999/11/17 19:40:47  evaleev
 
44
/* Made all the adjustments necessary to have direct UHF working. Still doesn't work though..
 
45
/*
 
46
/* Revision 1.2  1999/11/04 19:24:31  localpsi
 
47
/* STB (11/4/99) - Added the orb_mix feature which is equivalent to guess = mix
 
48
/* in G94 and also fixed restarting so that if you have different wavefuntions,
 
49
/* everything works.  Also if you specify no DOCC and SOCC and restart, if the
 
50
/* wavefunctions are different, it will guess again.
 
51
/*
 
52
/* Revision 1.1  1999/11/02 23:56:00  localpsi
 
53
/* Shawn Brown - (11/2/99) Modified to the code in a few major ways.
 
54
/*
 
55
/* 1.  Added the capability to do UHF.  All of the features available with the
 
56
/* other refrences have been added for UHF.
 
57
/*
 
58
/* 2.  For UHF, I had to alter the structure of file30. (See cleanup.c for a
 
59
/* map)  This entailed adding a pointer array right after the header in the SCF
 
60
/* section of file30 that pointed to all of the data for the SCF caclulation.
 
61
/* Functions were added to libfile30 to account for this and they are
 
62
/* incorporated in this code.
 
63
/*
 
64
/* 3.  Updated and fixed all of the problems associated with my previous
 
65
/* guessing code.  The code no longer uses OPENTYPE to specify the type of
 
66
/* occupation.  The keword REFERENCE and MULTP can now be used to indicate any
 
67
/* type of calculation.  (e.g. ROHF with MULTP of 1 is an open shell singlet
 
68
/* ROHF calculation)  This code was moved to occ_fun.c.  The code can also
 
69
/* guess at any multplicity in a highspin case, provided enough electrons.
 
70
/*
 
71
/* Revision 1.1.1.1  1999/04/12 16:59:28  evaleev
 
72
/* Added a version of CSCF that can work with CINTS.
 
73
/* -Ed
 
74
/*
 
75
 * Revision 1.5  1998/06/30  14:11:12  sbrown
 
76
 * *************************************************************
 
77
 * *Program Modification                                       *
 
78
 * *By: Shawn Brown                                            *
 
79
 * *Date: June 30, 1998                                        *
 
80
 * *Altered program to make a guess at occupations from the    *
 
81
 * *diagonalized core hamiltonian matrix.  Program can now     *
 
82
 * *make a guess at the beginning of the calculation or at     *
 
83
 * *or at every iteration.  Use the latter at your own risk.   *
 
84
 * *See man pages for details on new keywords.                 *
 
85
 * *************************************************************
 
86
 *
 
87
 * Revision 1.4  1995/07/21  17:37:15  psi
 
88
 * Made Jan 1st 1995 cscf the current accepted version of cscf.  Some
 
89
 * unidentified changes made after that date were causing problems.
 
90
 *
 
91
 * Revision 1.1  1991/06/15  20:22:40  seidl
 
92
 * Initial revision
 
93
 * */
 
94
 
 
95
#define EXTERN
 
96
#include "includes.h"
 
97
#include "common.h"
 
98
 
 
99
namespace psi { namespace cscf {
 
100
 
 
101
void uhf_iter()
 
102
 
 
103
{
 
104
  int i,j,l,m,t,ij;
 
105
  int nn,num_mo,newci;
 
106
  double cimax;
 
107
  double **scr;
 
108
  double **fock_c;
 
109
  double **fock_ct;
 
110
  double **ctrans;
 
111
  double tol = 1.0e-14;
 
112
  struct symm *s;
 
113
  struct spin *sp;
 
114
 
 
115
  diiser = 0.0;
 
116
  scr = block_matrix(nsfmax,nsfmax);
 
117
  fock_c = block_matrix(nsfmax,nsfmax);
 
118
  fock_ct = block_matrix(nsfmax,nsfmax);
 
119
  ctrans = block_matrix(nsfmax,nsfmax);
 
120
   
 
121
  /* and iterate */
 
122
 
 
123
  for (iter=0; iter < itmax ; ) {
 
124
    for(t = 0; t<2 ;t++){
 
125
      sp = &spin_info[t];
 
126
 
 
127
      if(print & 4) {
 
128
        for(m=0; m < num_ir ; m++) {
 
129
          if (nn=scf_info[m].num_so) {
 
130
            fprintf(outfile,
 
131
                    "\n%s gmat for irrep %s",sp->spinlabel,scf_info[m].irrep_label);
 
132
            print_array(sp->scf_spin[m].gmat,nn,outfile);
 
133
            if (ksdft) {
 
134
              fprintf(outfile,
 
135
                      "\n%s xcmat for irrep %s",sp->spinlabel,scf_info[m].irrep_label);
 
136
              print_array(sp->scf_spin[m].xcmat,nn,outfile);
 
137
            }
 
138
          }
 
139
        }
 
140
      }
 
141
           
 
142
      for (m=0; m < num_ir ; m++) {
 
143
        s = &scf_info[m];
 
144
               
 
145
        if (nn=s->num_so) {
 
146
 
 
147
          /*               
 
148
          if(!m && !t) {
 
149
            fprintf(outfile, "Fock matrix in top of uhf_iter:\n");
 
150
            print_array(sp->scf_spin[m].fock_pac, nn, outfile);
 
151
          }
 
152
          */
 
153
 
 
154
          /*  form fock matrix = h+g */
 
155
          add_arr(s->hmat,sp->scf_spin[m].gmat,
 
156
                  sp->scf_spin[m].fock_pac,ioff[nn]);
 
157
 
 
158
        }
 
159
      }
 
160
    }
 
161
 
 
162
    /*----------------------------------------------------
 
163
      In KS DFT case, Fock matrix doesn't include Fxc yet
 
164
      add them up only after the computation of energy
 
165
      ----------------------------------------------------*/
 
166
    ecalc(tol);
 
167
    if (ksdft) {
 
168
      /* now form f = h + g + fxc */
 
169
      /* it should be alright to use fock_pac as 2 arguments */
 
170
      for(t = 0; t<2 ;t++){
 
171
        sp = &spin_info[t];
 
172
        for (m=0; m < num_ir ; m++) {
 
173
          s = &(sp->scf_spin[m]);
 
174
          if (nn=scf_info[m].num_so)
 
175
            add_arr(s->fock_pac,s->xcmat,s->fock_pac,ioff[nn]);
 
176
          if(print & 4) {
 
177
            fprintf(outfile,"\n%s fock for irrep %s"
 
178
                    ,sp->spinlabel,scf_info[m].irrep_label);
 
179
            print_array(sp->scf_spin[m].fock_pac,nn,outfile);
 
180
          }
 
181
        }
 
182
      }
 
183
    }
 
184
       
 
185
    /* create new fock matrix in fock_pac or fock_eff */
 
186
    if(!diisflg) diis_uhf();
 
187
       
 
188
    for(t=0;t<2;t++){
 
189
      sp = &spin_info[t];
 
190
      for (m=0; m < num_ir ; m++) {
 
191
        s = &scf_info[m];
 
192
        if (nn=s->num_so) {
 
193
          num_mo = s->num_mo;
 
194
        
 
195
          /*
 
196
          if(!m && !t) {
 
197
            fprintf(outfile, "\nuhf_iter Fock matrix irrep %d spin %d iter %d\n", m, t, iter);
 
198
            print_array(sp->scf_spin[m].fock_pac, nn, outfile);
 
199
          }
 
200
          */
 
201
           
 
202
          /* transform fock_pac to mo basis */
 
203
          tri_to_sq(sp->scf_spin[m].fock_pac,fock_ct,nn);
 
204
          /*               mxmb(sp->scf_spin[m].cmat,nn,1
 
205
                           ,fock_ct,1,nn,scr,1,nn,nn,nn,nn);
 
206
                           mxmb(scr,1,nn,sp->scf_spin[m].cmat,1
 
207
                           ,nn,fock_c,1,nn,nn,nn,nn);*/
 
208
 
 
209
          mmult(sp->scf_spin[m].cmat,1,fock_ct,0,scr,0,num_mo,nn,nn,0);
 
210
          mmult(scr,0,sp->scf_spin[m].cmat,0,fock_c,0,num_mo,nn,num_mo,0);
 
211
 
 
212
          /*
 
213
          if(!m && !t) {
 
214
            fprintf(outfile, "\nMO uhf_iter Fock matrix irrep %d spin %d iter %d\n", m, t, iter);
 
215
            print_mat(fock_c, num_mo, num_mo, outfile);
 
216
          }
 
217
          */
 
218
                   
 
219
          /*  diagonalize fock_c to get ctrans */
 
220
          sq_rsp(num_mo,num_mo,fock_c,sp->scf_spin[m].fock_evals
 
221
                 ,1,ctrans,tol);
 
222
                   
 
223
          if(print & 4) {
 
224
            fprintf(outfile,"\n %s eigenvector for irrep %s\n"
 
225
                    ,sp->spinlabel,s->irrep_label);
 
226
            eivout(ctrans,sp->scf_spin[m].fock_evals,num_mo,num_mo,outfile);
 
227
          }
 
228
                   
 
229
          /*               mxmb(sp->scf_spin[m].cmat,1,nn,
 
230
                           ctrans,1,nn,scr,1,nn,nn,nn,nn);*/
 
231
          mmult(sp->scf_spin[m].cmat,0,ctrans,0,scr,0,nn,num_mo,num_mo,0);
 
232
                                   
 
233
          if(print & 4) {
 
234
            fprintf(outfile,"\n %s eigenvector after irrep %s\n",
 
235
                    sp->spinlabel,s->irrep_label);
 
236
            print_mat(scr,nn,num_mo,outfile);
 
237
          }
 
238
                   
 
239
          for (i=0; i < nn; i++)
 
240
            for (j=0; j < num_mo; j++)
 
241
              sp->scf_spin[m].cmat[i][j] = scr[i][j];
 
242
 
 
243
        }
 
244
      }
 
245
           
 
246
      if(converged) {
 
247
        free_block(scr);
 
248
        free_block(fock_c);
 
249
        free_block(fock_ct);
 
250
        free_block(ctrans);
 
251
        exit(PSI_RETURN_FAILURE);
 
252
        cleanup();
 
253
      }
 
254
    }
 
255
    schmit_uhf(1);
 
256
       
 
257
    if(print & 4) {
 
258
      for(j=0; j < 2; j++){
 
259
        for(i=0; i < num_ir ; i++) {
 
260
          s = &scf_info[i];
 
261
          if (nn=s->num_so) {
 
262
            num_mo = s->num_mo;
 
263
            fprintf(outfile,"\northogonalized mos irrep %s\n",
 
264
                    s->irrep_label);
 
265
            print_mat(spin_info[j].scf_spin[i].cmat,nn,num_mo,outfile);
 
266
          }
 
267
        }
 
268
      }
 
269
    }
 
270
       
 
271
    if(mixing && iter ==1)
 
272
      orb_mix();
 
273
 
 
274
    /* form new density matrix */
 
275
    dmatuhf();
 
276
       
 
277
    /* and form new fock matrix */
 
278
    if(iter < itmax) {
 
279
      if (!direct_scf)
 
280
        formg_open();
 
281
      else
 
282
        formg_direct();
 
283
    }
 
284
  }
 
285
}
 
286
 
 
287
}} // namespace psi::cscf