← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/cscf/scf_iter_2.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/cscf/scf_iter_2.cc
 
1
/*! \file
 
2
    \ingroup CSCF
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
/* $Log$
 
6
 * Revision 1.3  2003/08/09 17:39:56  crawdad
 
7
 * I added the ability to determine frozen core orbitals for UHF references to
 
8
 * cleanup.c.  I also commented out ip_cwk_clear and ip_cwk_add calls in
 
9
 * cleanup.c, guess.c, scf_input.c and scf_iter_2.c.  These calls were (1) poor
 
10
 * design and (2) interfering with default ip_tree behavior needed to simplify
 
11
 * the format of input.dat.
 
12
 * -TDC
 
13
 *
 
14
/* Revision 1.2  2002/03/25 02:17:36  janssen
 
15
/* Get rid of tmpl.  Use new naming scheme for libipv1 includes.
 
16
/*
 
17
/* Revision 1.1.1.1  2000/02/04 22:52:32  evaleev
 
18
/* Started PSI 3 repository
 
19
/*
 
20
/* Revision 1.2  1999/08/17 19:04:17  evaleev
 
21
/* Changed the default symmetric orthogonalization to the canonical
 
22
/* orthogonalization. Now, if near-linear dependencies in the basis are found,
 
23
/* eigenvectors of the overlap matrix with eigenvalues less than 1E-6 will be
 
24
/* left out. This will lead to num_mo != num_so, i.e. SCF eigenvector is no
 
25
/* longer a square matrix. Had to rework some routines in libfile30, and add some.
 
26
/* The progrem prints out a warning if near-linear dependencies are found. TRANSQT
 
27
/* and a whole bunch of other codes has to be fixed to work with such basis sets.
 
28
/*
 
29
/* Revision 1.1.1.1  1999/04/12 16:59:28  evaleev
 
30
/* Added a version of CSCF that can work with CINTS.
 
31
/* -Ed
 
32
 * */
 
33
 
 
34
#define EXTERN
 
35
#include "includes.h"
 
36
#include "common.h"
 
37
#include <libipv1/ip_lib.h>
 
