!
!     TOP_GRID: purpose genrate ELF on a 3-D grid parallel to the standard
!             axis defined in the MO calculation
!     
!     The wavefunction is read on a file generated by Gaussian92/94 
!     with the option out=wfn
!     here is an example of a Gaussain94 input
!
!----------------------------------------------------------------------------
!     # HF/6-31G out=wfn
!     
!     water molecule
!     
!     0  1
!     O  
!     H  1  0.96
!     H  1  0.96  2  104.5
!
!     h2o.wfn
!----------------------------------------------------------------------------
!
!     The input can be made either from the keybord or from an input file
!     Input description
!     1. wfn file name (a40) 
!     3. output of density grid (a1) : y output density /n otherwise
!     4. coordinates of the origin (free format) origin of the box defining the
!        grid
!     5. length of the edges of the box (free format}
!     6. number of intervals on each edge (free format}
!        ELF is computed by the formula
!        ELF = 1/(1*((D+e)/D0)**2)
!        In Becke and Edgecombe's work e=0.0
!        A recommanded value for the damping factor is 0.001
!     input example
!
!-----------------------------------------------------------------------------
!     h2o.wfn
!     y
!     -3.0 -3.0 -3.0
!      6.0  6.0  6.0
!      60   60   60
!-----------------------------------------------------------------------------
!
!     in output the program generates title_elf.sbf (and title_rho.btf if 
!     requested files) file(s) in which title is the first non-blank characters
!     of the title in gaussian input (here water)
!    
!     This files  are binary filies to be processed with TOP_BAS
!     and converted to stf to be visualized with the SciAn software 
!     SciAn is a scientific visualization and animation program for high   
!     performance graphic workstations developped by Eric Pepke, John Murray
!     Jim Lyons and Tzong-Yow Hwu of the Suprcomputer Computations Research
!     Institute of the Florida State University at Tallahassee.
!     SciAn is a free program you can get by anonymous ftp at
!     ftp.scri.fsu.edu 
!     login: anonymous
!     password: your internet address
!     the SciAn program and documentation are located in the subdirectories
!     of pub/SciAn
!     in pub/SciAn get the README file
!
!     References
!     A. Becke and K. E. Edgecombe, J. Chem. Phys., 92, 5397-5404 (1990)
!     B. Silvi and A. Savin, Nature, 371, 683-686 (1994)
!
!     please report bugs or problems via electronic mail at:
!     silvi@lct.jussieu.fr
      program top_grid
      use orbitals_module
      use grid_module
      use confs_module
      implicit none
      double precision thresh
      real t(2),elapsed
      character*1 molecule(1:40),blank
      character*40 filerho, fileelf, filelap, filegam, dumname
      character*40 filebas, filein
      integer intervalx, intervaly, intervalz,ibuf
      integer i,j,npmax
      logical if_file,rho_and_elf
      equivalence (molecule(1),dumname)
      common /files/ molecule,filein,filegam,fileelf,filebas,rho_and_elf
      data blank/' '/
!     input wavefunction
      call read_wfn
!     input grid specification (origin, edges, number of intervals)
!     create output file names
!
      i=1
      do while (molecule(i).eq.blank)
         do j=i,39
            molecule(j)=molecule(j+1)
         enddo
         i=i+1
      enddo
      do while (molecule(i).ne.blank)
         i=i+1
      enddo
      molecule(i)='_'
      molecule(i+1)='e'
      molecule(i+2)='l'
      molecule(i+3)='f'
      molecule(i+4)='.'
      molecule(i+5)='s'
      molecule(i+6)='b'
      molecule(i+7)='f'
      do j=i+8,40
         molecule(j)=blank
      enddo
      fileelf=dumname
      molecule(i+1)='g'
      molecule(i+2)='a'
      molecule(i+3)='m'
      filegam=dumname
      open(unit=2,file=fileelf,status='unknown',form='unformatted')
      molecule(i)='_'
      molecule(i+1)='r'
      molecule(i+2)='h'
      molecule(i+3)='o'
      molecule(i+4)='.'
      molecule(i+5)='s'
      molecule(i+6)='b'
      molecule(i+7)='f'
      filerho=dumname
      open(unit=3,file=filerho,status='unknown',form='unformatted')
      molecule(i)='_'
      molecule(i+1)='l'
      molecule(i+2)='a'
      molecule(i+3)='p'
      molecule(i+4)='.'
      molecule(i+5)='s'
      molecule(i+6)='b'
      molecule(i+7)='f'
      filelap=dumname
      open(unit=4,file=filelap,status='unknown',form='unformatted')
      nmola=nmol
      if(uhf) then
         nmola=n_alpha
      endif
      if(ifci) then
         nmola=nmol
      endif
      allocate(indij(nmola))
      do i=1,nmola
         indij(i)=(i-1)*(nmola+nmola-i)/2
      enddo
      inquire(file=filegam,exist=if_file)
      if(if_file) then
         write(*,'(2x,"correlated wave function coeff ",a40)') filegam
         open(unit=14,file=filegam,status='unknown',form='unformatted')
         call read_det_coeffs
         close(unit=14)
      else
         call build_det_coeffs
      endif
