read_qchem Subroutine

subroutine read_qchem(mol, unit, error)

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

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

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

   integer :: stat, pos, lnum, izp, iat
   integer :: charge, multiplicity
   type(token_type) :: token
   character(len=:), allocatable :: line
   real(wp) :: x, y, z
   character(len=symbol_length), allocatable :: sym(:)
   real(wp), allocatable :: xyz(:, :), abc(:, :), lattice(:, :)
   logical :: is_frac, periodic(3)

   iat = 0
   lnum = 0
   stat = 0

   do while(stat == 0)
      call next_line(unit, line, pos, lnum, stat)
      if (stat /= 0) exit

      call next_token(line, pos, token)
      if (token%first > len(line)) cycle
      if (to_lower(line(token%first:token%last)) == '$molecule') exit
   end do

   if (stat /= 0) then
      call io_error(error, "No atoms found", &
         & line, token_type(0, 0), filename(unit), lnum+1, "expected molecule block")
      return
   end if

   call next_line(unit, line, pos, lnum, stat)
   if (stat == 0) &
      call read_next_token(line, pos, token, charge, stat)
   if (stat == 0) &
      call read_next_token(line, pos, token, multiplicity, stat)
   if (stat /= 0) then
      call io_error(error, "Failed to read charge and multiplicity", &
         & line, token, filename(unit), lnum, "expected integer value")
      return
   end if

   allocate(sym(initial_size), source=repeat(' ', symbol_length))
   allocate(xyz(3, initial_size), source=0.0_wp)

   do while(stat == 0)
      call next_line(unit, line, pos, lnum, stat)
      if (stat /= 0) exit

      call next_token(line, pos, token)
      if (to_lower(line(token%first:token%last)) == '$end') exit

      if (iat >= size(sym)) call resize(sym)
      if (iat >= size(xyz, 2)) call resize(xyz)
      iat = iat + 1

      token%last = min(token%last, token%first + symbol_length - 1)
      sym(iat) = line(token%first:token%last)
      if (to_number(sym(iat)) == 0) then
         call read_token(line, token, izp, stat)
         sym(iat) = to_symbol(izp)
      end if
      if (stat /= 0) then
         call io_error(error, "Cannot map symbol to atomic number", &
            & line, token, filename(unit), lnum, "unknown element")
         return
      end if

      call read_next_token(line, pos, token, x, stat)
      if (stat == 0) &
         call read_next_token(line, pos, token, y, stat)
      if (stat == 0) &
         call read_next_token(line, pos, token, z, stat)
      if (stat /= 0) then
         call io_error(error, "Cannot read coordinates", &
            & line, token, filename(unit), lnum, "expected real value")
         return
      end if

      xyz(:, iat) = [x, y, z] * aatoau
   end do

   if (stat /= 0) then
      call io_error(error, "Failed to read molecule block", &
         & line, token_type(0, 0), filename(unit), lnum, "unexpected end of input")
      return
   end if

   call new(mol, sym(:iat), xyz, charge=real(charge, wp), uhf=multiplicity-1)
end subroutine read_qchem