reader for xmol and turbomole geometry files

coord.f

      subroutine readxyz(fname)
      implicit real*8(x-z)
      integer,parameter   :: idx = 42
      character(len=*)    :: fname
      character(len=100)  :: line
      real*8, allocatable :: xyz(:,:)
      integer,allocatable :: atom(:)
      logical             :: exist
      integer             :: c
      
      inquire(file=fname,exist=exist)
      if(.not.exist) stop 'coordinates not found'
      
      open(idx,file=fname)
! - - xmol format
      if(index(fname,'.xyz').ne.0) then
        read(idx,*) c
        read(idx,*)
        allocate( xyz(3,c) )
        allocate( atom(c) )
        do i = 1, c
          read(idx,*) at,x,y,z
          if(index(pse,elem).eq.0) stop 'unknown element!'
          xyz(1,i) = x*aata
          xyz(2,i) = y*aata
          xyz(3,i) = z*aata
          atom(i)  = at
c         write(*,*)  xyz(:,i),atom(i)
        enddo
      else
! - - - turbomole format
        do
          read(idx,'(a)') line
          if(index(line,'$coord').ne.0) then
            c = 0
            do
              read(idx,'(a)') line
              if(index(line,'$').ne.0) exit
              c = c+1
            enddo
          endif
          if(index(line,'$end').ne.0) exit
        enddo
        rewind(idx)
        do
          read(idx,'(a)') line
          if(index(line,'$coord').ne.0) then
            allocate( xyz(3,c) )
            allocate( atom(c) )
            do i = 1, c
              read(idx,'(a)') line
              read(line,*) x,y,z,at
              xyz(1,i) = x
              xyz(2,i) = y
              xyz(3,i) = z
              atom(i)  = at
c             write(*,*) xyz(:,i),atom(i)
            enddo
          endif
          if(index(line,'$end').ne.0) exit
        enddo
      endif
      close(idx)
      
      if(c.eq.0) stop 'coordinates not found!'
      write(*,'(2x''number of atoms: '',i26)') c
      
      return
      end subroutine readxyz