read_gaussian_external Subroutine

subroutine read_gaussian_external(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 :: stat, n, mode, chrg, spin, iat, ii, pos, lnum
   type(token_type) :: token, tnat
   character(len=:), allocatable :: line
   integer, allocatable :: at(:)
   real(wp), allocatable :: xyz(:,:)
   real(wp) :: coord(3), q

   lnum = 0
   call next_line(unit, line, pos, lnum, stat)
   if (stat == 0) then
      token = token_type(1, 10)
      tnat = token
      call read_token(line, token, n, stat)
   end if
   if (stat == 0) then
      token = token_type(11, 20)
      call read_token(line, token, mode, stat)
   end if
   if (stat == 0) then
      token = token_type(21, 30)
      call read_token(line, token, chrg, stat)
   end if
   if (stat == 0) then
      token = token_type(31, 40)
      call read_token(line, token, spin, stat)
   end if
   if (stat /= 0) then
      call io_error(error, "Could not read number of atoms", &
         & line, token, filename(unit), lnum, "expected integer value")
      return
   end if

   if (n <= 0) then
      call io_error(error, "Found no atoms, cannot work without atoms!", &
         & line, tnat, filename(unit), lnum, "expected positive integer")
      return
   end if

   allocate(xyz(3, n))
   allocate(at(n))

   ii = 0
   do while (ii < n)
      call next_line(unit, line, pos, lnum, stat)
      if (is_iostat_end(stat)) exit
      if (stat == 0) then
         token = token_type(1, 10)
         tnat = token
         call read_token(line, token, iat, stat)
      end if
      if (stat == 0) then
         token = token_type(11, 30)
         call read_token(line, token, coord(1), stat)
      end if
      if (stat == 0) then
         token = token_type(31, 50)
         call read_token(line, token, coord(2), stat)
      end if
      if (stat == 0) then
         token = token_type(51, 70)
         call read_token(line, token, coord(3), stat)
      end if
      if (stat == 0) then
         token = token_type(71, 90)
         call read_token(line, token, q, stat)
      end if
      if (stat /= 0) then
         call io_error(error, "Could not read geometry from Gaussian file", &
            & line, token, filename(unit), lnum, "unexpected value")
         return
      end if
      if (iat > 0) then
         ii = ii+1
         at(ii) = iat
         xyz(:, ii) = coord
      else
         call io_error(error, "Invalid atomic number", &
            & line, tnat, filename(unit), lnum, "expected positive integer")
         return
      end if
   end do

   call new(self, at, xyz, charge=real(chrg, wp), uhf=spin)

   if (ii /= n) then
      call fatal_error(error, "Atom number missmatch in Gaussian file")
      return
   end if

end subroutine read_gaussian_external