program mcdvdiv
      implicit none
      integer*4 ln,ig,i,nat0,ib,nat,it
      parameter (nat0=200)
      integer*4 ty(nat0)
      real*8 r(3*nat0),step,x,y,z
      character*180 s80
c
      write(6,600)
600   format(
     1' Devides gaussian output with numerical differentiation',/,
     1' to individual geometries',/,
     2'    GV.OUT - only one calculation (# .... #) possible.',/)
c
      ln=0
      ig=0
      ib=0
      step=0.001d0
      call inif(91,ig)
      open(9,file='GV.OUT',status='old')
11    read(9,900,end=88,err=88)s80
900   format(a180)
      ln=ln+1

      if(s80(27:44).eq.'Input orientation:')then
       write(6,900)s80
       write(6,*)'Line ',ln
       read(9,*)
       read(9,*)
       read(9,*)
       read(9,*)
       nat=0
200    read(9,*,err=87)it,it,x,x,y,z
       nat=nat+1
       if(nat.gt.nat0)then
        write(6,*)'too many atoms',nat,nat0
        stop
       endif
       ty(nat)=it
       r(1+3*(nat-1))=x
       r(2+3*(nat-1))=y
       r(3+3*(nat-1))=z
       goto 200
87     continue
       write(6,*)nat,' atoms'
       if(ig.eq.0)call wrgeo(r,nat,ig,ty,step)
       goto 11
      endif


      if(s80(2:16).eq.'Standard basis:'.and.ig.eq.0.and.ib.eq.0)then
       ib=ib+1
       open(80,file='basis')
       call ws80(80,s80)
6      read(9,900)s80
       if(s80(2:17).ne.'Integral buffers')then
        call ws80(80,s80)
        goto 6
       endif
       close(80)
       if(ig.eq.0)call wrbas
      endif

      call ws80(91,s80)

      if(s80(2:14).eq.'Nuclear step=')read(s80(15:23),*)step

      if(s80(2:18).eq.'Cartesian Forces:')then
       do 1 i=1,6
       read(9,900,end=88,err=88)s80
       if(s80(2:8).eq.'NDeriv=')then
        close(91)
        if(ig.eq.0)close(80)
        ig=ig+1
        call inif(91,ig)
        if(ig.gt.0)then
         call wrgeo(r,nat,ig,ty,step)
         call wrbas
        endif
        goto 11
       endif
1      continue
      endif

      goto 11
88    close(9)
      close(91)
      write(6,*)ln,' lines of GV.OUT'
      write(6,*)ig,' differentiation points'
      end

      subroutine wrbas
      character*180 s80
      open(80,file='basis')
5     read(80,900,end=44,err=44)s80
900   format(a180)
      call ws80(91,s80)
      goto 5
44    close(80)
      return
      end


      subroutine wrgeo(r,nat,ig,ty,step)
      implicit none
      integer*4 nat,ig,ty(*),ia,ar,xr,is
      real*8 r(*),x,y,z,si,step
      character*8 yy(2)
      data yy/'forward ','backward'/
      character*1 xy(3)
      data xy/'x','y','z'/
      write(91,1000)
1000  format(26x,'Input orientation:',/,
     11x,69(1h-),/,
     1' Center     Atomic      Atomic             Coordinates',
     1' (Angstroms)',/,
     1' Number     Number       Type             X           Y',
     1'      Z',/,
     11x,69(1h-))
      ar=(ig+5)/6
      xr=(ig-6*(ar-1)+1)/2
      is=ig-6*(ar-1)-2*(xr-1)
      if(is.eq.1)then
       si= step
      else
       si=-step
      endif
      if(ig.eq.0)then
       write(6,501)
501    format(' Zero geometry')
      else
       write(6,500)ar,xy(xr),yy(is)
500    format('IAt ',i3,1x,a1,2x,a8)
      endif
      do 1 ia=1,nat
      x=r(1+3*(ia-1))
      y=r(2+3*(ia-1))
      z=r(3+3*(ia-1))
      if(ia.eq.ar)then
       if(xr.eq.1)x=x+si
       if(xr.eq.2)y=y+si
       if(xr.eq.3)z=z+si
      endif
1     write(91,1001)ia,ty(ia),0,x,y,z
1001  format(i7,i11,i12,4x,3f12.6)
      write(91,1002)
1002  format(1x,69(1h-))
      return
      end

      subroutine inif(io,i)
      character*10 s10
      write(s10,10)i
10    format(i10)
      do 1 is=1,len(s10)
1     if(s10(is:is).ne.' ')goto 2
2     do 3 ie=len(s10),1,-1
3     if(s10(ie:ie).ne.' ')goto 4
4     call system('mkdir '//s10(is:ie))
      open(io,file=s10(is:ie)//'/GV.OUT')
      return
      end

      subroutine ws80(io,s)
      character*(*) s
      do 3 ie=len(s),1,-1
3     if(s(ie:ie).ne.' ')goto 4
4     write(io,901)(s(i:i),i=1,ie)
901   format(180a1)
      return
      end