      use orbitals_module
      implicit none
      double precision charge,pmatrix,enerhf,arg,dsq2
      integer ix(20),iy(20),iz(20)
      integer i,mu,k,nu,munu,j,munu_tot,newgto,irs
      integer deb,fin,atm_nb,deci,uni,nelec,nvac,natoms2
      integer, dimension (:), allocatable :: iatom,ipx,ipy,lc
      integer, dimension (:), allocatable :: idx2,idy2,idxy
      double precision, dimension (:,:),  allocatable :: cnew
      double precision, dimension (:),  allocatable :: ener,enernew
      double precision, dimension (:),  allocatable :: spinnew
      integer i1,i2,i3,i4,i5,k1,k2,k3,k4,k5,jmol,npx,npy,ndx2,ndy2,ndxy
      character*40 filein, format_elf, format_rho, fileout
      character*8 eof
      character*1 molecule(40)
      logical rho_and_elf,if_file
      common /files/ molecule, format_elf, format_rho, 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.
      dsq2=1.0d0/dsqrt(2.0d0)
      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(*,*) ' output wfn file '
      read(*,'(a40)') fileout
      open(unit=1,file=filein,status='old')
      read(1,'(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(2*nmol),occ(2,2*nmol),c(2*nmol,ngto),
     .         spinnew(3*nmol))
      allocate(cnew(3*nmol,ngto))
      allocate(itype(ngto),icent(ngto),ir(3,ngto),atom_name(natoms))
      allocate(lc(ngto),ipx(ngto),ipy(ngto))
      allocate(idx2(ngto),idy2(ngto),idxy(ngto))
      allocate(ener(nmol),enernew(2*nmol))
      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)
      read(1,'(20x,20i3)')(itype(i),i=1,ngto)
      read(1,'(10x,5d14.7)')(alfa(i),i=1,ngto)
      npx=0
      npy=0
      ndx2=0
      ndy2=0
      ndxy=0
      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))
         irs=ir(1,i)+ir(2,i)
         if(irs.eq.1) then
            if(ir(1,i).eq.1) then
               npx=npx+1
               ipx(npx)=i
            endif
            if(ir(2,i).eq.1) then
               npy=npy+1
               ipy(npy)=i
            endif
         else if (irs.eq.2) then
            if(ir(1,i).eq.2) then 
              ndx2=ndx2+1
              idx2(ndx2)=i
            endif
            if(ir(2,i).eq.2) then
              ndy2=ndy2+1
              idy2(ndy2)=i
            endif
            if((ir(1,i).eq.1).and.(ir(2,i).eq.1)) then
              ndxy=ndxy+1
              idxy(ndxy)=i
            endif
          endif
      enddo
      write(*,*) '  px basis functions ',npx,'  py ',npy
      write(*,*) ' dx2 basis functions ',ndx2,' dy2 ',ndy2,' dxy',ndxy
      do i=1, nmol
         read(1,'(2x,i5,7x,f12.8,10x,f12.8,14x,f12.7)')
     .            j,spin(i),occ(1,i),ener(i)
         read(1,'(5d16.8)')(c(i,j),j=1,ngto)
      enddo
      read(1,'(a8)') eof
      read(1,'(17x,f20.12,20x,f12.8)') enerhf,charge
      close(unit=1)
      occ(1,nmol)=0.5d0* occ(1,nmol)
      nmol=nmol+1
      occ(1,nmol)=occ(1,nmol-1)
      ener(nmol)=ener(nmol-1)
      spin(nmol)=spin(nmol-1)
      do j=1,ngto
         c(nmol,j)=c(nmol-1,j)
      enddo
      do j=1,npx
         k1=ipx(j)
         k2=ipy(j)
         c(nmol,k1)=c(nmol-1,k2)
         c(nmol,k2)=c(nmol-1,k1)
      enddo
      open(unit=2,file=fileout,status='unknown')
      write(2,'(40a1)') molecule
      write(2,"('GAUSSIAN',12x,i3,' MOL ORBITALS   ',i4,
     .         ' PRIMITIVES',6x,i3,' NUCLEI')") nmol,ngto,natoms
      do i=1,natoms
         write(2,"(a4,i4,4x,'(CENTRE  1)',1x,3f12.8,'  CHARGE =',
     .         f5.1)") atom_name(i),i,xat(i),yat(i),
     .    zat(i),znuc(i)
      enddo
      write(2,"('CENTRE ASSIGNMENTS  ',20i3)")(icent(i),i=1,ngto)
      write(2,"('TYPE ASSIGNMENTS    ',20i3)")(itype(i),i=1,ngto)
      write(2,"('EXPONENTS ',5e14.7)")(alfa(i),i=1,ngto)
      do i=1, nmol
         write(2,"('MO',i5,'     MO',f4.1,8x,'OCC NO =  ',
     .         f11.7,'  ORB. ENERGY =',f12.7)")i,spin(i),occ(1,i),
     .         ener(i)
         write(2,'(5d16.8)')(c(i,j),j=1,ngto)
      enddo
      write(2,"('END DATA')")
      write(2,"(' THE  HF ENERGY =  ',f18.12,' THE VIRIAL(-V/T)=  ',
     .          f12.8)") enerhf,charge
      close(unit=2)
      end
