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

« back to all changes in this revision

Viewing changes to src/bin/cscf/ecalc.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:
1
 
/* $Log$
2
 
 * Revision 1.10  2004/05/03 04:32:40  crawdad
3
 
 * Major mods based on merge with stable psi-3-2-1 release.  Note that this
4
 
 * version has not been fully tested and some scf-optn test cases do not run
5
 
 * correctly beccause of changes in mid-March 2004 to optking.
6
 
 * -TDC
7
 
 *
8
 
/* Revision 1.9.4.1  2004/04/10 19:41:32  crawdad
9
 
/* Fixed the DIIS code for UHF cases.  The new version uses the Pulay scheme of
10
 
/* building the error vector in the AO basis as FDS-SDF, followed by xformation
11
 
/* into the orthogonal AO basis.   This code converges faster for test cases
12
 
/* like cc8, but fails for linearly dependent basis sets for unknown reasons.
13
 
/* -TDC
14
 
/*
15
 
/* Revision 1.9  2003/04/14 17:25:47  sherrill
16
 
/* Change "total energy" to "SCF total energy" to make more explicit for
17
 
/* new users.  Yeah, this will probably break some test case perl scripts
18
 
/* temporarily :)
19
 
/*
20
 
/* Revision 1.8  2001/01/04 14:13:35  sbrown
21
 
/* Fixed the problem with iconv:  The new versions of linux had iconv already
22
 
/* assigned to something else so I changed all references of it to scf_conv.
23
 
/*
24
 
/* Revision 1.7  2000/12/05 19:40:03  sbrown
25
 
/* Added Unrestricted Kohn-Sham DFT.
26
 
/*
27
 
/* Revision 1.6  2000/09/02 20:48:51  evaleev
28
 
/* Print out one- and two-electron energies every iteration if iprint&2 .
29
 
/*
30
 
/* Revision 1.5  2000/08/23 17:15:16  sbrown
31
 
/* Added portions to separate out the correlation and exchange energy at the
32
 
/* end the calculation as well as do the consistency check on the integrated
33
 
/* density.
34
 
/*
35
 
/* Revision 1.4  2000/06/26 19:04:09  sbrown
36
 
/* Added DFT capapbilities to interface with cints using direct scf
37
 
/*
38
 
/* Revision 1.3  2000/06/22 22:15:00  evaleev
39
 
/* 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).
40
 
/*
41
 
/* Revision 1.2  2000/06/02 13:32:16  kenny
42
 
/*
43
 
/*
44
 
/* Added dynamic integral accuracy cutoffs for direct scf.  Added a few global
45
 
/* variables.  Added keyword 'dyn_acc'; true--use dynamic cutoffs.  Use of
46
 
/* 'dconv' and 'delta' to keep track of density convergence somewhat awkward,
47
 
/* but avoids problems when accuracy is switched and we have to wipe out density
48
 
/* matrices.  Also added error message and exit if direct rohf singlet is
49
 
/* attempted since it doesn't work.
50
 
/* --Joe Kenny
51
 
/*
52
 
/* Revision 1.1.1.1  2000/02/04 22:52:30  evaleev
53
 
/* Started PSI 3 repository
54
 
/*
55
 
/* Revision 1.3  1999/11/02 23:55:56  localpsi
56
 
/* Shawn Brown - (11/2/99) Modified to the code in a few major ways.
57
 
/*
58
 
/* 1.  Added the capability to do UHF.  All of the features available with the
59
 
/* other refrences have been added for UHF.
60
 
/*
61
 
/* 2.  For UHF, I had to alter the structure of file30. (See cleanup.c for a
62
 
/* map)  This entailed adding a pointer array right after the header in the SCF
63
 
/* section of file30 that pointed to all of the data for the SCF caclulation.
64
 
/* Functions were added to libfile30 to account for this and they are
65
 
/* incorporated in this code.
66
 
/*
67
 
/* 3.  Updated and fixed all of the problems associated with my previous
68
 
/* guessing code.  The code no longer uses OPENTYPE to specify the type of
69
 
/* occupation.  The keword REFERENCE and MULTP can now be used to indicate any
70
 
/* type of calculation.  (e.g. ROHF with MULTP of 1 is an open shell singlet
71
 
/* ROHF calculation)  This code was moved to occ_fun.c.  The code can also
72
 
/* guess at any multplicity in a highspin case, provided enough electrons.
73
 
/*
74
 
/* Revision 1.2  1999/08/17 19:04:14  evaleev
75
 
/* Changed the default symmetric orthogonalization to the canonical
76
 
/* orthogonalization. Now, if near-linear dependencies in the basis are found,
77
 
/* eigenvectors of the overlap matrix with eigenvalues less than 1E-6 will be
78
 
/* left out. This will lead to num_mo != num_so, i.e. SCF eigenvector is no
79
 
/* longer a square matrix. Had to rework some routines in libfile30, and add some.
80
 
/* The progrem prints out a warning if near-linear dependencies are found. TRANSQT
81
 
/* and a whole bunch of other codes has to be fixed to work with such basis sets.
82
 
/*
83
 
/* Revision 1.1.1.1  1999/04/12 16:59:26  evaleev
84
 
/* Added a version of CSCF that can work with CINTS.
85
 
/* -Ed
86
 
 * */
