      use orbitals_module
      implicit none
      double precision charge,pmatrix,ener
      integer ix(20),iy(20),iz(20)
      integer i,mu, nu, munu, j,munu_tot,newgto, k
      integer deb,fin,atm_nb,deci,uni,nelec,nvac,natoms2
      integer, dimension (:), allocatable :: iatom,icent2,itype2,lc
      double precision, dimension(:),allocatable :: alfa2,c2
      character*40 filein,filegam,fileout,filebas,fileelf
      character*1 molecule(40)
      character*4 string
      logical, dimension (:), allocatable :: if_cent
      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
      write(*,*) ' number of selected atoms '
      read(*,*) natoms2
      allocate(iatom(natoms2))
      read(*,*) (iatom(i),i=1,natoms2)
      write(*,*) ' output .wfn file '
      read(*,'(a40)') fileout
      open(unit=1,file=filein,status='old')
      open(unit=2,file=fileout,status='unknown')
      read(1,'(40a1)') molecule
      write(2,'(40a1)') molecule
      read(1,"('GAUSSIAN',12x,i3,' MOL ORBITALS   ',i4,
     .         'PRIMITIVES',7x,i3,'  NUCLEI')") 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))
      allocate(lc(ngto))
      do i=1,natoms
         read(1,'(a4,20x,3f12.8,10x,f6.2)') atom_name(i),xat(i),yat(i),
     .    zat(i),znuc(i)
      enddo
      read(1,'(20x,20i3)')(icent(i),i=1,ngto)
      newgto=0
      allocate(if_cent(ngto))
      do i=1,ngto
         if_cent(i)=.false.
         do j=1,natoms2
            if(icent(i).eq.iatom(j)) then
               newgto=newgto+1
               lc(i)=j
               if_cent(i)=.true.
            endif
         enddo
      enddo
      write(2,"('GAUSSIAN',12x,i3,' MOL ORBITALS   ',i4,
     .         ' PRIMITIVES',6x,i3,' NUCLEI')") nmol,newgto,natoms2
      do i=1,natoms2
         j=iatom(i)
         write(2,"(a4,i4,4x,'(CENTRE  1)',1x,3f12.8,'  CHARGE =',
     .         f5.1)") atom_name(j),i,xat(j),yat(j),
     .    zat(j),znuc(j)
      enddo
      allocate(icent2(newgto),itype2(newgto),alfa2(newgto))
      allocate(c2(newgto))
      read(1,'(20x,20i3)')(itype(i),i=1,ngto)
      read(1,'(10x,5d14.8)')(alfa(i),i=1,ngto)
      k=0
      do i=1,ngto
         if(if_cent(i)) then
            k=k+1
            icent2(k)=lc(i)
            itype2(k)=itype(i)
            alfa2(k)=alfa(i)
         endif
      enddo
      write(2,"('CENTRE ASSIGNMENTS  ',20i3)")(icent2(i),i=1,newgto)
      write(2,"('TYPE ASSIGNMENTS    ',20i3)")(itype2(i),i=1,newgto)
      write(2,"('EXPONENTS ',5e14.7)")(alfa2(i),i=1,newgto)
      do i=1, nmol
         read(1,'(2x,i5,7x,f12.8,10x,f12.8,14x,f12.7)') 
     .            j,spin(i),occ(1,i),ener
         write(2,"('MO',i5,'     MO',f4.1,8x,'OCC NO =  ',
     .         f11.7,'  ORB. ENERGY =',f12.7)"),j,spin(i),occ(1,i),
     .         ener   
         read(1,'(5d16.8)')(c(i,j),j=1,ngto)
         k=0
         do j=1,ngto
            if(if_cent(j)) then
               k=k+1
               c2(k)=c(i,j)
            endif
         enddo
         write(2,'(5d16.8)')(c2(j),j=1,newgto)
      enddo
      write(2,"('END DATA')")
      read(1,'(17x,f20.12,20x,f12.8)') ener,charge
      write(2,"(' THE  HF ENERGY =  ',f18.12,' THE VIRIAL(-V/T)=  ',
     .          f12.8)") ener,charge
      close(unit=1)
      close(unit=2)
      end
