← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/cscf/shalf.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/cscf/shalf.c
1
 
/* $Log$
2
 
 * Revision 1.5  2007/04/05 15:45:25  crawdad
3
 
 * Fixed a few memory leaks identified by valgrind. -TDC
4
 
 *
5
 
/* Revision 1.4  2004/05/03 04:32:40  crawdad
6
 
/* Major mods based on merge with stable psi-3-2-1 release.  Note that this
7
 
/* version has not been fully tested and some scf-optn test cases do not run
8
 
/* correctly beccause of changes in mid-March 2004 to optking.
9
 
/* -TDC
10
 
/*
11
 
/* Revision 1.3.6.4  2004/04/21 15:45:07  evaleev
12
 
/* Modified DIIS algorithm for RHF and ROHF to work in OSO basis rather than in
13
 
/* AO basis, to avoid difficulties of transforming between MO and AO bases
14
 
/* when linear dependencies are present.
15
 
/*
16
 
/* Revision 1.3.6.3  2004/04/10 19:41:32  crawdad
17
 
/* Fixed the DIIS code for UHF cases.  The new version uses the Pulay scheme of
18
 
/* building the error vector in the AO basis as FDS-SDF, followed by xformation
19
 
/* into the orthogonal AO basis.   This code converges faster for test cases
20
 
/* like cc8, but fails for linearly dependent basis sets for unknown reasons.
21
 
/* -TDC
22
 
/*
23
 
/* Revision 1.3.6.2  2004/04/09 00:16:29  evaleev
24
 
/* Added messages explaining why DGETRF and DGETRI most likely fail.
25
 
/*
26
 
/* Revision 1.3.6.1  2004/04/06 21:29:05  crawdad
27
 
/* Corrections to the RHF/ROHF DIIS algorithm, which was simply incorrect.
28
 
/* The backtransformation of the DIIS error vectors to the AO basis was not
29
 
/* mathematically right.
30
 
/* -TDC and EFV
31
 
/*
32
 
/* Revision 1.3  2002/11/24 22:52:17  crawdad
33
 
/* Merging the gbye-file30 branch into the main trunk.
34
 
/* -TDC
35
 
/*
36
 
/* Revision 1.2.6.1  2002/11/23 21:54:45  crawdad
37
 
/* Removal of mxcoef stuff for chkpt runs.
38
 
/* -TDC
39
 
/*
40
 
/* Revision 1.2  2000/10/13 19:51:22  evaleev
41
 
/* Cleaned up a lot of stuff in order to get CSCF working with the new "Mo-projection-capable" INPUT.
42
 
/*
43
 
/* Revision 1.1.1.1  2000/02/04 22:52:33  evaleev
44
 
/* Started PSI 3 repository
45
 
/*
46
 
/* Revision 1.3  1999/08/17 19:04:18  evaleev
47
 
/* Changed the default symmetric orthogonalization to the canonical
48
 
/* orthogonalization. Now, if near-linear dependencies in the basis are found,
49
 
/* eigenvectors of the overlap matrix with eigenvalues less than 1E-6 will be
50
 
/* left out. This will lead to num_mo != num_so, i.e. SCF eigenvector is no
51
 
/* longer a square matrix. Had to rework some routines in libfile30, and add some.
52
 
/* The progrem prints out a warning if near-linear dependencies are found. TRANSQT
53
 
/* and a whole bunch of other codes has to be fixed to work with such basis sets.
54
 
/*
55
 
/* Revision 1.2  1999/08/11 18:39:03  evaleev
56
 
/* Added some checks on the lowest eigenvalue of the overlap matrix.
57
 
/*
58
 
/* Revision 1.1.1.1  1999/04/12 16:59:28  evaleev
59
 
/* Added a version of CSCF that can work with CINTS.
60
 
/* -Ed
61
 
 * */
62
 
 
63
 
static char *rcsid = "$Id: shalf.c 3324 2007-04-05 15:45:25Z crawdad $";
64
 
 
65
 
/* construct S-1/2 matrix 'sahalf' using CANONICAL orthogonalization  */
66
 
 
67
 
#define EXTERN
68
 
#include "includes.h"
69
 
#include "common.h"
70
 
 
71
 
void shalf(void)
72
 
