      subroutine read_elf
      use grid_module
      use confs_module
      use orbitals_module
      use basins_module
      implicit none
      integer i,j,ij,kl,npyz,imol,jmol
      integer*2 count
      integer low_limit
      integer npx,npy,npz
      integer, dimension (:), allocatable :: atom_code
      integer new_attr
      real t(2),elapsed,spent,thresh
      real elf_min
      double precision eps,tota
c     double precision origin(3), edge(3),step(3),dv
      double precision dv
      character*1 molecule(40),blank
      character*40 filebas,filein,fileelf,filegam,dumname
      logical func_type,if_file,rho_and_elf
      common /files/molecule,filein,filegam,fileelf,filebas,rho_and_elf
      common/ftype/func_type
      common/valence/thresh,low_limit
      equivalence (molecule(1),dumname)
      data blank/' '/
      allocate(atom_code(natoms+1))
      count=0
      new_attr=0
      eps=1.0d-5
      nmola=nmol
      if(uhf) then
         nmola=n_alpha
      endif
      if(ifci) then
         nmola=nmol
      endif
      deallocate(indij)
      allocate(indij(nmola*(nmola+1)/2))
      do i=1,nmola
         indij(i)=(i-1)*(nmola+nmola-i)/2
      enddo
*
*     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)='_'
      if(func_type) then
         molecule(i+1)='e'
         molecule(i+2)='l'
         molecule(i+3)='f'
      else
         molecule(i+1)='r'
         molecule(i+2)='h'
         molecule(i+3)='o'
      endif
      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
      molecule(i)='_'
      if(func_type) then
         molecule(i+1)='e'
      else
         molecule(i+1)='r'
      endif
      molecule(i+2)='b'
      molecule(i+3)='a'
      molecule(i+4)='s'
      molecule(i+5)='.'
      molecule(i+6)='s'
      molecule(i+7)='b'
      molecule(i+8)='f'
      filebas=dumname
      do i=1,natoms
         atom_code(i)=0
      enddo
      atom_name(natoms+1)='    '
      inquire(file=fileelf,exist=if_file)
         if(if_file) then
         call read_sbf
         origin(1)=xmin
         origin(2)=ymin
         origin(3)=zmin
         edge(1)=xmax-origin(1)
         edge(2)=ymax-origin(2)
         edge(3)=zmax-origin(3)
         npx=np(1)-1
         npy=np(2)-1
         npz=np(3)-1
         if(npx.eq.0) then
            step(1)=0.0d0
         else
            step(1)=edge(1)/npx
         endif
         if(npy.eq.0) then
            step(2)=0.0d0
         else
            step(2)=edge(2)/npy
         endif
         if(npz.eq.0) then
            step(3)=0.0d0
         else
            step(3)=edge(3)/npz
         endif
         dv=step(1)*step(2)*step(3)
         allocate(fx(np(1),ngto),dfx(np(1),ngto),d2fx(np(1),ngto))
         allocate(fy(np(2),ngto),dfy(np(2),ngto),d2fy(np(2),ngto))
         allocate(fz(np(3),ngto),dfz(np(3),ngto),d2fz(np(3),ngto))
         call generf2
         write(*,'(/)')
         write(*,'(2x,''origin      '',3f12.4)') origin
         write(*,'(2x,''edge        '',3f12.4)') edge  
         write(*,'(2x,''step        '',3f12.4)') step  
         write(*,'(2x,''grid points '',3i12)') npx,npy,npz
         npts=np(1)*np(2)*np(3)
         write(*,'(/,2x,''total number of grid points: '',i9,/)') npts
         call xtime(t)
         spent=t(1)+t(2)
*
*     sort elf array -----> indices in indx
*
         if(func_type) then
            thresh=1.0e-5
         else
            thresh=1.0e-4
         endif
         call elfsort(thresh)
         elf_min=0.0
         i=npts+1
         do while(elf_min.lt.thresh)
            i=i-1
            j=indx(i)
            elf_min=elf(j)
         enddo
         low_limit=i
         npyz=np(2)*np(3)