38
 
 
39
namespace psi { namespace cscf {
 
40
 
 
41
/* scf procedure for two-configuration scf */
 
42
 
 
43
void scf_iter_2()
 
44
 
 
45
{
 
46
   int i,j,k,m,ij;
 
47
   int joff;
 
48
   int nn,num_mo;
 
49
   int iju=0;
 
50
   int newci=1;
 
51
   int opblk=0;
 
52
   int *optest;
 
53
   int errcod;
 
54
   double ciconv = 10.0e-9;  // FAE it was 10.0e-9
 
55
   double c1i,c1ii;
 
56
   double eci1,eci2,eci3;
 
57
   double ci1,ci2,ci1old,term,hoff,amax,cn,sn;
 
58
   double cimax,cimin,incr;
 
59
   double occi, occj, occ0, den;
 
60
   double **scr;
 
61
   double *c1, *c2;
 
62
   double **fock_c, **fock_o;
 
63
   double **fock_ct,**fock_eff;
 
64
   double **ctrans;
 
65
   double tol = 1.0e-16;
 
66
   struct symm *s;
 
67
 
 
68
   diiser = 0.0;
 
69
   /* optest = (int *) init_array(ioff[nbasis]/2); */
 
70
   optest = (int *) init_int_array(ioff[nbasis]);
 
71
   scr = (double **) init_matrix(nsfmax,nsfmax);
 
72
   fock_c = (double **) init_matrix(nsfmax,nsfmax);
 
73
   fock_ct = (double **) init_matrix(nsfmax,nsfmax);
 
74
   ctrans = (double **) init_matrix(nsfmax,nsfmax);
 
75
   c1 = (double *) init_array(num_ir);
 
76
   c2 = (double *) init_array(num_ir);
 
77
   fock_o = (double **) init_matrix(nsfmax,nsfmax);
 
78
 
 
79
   for (i=0; i < num_ir ; i++) {
 
80
                 s = &scf_info[i];
 
81
                 if (nn=s->num_so) {
 
82
                   num_mo = s->num_mo;
 
83
                   s->ucmat = block_matrix(num_mo, num_mo);
 
84
                   for(j=0; j<num_mo; j++)
 
85
                     s->ucmat[j][j] = 1.0;
 
86
     }
 
87
   }
 
88
 
 
89
         
 
90
   /*
 
91
   ip_cwk_clear();
 
92
   ip_cwk_add(":DEFAULT");
 
93
   ip_cwk_add(":SCF");
 
94
   */
 
95
 
 
96
   incr=0.25;
 
97
   errcod = ip_data("INCR","%lf",&incr,0);
 
98
   fprintf(outfile,"\n  tcscf increment = %f\n",incr);
 
99
 
 
100
 
 
101
/* find open shells */
 
102
 
 
103
   for (i=0; i < num_ir ; i++) {
 
104
      iju += scf_info[i].num_mo;
 
105
      if(scf_info[i].nopen) {
 
106
         iju = ioff[iju];
 
107
         break;
 
108
         }
 
109
      }
 
110
           
 
111
   for (i=0; i < num_ir ; i++) {
 
112
      s= &scf_info[i];
 
113
      if(nn=s->num_so) {
 
114
         num_mo = s->num_mo;
 
115
         for (j=0; j < num_mo ; j++) {
 
116
            if (s->occ_num[j] == 1.0) {
 
117
               if (!opblk) {
 
118
                  opblk++;
 
119
                  opshl1 = j;
 
120
                  opblk1 = i;
 
121
                  }
 
122
               else {
 
123
                  opshl2 = j;
 
124
                  opblk2 = i;
 
125
                  }
 
126
               }
 
127
            }
 
128
         }
 
129
      }
 
130
 
 
131
/* set up array of flags indicating open shells */
 
132
 
 
133
   for (k=joff=0; k < num_ir ; k++) {
 
134
      s = &scf_info[k];
 
135
      if (nn=s->num_so) {
 
136
         for (i=0; i < nn ; i++)
 
137
            for (j=0; j <= i ; j++)
 
138
               if(s->nopen) optest[ioff[i+joff]+j+joff] = 1;
 
139
         joff += nn;
 
140
         }
 
141
      }
 
142
 
 
143
/* set up c1 and c2  */
 
144
   for (i=0; i < num_ir ; i++) {
 
145
      s = &scf_info[i];
 
146
      if (nn=s->num_so) {
 
147
         num_mo = s->num_mo;
 
148
         for (j=0; j < num_mo ; j++) {
 
149
            if(s->occ_num[j]==2.0) c1[i]=2.0;
 
150
            if(s->occ_num[j]==1.0) c2[i]=1.0;
 
151
            }
 
152
         if(i == opblk1) c1i=c1[i];
 
153
         if(i == opblk2) c1ii=c1[i];
 
154
         den = c1[i]-c2[i];
 
155
         if(den != 0.0) {
 
156
            c1[i] /= den;
 
157
            c2[i] /= den;
 
158
            }
 
159
         }
 
160
      }
 
161
   ci1    =  1.0;
 
162
   ci1old =  0.0;  // FAE Give it some arbitrary value
 
163
 
 
164
/* and iterate */
 
165
   
 
166
   if (second_root) 
 
167
      fprintf(outfile, "Warning: Finding second root of TCSCF\n\n");
 
168
 
 
169
   for (iter=0; iter < itmax ; ) {
 
170
      if(iter==stop_lshift){
 
171
        lshift = 1.0;
 
172
        fprintf(outfile,"Setting level shift to %lf\n",lshift);
 
173
      }
 
174
      schmit(1);
 
175
 
 
176
      if(!newci || fabs(ci1old-ci1) < ciconv && iter) goto L2;
 
177
 
 
178
/* do some funky stuff to get coeffs ? */
 
179
 
 
180
      dmat_2(opblk1);
 
181
      dmat_2(opblk2);
 
182
 
 
183
      for (i=0; i < num_ir ; i++) {
 
184
         double d1,d2;
 
185
         s = &scf_info[i];
 
186
         if (nn=s->num_so) {
 
187
            for (j=0; j < ioff[nn] ; j++) {
 
188
               d1 = s->pmat2[j];
 
189
               d2 = s->pmato2[j];
 
190
               s->pmat2[j] = d1+d2;
 
191
               s->pmato2[j] = d1-d2;
 
192
               }
 
193
            }
 
194
         }
 
195
 
 
196
      formg_open();
 
197
 
 
198
      eci1 = eci2 = eci3 = 0.0;
 
199
      for (i=0; i < num_ir ; i++) {
 
200
         s = &scf_info[i];
 
201
         if (nn=s->num_so) {
 
202
            for (j=0; j < ioff[nn] ; j++) {
 
203
               eci1 += s->pmato2[j]*s->hmat[j];
 
204
               eci2 += s->pmato2[j]*s->dpmat[j];
 
205
               eci3 += s->pmato2[j]*s->dpmato[j];
 
206
               }
 
207
            }
 
208
         }
 
209
 
 
210
      term = 0.5*(eci1+0.5*eci2);
 
211
      hoff = -0.25*eci3;
 
212
      amax = term*term + hoff*hoff;
 
213
      amax = sqrt(amax);
 
214
      if(term < 0.0) amax = -(amax);
 
215
      cn = (amax+term)/(2*amax);
 
216
      cn = sqrt(cn);
 
217
      sn = hoff/(cn*(amax+amax));
 
218
      ci1old = ci1;
 
219
 
 
220
      if (!second_root) {
 
221
         if(term >= 0.0) {
 
222
            ci1 = cn;
 
223
            ci2 = -sn;
 
224
            }
 
225
         else {
 
226
            ci1=sn;
 
227
            ci2=cn;
 
228
            }
 
229
         }
 
230
      else {
 
231
         if(term >= 0.0) {
 
232
            ci1 = sn;
 
233
            ci2 = cn;
 
234
            }
 
235
         else {
 
236
            ci1=cn;
 
237
            ci2=-sn;
 
238
            }
 
239
      }
 
240
 
 
241
      save_ci1 = ci1;
 
242
      save_ci2 = ci2;
 
243
 
 
244
      scf_info[opblk1].occ_num[opshl1] = ci1*ci1*2.0;
 
245
      scf_info[opblk2].occ_num[opshl2] = ci2*ci2*2.0;
 
246
      den = c1i-scf_info[opblk1].occ_num[opshl1];
 
247
      c1[opblk1] = c1i/den;
 
248
      c2[opblk1] = scf_info[opblk1].occ_num[opshl1]/den;
 
249
      den = c1ii-scf_info[opblk2].occ_num[opshl2];
 
250
      c1[opblk2] = c1ii/den;
 
251
      c2[opblk2] = scf_info[opblk2].occ_num[opshl2]/den;
 
252
      alpha1 = 1.0 - 2.0/scf_info[opblk1].occ_num[opshl1];
 
253
      alpha2 = -1.0/(ci1*ci2);
 
254
      alpha3 = 1.0 - 2.0/scf_info[opblk2].occ_num[opshl2];
 
255
 
 
256
      cimax = MAX0(ci1,ci2);
 
257
      cimin = MIN0(ci1,ci2);
 
258
      cimax = 1.0/(cimax*cimax+dampsv);
 
259
      if(fabs(cimin) >= 0.1) cimax = 1.0; 
 
260
 
 
261
      fprintf(outfile,"  ci coeffs %14.10f %14.10f\n",ci1,ci2);
 
262
 
 
263
L2:
 
264
 
 
265
/* form density matrix */
 
266
      
 
267
      dmat();
 
268
 
 
269
/* form g matrix */
 
270
 
 
271
      formg_two(iju,optest);
 
272
 
 
273
      for (m=0; m < num_ir ; m++) {
 
274
         s = &scf_info[m];
 
275
         if (nn = s->num_so) {
 
276
 
 
277
          /*  form fock matrix = h+g */
 
278
            add_arr(s->hmat,s->gmat,s->fock_pac,ioff[nn]);
 
279
 
 
280
          /* form fock_open = h+g-q */
 
281
 
 
282
            for (i=0; i < ioff[nn] ; i++)
 
283
               s->fock_open[i]=s->fock_pac[i]-s->gmato[i];
 
284
            }
 
285
         }
 
286
 
 
287
      newci = ecalc(incr);
 
288
 
 
289
   /* create new fock matrix in fock_eff */
 
290
      if(!diisflg) diis(scr,fock_c,fock_ct,c1,c2,cimax,newci);
 
291
 
 
292
      for (m=0; m < num_ir ; m++) {
 
293
         s = &scf_info[m];
 
294
         if (nn=s->num_so) {
 
295
            num_mo = s->num_mo;
 
296
         /* transform fock_pac to mo basis */
 
297
            tri_to_sq(s->fock_pac,fock_ct,nn);
 
298
/*            mxmb(s->cmat,nn,1,fock_ct,1,nn,scr,1,nn,nn,nn,nn);
 
299
            mxmb(scr,1,nn,s->cmat,1,nn,fock_c,1,nn,nn,nn,nn);*/
 
300
            mmult(s->cmat,1,fock_ct,0,scr,0,num_mo,nn,nn,0);
 
301
            mmult(scr,0,s->cmat,0,fock_c,0,num_mo,nn,num_mo,0);
 
302
 
 
303
         /* transform fock_open to mo basis */
 
304
            tri_to_sq(s->fock_open,fock_ct,nn);
 
305
/*            mxmb(s->cmat,nn,1,fock_ct,1,nn,scr,1,nn,nn,nn,nn);
 
306
            mxmb(scr,1,nn,s->cmat,1,nn,fock_o,1,nn,nn,nn,nn);*/
 
307
            mmult(s->cmat,1,fock_ct,0,scr,0,num_mo,nn,nn,0);
 
308
            mmult(scr,0,s->cmat,0,fock_o,0,num_mo,nn,num_mo,0);
 
309
 
 
310
         /* form effective fock matrix in mo basis */
 
311
 
 
312
            occ0 = s->occ_num[0];
 
313
            for (i=ij=0; i < num_mo; i++) {
 
314
               for (j=0; j <= i; j++,ij++) {
 
315
                  occi = s->occ_num[i];
 
316
                  occj = s->occ_num[j];
 
317
 
 
318
              /* default: Guest & Saunders general form */
 
319
                  if(iter < itmax-1 && !converged && !fock_typ) {
 
320
                     if(occi == occj)
 
321
                       s->fock_eff[ij] = fock_c[i][j];
 
322
                     else if(occi)
 
323
                       s->fock_eff[ij] = 
 
324
                         cimax*(c1[m]*fock_c[i][j]-c2[m]*fock_o[i][j]);
 
325
                     else {
 
326
                       if(occj==2.0)
 
327
                         s->fock_eff[ij] = cimax*fock_c[i][j];
 
328
                       else
 
329
                         s->fock_eff[ij] = cimax*fock_o[i][j];
 
330
                       }
 
331
                     }
 
332
 
 
333
               /* Guest & Saunders' form for high spin */
 
334
                  else if(iter < itmax-1 && !converged && fock_typ == 1) {
 
335
                     if (occi == occj || occi)
 
336
                        s->fock_eff[ij] = c1[m]*fock_c[i][j]-c2[m]*fock_o[i][j];
 
337
                     else if (occj == 2.0) s->fock_eff[ij] = fock_c[i][j];
 
338
                     else s->fock_eff[ij] = fock_o[i][j];
 
339
                     }
 
340
 
 
341
               /* test form (fo fo fo) */
 
342
                  else if(iter < itmax-1 && !converged && fock_typ == 2) {
 
343
                     if (occi == occj) s->fock_eff[ij] = fock_o[i][j];
 
344
                     else if(occi)
 
345
                        s->fock_eff[ij] = c1[m]*fock_c[i][j]-c2[m]*fock_o[i][j];
 
346
                     else if(occj == 2.0)
 
347
                        s->fock_eff[ij] = fock_c[i][j];
 
348
                     else s->fock_eff[ij] = fock_o[i][j];
 
349
                     }
 
350
 
 
351
            /* test form a*(fc fc fc) */
 
352
                  else if(iter < itmax-1 && !converged && fock_typ == 3) {
 
353
                     if (occi == occj) s->fock_eff[ij] = dampd*fock_c[i][j];
 
354
                     else if(occi)
 
355
                        s->fock_eff[ij] =
 
356
                          dampo*(c1[m]*fock_c[i][j]-c2[m]*fock_o[i][j]);
 
357
                     else if(occj == 2.0)
 
358
                        s->fock_eff[ij] = dampo*fock_c[i][j];
 
359
                     else s->fock_eff[ij] = dampo*fock_o[i][j];
 
360
                     }
 
361
 
 
362
            /* test form a*(fo fo fo) */
 
363
                  else if(iter < itmax-1 && !converged && fock_typ == 4) {
 
364
                     if (occi == occj) s->fock_eff[ij] = dampd*fock_o[i][j];
 
365
                     else if(occi)
 
366
                        s->fock_eff[ij] = 
 
367
                          dampo*(c1[m]*fock_c[i][j]-c2[m]*fock_o[i][j]);
 
368
                     else if(occj == 2.0)
 
369
                        s->fock_eff[ij] = dampo*fock_c[i][j];
 
370
                     else s->fock_eff[ij] = dampo*fock_o[i][j];
 
371
                     }
 
372
 
 
373
            /* test form a*(2fc-fo 2fc-fo 2fc-fo) */
 
374
                  else if(iter < itmax-1 && !converged && fock_typ == 5) {
 
375
                     if (occi == occj)
 
376
                        s->fock_eff[ij] =
 
377
                          dampd*(c1[m]*fock_c[i][j]-c2[m]*fock_o[i][j]);
 
378
                     else if(occi)
 
379
                        s->fock_eff[ij] = 
 
380
                          dampo*(c1[m]*fock_c[i][j]-c2[m]*fock_o[i][j]);
 
381
                     else if(occj == 2.0)
 
382
                        s->fock_eff[ij] = dampo*fock_c[i][j];
 
383
                     else s->fock_eff[ij] = dampo*fock_o[i][j];
 
384
                     }
 
385
 
 
386
            /* form for converged wavefunction */
 
387
                  else {
 
388
                     if (occi == 2.0) s->fock_eff[ij]=fock_c[i][j];
 
389
                     else if(occj != 2.0)
 
390
                       s->fock_eff[ij]=fock_o[i][j];
 
391
                     else {
 
392
                       if(occi)
 
393
                         s->fock_eff[ij] = 
 
394
                           c1[m]*fock_c[i][j]-c2[m]*fock_o[i][j];
 
395
                       else s->fock_eff[ij]=fock_c[i][j];
 
396
                       }
 
397
                     }
 
398
 
 
399
                  if(j==i) {
 
400
                     if (occi == occ0 && occi) s->fock_eff[ij] -= lshift;
 
401
                     else if (occi) s->fock_eff[ij] -= 0.5*lshift;
 
402
                     }
 
403
                  }
 
404
               }
 
405
            }
 
406
         }
 
407
 
 
408
      for (m=0; m < num_ir ; m++) {
 
409
         s = &scf_info[m];
 
410
         if (nn=s->num_so) {
 
411
            num_mo = s->num_mo;
 
412
            occ0 = s->occ_num[0];
 
413
 
 
414
            rsp(num_mo,num_mo,ioff[num_mo],s->fock_eff,s->fock_evals,1,ctrans,tol);
 
415
 
 
416
/*            mxmb(s->cmat,1,nn,ctrans,1,nn,scr,1,nn,nn,nn,nn);*/
 
417
            mmult(s->cmat,0,ctrans,0,scr,0,nn,num_mo,num_mo,0);
 
418
 
 
419
            for (i=0; i < num_mo; i++) {
 
420
               occi = s->occ_num[i];
 
421
               if (occi == occ0 && occi) s->fock_evals[i] += lshift;
 
422
               else if (occi) s->fock_evals[i] += 0.5*lshift;
 
423
 
 
424
               }
 
425
            for(i=0; i < nn; i++)
 
426
              for (j=0; j < num_mo; j++)
 
427
                s->cmat[i][j] = scr[i][j];
 
428
            }
 
429
         }
 
430
 
 
431
      if(converged) {
 
432
         free_matrix(scr,nsfmax);
 
433
         free_matrix(fock_c,nsfmax);
 
434
         free_matrix(fock_ct,nsfmax);
 
435
         free_matrix(ctrans,nsfmax);
 
436
         free_matrix(fock_o,nsfmax);
 
437
         free(c1);
 
438
         free(c2);
 
439
         free(optest);
 
440
         /* Clean up */
 
441
                                 for(i=0; i < num_ir ; i++) {
 
442
                                   s = &scf_info[i];
 
443
                                   if (nn=s->num_so)
 
444
                                   free_block(s->ucmat);
 
445
                                 }
 
446
         cleanup();
 
447
         }
 
448
      }
 
449
   }
 
450
 
 
451
}} // namespace psi::cscf