read_xyz Subroutine

subroutine read_xyz(self, unit, error)

   !> Instance of the molecular structure data
   type(structure_type), intent(out) :: self

   !> File handle
   integer, intent(in) :: unit

   !> Error handling
   type(error_type), allocatable, intent(out) :: error

   integer :: ii, n, iat, stat, pos, lnum
   integer :: species_col, z_col, pos_col, ncols
   real(wp) :: x, y, z, conv, lattice(3, 3)
   real(wp), allocatable :: xyz(:, :)
   logical :: extended, have_lattice, have_pbc, periodic(3)
   type(token_type) :: token, tsym, tnat
   character(len=symbol_length) :: chdum
   character(len=symbol_length), allocatable :: sym(:)
   character(len=:), allocatable :: line, comment, fline, properties, ext_comment

   conv = aatoau
   lnum = 0

   call next_line(unit, fline, pos, lnum, stat)
   call read_next_token(fline, pos, tnat, n, stat)
   if (stat /= 0) then
      call io_error(error, "Could not read number of atoms", &
         & fline, tnat, filename(unit), lnum, "expected integer value")
      return
   end if

   if (n<1) then
      call io_error(error, "Impossible number of atoms provided", &
         & fline, tnat, filename(unit), lnum, "expected positive integer value")
      return
   end if

   allocate(sym(n))
   allocate(xyz(3, n))

   ! next record is either a plain XYZ comment or an Extended XYZ header
   call next_line(unit, comment, pos, lnum, stat)
   if (stat /= 0) then
      call io_error(error, "Unexpected end of file", &
         & "", token_type(0, 0), filename(unit), lnum+1, "expected value")
      return
   end if

   call parse_extxyz_header(comment, extended, properties, lattice, have_lattice, &
      & periodic, have_pbc, ext_comment, stat)
   if (stat /= 0) then
      call fatal_error(error, "Could not parse Extended XYZ header in '"//filename(unit)//"'")
      return
   end if

   species_col = 0
   z_col = 0
   pos_col = 0
   ncols = 4
   if (extended) then
      call parse_properties(properties, species_col, z_col, pos_col, ncols, stat)
      if (stat /= 0) then
         call fatal_error(error, "Invalid Properties specification in Extended XYZ file '"// &
            & filename(unit)//"'")
         return
      end if
   end if

   ii = 0
   do while (ii < n)
      call next_line(unit, line, pos, lnum, stat)
      if (is_iostat_end(stat)) exit
      if (stat /= 0) then
         call io_error(error, "Could not read geometry from xyz file", &
            & "", token_type(0, 0), filename(unit), lnum+1, "expected value")
         return
      end if

      if (extended) then
         call read_extxyz_atom(line, species_col, z_col, pos_col, ncols, &
            & chdum, x, y, z, token, stat)
         if (stat /= 0) then
            call io_error(error, "Could not parse atom data from Extended XYZ file", &
               & line, token, filename(unit), lnum, "unexpected value")
            return
         end if
         iat = to_number(chdum)
      else
         call next_token(line, pos, tsym)
         if (stat == 0) then
            call read_next_token(line, pos, token, x, stat)
         end if
         if (stat == 0) then
            call read_next_token(line, pos, token, y, stat)
         end if
         if (stat == 0) then
            call read_next_token(line, pos, token, z, stat)
         end if
         if (stat /= 0) then
            call io_error(error, "Could not parse coordinates from xyz file", &
               & line, token, filename(unit), lnum, "expected real value")
            return
         end if

         ! Adjust the token length to faithfully report the used chars in case of an error
         tsym%last = min(tsym%last, tsym%first + symbol_length - 1)
         chdum = line(tsym%first:tsym%last)
         iat = to_number(chdum)
         if (iat <= 0) then
            read(chdum, *, iostat=stat) iat
            if (stat == 0) then
               chdum = to_symbol(iat)
            else
               iat = 0
            end if
         end if
      end if

      if (iat > 0) then
         ii = ii+1
         sym(ii) = trim(chdum)
         xyz(:, ii) = [x, y, z]*conv
      else
         if (extended) then
            call io_error(error, "Cannot map symbol to atomic number", &
               & line, token, filename(unit), lnum, "unknown element")
         else
            call io_error(error, "Cannot map symbol to atomic number", &
               & line, tsym, filename(unit), lnum, "unknown element")
         end if
         return
      end if
   end do

   if (ii /= n) then
      call io_error(error, "Atom number missmatch in xyz file", &
         & fline, tnat, filename(unit), 1, "found "//to_string(ii)//" atoms in input")
      return
   end if

   if (extended) then
      if (have_lattice) then
         call new(self, sym, xyz, lattice=lattice, periodic=periodic)
      else
         call new(self, sym, xyz, periodic=periodic)
      end if
      if (allocated(ext_comment)) self%comment = ext_comment
   else
      call new(self, sym, xyz)
      if (len(comment) > 0) self%comment = comment
   end if

end subroutine read_xyz