      use orbitals_module
      implicit none
      double precision rab,rcov(103),ratoms
      integer i, j
      integer deb,fin,atm_nb,deci,uni,iline
      character*40 filine,dumname,filein,filegam,fileelf,filebas
      character*1 molecule(40), blank
      character*4 string
      logical rho_and_elf,if_file,clean
      equivalence (molecule(1),dumname)
      common /files/ molecule,filein,filegam,fileelf,filebas,rho_and_elf
      data blank/' '/
      data rcov/0.7d0,1.9D1,2.3d0,2.1d0,1.7d0,1.5d0,1.3d0,1.3d0,1.3d0,
     .                 1.9d1,3.2d0,3.0d0,2.8d0,2.5d0,2.3d0,2.1d0,1.9d0,
     .                1.9d1,3.8d0,4.2d0,10*3.8d0,2*2.6d0,2*2.5d0,2.3d0,
     .                68*4.3d0/      
      call read_wfn
*
      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)='l'
      molecule(i+2)='i'
      molecule(i+3)='n'
      molecule(i+4)='e'
      molecule(i+5)='.'
      molecule(i+6)='a'
      molecule(i+7)='m'
      do j=i+8,40
         molecule(j)=blank
      enddo
      filine=dumname
      write(*,'(a40)') filine
      open(unit=2,file=filine,status='unknown')
      write(2,'("#AmiraMesh ASCII 1.0 ")')
      iline=0
      do i=1,natoms-1
         do j=i+1, natoms
            rab=(xat(j)-xat(i))**2+(yat(j)-yat(i))**2
     .          +(zat(j)-zat(i))**2
            rab=dsqrt(rab)
            ratoms=rcov(nat(i))+rcov(nat(j))
            if(rab.le.ratoms) then
               iline=iline+1
            endif
          enddo
       enddo
       write(2,'("define Lines ",i3)') 3*iline
       write(2,'("define Vertices ",i3)') natoms
       write(2,'("Parameters {")')
       write(2,'(''      ContentType "HxLineSet"'')')
       write(2,'("}")')
       write(2,'("Vertices { float[3] Coordinates } = @1")')
       write(2,'("Vertices { float Data} = @2")')
       write(2,'("Lines { int LineIdx } = @3")')
       write(2,'("@1")')
       do i=1,natoms
          write(2,'(3f12.4)') zat(i),yat(i),xat(i)
       enddo
       write(2,'("@2")')
       write(2,'(20i4)') (nat(i),i=1,natoms)
       write(2,'("@3")')
      do i=1,natoms-1
         do j=i+1, natoms
            rab=(xat(j)-xat(i))**2+(yat(j)-yat(i))**2
     .          +(zat(j)-zat(i))**2
            rab=dsqrt(rab)
            ratoms=rcov(nat(i))+rcov(nat(j))
            if(rab.le.ratoms) then
               write(2,'(3i4)') i-1,j-1,-1
            endif
          enddo
       enddo
       close(unit=2)
       end
