← Back to branch summary

~ubuntu-branches/ubuntu/utopic/nwchem/utopic

« back to all changes in this revision

Viewing changes to src/tools/ga-5-2/tcgmsg/examples/scf.f

  • Committer: Package Import Robot
  • Author(s): Michael Banck, Daniel Leidert, Andreas Tille, Michael Banck
  • Date: 2013-07-04 12:14:55 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20130704121455-5tvsx2qabor3nrui
Tags: 6.3-1
http://bugs.debian.org/696361
* New upstream release.
* Fixes anisotropic properties (Closes: #696361).
* New features include:
  + Multi-reference coupled cluster (MRCC) approaches
  + Hybrid DFT calculations with short-range HF 
  + New density-functionals including Minnesota (M08, M11) and HSE hybrid
    functionals
  + X-ray absorption spectroscopy (XAS) with TDDFT
  + Analytical gradients for the COSMO solvation model
  + Transition densities from TDDFT 
  + DFT+U and Electron-Transfer (ET) methods for plane wave calculations
  + Exploitation of space group symmetry in plane wave geometry optimizations
  + Local density of states (LDOS) collective variable added to Metadynamics
  + Various new XC functionals added for plane wave calculations, including
    hybrid and range-corrected ones
  + Electric field gradients with relativistic corrections 
  + Nudged Elastic Band optimization method
  + Updated basis sets and ECPs 

[ Daniel Leidert ]
* debian/watch: Fixed.

[ Andreas Tille ]
* debian/upstream: References

[ Michael Banck ]
* debian/upstream (Name): New field.
* debian/patches/02_makefile_flags.patch: Refreshed.
* debian/patches/06_statfs_kfreebsd.patch: Likewise.
* debian/patches/07_ga_target_force_linux.patch: Likewise.
* debian/patches/05_avoid_inline_assembler.patch: Removed, no longer needed.
* debian/patches/09_backported_6.1.1_fixes.patch: Likewise.
* debian/control (Build-Depends): Added gfortran-4.7 and gcc-4.7.
* debian/patches/10_force_gcc-4.7.patch: New patch, explicitly sets
  gfortran-4.7 and gcc-4.7, fixes test suite hang with gcc-4.8 (Closes:
  #701328, #713262).
* debian/testsuite: Added tests for COSMO analytical gradients and MRCC.
* debian/rules (MRCC_METHODS): New variable, required to enable MRCC methods.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/tools/ga-5-2/tcgmsg/examples/scf.f
 
1
      program scf
 
2
C$Id: scf.f,v 1.3 1997-03-04 06:17:21 d3e129 Exp $
 
3
      implicit double precision (a-h, o-z)
 
4
      include 'cscf.h'
 
5
      include 'msgtypesf.h'
 
6
      dimension orbs(nbfn*nbfn), dens(nnbfn), fock(nnbfn),
 
7
     $     work(nbfn*nbfn), evals(nbfn)
 
8
      data tinit, tonel, ttwoel, tdiag, tdens, tprint /6*0.0d0/
 
9
      data eone, etwo, energy, deltad /4*0.0d0/
 
10
c     
 
11
c     initalize the parallel message passing environment
 
12
c     
 
13
      call pbeginf
 
14
*      call ieeetrap
 
15
      me = nodeid()
 
16
      nproc = nnodes()
 
17
c     
 
18
c     initialize a bunch of stuff and initial density matrix
 
19
c     
 
20
      rjunk = timer()
 
21
      call ininrm
 
22
      call denges(dens, work)
 
23
      tinit = timer()
 
24
c
 
25
c     make initial orthogonal orbital set for jacobi diagonalizer
 
26
c
 
27
      call makeob(orbs, work)
 
28
c     
 
29
c     iterate
 
30
c     
 
31
      do 10 iter = 1, mxiter
 
32
c     
 
33
c     make info for sparsity test ... redone every iter to save space
 
34
c     
 
35
         call makesz(work, schwmax)
 
36
c     
 
37
c     make the one particle contribution to the fock matrix (in fock)
 
38
c     and the partial contribution to the energy
 
39
c     
 
40
         call oneel(dens, work, schwmax, fock, eone)
 
41
         call dgop(1, eone, 1, '+')
 
42
         tonel = tonel + timer()
 
43
c     
 
44
c     compute the two particle contribution and then add up the
 
45
c     contributions from each process with dgop
 
46
c     
 
47
         call twoel(dens, work, schwmax, fock, etwo)
 
48
         call dgop(2, etwo, 1, '+')
 
49
         call dgop(3, fock, nnbfn, '+')
 
50
         ttwoel = ttwoel + timer()
 
51
c     
 
52
c     only process 0 diagonalizes and updates the density
 
53
c     
 
54
         if (me.eq.0) then
 
55
c     
 
56
c     diagonalize the fock matrix ...
 
57
c     
 
58
            call diagon(fock, orbs, evals, work, tester, iter)
 
59
      call flush(6)
 
60
            tdiag = tdiag + timer()
 
61
c     
 
62
c     make the new density matrix in work from orbitals in orbs,
 
63
c     compute the norm of the change in the density matrix and
 
64
c     then update the density matrix in dens with damping.
 
65
c     
 
66
            call makden(orbs, work)
 
67
            deltad = dendif(dens, work)
 
68
            if (iter.eq.1) then
 
69
               scale = 0.0d0
 
70
            else if (iter .le. 5) then
 
71
               if (nbfn .gt. 60) then
 
72
                  scale = 0.5d0
 
73
               else
 
74
                  scale = 0.0d0
 
75
               endif
 
76
            else
 
77
               scale = 0.0d0
 
78
            endif
 
79
            call damp(scale, dens, work)
 
80
            tdens = tdens + timer()
 
81
c     
 
82
c     add up energy and print out convergence information
 
83
c     
 
84
            energy = enrep + eone + etwo
 
85
            call prnout(iter, energy, deltad, tester)
 
86
            tprint = tprint + timer()
 
87
         endif
 
88
c     
 
89
c     brodcast new density matrix and deltad to everyone
 
90
c     
 
91
         call brdcst(4+MSGDBL, dens, mdtob(nnbfn), 0)
 
92
         call brdcst(5+MSGDBL, deltad, mdtob(1), 0)
 
93
c     
 
94
c     if converged then exit iteration loop
 
95
c     
 
96
         if (deltad .lt. tol) goto 20
 
97
 10   continue
 
98
      if(me.eq.0)
 
99
     $     write(6,*) ' SCF failed to converge in ', mxiter, ' iters'
 
100
c     
 
101
c     finished ... print out eigenvalues and occupied orbitals
 
102
c     
 
103
 20   continue
 
104
      call igop(6, icut1, 1, '+')
 
105
      call igop(7, icut2, 1, '+')
 
106
      call igop(8, icut3, 1, '+')
 
107
      if (me.eq.0) then
 
108
c     
 
109
c     print out timing information
 
110
c     
 
111
         call prnfin(energy, evals, orbs)
 
112
         write(6,1) tinit, tonel, ttwoel, tdiag, tdens, tprint,
 
113
     $        nproc
 
114
 1       format(/'   init   onel  twoel   diag   dens   print  ncpu'/
 
115
     $        '  ------ ------ ------ ------ ------ ------ ------'/
 
116
     $        1x, 6f7.2, i7/)
 
117
c     
 
118
c     print out information on # integrals evaulated each iteration
 
119
c     
 
120
         nints = nnbfn*(nnbfn+1)/2
 
121
         frac  = dble(icut3)/dble(nints)
 
122
         write(6,2) icut1, icut2, icut3, nints, frac
 
123
 2       format(/'       No. of integrals screened or computed '
 
124
     $        /'       -------------------------------------'/
 
125
     $        /1x,'   #ij test   #kl test   #compute     #total',
 
126
     $        '  fraction',
 
127
     $        /1x,'  ---------  ---------  ---------  ---------',
 
128
     $        '  --------',
 
129
     $        /1x,4(2x,i9),f9.3)
 
130
         call stats
 
131
      endif
 
132
c     
 
133
      call pend
 
134
      call fexit
 
135
c     
 
136
      end
 
137
      subroutine makesz(schwarz, schwmax)
 
138
      implicit double precision (a-h, o-z)
 
139
      include 'cscf.h'
 
140
      include 'msgtypesf.h'
 
141
      dimension schwarz(nnbfn)
 
142
c
 
143
c     schwarz(ij) = (ij|ij) for sparsity test
 
144
c
 
145
      icut1 = 0
 
146
      icut2 = 0
 
147
      icut3 = 0
 
148
      ij = 0
 
149
      schwmax = 0.0d0
 
150
      call dfill(nnbfn, 0.0d0, schwarz, 1)
 
151
c
 
152
      me = nodeid()
 
153
      nproc = nnodes()
 
154
      do 10 i = 1, nbfn
 
155
         do 20 j = 1, i
 
156
            ij = ij + 1
 
157
            if (mod(ij,nproc).eq.me) then
 
158
               call g(gg, i, j, i, j)
 
159
               schwarz(ij) = sqrt(gg)
 
160
               schwmax = max(schwmax, schwarz(ij))
 
161
            endif
 
162
 20      continue
 
163
 10   continue
 
164
c
 
165
      call dgop(101+MSGDBL, schwarz, nnbfn, '+')
 
166
      call dgop(102+MSGDBL, schwmax, 1, 'max')
 
167
c
 
168
      end
 
169
      subroutine ininrm
 
170
      implicit double precision (a-h, o-z)
 
171
      include 'cscf.h'
 
172
c
 
173
c     write a little welcome message
 
174
c
 
175
      if (nodeid().eq.0) write(6,1) natom, nocc, nbfn, tol
 
176
1     format(/' Example Direct Self Consistent Field Program '/
 
177
     $        ' -------------------------------------------- '//
 
178
     $        ' no. of atoms ............... ',i3/
 
179
     $        ' no. of occupied orbitals ... ',i3/
 
180
     $        ' no. of basis functions ..... ',i3/
 
181
     $        ' convergence threshold ...... ',d9.2//)
 
182
c
 
183
c     generate normalisation coefficients for the basis functions
 
184
c     and the index array iky
 
185
c
 
186
      do 10 i = 1, nbfn
 
187
         iky(i) = i*(i-1)/2
 
188
 10   continue
 
189
c
 
190
      do 20 i = 1, nbfn
 
191
         rnorm(i) = (expnt(i)*2.0d0/pi)**0.75d0
 
192
 20   continue
 
193
c
 
194
c     initialize common for computing f0
 
195
c
 
196
      call setfm
 
197
c
 
198
      end
 
199
      double precision function h(i,j)
 
200
      implicit double precision (a-h, o-z)
 
201
      include 'cscf.h'
 
202
cvd$r novector
 
203
cvd$r noconcur
 
204
c
 
205
c     generate the one particle hamiltonian matrix element
 
206
c     over the normalized primitive 1s functions i and j
 
207
c
 
208
      f0val = 0.0d0
 
209
      sum = 0.0d0
 
210
      rab2 = (x(i)-x(j))**2 + (y(i)-y(j))**2 + (z(i)-z(j))**2
 
211
      facij = expnt(i)*expnt(j)/(expnt(i)+expnt(j))
 
212
      expij = exprjh(-facij*rab2)
 
213
      repij = (2.0d0*pi/(expnt(i)+expnt(j))) * expij
 
214
c
 
215
c     first do the nuclear attraction integrals
 
216
c
 
217
      do 10 iat = 1, natom
 
218
         xp = (x(i)*expnt(i) + x(j)*expnt(j))/(expnt(i)+expnt(j))
 
219
         yp = (y(i)*expnt(i) + y(j)*expnt(j))/(expnt(i)+expnt(j))
 
220
         zp = (z(i)*expnt(i) + z(j)*expnt(j))/(expnt(i)+expnt(j))
 
221
         rpc2 = (xp-ax(iat))**2 + (yp-ay(iat))**2 + (zp-az(iat))**2
 
222
c
 
223
         call f0(f0val, (expnt(i)+expnt(j))*rpc2)
 
224
         sum = sum - repij * q(iat) * f0val
 
225
 10   continue
 
226
c
 
227
c     add on the kinetic energy term
 
228
c
 
229
      sum = sum + facij*(3.0d0-2.0d0*facij*rab2) *
 
230
     $     (pi/(expnt(i)+expnt(j)))**1.5d0 * expij
 
231
c
 
232
c     finally multiply by the normalization constants
 
233
c
 
234
      h = sum * rnorm(i) * rnorm(j)
 
235
c
 
236
      end
 
237
      double precision function s(i,j)
 
238
      implicit double precision (a-h, o-z)
 
239
      include 'cscf.h'
 
240
c
 
241
c     generate the overlap matrix element between the normalized
 
242
c     primitve gaussian 1s functions i and j
 
243
c
 
244
      rab2 = (x(i)-x(j))**2 + (y(i)-y(j))**2 + (z(i)-z(j))**2
 
245
      facij = expnt(i)*expnt(j)/(expnt(i)+expnt(j))
 
246
      s = (pi/(expnt(i)+expnt(j)))**1.5d0 * exprjh(-facij*rab2) *
 
247
     $     rnorm(i)*rnorm(j)
 
248
c
 
249
      end
 
250
      subroutine makden(orbs, dens)
 
251
      implicit double precision (a-h, o-z)
 
252
      include 'cscf.h'
 
253
      dimension orbs(nbfn, nbfn), dens(nnbfn)
 
254
c
 
255
c     generate density matrix from orbitals in orbs. the first
 
256
c     nocc orbitals are doubly occupied. Note that the diagonal
 
257
c     elements are scaled by 0.5
 
258
c
 
259
      ij = 0
 
260
      do 10 i = 1, nbfn
 
261
         do 20 j = 1, i
 
262
            p = 0.0d0
 
263
            do 30 k = 1, nocc
 
264
               p = p + orbs(i,k)*orbs(j,k)
 
265
 30         continue
 
266
            ij = ij + 1
 
267
            dens(ij) = 2.0d0 * p
 
268
 20      continue
 
269
         dens(ij) = dens(ij)*0.5d0
 
270
 10   continue
 
271
c
 
272
      end
 
273
      subroutine oneel(dens, schwarz, schwmax, fock, eone)
 
274
      implicit double precision (a-h, o-z)
 
275
      include 'cscf.h'
 
276
      dimension dens(nnbfn), fock(nnbfn), schwarz(nnbfn)
 
277
c
 
278
c     fill in the one-electron part of the fock matrix and
 
279
c     compute the one-electron energy contribution
 
280
c
 
281
c     simple structure to share out the work between processes
 
282
c
 
283
      me = nodeid()
 
284
      nproc = nnodes()
 
285
c
 
286
      call dfill(nnbfn, 0.0d0, fock, 1)
 
287
      do 10 i = me+1, nbfn, nproc
 
288
         do 20 j = 1,i
 
289
            ij = iky(i) + j
 
290
            if (schwarz(ij)*schwmax.gt.tol2e) fock(ij) = h(i,j)
 
291
 20      continue
 
292
 10   continue
 
293
      eone = ddot(nnbfn, fock, 1, dens, 1)
 
294
c
 
295
      end
 
296
      integer function nxtask(nproc)
 
297
      parameter (ichunk = 10)
 
298
      save icount, nleft
 
299
      data nleft, icount /0, 0/
 
300
c
 
301
c     wrapper round nxtval() to increase granularity
 
302
c     and thus reduce no. of requests to shared counter
 
303
c
 
304
      if(nproc.gt.0) then
 
305
         if(nleft.eq.0) then
 
306
            icount = nxtval(nproc) * ichunk
 
307
            nleft = ichunk
 
308
         endif
 
309
         nxtask = icount
 
310
         icount = icount + 1
 
311
         nleft = nleft -1
 
312
      else
 
313
          nleft = 0
 
314
          nxtask = 0
 
315
          junk = nxtval(nproc)
 
316
      endif
 
317
c
 
318
c     following does dumb static load balancing
 
319
c
 
320
c$$$      if(nproc.gt.0) then
 
321
c$$$         if (nleft .eq. 0) then
 
322
c$$$            icount = nodeid()
 
323
c$$$            nleft = 1
 
324
c$$$         endif
 
325
c$$$         nxtask = icount
 
326
c$$$         icount = icount + nnodes()
 
327
c$$$      else
 
328
c$$$          nleft = 0
 
329
c$$$          nxtask = 0
 
330
c$$$      endif
 
331
      end
 
332
      subroutine twoel(dens, schwarz, schwmax, fock, etwo)
 
333
      implicit double precision (a-h, o-z)
 
334
      include 'cscf.h'
 
335
      dimension dens(nnbfn), fock(nnbfn), schwarz(nnbfn)
 
336
c     
 
337
c     add in the two-electron contribution to the fock matrix
 
338
c     
 
339
      nproc = nnodes()
 
340
c     
 
341
      gg = 0.0d0
 
342
c     
 
343
      next = nnbfn - nxtask(nproc)
 
344
      do 10 i = nbfn, 1, -1
 
345
         do 20 j = i, 1, -1
 
346
            ij = iky(i) + j
 
347
            if (ij .eq. next) then
 
348
               if (schwarz(ij)*schwmax .lt. tol2e) then
 
349
                  icut1 = icut1 + ij
 
350
               else
 
351
                  do 30 k = 1, i
 
352
                     lhi = k
 
353
                     if (k.eq.i) lhi = j
 
354
                     do 40 l = 1, lhi
 
355
                        kl = iky(k) + l
 
356
                        if (schwarz(ij)*schwarz(kl).lt.tol2e) then
 
357
                           icut2 = icut2 + 1
 
358
                        else
 
359
                           icut3 = icut3 + 1
 
360
c     
 
361
c     compute value of integral (ij|kl) and add into fock matrix
 
362
c     
 
363
                           call g(gg, i, j, k, l)
 
364
                           call addin(gg*0.5d0, i, j, k, l, fock,
 
365
     $                          dens, iky)
 
366
                        endif
 
367
 40                  continue
 
368
 30               continue
 
369
               endif
 
370
               next = nnbfn - nxtask(nproc)
 
371
            endif
 
372
 20      continue
 
373
 10   continue
 
374
c     
 
375
      etwo = ddot(nnbfn, fock, 1, dens, 1)
 
376
c     
 
377
c     wait for all processes to finish work
 
378
c     
 
379
      junk = nxtask(-nproc)
 
380
c     
 
381
      end
 
382
      subroutine damp(fac, dens, work)
 
383
      implicit double precision (a-h, o-z)
 
384
      include 'cscf.h'
 
385
      dimension dens(nnbfn), work(nnbfn)
 
386
c
 
387
      ofac = 1.0d0 - fac
 
388
      do 10 i = 1, nnbfn
 
389
         dens(i) = fac*dens(i) + ofac*work(i)
 
390
10    continue
 
391
c
 
392
      end
 
393
      subroutine prnout(iter, energy, deltad, tester)
 
394
      implicit double precision (a-h, o-z)
 
395
      include 'cscf.h'
 
396
c
 
397
c     printout results of each iteration
 
398
c
 
399
      write(6,1) iter, energy, deltad, tester
 
400
      call flush(6)
 
401
1     format(' iter=',i3,', energy=',f13.8,', deltad=',d9.2,
 
402
     $     ', deltaf=',d9.2)
 
403
c
 
404
      end
 
405
      double precision function dendif(dens, work)
 
406
      implicit double precision (a-h, o-z)
 
407
      include 'cscf.h'
 
408
      dimension dens(nnbfn), work(nnbfn)
 
409
c
 
410
c     compute largest change in density matrix elements
 
411
c
 
412
      denmax = 0.0d0
 
413
      do 10 i = 1, nnbfn
 
414
         denmax = max(denmax, abs(work(i)-dens(i)))
 
415
10    continue
 
416
      dendif = denmax
 
417
c
 
418
      end
 
419
      subroutine prnfin(energy, evals, orbs)
 
420
      implicit double precision (a-h, o-z)
 
421
      include 'cscf.h'
 
422
      dimension evals(nbfn), orbs(nbfn, nbfn)
 
423
c
 
424
c     printout final results
 
425
c
 
426
      write(6,1) energy
 
427
 1    format(//' final energy = ',f16.11//' eigenvalues')
 
428
      call output(evals, 1, min(nbfn,nocc+5), 1, 1, nbfn, 1, 1)
 
429
      write(6,2)
 
430
2     format(//' eigenvectors ')
 
431
      call output(orbs, 1, nbfn, 1, nocc, nbfn, nbfn, 1)
 
432
c
 
433
      end
 
434
      subroutine g(value,i,j,k,l)
 
435
      implicit double precision (a-h, o-z)
 
436
      include 'cscf.h'
 
437
c
 
438
c     compute the two electon integral (ij|kl) over normalized
 
439
c     primitive 1s gaussians
 
440
c
 
441
      f0val = 0.0d0
 
442
      rab2 = (x(i)-x(j))**2 + (y(i)-y(j))**2 + (z(i)-z(j))**2
 
443
      rcd2 = (x(k)-x(l))**2 + (y(k)-y(l))**2 + (z(k)-z(l))**2
 
444
      facij = expnt(i)*expnt(j)/(expnt(i)+expnt(j))
 
445
      fackl = expnt(k)*expnt(l)/(expnt(k)+expnt(l))
 
446
      exijkl = exprjh(- facij*rab2 - fackl*rcd2)
 
447
      denom = (expnt(i)+expnt(j))*(expnt(k)+expnt(l)) *
 
448
     $        sqrt(expnt(i)+expnt(j)+expnt(k)+expnt(l))
 
449
      fac = (expnt(i)+expnt(j))*(expnt(k)+expnt(l)) /
 
450
     $        (expnt(i)+expnt(j)+expnt(k)+expnt(l))
 
451
c
 
452
      xp = (x(i)*expnt(i) + x(j)*expnt(j))/(expnt(i)+expnt(j))
 
453
      yp = (y(i)*expnt(i) + y(j)*expnt(j))/(expnt(i)+expnt(j))
 
454
      zp = (z(i)*expnt(i) + z(j)*expnt(j))/(expnt(i)+expnt(j))
 
455
      xq = (x(k)*expnt(k) + x(l)*expnt(l))/(expnt(k)+expnt(l))
 
456
      yq = (y(k)*expnt(k) + y(l)*expnt(l))/(expnt(k)+expnt(l))
 
457
      zq = (z(k)*expnt(k) + z(l)*expnt(l))/(expnt(k)+expnt(l))
 
458
      rpq2 = (xp-xq)**2 + (yp-yq)**2 + (zp-zq)**2
 
459
c
 
460
      call f0(f0val, fac*rpq2)
 
461
      value = (2.0d0 * pi**2.5d0 / denom) * exijkl * f0val *
 
462
     $    rnorm(i)*rnorm(j)*rnorm(k)*rnorm(l)
 
463
c
 
464
      end
 
465
      subroutine diagon(fock, orbs, evals, work, tester, iter)
 
466
      implicit double precision (a-h, o-z)
 
467
      include 'cscf.h'
 
468
      dimension fock(nnbfn), orbs(nbfn, nbfn), evals(nbfn),
 
469
     $     work(nbfn*nbfn)
 
470
      dimension ilifq(nbfn), work2(nbfn*nbfn)
 
471
c
 
472
      do 10 i = 1, nbfn
 
473
         ilifq(i) = (i-1)*nbfn
 
474
 10   continue
 
475
c
 
476
c     transform fock matrix from AO to MO basis
 
477
c
 
478
      call zsqua(nbfn, fock, work)
 
479
      call dgemm('n', 'n', nbfn, nbfn, nbfn, 1.0d0, 
 
480
     $     work, nbfn, orbs, nbfn, 0.0d0, work2, nbfn)
 
481
      call dgemm('t', 'n', nbfn, nbfn, nbfn, 1.0d0, 
 
482
     $     orbs, nbfn, work2, nbfn, 0.0d0, work, nbfn)
 
483
      call zfold(nbfn, fock, work)
 
484
c
 
485
      iop1 = 1
 
486
      iop2 = 2
 
487
      tester = 0.0d0
 
488
      do 20 i = 2, nbfn
 
489
         do 30 j = 1, i-1
 
490
            tester = max(tester, abs(fock(iky(i)+j)))
 
491
 30      continue
 
492
 20   continue
 
493
c
 
494
      if (tester.gt.0.3d0) then
 
495
        shift = 0.3d0
 
496
      else
 
497
        if (nbfn .gt. 60) then
 
498
           shift = 0.1d0
 
499
        else
 
500
           shift = 0.0d0
 
501
        endif
 
502
      endif    
 
503
      if (iter.ge.2) then
 
504
         do 40 i = nocc+1, nbfn
 
505
            fock(iky(i)+i) = fock(iky(i)+i) + shift
 
506
 40      continue
 
507
      endif
 
508
      thresh = min(0.0001d0,max(1.0d-12,tester*0.01d0))
 
509
      call jacobi(fock, iky, nbfn, orbs, ilifq, nbfn, evals, iop1,
 
510
     *     iop2,thresh)
 
511
      if (iter.ge.2) then
 
512
         do 50 i = nocc+1, nbfn
 
513
            evals(i) = evals(i) - shift
 
514
 50      continue
 
515
      endif
 
516
c
 
517
      end
 
518
      subroutine makeob(orbs, work)
 
519
      implicit double precision (a-h, o-z)
 
520
      include 'cscf.h'
 
521
      include 'msgtypesf.h'
 
522
      dimension orbs(nbfn,nbfn), work(nbfn,nbfn)
 
523
      dimension tmp1(nbfn), tmp2(nbfn)
 
524
c
 
525
c     generate set of orthonormal vectors by orthoging twice
 
526
c     a set of random vectors ... don't do this at home!
 
527
c     ... should really diagonalize the overlap to get sym adaption
 
528
c
 
529
      if (nodeid().eq.0) then
 
530
         call srand48(12345)
 
531
         do 10 i = 1, nbfn
 
532
            do 20 j = 1, nbfn
 
533
               work(j,i) = s(i,j)
 
534
               orbs(j,i) = drand48(0)
 
535
 20         continue
 
536
 10      continue
 
537
         call orthv2(nbfn, orbs, work, tmp1, tmp2)
 
538
         call orthv2(nbfn, orbs, work, tmp1, tmp2)
 
539
      endif
 
540
      call brdcst(99+MSGDBL, orbs, mdtob(nbfn*nbfn), 0)
 
541
c
 
542
      end
 
543
      subroutine denges(dens, work)
 
544
      implicit double precision (a-h, o-z)
 
545
      include 'cscf.h'
 
546
      dimension dens(nnbfn), work(nbfn, nbfn)
 
547
c
 
548
c     Form guess density from superposition of atomic densities in the AO
 
549
c     basis set ... instead of doing the atomic SCF hardwire for this
 
550
c     small basis set for the Be atom.
 
551
c
 
552
      dimension atdens(15,15)
 
553
      data atdens/
 
554
     $     0.000002,0.000027,0.000129,0.000428,0.000950,0.001180,
 
555
     $     0.000457,-0.000270,-0.000271,0.000004,0.000004,0.000004,
 
556
     $     0.000004,0.000004,0.000004,0.000027,0.000102,0.000987,
 
557
     $     0.003269,0.007254,0.009007,0.003492,-0.002099,-0.002108,
 
558
     $     0.000035,0.000035,0.000035,0.000035,0.000035,0.000035,
 
559
     $     0.000129,0.000987,0.002381,0.015766,0.034988,0.043433,
 
560
     $     0.016835,-0.010038,-0.010082,0.000166,0.000166,0.000166,
 
561
     $     0.000166,0.000166,0.000166,0.000428,0.003269,0.015766,
 
562
     $     0.026100,0.115858,0.144064,0.055967,-0.035878,-0.035990,
 
563
     $     0.000584,0.000584,0.000584,0.000584,0.000584,0.000584,
 
564
     $     0.000950,0.007254,0.034988,0.115858,0.128586,0.320120,
 
565
     $     0.124539,-0.083334,-0.083536,0.001346,0.001346,0.001346,
 
566
     $     0.001346,0.001346,0.001346,0.001180,0.009007,0.043433,
 
567
     $     0.144064,0.320120,0.201952,0.159935,-0.162762,-0.162267,
 
568
     $     0.002471,0.002471,0.002471,0.002471,0.002471,0.002471,
 
569
     $     0.000457,0.003492,0.016835,0.055967,0.124539,0.159935,
 
570
     $     0.032378,-0.093780,-0.093202,0.001372,0.001372,0.001372,
 
571
     $     0.001372,0.001372,0.001372,-0.000270,-0.002099,-0.010038,
 
572
     $     -0.035878,-0.083334,-0.162762,-0.093780,0.334488,0.660918,
 
573
     $     -0.009090,-0.009090,-0.009090,-0.009090,-0.009090,-0.009090,
 
574
     $     -0.000271,-0.002108,-0.010082,-0.035990,-0.083536,-0.162267,
 
575
     $     -0.093202,0.660918,0.326482,-0.008982,-0.008982,-0.008981,
 
576
     $     -0.008981,-0.008981,-0.008982,0.000004,0.000035,0.000166,
 
577
     $     0.000584,0.001346,0.002471,0.001372,-0.009090,-0.008982,
 
578
     $     0.000062,0.000124,0.000124,0.000124,0.000124,0.000124,
 
579
     $     0.000004,0.000035,0.000166,0.000584,0.001346,0.002471,
 
580
     $     0.001372,-0.009090,-0.008982,0.000124,0.000062,0.000124,
 
581
     $     0.000124,0.000124,0.000124,0.000004,0.000035,0.000166,
 
582
     $     0.000584,0.001346,0.002471,0.001372,-0.009090,-0.008981,
 
583
     $     0.000124,0.000124,0.000062,0.000124,0.000124,0.000124,
 
584
     $     0.000004,0.000035,0.000166,0.000584,0.001346,0.002471,
 
585
     $     0.001372,-0.009090,-0.008981,0.000124,0.000124,0.000124,
 
586
     $     0.000062,0.000124,0.000124,0.000004,0.000035,0.000166,
 
587
     $     0.000584,0.001346,0.002471,0.001372,-0.009090,-0.008981,
 
588
     $     0.000124,0.000124,0.000124,0.000124,0.000062,0.000124,
 
589
     $     0.000004,0.000035,0.000166,0.000584,0.001346,0.002471,
 
590
     $     0.001372,-0.009090,-0.008982,0.000124,0.000124,0.000124,
 
591
     $     0.000124,0.000124,0.000062/
 
592
c
 
593
      call dfill(nbfn*nbfn,0.0d0,work,1)
 
594
      do 10 iat = 1,natom
 
595
         ioff = (iat-1)*15
 
596
         do 20 i = 1,15
 
597
            do 30 j = 1,15
 
598
               work(ioff+j,ioff+i) = atdens(j,i)*0.5d0
 
599
 30         continue
 
600
 20      continue
 
601
 10   continue
 
602
c
 
603
      call zfold(nbfn,dens,work)
 
604
c
 
605
      end