      subroutine read_wfn
      use orbitals_module
      use confs_module
      use grid_module
      implicit none
      double precision charge,pmatrix,occtest
      integer ix(20),iy(20),iz(20)
      integer i,mu, nu, munu, j,munu_tot,nofirst
      integer deb,fin,atm_nb,deci,uni,nelec,nvac
      character*40 filebas, fileelf,filein,filegam
      character*1 molecule(40)
      character*4 string
      logical rho_and_elf,if_file
      common /files/ molecule,filein,filegam,fileelf,filebas,rho_and_elf
      data ix/0,1,0,0,2,0,0,1,1,0,3,0,0,2,2,0,1,1,0,1/
      data iy/0,0,1,0,0,2,0,1,0,1,0,3,0,1,0,2,2,0,1,1/
      data iz/0,0,0,1,0,0,2,0,1,1,0,0,3,0,1,1,0,2,2,1/
      open_shell=.false.
      uhf=.false.
      ifno=.false.
      ifci=.false.
      fukui=.false.
      write(*,*) ' input .wfn file '
      read(*,'(a40)') filein
      inquire(file=filein,exist=if_file)
      if(.not.if_file) then
        write(*,'(2x,"requested file does not exist:",a40)') filein
        stop
      endif
      open(unit=1,file=filein,status='old')
      read(1,'(40a1)') molecule
      read(1,'(20x,i3,16x,i4,17x,i3)') nmol,ngto,natoms
      write(*,'(i4,'' molecular orbitals'')') nmol
      write(*,'(i4,'' primitive functions'')') ngto
      write(*,'(i4,'' atomic centres'')') natoms
      allocate(p(ngto*(ngto+1)/2), alfa(ngto), rc(3,ngto))
      allocate(xat(natoms),yat(natoms),zat(natoms),znuc(natoms),
     .                     nat(natoms))
      allocate(spin(nmol),occ(2,nmol), c(nmol,ngto))
      allocate(itype(ngto),icent(ngto),ir(3,ngto),atom_name(natoms+1))
      do i=1,natoms
         read(1,'(a4,20x,3f12.8,10x,f6.2)') atom_name(i),xat(i),yat(i),
     .                                      zat(i),znuc(i)
         nat(i)=int(znuc(i))
         call sizstr(atom_name(i),deb,fin)
         atom_name(i)=atom_name(i)(deb:fin)
      enddo