{
73
 
  int i,nn,num_mo;
74
 
  int ii,jj,kk;
75
 
  double *eig_vals, **eig_vecs;
76
 
  double **shalf;
77
 
  double tol = 1.0e-20;
78
 
  double min_eval = 100000.0;
79
 
  struct symm *s;
80
 
  int info, lwork, *ipiv;
81
 
  double *work, **P, **T, **X;
82
 
     
83
 
  eig_vals = (double *) init_array(nsfmax);
84
 
  eig_vecs = (double **) init_matrix(nsfmax,nsfmax);
85
 
  shalf = block_matrix(nsfmax,nsfmax);
86
 
 
87
 
#if !USE_LIBCHKPT
88
 
  mxcoef = 0;
89
 
#endif
90
 
  mxcoef2 = 0;
91
 
  nmo = 0;
92
 
     
93
 
  /*  diagonalize smat to get eigenvalues and eigenvectors  */
94
 
 
95
 
  for (i=0; i < num_ir ; i++) {
96
 
    s = &scf_info[i];
97
 
    if (nn=s->num_so) {
98
 
 
99
 
      rsp(nn,nn,ioff[nn],s->smat,eig_vals,1,eig_vecs,tol);
100
 
 
101
 
      /*--- Find the lowest eigenvalue ---*/
102
 
      for(ii=0; ii < nn ; ii++)
103
 
        if (eig_vals[ii] < min_eval)
104
 
          min_eval = eig_vals[ii];
105
 
 
106
 
      if(print & 64) {
107
 
        fprintf(outfile,"\noverlap eigenstuff\n");
108
 
        eivout(eig_vecs,eig_vals,nn,nn,outfile);
109
 
      }
110
 
 
111
 
      /*--- Go through the eigenvalues and "throw away"
112
 
        the dangerously small ones ---*/
113
 
      num_mo = 0;
114
 
      for(ii=0; ii < nn ; ii++)
115
 
        if (eig_vals[ii] >= LINDEP_CUTOFF) {
116
 
          eig_vals[ii] = 1.0/sqrt(eig_vals[ii]);
117
 
          num_mo++;
118
 
        }
119
 
      s->num_mo = num_mo;
120
 
#if !USE_LIBCHKPT
121
 
      mxcoef += num_mo * nn;
122
 
#endif
123
 
      mxcoef2 += ioff[nn];
124
 
      nmo += num_mo;
125
 
      if (num_mo < nn)
126
 
        fprintf(outfile,"\n  In symblk %d %d eigenvectors of S with eigenvalues < %lf are thrown away",
127
 
                i,nn-num_mo,LINDEP_CUTOFF);
128
 
 
129
 
      /* form 'sahalf' matrix sahalf = U*(s-1/2)  */
130
 
 
131
 
      for (ii=0; ii < nn ; ii++)
132
 
        for (jj=0; jj < num_mo ; jj++)
133
 
          s->sahalf[ii][jj] += eig_vecs[ii][jj+nn-num_mo]*eig_vals[jj+nn-num_mo];
134
 
 
135
 
      /* form 'shalf' matrix shalf = U*(s^1/2)  */
136
 
 
137
 
      for (ii=0; ii < nn ; ii++)
138
 
        for (jj=0; jj < num_mo ; jj++)
139
 
          shalf[ii][jj] = eig_vecs[ii][jj+nn-num_mo]*(1.0/eig_vals[jj+nn-num_mo]);
140
 
 
141
 
      /*
142
 
      fprintf(outfile, "S^-1/2 Original:\n");
143
 
      print_mat(s->sahalf, nn, num_mo, outfile);
144
 
      */
145
 
 
146
 
      /* new transformation matrix for the DIIS error vectors: P^-1 = S^{1/2} * (S^{1/2})^+ */
147
 
      C_DGEMM('n', 't', nn, nn, num_mo, 1.0, shalf[0], nsfmax, shalf[0], nsfmax,
148
 
              0.0, s->pinv[0], nn);
149
 
 
150
 
    }
151
 
  }
152
 
 
153
 
  free(eig_vals);
154
 
  free_matrix(eig_vecs,nsfmax);
155
 
 
156
 
  fprintf(outfile,"\n  The lowest eigenvalue of the overlap matrix was %e\n\n",
157
 
          min_eval);
158
 
   
159
 
}