← Back to branch summary

~ubuntu-branches/ubuntu/trusty/nwchem/trusty-proposed

« back to all changes in this revision

Viewing changes to src/tools/ga-5-1/armci/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-1/armci/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