!     form atom labels
      do i=1,natoms-1
        atm_nb=1
        do j=i+1,natoms
          if(atom_name(j).eq.atom_name(i)) atm_nb=atm_nb+1
        enddo
        if(atm_nb.gt.1) then
          atm_nb=1
          call sizstr(atom_name(i),deb,fin)
          string=atom_name(i)
          do j=1,natoms
            if(atom_name(j).eq.string) then
              if(atm_nb.ge.10) then
                 deci=atm_nb/10
                 uni=atm_nb-deci*10
                 atom_name(j)=atom_name(j)(deb:fin)//char(ichar('0')+
     .                        deci)//char(ichar('0')+uni)
               else
                 atom_name(j)=atom_name(j)(deb:fin)//char(ichar('0')+
     .                        atm_nb)
               endif
               atm_nb=atm_nb+1
            endif
          enddo
        endif
      enddo
      read(1,'(20x,20i3)')(icent(i),i=1,ngto)
      read(1,'(20x,20i3)')(itype(i),i=1,ngto)
      read(1,'(10x,5d14.8)')(alfa(i),i=1,ngto)
      do i=1,ngto
         rc(1,i)=xat(icent(i))
         rc(2,i)=yat(icent(i))
         rc(3,i)=zat(icent(i))
         ir(1,i)=ix(itype(i))
         ir(2,i)=iy(itype(i))
         ir(3,i)=iz(itype(i))
      enddo
      uhf=.false.
      nvac=0
      n_alpha=0
      n_beta=0
      nofirst=0
      occtest=1.0d6
      do i=1, nmol
         read(1,'(2x,i5,7x,f12.8,10x,f12.8)') j,spin(i),occ(1,i)
         if(occ(1,i).gt.occtest) then
            open_shell=.true.
            ifno=.true.
            uhf=.true.
            nofirst=i
         endif
         occtest=occ(1,i)
         if(i.eq.1) then
            charge=occ(1,i)
         endif
         if(occ(1,i).ne.0.0d0) then
            if(occ(1,i).eq.2.0d0) then
               n_alpha=n_alpha+1
               n_beta=n_beta+1
               spin(i)=0.0d0
               occ(1,i)=1.0d0
               occ(2,i)=1.0d0
            else if (occ(1,i).eq.1.0d0) then
               if(j.eq.0) then
                  n_alpha=n_alpha+1
                  n_beta=n_beta+1
                  spin(i)=0.0d0
                  occ(1,i)=0.5d0
                  occ(2,i)=0.5d0
               else if(i.eq.j) then
                  if(nvac.eq.0) then
                     n_alpha=n_alpha+1
                     spin(i)=1.0d0
                     occ(1,i)=1.0d0
                     occ(2,i)=0.0d0
                  else
                     n_beta=n_beta+1
                     spin(i)=-1.0d0
                     occ(1,i)=0.0d0
                     occ(2,i)=1.0d0
                  endif
               else
                  n_beta=n_beta+1
                  spin(i)=-1.0d0
                  occ(1,i)=0.0d0
                  occ(2,i)=1.0d0
               endif
            else if(spin(i).lt.0.0d0) then
               n_beta=n_beta+1
               occ(2,i)=occ(1,i)
            else if(spin(i).eq.0.0d0) then
               n_alpha=n_alpha+1
               n_beta=n_beta+1
               spin(i)=0.0d0
               occ(1,i)=0.5d0*occ(1,i)
               occ(2,i)=occ(1,i)
            else
               n_alpha=n_alpha+1
            endif
         else
            nvac=nvac+1
            occ(1,i)=0.0d0
            occ(2,i)=0.0d0
         endif
         read(1,'(5d16.8)')(c(i,j),j=1,ngto)
      enddo
      if(ifno) then
         n_alpha=nofirst-1
         n_beta=nmol-n_alpha
         do i=1,nofirst-1
            occ(1,i)=2.0d0*occ(1,i)
            if(occ(1,i).eq.2.0d0) occ(1,i)=1.0d0
            occ(2,i)=0.0d0
            spin(i)=1.0d0
         enddo
         do i=nofirst,nmol
            occ(2,i)=2.0d0*occ(1,i)
            if(occ(2,i).eq.2.0d0) occ(2,i)=1.0d0
            occ(1,i)=0.0d0
            spin(i)=-1.0d0
         enddo
      endif
      if(charge.eq.1.0d0) then
         uhf=.true.
         open_shell=.true.
      else if(n_alpha.ne.n_beta) then
         open_shell=.true.
      endif
      charge=0.0d0
      do i=1,nmol
         charge=charge+dabs(occ(1,i)+occ(2,i))
      enddo
      nelec=nint(charge)
      if(nelec.ne.(n_alpha+n_beta)) then
         ifno=.true.
         ifci=.true.
      endif
      if(nmol.gt.(n_alpha+n_beta)) then
         ifci=.true.
      endif
      if(open_shell) then
         if(uhf) then
            if(ifno) then
               write(*,*) ' Natural orbitals open shell'
            else if(nvac.eq.0) then
               write(*,*) ' UHF open shell SCF wave function '
            else
               write(*,*) ' UHF open shell CI wave function '
            endif
         else
            if(ifno) then
               write(*,*) ' Natural orbitals open shell'
            else if (nvac.eq.0) then
               write(*,*) ' ROHF open shell SCF wave function '
            else
               write(*,*) ' ROHF open shell CI wave function '
            endif
         endif
      else
         if(ifno) then
            write(*,*) ' Natural orbitals closed shell'
         else if (nvac.eq.0) then
            write(*,*) ' Closed shell SCF wave function '
         else if (nvac.eq.1) then
            write(*,*) ' Closed shell CI or SCF wave function '
            fukui=.true.
         else
            write(*,*) ' Closed shell CI wave function '
         endif   
      endif
!**   update n_alpha and n_beta (CI wavefunction)
      if(ifno) then
         ifci=.false.
      endif
      if(ifci) then
         if(uhf) then
            i=n_alpha+1
            do while (occ(2,i).eq.0.0d0) 
               i=i+1
            enddo
            n_alpha=i-1
            n_beta=nmol-n_alpha
         else
            n_alpha=nmol
            n_beta=n_alpha
         endif
       endif
       write(*,*) ' alpha MOs ',n_alpha,' beta MOs ',n_beta
      munu=0
      do mu=1,ngto
         do nu=mu,ngto
            pmatrix=0.0d0
            do i=1,nmol
               pmatrix=pmatrix+(occ(1,i)+occ(2,i))*c(i,mu)*c(i,nu)
            enddo
            if(mu.ne.nu) then
               pmatrix=2.0d0*pmatrix
            endif
            munu=munu+1
            p(munu)=pmatrix
         enddo
      enddo
      munu_tot=munu
      write(*,'(i10,'' non zero density matrix elements among'',i10)') 
     .           munu_tot, ngto*(ngto+1)/2
      close(unit=1)
      return
      end