!
!     find indicative limits
!     
      xmin=xat(1)
      ymin=yat(1)
      zmin=zat(1)
      xmax=xat(1)
      ymax=yat(1)
      zmax=zat(1)
      if(natoms.gt.1) then
         do i=1,natoms
            if(xat(i).lt.xmin) then
               xmin=xat(i)
            endif
            if(xat(i).gt.xmax) then
               xmax=xat(i)
            endif
            if(yat(i).lt.ymin) then
               ymin=yat(i)
            endif
            if(yat(i).gt.ymax) then
               ymax=yat(i)
            endif
            if(zat(i).lt.zmin) then
               zmin=zat(i)
            endif
            if(zat(i).gt.zmax) then
               zmax=zat(i)
            endif
         enddo
      endif
      write(*,*) ' indicative origin of the grid '
      write(*,'(4x,3f12.6)') xmin-5.0d0,ymin-5.0d0,zmin-5.0d0
      write(*,*) ' indicative edges along x,y,z '
      write(*,'(4x,3f12.6)') xmax+10.0d0-xmin,ymax+10.0d0-ymin,
     .                       zmax+10.0d0-zmin
*      write(*,*) ' input the origin of the grid x,y,z '
*      read(*,*) origin(1), origin(2), origin(3)
		origin(1) =  xmin-5.0d0
		origin(2) =  ymin-5.0d0
		origin(3) = zmin-5.0d0
*      write(*,*) ' input edges along x,y,z '
*      read(*,*) edge(1), edge(2), edge(3)
		edge(1) = xmax+10.0d0-xmin
		edge(2) = ymax+10.0d0-ymin
		edge(3) = zmax+10.0d0-zmin
      write(*,*) ' input the number of intervals along x,y,z '
      read(*,*) intervalx, intervaly, intervalz
      np(1)=intervalx+1
      np(2)=intervaly+1
      np(3)=intervalz+1
      write(*,'(i4,'' grid points along x axis'')') np(1)
      write(*,'(i4,'' grid points along y axis'')') np(2)
      write(*,'(i4,'' grid points along z axis'')') np(3)
      npmax=max(np(1),np(2),np(3))
      write(2) np(1),np(2),np(3)
      write(2) origin(1),origin(1)+edge(1),origin(2),origin(2)+
     .         edge(2),origin(3),origin(3)+edge(3) 
      write(3) np(1),np(2),np(3)
      write(3) origin(1),origin(1)+edge(1),origin(2),origin(2)+
     .         edge(2),origin(3),origin(3)+edge(3)
      write(4) np(1),np(2),np(3)
      write(4) origin(1),origin(1)+edge(1),origin(2),origin(2)+
     .         edge(2),origin(3),origin(3)+edge(3)
      allocate(fx(np(1),ngto),dfx(np(1),ngto),dfx2(np(1),ngto))
      allocate(fy(np(2),ngto),dfy(np(2),ngto),dfy2(np(2),ngto))
      allocate(fz(np(3),ngto),dfz(np(3),ngto),dfz2(np(3),ngto))
      allocate(buffer_rho(np(1)*np(2)*np(3)))
      allocate(buffer_elf(np(1)*np(2)*np(3)))
      allocate(buffer_lap(np(1)*np(2)*np(3)))
      do i=1,3
         if(edge(i).ne.0.0d0) then
            step(i)=edge(i)/(np(i)-1)
         else
            step(i)=0.0d0
         endif
      enddo
      call xtime(t)
      elapsed=t(1)+t(2)
      call generf
      thresh=1.0d-5
      call elfcalc(thresh)
      ibuf=np(1)*np(2)*np(3)
      if(ibuf.ne.0) then
         write(2)(buffer_elf(i),i=1,ibuf)
         write(3)(buffer_rho(i),i=1,ibuf)
         write(4)(buffer_lap(i),i=1,ibuf)
      endif
      call xtime(t)
      elapsed=t(1)+t(2)-elapsed
      write(*,'(2x,'' elapsed time '',f9.2, ''sec'',/)') elapsed
      write(*,'(i9,'' elf values have been written on file: '',a40)') 
     .       np(1)*np(2)*np(3),fileelf
      close(unit=2)
      write(*,'(i9,'' values of the density have been written on file: '
     .'        ,a40)') np(1)*np(2)*np(3),filerho
      close(unit=3)
      write(*,'(i9,'' values of the laplacian have been written on file:
     . ''           ,a40)') np(1)*np(2)*np(3),filelap
      close(unit=4)
      stop 'normal termination '
      end
