1C
2c 11 jul 01: header fixed; send to Gijs
3c  7 jun 01: determine runtype, scftype, dft, mp2 from output
4c  6 jun 01: uhf fixed
5c  4 jun 01: prepares MOLDEN Format file from CADPAC output for:
6c  { rhf | uhf }
7c  { ene | opt | hess }
8c  { dft | mp2 }
9
10c
11c Mariusz Klobukowski
12c cad2mol is a minor modification of Chris Lovallo's hon2mol
13c
14C This program will take a CADPAC output file and convert it to
15C Molden Format for use with Molden (as filename.mdn). The program
16C should work for SP energy, geometry optimizations, and Hessians for
17C any single-determinant wavefunction
18
19C For the orbitals to be read correctly, NOPRINT OCCVECTORS must not
20c be used
21
22      program cad2mol
23      implicit double precision(a-e,g-h,o-z)
24      implicit logical(f)
25      character filenm*70
26      parameter(inp=1,iout=2,itmp=3,maxbf=2000)
27      common /flags/ flgrhf,flguhf,flgopn,flgdft,flggvb,flgmp2,flgmp4
28      common /flags2/ flgmc
29      common /rundat/ irun,iwfn,nbf,ndoc,nsoc,maxit,nat,norb
30      dimension occ(maxbf)
31
32C Open necessary files (input, output, and temp)
33
34      call open_files(inp,iout,itmp,filenm)
35
36C Print title onto files
37
38      call write_files(iout,itmp,'[Molden Format]')
39
40C Read and setup parameters of run
41
42      call find_flags(inp,occ)
43      rewind(inp)
44
45C Find and write coordinates and basis set
46
47      call find_coord(inp,iout,itmp)
48
49C Find and write eigenvectors and eigenvalues
50
51      call find_eigen(inp,iout,occ)
52
53C Next (and last for SP energy) is SCF energy convergence
54
55      rewind(inp)
56      call find_conv1(inp,iout)
57
58C End if the run is SP energy, skip to frequencies if run is Hessian
59
60      if (irun .eq. 0) goto 9999
61      if (irun .eq. 2) goto 10
62
63C Find and write out geometry optimization data
64
65      call find_gopt(inp,iout)
66      if (irun .eq. 3) goto 10
67      goto 9999
68
69C Find and write frequencies and normal modes (for Hessian runs)
70
71 10   call find_freq(inp,iout,itmp)
72
73C Close files
74
75 9999 call close_files(inp,iout,itmp)
76
77C That's it! It's over! Bye now!
78
79      end
80C-----------------------------------------------------------------------
81      subroutine open_files(inp,iout,itmp,filenm)
82      character filenm*70
83
84C Open input file (CADPAC output file)
85
86      write(6,*) 'Opening necessary files...'
87      call getarg(1,filenm)
88      if(filenm.eq.' ') then
89        write(6,fmt='('' Enter the input file name (with .out): '',$)')
90        read (5,fmt='(a70)') filenm
91      endif
92      open(unit=inp,file=filenm,access='SEQUENTIAL',
93     1     status='OLD',form='FORMATTED')
94
95C Open output and temp files
96
97      do i=1,67
98        if (filenm(i:i+3) .eq. '.out') filenm(i:i+3)='.mdn'
99      enddo
100      open(unit=iout,file=filenm,access='SEQUENTIAL',
101     1     status='UNKNOWN',form='FORMATTED')
102      open(unit=itmp,access='SEQUENTIAL',
103     1     status='UNKNOWN',name='TEMP',form='FORMATTED')
104c    1     status='SCRATCH',form='FORMATTED')
105      return
106      end
107C-----------------------------------------------------------------------
108      subroutine find_flags(inp,occ)
109      implicit double precision(a-e,g-h,o-z)
110      implicit logical(f)
111      character line*130
112c     character kyword*26
113      parameter(maxbf=2000)
114      common /flags/ flgrhf,flguhf,flgopn,flgdft,flggvb,flgmp2,flgmp4
115      common /flags2/ flgmc
116      common /rundat/ irun,iwfn,nbf,ndoc,nsoc,maxit,nat,norb
117      dimension occ(*)
118
119C This subroutine searches the CADPAC file for various data pertaining
120C to the run (type of wavefunction, number of basis functions and
121C occupied orbitals, etc.)
122
123      write(6,*) 'Beginning initialization...'
124
125C Initialize flags
126
127      flgrhf=.false.
128      flguhf=.false.
129      flgopn=.false.
130      flgdft=.false.
131      flggvb=.false.
132      flgmp2=.false.
133      flgmp4=.false.
134      flgmc= .false.
135
136C Search for type of run and properties - decode from command line
137C Determine type of wavefunction - decode from command line
138c cad2mol job.out
139c  { rhf | uhf | rohf | gvb | mp2 | mp4 }
140c  { ene | opt | hess }
141c  { dft }
142
143*     inparg=iargc()
144*     write(6,'('' inparg:'',i3)') inparg
145
146      irun = -1
147      iwfn = -1
148*     do narg=2,inparg
149*       call getarg(narg,kyword)
150*       write(6,'(''   narg:'',i3,'' kyword:'',a)') narg,kyword
151*       call up_case(kyword,26)
152*       call last_char(kyword,26,loc)
153*       if     (kyword(1:loc).eq.'RHF') then
154*         flgrhf=.true.
155*         iwfn=0
156*       else if(kyword(1:loc).eq.'ROHF') then
157*         flguhf=.true.
158*         flgopn=.true.
159*         iwfn=0
160*       else if(kyword(1:loc).eq.'GVB') then
161*         flguhf=.true.
162*         flgopn=.true.
163*         iwfn=0
164*       else if(kyword(1:loc).eq.'UHF') then
165*         flguhf=.true.
166*         flgopn=.true.
167*         iwfn=0
168*       else if(kyword(1:loc).eq.'MP2') then
169*         iwfn=5
170*         flgmp2=.true.
171*       else if(kyword(1:loc).eq.'MP4') then
172*         iwfn=6
173*         flgmp4=.true.
174*       else if(kyword(1:loc).eq.'CAS') then
175*         iwfn=1
176*       else if(kyword(1:loc).eq.'DFT') then
177*         flgdft=.true.
178*       else if(kyword(1:loc).eq.'ENE') then
179*         irun=0
180*       else if(kyword(1:loc).eq.'OPT') then
181*         irun=1
182*       else if(kyword(1:loc).eq.'HESS') then
183*         irun=2
184*       else if(kyword(1:loc).eq.'OPHS') then
185*         irun=3
186*       else
187*         write(6,*) 'wrong parameter on command line'
188*         call error
189*       endif
190*     enddo
191
192      rewind(inp)
193      ifound=0
194      call search(inp,line,'RUN TYPE = ENERGY',ifound,0)
195      if(ifound.eq.1) then
196        write(6,*) 'RUN TYPE = ENERGY'
197        irun=0
198      endif
199      rewind(inp)
200      ifound=0
201      call search(inp,line,'RUN TYPE = OPTIMIZE',ifound,0)
202      if(ifound.eq.1) then
203        write(6,*) 'RUN TYPE = OPTIMIZE'
204        irun=1
205      endif
206      rewind(inp)
207      ifound=0
208      call search(inp,line,'RUN TYPE = SECDER',ifound,0)
209      if(ifound.eq.1.and.irun.eq.-1) then
210        irun=2
211        write(6,*) 'RUN TYPE = SECDER -- only'
212      endif
213      if(ifound.eq.1.and.irun.eq.1) then
214        irun=3
215        write(6,*) 'RUN TYPE = SECDER -- after OPTIMIZE'
216      endif
217      rewind(inp)
218      ifound=0
219      call search(inp,line,'RUN TYPE = FORCE',ifound,0)
220      if(ifound.eq.1.and.irun.eq.-1) then
221        irun=2
222        write(6,*) 'RUN TYPE = FORCE -- only'
223      endif
224      if(ifound.eq.1.and.irun.eq.1) then
225        irun=3
226        write(6,*) 'RUN TYPE = FORCE -- after OPTIMIZE'
227      endif
228      rewind(inp)
229      ifound=0
230      call search(inp,line,'SCF TYPE = CLOSED',ifound,0)
231      if(ifound.eq.1) then
232        flgrhf=.true.
233        iwfn=0
234        write(6,*) 'SCF TYPE = CLOSED'
235      endif
236      rewind(inp)
237      ifound=0
238      call search(inp,line,'UHF Calculation',ifound,0)
239      if(ifound.eq.1) then
240        iwfn=0
241        flguhf=.true.
242        flgopn=.true.
243        write(6,*) 'UHF Calculation'
244      endif
245      rewind(inp)
246      ifound=0
247      call search(inp,line,'CLOSED-SHELL KOHN-SHAM',ifound,0)
248      if(ifound.eq.1) then
249        flgdft=.true.
250        write(6,*) 'CLOSED-SHELL KOHN-SHAM'
251      endif
252      rewind(inp)
253      ifound=0
254      call search(inp,line,'Unrestricted KOHN-SHAM',ifound,0)
255      if(ifound.eq.1) then
256        iwfn=0
257        flgdft=.true.
258        flguhf=.true.
259        flgopn=.true.
260        write(6,*) 'Unrestricted KOHN-SHAM'
261      endif
262      rewind(inp)
263      ifound=0
264      call search(inp,line,'RHF MP2 Calculation',ifound,0)
265      if(ifound.eq.1) then
266        iwfn=5
267        flgmp2=.true.
268        write(6,*) 'RHF MP2 Calculation'
269      endif
270      rewind(inp)
271      ifound=0
272      call search(inp,line,'UHF MP2 Calculation',ifound,0)
273      if(ifound.eq.1) then
274        iwfn=5
275        flgmp2=.true.
276        flguhf=.true.
277        flgopn=.true.
278        write(6,*) 'UHF MP2 Calculation'
279      endif
280
281      write(6,fmt='('' irun   ='',i2)') irun
282      write(6,fmt='('' iwfn   ='',i2)') iwfn
283      write(6,fmt='('' flgrhf ='',L2)') flgrhf
284      write(6,fmt='('' flguhf ='',L2)') flguhf
285      write(6,fmt='('' flgopn ='',L2)') flgopn
286      write(6,fmt='('' flgdft ='',L2)') flgdft
287      write(6,fmt='('' flgmp2 ='',L2)') flgmp2
288
289      if (irun.ne.0. and.
290     1    irun.ne.1 .and.
291     1    irun.ne.2 .and.
292     1    irun.ne.3) then
293        write(6,*)  'Only RUNTYPs=0-3 (SP energy, geom opt, and Hessian)
294     1 are currently implemented!'
295        call error
296      endif
297
298      if (iwfn .ne. 0 .and. iwfn .ne. 5 .and. iwfn .ne. 1) then
299       write(6,*)   'Only single-determinant and MCSCF wavefunctions are
300     1 currently implemented!'
301       call error
302      endif
303
304      rewind (inp)
305
306C Find number of basis functions
307
308      call search(inp,line,'Number of Basis Functions',ifound,-1)
309      read(unit=line,fmt='(46x,i5)') nbf
310      write(*,*) 'nbf  =',nbf
311      if (nbf .ge. maxbf) then
312        write(6,*)   'Error code 5! Number of basis functions is greater
313     1 than parameter MAXBF!'
314        write(6,*) 'Check parameter MAXBF and run program again.'
315        call error
316      endif
317
318
319C Determine type of run and orbital occupancies
320
321      if (iwfn .eq. 1 .or. iwfn .eq. 10 .or. iwfn .eq. 12) flgmc=.true.
322      if (flgmc .or. iwfn .eq. 2) then
323        call search(inp,line,'NUMBER OF ORBITALS',ifound,-1)
324        read(unit=line,fmt='(27x,i5)') norb
325        call search(inp,line,'# OF CORE   ORBITALS',ifound,-1)
326        read(unit=line,fmt='(26x,i5)') ndoc
327        goto 20
328      endif
329
330C The program does not search for MP2/4 lines, it sets the
331C flags according to iwfn
332
333 10   if (iwfn .eq. 5) flgmp2=.true.
334      if (iwfn .eq. 6) flgmp4=.true.
335      if(flguhf) then
336        call search(inp,line,'occupied alpha orbitals',ifound,-1)
337        read(unit=line,fmt='(46x,i5)') ndoc
338      else
339        call search(inp,line,
340     1       'Number of doubly occupied orbitals',ifound,-1)
341        read(unit=line,fmt='(46x,i5)') ndoc
342      endif
343      write(*,*) 'ndoc =',ndoc
344      if (.not. flgopn) goto 20
345      if (flguhf) then
346        call search(inp,line,'occupied beta orbitals',ifound,-1)
347        read(unit=line,fmt='(46x,i5)') nsoc
348        write(*,*) 'nsoc=',nsoc
349        goto 20
350      endif
351
352 20   if (.not. flguhf) then
353        do i=1,ndoc
354          occ(i)=2.0
355        enddo
356        if (nsoc .gt. 0) then
357          do i=ndoc+1,ndoc+1+nsoc
358            occ(i)=1.0
359          enddo
360          do i=ndoc+1+nsoc,nbf
361            occ(i)=0.0
362          enddo
363        else
364          do i=ndoc+1,nbf
365            occ(i)=0.0
366          enddo
367        endif
368        goto 30
369      else
370        do i=1,ndoc
371          occ(i)=1.0
372        enddo
373        do i=ndoc+1,nbf
374          occ(i)=0.0
375        enddo
376      endif
377
378C Determine maximum number of iterations (SCF or CI)
379
380 30   if (iwfn .eq. 2) then
381        call search(inp,line,'MAX. NUMB. OF ITERATIONS',ifound,-1)
382        read(unit=line,fmt='(29x,i6)') maxit
383      else
384        if(flguhf) then
385         call search(inp,line,'Maximum Number of Iterations',ifound,-1)
386        else
387         call search(inp,line,'Maximum number of iterations',ifound,-1)
388        endif
389        read(unit=line,fmt='(32x,i6)') maxit
390        write(*,*) 'maxit=',maxit
391      endif
392      maxit=maxit+1
393
394      write(6,*) 'Initialization done!'
395      return
396      end
397C-----------------------------------------------------------------------
398      subroutine find_coord(inp,iout,itmp)
399      implicit double precision(a-e,g-h,o-z)
400      implicit logical(f)
401      common /rundat/ irun,iwfn,nbf,ndoc,nsoc,maxit,nat,norb
402      character line*130,lintmp*130,atmnam*8
403      parameter(maxat=100)
404      dimension atmnam(maxat)
405
406C This subroutine finds and writes out the coordinates of the atoms
407C It searches for CHARGE, and then skips the --- and blank lines
408C It then calls find_bas to write out the basis set
409
410      write(6,*) 'Reading coordinates...'
411      call write_files(iout,itmp,'[Atoms] AU')
412      call search(inp,line,'Atom     Nuclear',ifound,-1)
413      call skip(inp,2)
414      i=1
415 10   call read_line(inp,line)
416      read(unit=line,fmt='(4x,a8)') atmnam(i)
417      if (atmnam(i) .eq. '--------') goto 20
418      read(unit=line,fmt='(4x,a8,2x,f5.1,3f15.8)') atmnam(i),chrg
419     1 ,x,y,z
420      ichrg=int(chrg)
421      write(*,fmt='(1x,a8,i3,i3,3(2x,f12.8))') atmnam(i),i,ichrg,x,y,z
422      write(unit=lintmp,fmt='(1x,a8,i3,i3,3(2x,f12.7))') atmnam(i),i
423     1      ,ichrg,x,y,z
424      call write_files(iout,itmp,lintmp)
425      i=i+1
426      goto 10
427 20   nat=i-1
428      write(*,*) 'nat=',nat
429      write(6,*) 'Done writing coordinates!'
430      call find_bas(inp,iout,itmp,atmnam)
431      return
432      end
433C-----------------------------------------------------------------------
434      subroutine find_bas(inp,iout,itmp,atmnam)
435      implicit double precision(a-e,g-h,o-z)
436      implicit logical(f)
437      common /rundat/ irun,iwfn,nbf,ndoc,nsoc,maxit,nat,norb
438      character line*130,lintmp*130,temp8*8,temp4*4,atmnam*8,shltyp*4
439      parameter(mxprim=30)
440      dimension atmnam(*),coef1(mxprim),coef2(mxprim),exp(mxprim)
441
442C This subroutine reads and writes out the basis set, including looping
443C over the non-symmetry unique atoms
444
445      write(6,*) 'Reading basis set...'
446      call skip(inp,6)
447      call write_files(iout,itmp,'[GTO]')
448      do j=1,nat
449*     write(*,*) j,atmnam(j)
450c ... reset cgtf count for each new atom
451        ikount=0
452        write(unit=lintmp,fmt='(i3,a5)') j, ' 0   '
453        call write_files(iout,itmp,lintmp)
454        call read_line(inp,line)
455*       write(*,*) '>'//line(1:76)
456        read(unit=line,fmt='(1x,a8)') temp8
457*       write(*,*) 'temp8'//temp8
458        if (temp8 .ne. atmnam(j)) then
459C Either an error has occurred, or this is not a symmetry unique atom.
460C Checking for an atom with the same name and copying basis set...
461          backspace(inp)
462          temp8=atmnam(j)
463          do k=1,j-1
464            if (temp8 .eq. atmnam(k)) then
465              rewind(itmp)
466              call search(itmp,line,atmnam(k),ifound,-1)
467              read(unit=line,fmt='(10x,i3)') n
468              write(unit=temp8,fmt='(i3,a5)') n, ' 0   '
469              call search(itmp,line,temp8,ifound,-1)
470              write(unit=temp8,fmt='(i3,a5)') n+1, ' 0   '
471 10           call read_line(itmp,line)
472*             write(*,*) '>'//line(1:76)
473              if (line(1:8) .ne. temp8) then
474                call write_files(iout,0,line)
475                goto 10
476              endif
477              call ffwrd(itmp)
478              goto 40
479            endif
480          enddo
481
482C Nope! There was an error! Stopping...
483
484          write(6,*) 'Error code 1 reading basis set!'
485          call error
486        endif
487
488        l=1
489 20     call read_line(inp,line)
490*       write(*,*) '>'//line(1:76)
491        if(line(1:6).eq.' Model') call read_line(inp,line)
492        read(unit=line,fmt='(12x,i2,5x,a4,6x,f16.6,f15.6)')
493     1     jkount,temp4,exp(l),coef1(l)
494*       write(*,*) 'l ikount jkount',l,ikount,jkount,exp(l)
495        if(jkount.eq.0) then
496*          backspace(inp)
497           ikount=jkount
498           goto 30
499        endif
500        if(ikount.gt.0 .and. jkount.ne.ikount) then
501          ikount=jkount
502          backspace(inp)
503          goto 30
504        else
505          ikount=jkount
506        endif
507        shltyp=temp4
508        if (shltyp(1:1) .eq. 'L') then
509          read(unit=line,fmt='(19x,a4,6x,f16.6,2f15.6)')
510     1         shltyp,exp(l),coef1(l),coef2(l)
511        endif
512        l=l+1
513        goto 20
514 30     nprim=l-1
515*       write(*,*) 'nprim exp',nprim,exp(nprim)
516
517C Write basis set
518
519        if (shltyp(1:1) .ne. 'L') then
520          write(unit=lintmp,fmt='(a2,i5,a5)') shltyp(1:1),nprim, ' 1.00'
521          call write_files(iout,itmp,lintmp)
522          do l=1,nprim
523            write(unit=lintmp,fmt='(2(1x,d17.10))') exp(l),coef1(l)
524            call write_files(iout,itmp,lintmp)
525          enddo
526        elseif (shltyp(1:1) .eq. 'L') then
527          write(unit=lintmp,fmt='(1x,a2,i4,a5)') 'SP',nprim,
528     1          ' 1.00'
529          call write_files(iout,itmp,lintmp)
530          do l=1,nprim
531            write(unit=lintmp,fmt='(3(1x,d17.10))') exp(l),coef1(l)
532     1            ,coef2(l)
533            call write_files(iout,itmp,lintmp)
534          enddo
535        else
536          write(6,*) 'Error code 3 writing basis set!'
537          call error
538        endif
539        if(jkount.eq.0) goto 40
540        l=1
541        goto 20
542 40     call blank_line(iout)
543      enddo
544      call blank_line(iout)
545      write(6,*) 'Done writing basis set!'
546      return
547      end
548C-----------------------------------------------------------------------
549      subroutine find_eigen(inp,iout,occ)
550      implicit double precision(a-e,g-h,o-z)
551      implicit logical(f)
552      common /flags/ flgrhf,flguhf,flgopn,flgdft,flggvb,flgmp2,flgmp4
553      common /flags2/ flgmc
554      common /rundat/ irun,iwfn,nbf,ndoc,nsoc,maxit,nat,norb
555      character line*130,lintmp*130
556      parameter(maxbf=2000)
557      parameter (zero=0.0d+00, one=1.0d+00)
558      dimension occ(*),orbs(maxbf,maxbf),eval(maxbf)
559
560C This subroutine finds and writes out the molecular orbitals and
561C orbital energies, along with their occupancies
562
563      write(6,*) 'Reading orbitals and eigenvalues...'
564      call write_files(iout,0,'[MO]')
565c ... 1st MO output
566      first_SCF=.true.
567      if(flguhf) then
568        call search(inp,line,'Alpha-spin Orbitals',ifound,-1)
569      else
570        call search(inp,line,'Molecular Orbitals',ifound,-1)
571      endif
572c ... 2nd MO output - used if opt+hess in the same run
573      if (irun .eq. 1 .or. irun.eq.2 .or. irun.eq.3) then
574        if(flguhf) then
575          call search(inp,line,'Alpha-spin orbitals',ifound,0)
576        else
577          call search(inp,line,'Molecular orbitals',ifound,0)
578        endif
579      endif
580      first_SCF=.false.
581      if(ifound.eq.0) then
582        write(*,*) '... only hessian'
583        first_SCF=.true.
584        rewind(inp)
585        if(flguhf) then
586          call search(inp,line,'Alpha-spin Orbitals',ifound,-1)
587        else
588          call search(inp,line,'Molecular Orbitals',ifound,-1)
589        endif
590      endif
591c ... locate the 1st eigenvalue
592      call search(inp,line,'     1    ',ifound,-1)
593      backspace(inp)
594c ... read the number of MOs printed
595      found=.false.
596      nmos=0
597      do while (.not.found)
598        call read_line(inp,line)
599        if(line(1:20).ne.'                    ' .and.
600     1     line(1:13).ne.' Eigenvectors') then
601          nmos=nmos+1
602        else
603          found=.true.
604          backspace(inp)
605        endif
606      enddo
607      write(*,*) 'nmos=',nmos
608      write(*,*) 'first_SCF ',first_SCF
609      if(.not.flguhf) then
610        if(irun.eq.0.or.first_SCF) call skip(inp,2)
611      endif
612      if(flguhf.and.first_SCF) call skip(inp,1)
613      if (flggvb) call skip(inp,1)
614      j=1
615 10   k=j+4
616      call skip(inp,2)
617      call read_line(inp,line)
618      read(unit=line,fmt='(12x,5(f13.6))') (eval(i), i=j,k)
619      write(*,fmt='(5(i4,f13.6))') (i,eval(i), i=j,k)
620      call skip(inp,1)
621      do l=1,nbf
622        call read_line(inp,line)
623        read(unit=line,fmt='(12x,5(f13.6))') (orbs(i,l), i=j,k)
624      enddo
625      if (k .ge. nmos) then
626        do i=1,nmos
627          write(unit=lintmp,fmt='(a5,2x,f9.4)') ' Ene=',eval(i)
628          call write_files(iout,0,lintmp)
629          call write_files(iout,0,' Spin= Alpha')
630          if(flguhf.and.i.gt.ndoc) occ(i)=zero
631          write(unit=lintmp,fmt='(a7,f11.6)') ' Occup=',occ(i)
632          call write_files(iout,0,lintmp)
633          do j=1,nbf
634            write(unit=lintmp,fmt='(i4,f11.6)') j,orbs(i,j)
635            call write_files(iout,0,lintmp)
636          enddo
637        enddo
638        if (.not. flguhf) goto 20
639
640C Now read beta orbitals (for UHF runs only)
641
642        write(6,*) 'Done writing alpha orbitals...'
643        call find_beta(inp,iout,occ)
644        write(6,*) 'Done writing beta  orbitals...'
645        goto 20
646      endif
647      j=j+5
648      goto 10
649 20   write(6,*) 'Done writing orbitals and eigenvalues!'
650      return
651      end
652C-----------------------------------------------------------------------
653      subroutine find_beta(inp,iout,occ)
654
655C This subroutine finds and write the beta molecular orbitals and
656C orbital energies for a UHF wavefunction
657
658      implicit double precision(a-e,g-h,o-z)
659      implicit logical(f)
660      common /flags/ flgrhf,flguhf,flgopn,flgdft,flggvb,flgmp2,flgmp4
661      common /flags2/ flgmc
662      common /rundat/ irun,iwfn,nbf,ndoc,nsoc,maxit,nat,norb
663      character line*130,lintmp*130
664      parameter(maxbf=2000)
665      dimension occ(*),orbs(maxbf,maxbf),eval(maxbf)
666
667c ... locate start of current beta's
668      call search(inp,line,'Beta-spin',ifound,-1)
669c ... locate 1st beta MO
670      call search(inp,line,'     1     ',ifound,-1)
671      backspace(inp)
672c ... read the number of MOs printed
673      found=.false.
674      nmos=0
675      do while (.not.found)
676        call read_line(inp,line)
677        if(line(1:20).ne.'                    '.and.
678     1     line(1:13).ne.' Eigenvectors') then
679          nmos=nmos+1
680        else
681          found=.true.
682          backspace(inp)
683        endif
684      enddo
685      write(*,*) 'nmos=',nmos
686      call read_line(inp,line)
687      if(line(1:13).ne.' Eigenvectors') backspace(inp)
688      j=1
689 10   k=j+4
690      call skip(inp,2)
691      call read_line(inp,line)
692      read(unit=line,fmt='(12x,5(f13.6))') (eval(i), i=j,k)
693      write(*,fmt='(5(i4,f13.6))') (i,eval(i), i=j,k)
694      call skip(inp,1)
695      do l=1,nbf
696        call read_line(inp,line)
697        read(unit=line,fmt='(12x,5(f13.6))') (orbs(i,l), i=j,k)
698      enddo
699      if (k .ge. nmos) then
700        do i=1,nmos
701          write(unit=lintmp,fmt='(a5,2x,f9.4)') ' Ene=',eval(i)
702          call write_files(iout,0,lintmp)
703          call write_files(iout,0,' Spin= Beta')
704          if(flguhf.and.i.gt.nsoc) occ(i)=zero
705          write(unit=lintmp,fmt='(a7,f11.6)') ' Occup=',occ(i)
706          call write_files(iout,0,lintmp)
707          do j=1,nbf
708            write(unit=lintmp,fmt='(i4,f11.6)') j,orbs(i,j)
709            call write_files(iout,0,lintmp)
710          enddo
711        enddo
712        goto 20
713      endif
714      j=j+5
715      goto 10
716 20   return
717      end
718C----------------------------------------------------------------------
719      subroutine find_eigmc(inp,iout,occ)
720      implicit double precision(a-e,g-h,o-z)
721      implicit logical(f)
722      common /flags/ flgrhf,flguhf,flgopn,flgdft,flggvb,flgmp2,flgmp4
723      common /flags2/ flgmc
724      common /rundat/ irun,iwfn,nbf,ndoc,nsoc,maxit,nat,norb
725      character line*130,lintmp*130
726      parameter(maxbf=2000)
727      dimension occ(*),orbs(maxbf,maxbf),eval(maxbf)
728
729C This subroutine finds and writes out the natural orbitals and their
730C occupancies (energies are all set to 0 if the Fock energies are not
731C printed)
732
733      do i=1,nbf
734        eval(i)=0.0
735      enddo
736      rewind(inp)
737      call search(inp,line,'EIGENVALUES OF -MC-',ifound,0)
738      if (ifound .eq. 1) then
739        call skip(inp,2)
740        j=1
741 10     k=j+5
742        call read_line(inp,line)
743        read(unit=line,fmt='(6(10x,f11.4))') (eval(i), i=j,k)
744        j=j+6
745        if (k .lt. nbf) goto 10
746      endif
747      rewind(inp)
748      call write_files(iout,0,'[MO]')
749      if (irun .ne. 1) goto 30
750      n=0
751 20   call search(inp,line,'CAN.CORE + ) NAT.VAL',ifound,0)
752      if (ifound .eq. 1) then
753        n=n+1
754        goto 20
755      endif
756      rewind(inp)
757      do i=1,n
758        call search(inp,line,'CAN.CORE + ) NAT.VAL',ifound,-1)
759      enddo
760      goto 40
761 30   call search(inp,line,'CAN.CORE + ) NAT.VAL',ifound,-1)
762 40   call skip(inp,9)
763      call read_line(inp,line)
764      j=1
765 50   k=j+6
766      read(unit=line,fmt='(15x,7(f15.10))') (occ(i), i=j,k)
767      call skip(inp,2)
768      do l=1,nbf
769        call read_line(inp,line)
770        read(unit=line,fmt='(15x,7(f15.10))') (orbs(i,l), i=j,k)
771      enddo
772      if (k .ge. norb) then
773        do i=1,norb
774          write(unit=lintmp,fmt='(a5,2x,f9.4)') ' Ene=',eval(i)
775          call write_files(iout,0,lintmp)
776          call write_files(iout,0,' Spin= Alpha')
777          write(unit=lintmp,fmt='(a7,f11.6)') ' Occup=',occ(i)
778          call write_files(iout,0,lintmp)
779          do j=1,nbf
780            write(unit=lintmp,fmt='(i4,f11.6)') j,orbs(i,j)
781            call write_files(iout,0,lintmp)
782          enddo
783        enddo
784        goto 60
785      endif
786      j=j+7
787      call skip(inp,8)
788      call read_line(inp,line)
789      goto 50
790 60   return
791      end
792C-----------------------------------------------------------------------
793      subroutine find_eigci(inp,iout,occ)
794      implicit double precision(a-e,g-h,o-z)
795      implicit logical(f)
796      common /flags/ flgrhf,flguhf,flgopn,flgdft,flggvb,flgmp2,flgmp4
797      common /flags2/ flgmc
798      common /rundat/ irun,iwfn,nbf,ndoc,nsoc,maxit,nat,norb
799      character line*130,lintmp*130
800      parameter(maxbf=2000)
801      dimension occ(*),orbs(maxbf,maxbf),eval(maxbf)
802
803C This subroutine finds and writes out the natural orbitals and their
804C occupancies (energies are all set to 0)
805
806      do i=1,nbf
807        eval(i)=0.0
808      enddo
809      rewind(inp)
810      call write_files(iout,0,'[MO]')
811      if (irun .ne. 1) goto 30
812      n=0
813 20   call search(inp,line,'NATURAL ORBITALS IN ATOMIC',ifound,0)
814      if (ifound .eq. 1) then
815        n=n+1
816        goto 20
817      endif
818      rewind(inp)
819      do i=1,n
820        call search(inp,line,'NATURAL ORBITALS IN ATOMIC',ifound,-1)
821      enddo
822      goto 40
823 30   call search(inp,line,'NATURAL ORBITALS IN ATOMIC',ifound,-1)
824 40   call skip(inp,6)
825      call read_line(inp,line)
826      j=1
827 50   k=j+6
828      read(unit=line,fmt='(15x,7(f15.10))') (occ(i), i=j,k)
829      call skip(inp,2)
830      do l=1,nbf
831        call read_line(inp,line)
832        read(unit=line,fmt='(15x,7(f15.10))') (orbs(i,l), i=j,k)
833      enddo
834      if (k .ge. norb) then
835        do i=1,norb
836          write(unit=lintmp,fmt='(a5,2x,f9.4)') ' Ene=',eval(i)
837          call write_files(iout,0,lintmp)
838          call write_files(iout,0,' Spin= Alpha')
839          write(unit=lintmp,fmt='(a7,f11.6)') ' Occup=',occ(i)
840          call write_files(iout,0,lintmp)
841          do j=1,nbf
842            write(unit=lintmp,fmt='(i4,f11.6)') j,orbs(i,j)
843            call write_files(iout,0,lintmp)
844          enddo
845        enddo
846        goto 60
847      endif
848      j=j+7
849      call skip(inp,5)
850      call read_line(inp,line)
851      goto 50
852 60   return
853      end
854C-----------------------------------------------------------------------
855      subroutine find_conv1(inp,iout)
856      implicit double precision(a-e,g-h,o-z)
857      implicit logical(f)
858      common /flags/ flgrhf,flguhf,flgopn,flgdft,flggvb,flgmp2,flgmp4
859      common /flags2/ flgmc
860      common /rundat/ irun,iwfn,nbf,ndoc,nsoc,maxit,nat,norb
861      character line*130,lintmp*130,temp4*4
862      parameter(maxits=101)
863      dimension etot(maxits)
864
865C This subroutine finds and writes out the SCF energy convergence
866C (the first one)
867
868      write(6,*) 'Reading energy convergence...'
869      if(flgdft) then
870        write(6,*) '... DFT calculations - no energy convergence'
871        return
872      endif
873      if (maxit .ge. maxits) then
874        write(6,*)   'Error code 5! Number of iterations is greater than
875     1 parameter MAXITS!'
876        write(6,*) 'Check parameter MAXITS and run program again.'
877        call error
878      endif
879      if (flgmc) then
880        call search(inp,line,' CYCLE  TOTAL ENERGY',ifound,-1)
881      else if (iwfn .eq. 2) then
882        call search(inp,line,' ITER.    MAX.DEV.',ifound,-1)
883        call skip(inp,1)
884      else
885*       call search(inp,line,'Cycle      Total energy',ifound,-1)
886        call search(inp,line,'Cycle     ',ifound,-1)
887        if(.not.flguhf) call skip(inp,1)
888      endif
889      i=0
890      finish=.false.
891      do while(.not.finish)
892        call read_line(inp,line)
893        if (iwfn .eq. 2) then
894          read(unit=line,fmt='(a4,23x,f16.9)') temp4,etot(i)
895        else
896          if(line(1:5).ne.'     '.and.line(1:4).ne.' Den'
897     1                           .and.line(1:5).ne.' Last'
898     1                           .and.line(1:5).ne.' HOMO') then
899            i=i+1
900            read(unit=line,fmt='(5x,f18.9)') etot(i)
901*           write(*,*) 'i e',i,etot(i)
902          endif
903        endif
904        if (line(1:4) .eq. ' Den') then
905          finish=.true.
906        endif
907      enddo
908      ncyc=i
909 10   call write_files(iout,0,'[SCFCONV]')
910      write(unit=lintmp,fmt='(a22,i5)') 'scf-first    1 THROUGH',ncyc
911      call write_files(iout,0,lintmp)
912      do i=1,ncyc
913       write(unit=lintmp,fmt='(f14.6)') etot(i)
914        call write_files(iout,0,lintmp)
915      enddo
916      write(6,*) 'Finished writing energy convergence!'
917      return
918      end
919C-----------------------------------------------------------------------
920      subroutine find_gopt(inp,iout)
921
922C This subroutine finds and writes the data pertaining to the geometry
923C optimization process.
924
925      implicit double precision(a-e,g-h,o-z)
926      implicit logical(f)
927      common /flags/ flgrhf,flguhf,flgopn,flgdft,flggvb,flgmp2,flgmp4
928      common /flags2/ flgmc
929      common /rundat/ irun,iwfn,nbf,ndoc,nsoc,maxit,nat,norb
930      character line*130,lintmp*130,temp4*4,temp8*8,atmnam*8
931      parameter(maxstp=200,bohr2a=5.29177249D-01,maxits=101,maxat=100)
932      dimension estep(maxstp),xmxfor(maxstp),rmsfor(maxstp)
933      dimension ecyc(maxits),atmnam(maxat),geom(3,maxat,maxstp)
934*     dimension grdint(5*maxat)
935
936      write(6,*) 'Reading geometry optimization information...'
937
938c ... DFT calculations do not show SCF convergence
939c ... irrelevant for MP2
940      if(flgdft .or. flgmp2) goto 25
941      rewind(inp)
942      nstep=0
943 10   if (flgmc) then
944        call search(inp,line,' CYCLE  TOTAL ENERGY',ifound,0)
945      else if (iwfn .eq. 0) then
946        if(flguhf) then
947          call search(inp,line,' UHF Calculation',ifound,0)
948        else
949          call search(inp,line,' CLOSED-SHELL SCF CALCULATION',ifound,0)
950        endif
951      else
952        call search(inp,line,' Cycle     ',ifound,0)
953      endif
954      if (ifound .eq. 1) then
955        nstep=nstep+1
956        goto 10
957      endif
958      write(*,*) 'nstep=',nstep
959
960c ... position to the last SCF
961      rewind(inp)
962      do i=1,nstep
963        if (flgmc) then
964          call search(inp,line,' CYCLE  TOTAL ENERGY',ifound,-1)
965        else if (iwfn .eq. 0) then
966          if(flguhf) then
967          call search(inp,line,' UHF Calculation',ifound,0)
968          else
969          call search(inp,line,' CLOSED-SHELL SCF CALCULATION',ifound,0)
970          endif
971          call skip(inp,1)
972        else
973          call search(inp,line,' CYCLE   TOTAL ENERGY',ifound,-1)
974          call skip(inp,1)
975        endif
976      enddo
977c ... located last SCF; find 1st scf cycle
978      call search(inp,line,'  1  ',ifound,-1)
979      backspace(inp)
980      i=0
981      finish=.false.
982      do while(.not.finish)
983        call read_line(inp,line)
984        if (iwfn .eq. 2) then
985          read(unit=line,fmt='(a4,23x,f16.9)') temp4,ecyc(i)
986        else
987          if(line(1:5).ne.'     '.and.line(1:4).ne.' Den'
988     1                           .and.line(1:5).ne.' Last') then
989            i=i+1
990            read(unit=line,fmt='(5x,f18.9)') ecyc(i)
991            write(*,*) 'i e',i,ecyc(i)
992          endif
993        endif
994        if (line(1:4) .eq. ' Den') then
995          finish=.true.
996        endif
997      enddo
998      ncyc=i
999      write(*,*) 'ncyc=',ncyc
1000c ... write scf convergence - last SCF
1001 20   write(unit=lintmp,fmt='(a21,i5)') 'scf-last    1 THROUGH',ncyc
1002      call write_files(iout,0,lintmp)
1003      do i=1,ncyc
1004        write(unit=lintmp,fmt='(f14.6)') ecyc(i)
1005        call write_files(iout,0,lintmp)
1006      enddo
1007
1008  25  continue
1009
1010C The program only uses those steps that have a progress line at the end
1011C of them (one with NSERCH, etc.), so nstep above may not be the nstep
1012C used later (see line 30)
1013
1014      rewind(inp)
1015      nstep=0
1016 11   continue
1017        call search(inp,line,' RMS of gradient    ',ifound,0)
1018        if (ifound .eq. 1) then
1019          nstep=nstep+1
1020          goto 11
1021        endif
1022      write(*,*) 'n RMS=',nstep
1023
1024      rewind(inp)
1025      do i=1,nstep
1026        call search(inp,line,' RMS of gradient',ifound,0)
1027        if (ifound .eq. 0)  goto 30
1028        call backs(inp,line,' Molecular Geometry',ifound,0)
1029        call skip(inp,2)
1030        do j=1,nat
1031          call read_line(inp,line)
1032          read(unit=line,fmt='(6x,a8,2x,3f20.10)') atmnam(j),
1033     1         geom(1,j,i),geom(2,j,i), geom(3,j,i)
1034          write(*,*) j,geom(1,j,i),geom(2,j,i), geom(3,j,i)
1035          do k=1,3
1036            geom(k,j,i) = geom(k,j,i)*bohr2a
1037          enddo
1038        enddo
1039        call search(inp,line,' Current energy    ',ifound,-1)
1040        read(unit=line,fmt='(32x,f19.10)') estep(i)
1041        write(*,*) 'estep(i)',i,estep(i)
1042        call skip(inp,2)
1043        call read_line(inp,line)
1044        read(unit=line,fmt='(32x,f13.6)') rmsfor(i)
1045        write(*,*) 'rmsfor(i)',rmsfor(i)
1046        call read_line(inp,line)
1047        read(unit=line,fmt='(32x,f13.6)') xmxfor(i)
1048        write(*,*) 'xmxfor(i)',xmxfor(i)
1049      enddo
1050      nstep=i
1051      write(*,*) 'nstep',nstep
1052
1053C Now write everything out
1054
1055 30   nstep=i-1
1056      call write_files(iout,0,'[GEOCONV]')
1057      call write_files(iout,0,'energy')
1058      do i=1,nstep
1059        write(unit=lintmp,fmt='(1x,f13.6)') estep(i)
1060        call write_files(iout,0,lintmp)
1061      enddo
1062      call write_files(iout,0,'max-force')
1063      do i=1,nstep
1064        write(unit=lintmp,fmt='(4x,f10.6)') xmxfor(i)
1065        call write_files(iout,0,lintmp)
1066      enddo
1067      call write_files(iout,0,'rms-force')
1068      do i=1,nstep
1069        write(unit=lintmp,fmt='(4x,f10.6)') rmsfor(i)
1070        call write_files(iout,0,lintmp)
1071      enddo
1072
1073C Note that the program prints out only the Cartesian geometries
1074C (no matter what coordinates were used in the optimization)
1075
1076      call write_files(iout,0,'[GEOMETRIES] XYZ')
1077      write(unit=temp8,fmt='(i5)') nat
1078      do i=1,nstep
1079        write(unit=lintmp,fmt='(a9,f13.6)') 'scf done:', estep(i)
1080        call write_files(iout,0,temp8)
1081        call write_files(iout,0,lintmp)
1082        do j=1,nat
1083          write(unit=lintmp,fmt='(1x,a8,3(f11.6,2x))') atmnam(j),
1084     1          geom(1,j,i),geom(2,j,i), geom(3,j,i)
1085          call write_files(iout,0,lintmp)
1086        enddo
1087      enddo
1088      write(6,*) 'Finished writing geometry optimization information!'
1089      return
1090      end
1091C-----------------------------------------------------------------------
1092      subroutine find_freq(inp,iout,itmp)
1093
1094C This subroutine finds and writes the frequencies and normal modes of
1095C vibration.
1096
1097      implicit double precision(a-e,g-h,o-z)
1098      implicit logical(f)
1099      common /flags/ flgrhf,flguhf,flgopn,flgdft,flggvb,flgmp2,flgmp4
1100      common /flags2/ flgmc
1101      common /rundat/ irun,iwfn,nbf,ndoc,nsoc,maxit,nat,norb
1102      character line*130,lintmp*130
1103      parameter(maxat=100)
1104      dimension efreq(3*maxat),vibcor(3*maxat,3*maxat)
1105
1106      write(6,*) 'Reading vibrational frequencies and displacements...'
1107      call write_files(iout,0,'[FREQ]')
1108      rewind(inp)
1109      call search(inp,line,'Transformation matrix from',ifound,-1)
1110      j=1
1111 10   k=j+5
1112      call skip(inp,3)
1113      call read_line(inp,line)
1114      read(unit=line,fmt='(17x,6f10.3)') (efreq(i), i=j,k)
1115      write(*,fmt='(6(i3,f9.3))') (i,efreq(i), i=j,k)
1116      call skip(inp,2)
1117      do l=1,3*nat
1118        call read_line(inp,line)
1119        read(unit=line,fmt='(17x,6f10.5)') (vibcor(i,l),i=j,k)
1120*       write(*,fmt='(6f10.5)') (vibcor(i,l),i=j,k)
1121      enddo
1122      if (k .ge. 3*nat) then
1123        do i=1,3*nat
1124          write(unit=lintmp,fmt='(f10.4)') efreq(i)
1125          call write_files(iout,0,lintmp)
1126        enddo
1127        call write_files(iout,0,'[FR-COORD]')
1128        rewind(itmp)
1129        call skip(itmp,2)
1130        do i=1,nat
1131          call read_line(itmp,line)
1132          read(unit=line,fmt='(1x,a8,6x,3(2x,f12.7))') atmnam,x,y,z
1133          write(unit=lintmp,fmt='(a8,3(1x,f12.6))') atmnam,x,y,z
1134          call write_files(iout,0,lintmp)
1135        enddo
1136        call write_files(iout,0,'[FR-NORM-COORD]')
1137        do i=1,3*nat
1138          write(unit=lintmp,fmt='(a9,i5)') 'vibration',i
1139          call write_files(iout,0,lintmp)
1140          do j=1,3*nat,3
1141            write(unit=lintmp,fmt='(3(f12.6,1x))') vibcor(i,j),
1142     1            vibcor(i,j+1),vibcor(i,j+2)
1143            call write_files(iout,0,lintmp)
1144          enddo
1145        enddo
1146      goto 20
1147      endif
1148      j=j+6
1149      goto 10
1150 20   write(6,*)           'Finished writing vibrational frequencies and
1151     1 displacements!'
1152      return
1153      end
1154C----------------------------------------------------------------------
1155      subroutine close_files(inp,iout,itmp)
1156
1157C Closes all open files (temp file is deleted)
1158
1159      write(6,*) 'Closing files...'
1160      close(unit=inp, err=10, status='KEEP')
1161      close(unit=iout, err=20, status='KEEP')
1162      close(unit=itmp, err=30, status='DELETE')
1163      goto 40
1164 10   write(6,*) 'An error occurred while closing the input file!'
1165      call error
1166 20   write(6,*) 'An error occurred while closing the output file!'
1167      call error
1168 30   write(6,*) 'An error occurred while closing the temp file!'
1169      call error
1170 40   write(6,*) 'All done! Ending program...'
1171      return
1172      end
1173C-----------------------------------------------------------------------
1174      subroutine blank_line(iout)
1175
1176C Writes a blank line to unit=iout.
1177
1178      write(iout,fmt='(a)') ' '
1179      return
1180      end
1181C----------------------------------------------------------------------
1182      subroutine error
1183      write(6,*) 'An error has occurred while processing the file.'
1184      write(6,*) 'Aborting run...'
1185      stop
1186      end
1187C-----------------------------------------------------------------------
1188      subroutine ffwrd(iunit)
1189
1190C Fast-forward unit=iunit to the end of the file
1191
1192 10   read(iunit,fmt='(a130)',end=20) line
1193      goto 10
1194 20   return
1195      end
1196C----------------------------------------------------------------------
1197      subroutine last_char(line,length,loc)
1198      character*(*) line
1199
1200C  Given string LINE of length LENGTH locate and return in LOC
1201C  the position of the last character.
1202
1203      loc=length
1204      do i=length,1,-1
1205        if(line(i:i).ne.' ') then
1206          loc=i
1207          return
1208        endif
1209      enddo
1210      return
1211      end
1212C-----------------------------------------------------------------------
1213      subroutine read_line(iunit,line)
1214      character line*130
1215      read(iunit,fmt='(a130)') line
1216      return
1217      end
1218C-----------------------------------------------------------------------
1219      subroutine search(iunit,line,text,ifound,idbug)
1220      character*(*) text
1221      character line*130
1222
1223C Scan unit IUNIT for the next occurence of TEXT.
1224C If idbug > 0 write line to unit idbug (probably 6).
1225C If idbug = -1 then stop program if TEXT not found
1226
1227      ifound=0
1228      l=len(text)
1229 10   read(iunit,fmt='(a130)',end=20) line
1230      if (idbug .gt. 0) then
1231        call write_files(idbug,0,line)
1232      endif
1233      do i=1,130-l+1
1234        if (line(i:i+l-1) .eq. text) then
1235          ifound=1
1236          return
1237        endif
1238      enddo
1239      goto 10
1240 20   if (idbug .eq. 1) then
1241        write(6,*) 'Cannot find ', text(1:l+1)
1242        call error
1243      endif
1244      return
1245      end
1246C-----------------------------------------------------------------------
1247      subroutine skip(iunit,nlines)
1248      character char*1
1249
1250C Skip next NLINES on unit IUNIT
1251
1252      do i=1,nlines
1253        read(iunit,fmt='(1a1)') char
1254      enddo
1255      return
1256      end
1257C-----------------------------------------------------------------------
1258      subroutine write_files(iout1,iout2,text)
1259      character*(*) text
1260
1261C Writes TEXT to files iout1 and iout2 (usually output and temp files)
1262C or just output file (iout2=0)
1263
1264      loc1=len(text)
1265      call last_char(text,loc1,loc)
1266      write(iout1,'(a)') text(1:loc)
1267      if (iout2 .eq. 0) goto 10
1268      write(iout2,'(a)') text(1:loc)
1269
1270 10   return
1271      end
1272C----------------------------------------------------------------------
1273
1274      subroutine up_case(string,lenstr)
1275c convert all lower case characters to upper case
1276C
1277      CHARACTER*(*) STRING
1278      CHARACTER*26 UCASE,LCASE
1279C
1280      DATA UCASE /'ABCDEFGHIJKLMNOPQRSTUVWXYZ'/
1281      DATA LCASE /'abcdefghijklmnopqrstuvwxyz'/
1282C
1283C        CONVERTS LOWER CASE TO UPPER CASE IN THE GIVEN STRING.
1284C
1285      DO 100 I=1,LENSTR
1286         IC = INDEX(LCASE,STRING(I:I))
1287         IF (IC.GT.0) STRING(I:I) = UCASE(IC:IC)
1288  100 CONTINUE
1289      RETURN
1290      END
1291C-----------------------------------------------------------------------
1292      subroutine backs(iunit,line,text,ifound,idbug)
1293      character*(*) text
1294      character line*130
1295
1296C Scan unit IUNIT backwards for the occurence of TEXT.
1297C If idbug > 0 write line to unit idbug (probably 6).
1298C If idbug = -1 then stop program if TEXT not found
1299
1300      ifound=0
1301      l=len(text)
1302 10   read(iunit,fmt='(a130)',end=20) line
1303      if (idbug .gt. 0) then
1304        call write_files(idbug,0,line)
1305      endif
1306      do i=1,130-l+1
1307        if (line(i:i+l-1) .eq. text) then
1308          ifound=1
1309          return
1310        endif
1311      enddo
1312      backspace(iunit)
1313      backspace(iunit)
1314      goto 10
1315 20   if (idbug .eq. 1) then
1316        write(6,*) 'Cannot find ', text(1:l+1)
1317        call error
1318      endif
1319      return
1320      end
1321