87
 
 
88
 
static char *rcsid = "$Id: ecalc.c 2455 2004-05-03 04:32:41Z crawdad $";
89
 
 
90
 
#define EXTERN
91
 
#include "includes.h"
92
 
#include "common.h"
93
 
 
94
 
static double twocut=1.0;
95
 
static double eelec;       /*--- elec. energy from the previous iteration ---*/
96
 
double dconv;
97
 
 
98
 
int ecalc(incr)
99
 
   double incr;
100
 
 
101
 
{
102
 
  int i,j,k,ij,nn;
103
 
  double edif;
104
 
  double plimit = pow(10.0,(double) -scf_conv);
105
 
  double neelec = 0.0;
106
 
  double oe_energy, te_energy, dtmp, dtmp1;
107
 
  double cinext;
108
 
  struct symm *s;
109
 
 
110
 
  delta=0.0;
111
 
  oe_energy = te_energy = 0.0;
112
 
  for (k=0; k < num_ir ; k++) {
113
 
    s = &scf_info[k];
114
 
    if (nn=s->num_so) {
115
 
 
116
 
      for (i=ij=0; i < nn ; i++) {
117
 
        for (j = 0 ; j <= i ; j++,ij++) {
118
 
          oe_energy += 0.5*s->pmat[ij]*s->hmat[ij];
119
 
          if(uhf) {
120
 
            te_energy += 0.5*((spin_info[0].scf_spin[k].pmat[ij]
121
 
                               *spin_info[0].scf_spin[k].fock_pac[ij])
122
 
                              +(spin_info[1].scf_spin[k].pmat[ij]
123
 
                                *spin_info[1].scf_spin[k].fock_pac[ij]));
124
 
          }
125
 
          else if(!iopen) {
126
 
            te_energy += 0.5*s->pmat[ij]*s->fock_pac[ij];
127
 
          }
128
 
          else {
129
 
            te_energy += 0.5*s->pmat[ij]*s->fock_pac[ij]
130
 
              - 0.5*s->pmato[ij]*s->gmato[ij];
131
 
          }
132
 
        }
133
 
      }
134
 
      if (iter) {
135
 
        if(uhf){
136
 
          for (i = 0; i < ioff[nn] ; i++) {
137
 
            dtmp = spin_info[0].scf_spin[k].dpmat[i];
138
 
            dtmp1 = spin_info[1].scf_spin[k].dpmat[i];
139
 
            delta += dtmp*dtmp;
140
 
            delta += dtmp1*dtmp1;
141
 
          }
142
 
        }
143
 
        else {
144
 
          for (i = 0; i < ioff[nn] ; i++) {
145
 
            dtmp = s->dpmat[i];
146
 
            delta += dtmp*dtmp;
147
 
          }
148
 
        }
149
 
      }
150
 
    }
151
 
  }
152
 
  neelec = oe_energy + te_energy;
153
 
 
154
 
  /*JPK(6/1/00) dynamic integral accuracy modifications*/
155
 
  dconv = sqrt(delta)/mxcoef2;
156
 
  delta = dconv;
157
 
  if(acc_switch==1 || iter==0) {
158
 
    delta=1.0;
159
 
    acc_switch=0;
160
 
  }                            
161
 
  coulomb_energy = neelec;
162
 
  if (ksdft){      
163
 
    neelec += exc;
164
 
    /*printf("XC_energy = %10.10lf",exc);*/
165
 
  }
166
 
   
167
 
  etot = repnuc + neelec;
168
 
  edif =  eelec - neelec;
169
 
  ediff = edif;
170
 
 
171
 
  if (!iter) fprintf(outfile,"\n  iter       total energy        delta E         delta P          diiser\n");
172
 
  fprintf(outfile, "%5d %20.10f %15.6e %15.6e %15.6e\n", 
173
 
          ++iter, etot, edif, dconv, diiser);
174
 
  if (print & 2) {
175
 
    fprintf(outfile, "one-electron energy = %25.15f\n", oe_energy);
176
 
    fprintf(outfile, "two-electron energy = %25.15f\n", te_energy);
177
 
    fprintf(outfile, "coulomb energy      = %25.15f\n",coulomb_energy);
178
 
    fprintf(outfile, "SCF total energy    = %25.15f\n", etot);
179
 
  }
180
 
  fflush(outfile);
181
 
  diiser=0.0;
182
 
 
183
 
  if ( delta < plimit && iter > 1) {
184
 
    converged=1;
185
 
    if(!iopen || iopen && fock_typ >= 2) cleanup();
186
 
  }
187
 
 
188
 
  eelec = neelec;
189
 
 
190
 
  cinext = pow(10.0,-twocut);
191
 
  if (delta < cinext && delta && !converged) {
192
 
    twocut += incr;
193
 
    return(1);
194
 
  }
195
 
  else return(0);
196
 
}