*
*     assign attractor grid points
*
         do i=1,npts
            attract_code(i)=0
         enddo
         call xtime(t)
         elapsed=t(1)+t(2)-spent
         spent=t(1)+t(2)
         write(*,'(/,2x,''elapsed time in sort procedure'',f9.2,
     .               '' sec.'')') elapsed
      else
         write(*,'(2x,"requested file does not exist:",a40)') fileelf
         stop
      endif
*     input 1st density matrix coefficients
*     read(unit) na,ia,pa
      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')
      read(14) na
      allocate(ia(na),pa(na))
      read(14)(pa(i),i=1,na)
      read(14)(ia(i),i=1,na)
      read(14) nb
      nmola=nmol
      allocate(ib(nb),pb(nb))
      if(nb.ne.0) then
         nmola=nmol/2
         read(14)(pb(i),i=1,nb)
         read(14)(ib(i),i=1,nb)
      endif
*     input 2nd density matrix coefficients
      read(14) naa
      allocate(iaa(naa),jaa(naa),kaa(naa),laa(naa),paa(naa))
      read(14)(paa(i),i=1,naa)
      read(14)(iaa(i),i=1,naa)
      read(14)(jaa(i),i=1,naa)
      read(14)(kaa(i),i=1,naa)
      read(14)(laa(i),i=1,naa)
      do i=1,naa
         ij=indij(iaa(i))+jaa(i)
         kl=indij(kaa(i))+laa(i)
         if(ij.ne.kl) then
            paa(i)=paa(i)+paa(i)
         endif
      enddo
      read(14) nbb
      allocate(ibb(nbb),jbb(nbb),kbb(nbb),lbb(nbb),pbb(nbb))
      if(nbb.ne.0) then
         read(14)(pbb(i),i=1,nbb)
         read(14)(ibb(i),i=1,nbb)
         read(14)(jbb(i),i=1,nbb)
         read(14)(kbb(i),i=1,nbb)
         read(14)(lbb(i),i=1,nbb)
         do i=1,nbb
            ij=indij(ibb(i))+jbb(i)
            kl=indij(kbb(i))+lbb(i)
            if(ij.ne.kl) then
               pbb(i)=pbb(i)+pbb(i)
            endif
         enddo
      endif
      close(unit=14)
      else if(open_shell) then
         allocate(iaa(nmol*(nmol+1)/2),jaa(nmol*(nmol+1)/2))
         allocate(kaa(nmol*(nmol+1)/2),laa(nmol*(nmol+1)/2))
         allocate(ibb(nmol*(nmol+1)/2),jbb(nmol*(nmol+1)/2))
         allocate(kbb(nmol*(nmol+1)/2),lbb(nmol*(nmol+1)/2))
         allocate(paa(nmol*(nmol+1)/2),pbb(nmol*(nmol+1)/2))
         if(uhf) then
            na=0
            nb=0
            naa=0
            nbb=0
            nmola=n_alpha
            allocate(ia(n_alpha),pa(n_alpha),ib(n_beta),pb(n_beta))
            do imol=1,n_alpha
               if(occ(1,imol).gt.0.0d0) then
                  na=na+1
                  pa(na)=1.0d0
                  if(ifno) then
                     pa(na)=occ(1,imol)
                  endif
                  ia(na)=indij(imol)+imol
                  do jmol=imol,nmol
                     if(occ(1,jmol).gt.0.0d0) then
                        naa=naa+1
                        paa(naa)=1.0d0
                        iaa(naa)=imol
                        jaa(naa)=jmol
                        kaa(naa)=imol
                        laa(naa)=jmol
                        if(ifno) then
                           tota=occ(1,jmol)*occ(1,imol)
                           paa(naa)=tota
                        endif
                     endif
                  enddo
               endif
            enddo
            do imol=n_alpha+1,n_alpha+n_beta
               if(occ(2,imol).gt.0.0d0) then
                  nb=nb+1
                  pb(nb)=1.0d0
                  if(ifno) then
                     pb(nb)=occ(2,imol)
                  endif
                  ib(nb)=indij(imol-n_alpha)+imol-n_alpha
                  do jmol=imol,n_alpha+n_beta
                     if(occ(2,jmol).gt.0.0d0) then
                        nbb=nbb+1
                        pbb(nbb)=1.0d0
                        ibb(nbb)=imol-n_alpha
                        jbb(nbb)=jmol-n_alpha
                        kbb(nbb)=imol-n_alpha
                        lbb(nbb)=jmol-n_alpha
                        if(ifno) then
                           tota=occ(2,imol)*occ(2,jmol)
                           pbb(nbb)=tota
                        endif
                     endif
                  enddo
               endif
            enddo
         else
            na=0
            nb=0
            naa=0
            nbb=0
            nmola=nmol
            allocate(ia(nmol),pa(nmol),ib(nmol),pb(nmol))
            do imol=1,nmol
               if(occ(1,imol).gt.0.0d0) then
                  na=na+1
                  pa(na)=1.0d0
                  ia(na)=indij(imol)+imol
                  do jmol=imol,nmol
                     if(occ(1,jmol).gt.0.0d0) then
                        naa=naa+1
                        paa(naa)=1.0d0
                        iaa(naa)=imol
                        jaa(naa)=jmol
                        kaa(naa)=imol
                        laa(naa)=jmol
                     endif
                  enddo
               endif
            enddo
            do imol=1,nmol
               if(occ(2,imol).gt.0.0d0) then
                  nb=nb+1
                  pb(nb)=1.0d0
                  ib(nb)=indij(imol)+imol
                  do jmol=imol,nmol
                     if(occ(2,jmol).gt.0.0d0) then
                        nbb=nbb+1
                        pbb(nbb)=1.0d0
                        ibb(nbb)=imol
                        jbb(nbb)=jmol
                        kbb(nbb)=imol
                        lbb(nbb)=jmol
                     endif
                  enddo
               endif
            enddo
         endif
      else
        if(ifno) then
        allocate(ia(nmol),pa(nmol),iaa(nmol*(nmol+1)))
        allocate(jaa(nmol*(nmol+1)),kaa(nmol*(nmol+1)))
        allocate(laa(nmol*(nmol+1)),paa(nmol*(nmol+1)))
        allocate(ib(nmol),pb(nmol),ibb(nmol*(nmol+1)))
        allocate(jbb(nmol*(nmol+1)),kbb(nmol*(nmol+1)))
        allocate(lbb(nmol*(nmol+1)),pbb(nmol*(nmol+1)))
        else
        allocate(ia(nmol),pa(nmol),iaa(nmol*(nmol+1)/2))
        allocate(jaa(nmol*(nmol+1)/2),kaa(nmol*(nmol+1)/2))
        allocate(laa(nmol*(nmol+1)/2),paa(nmol*(nmol+1)/2))
        allocate(ib(nmol),pb(nmol),ibb(nmol*(nmol+1)/2))
        allocate(jbb(nmol*(nmol+1)/2),kbb(nmol*(nmol+1)/2))
        allocate(lbb(nmol*(nmol+1)/2),pbb(nmol*(nmol+1)/2))
        endif
        na=0
        nb=0
        naa=0
        nbb=0
        nmola=nmol
        do imol=1,nmol
           if(occ(1,imol).gt.0.0d0) then
              na=na+1
              pa(na)=1.0d0
              if(ifno) then
                 pa(na)=occ(1,imol)
              endif
              ia(na)=indij(imol)+imol
              do jmol=imol,nmol
                 if(occ(1,jmol).gt.0.0d0) then
                    naa=naa+1
                    paa(naa)=1.0d0
                    iaa(naa)=imol
                    jaa(naa)=jmol
                    kaa(naa)=imol
                    laa(naa)=jmol
                    if(ifno) then
                       tota=occ(1,jmol)*occ(1,imol)
                       paa(naa)=tota
                    endif
                 endif
              enddo
           endif
        enddo
      endif
      deallocate(atom_code)
      return
      end
