read_genformat Subroutine

subroutine read_genformat(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

   character(len=:), allocatable :: line
   integer :: natoms, nspecies, iatom, dummy, isp, ilat, stat, istart, iend
   logical :: cartesian, periodic(3)
   real(wp) :: coord(3), origin(3)
   character(len=1) :: variant
   type(token_type) :: token
   character(len=symbol_length), allocatable :: species(:), sym(:)
   real(wp), allocatable :: xyz(:, :), abc(:, :), lattice(:, :)
   type(structure_info) :: info
   integer :: pos, lnum

   lnum = 0
   call advance_line(unit, line, pos, lnum, stat)
   call read_next_token(line, pos, token, natoms, stat)
   if (stat /= 0 .or. natoms < 1) then
      call io_error(error, "Could not read number of atoms", &
         & line, token, filename(unit), lnum, "expected integer value")
      return
   end if

   allocate(species(natoms))
   allocate(sym(natoms))
   allocate(xyz(3, natoms))
   allocate(abc(3, natoms))

   call next_token(line, pos, token)
   select case(line(token%first:token%last))
   case('c', 'C')
      cartesian = .true.
      periodic = .false.
   case('s', 'S')
      cartesian = .true.
      periodic = .true.
      allocate(lattice(3, 3), source=0.0_wp)
   case('f', 'F')
      cartesian = .false.
      periodic = .true.
      allocate(lattice(3, 3), source=0.0_wp)
   case('h', 'H')
      cartesian = .true.
      periodic = [.false., .false., .true.]
      allocate(lattice(3, 1), source=0.0_wp)
   case default
      call io_error(error, "Invalid input version found", &
         & line, token, filename(unit), lnum, "unknown identifier")
      return
   end select

   call advance_line(unit, line, pos, lnum, stat)
   isp = 0
   do while(pos < len(line))
      call next_token(line, pos, token)
      isp = isp + 1
      token%last = min(token%last, token%first + symbol_length - 1)
      species(isp) = line(token%first:token%last)
      if (to_number(species(isp)) == 0) then
         call io_error(error, "Cannot map symbol to atomic number", &
            & line, token, filename(unit), lnum, "unknown element")
         return
      end if
   end do
   nspecies = isp

   do iatom = 1, natoms
      token = token_type(0, 0)
      call advance_line(unit, line, pos, lnum, stat)
      if (stat == 0) &
         call read_next_token(line, pos, token, dummy, stat)
      if (stat == 0) &
         call read_next_token(line, pos, token, isp, stat)
      if (stat == 0) &
         call read_next_token(line, pos, token, coord(1), stat)
      if (stat == 0) &
         call read_next_token(line, pos, token, coord(2), stat)
      if (stat == 0) &
         call read_next_token(line, pos, token, coord(3), stat)
      if (stat /= 0) then
         call io_error(error, "Cannot read coordinates", &
            & line, token, filename(unit), lnum, "unexpected value")
         return
      end if
      sym(iatom) = species(isp)
      if (cartesian) then
         xyz(:, iatom) = coord * aatoau
      else
         abc(:, iatom) = coord
      end if
   end do

   if (any(periodic)) then
      call advance_line(unit, line, pos, lnum, stat)
      if (stat /= 0) then
         call io_error(error, "Unexpected end of file", &
            & line, token_type(0, 0), filename(unit), lnum, "missing lattice information")
         return
      end if
      if (stat == 0) &
         call read_next_token(line, pos, token, origin(1), stat)
      if (stat == 0) &
         call read_next_token(line, pos, token, origin(2), stat)
      if (stat == 0) &
         call read_next_token(line, pos, token, origin(3), stat)
      if (stat /= 0) then
         call io_error(error, "Cannot read origin", &
            & line, token, filename(unit), lnum, "expected real value")
         return
         end if
   end if

   if (all(periodic)) then
      do ilat = 1, 3
         call advance_line(unit, line, pos, lnum, stat)
         if (stat == 0) &
            call read_next_token(line, pos, token, coord(1), stat)
         if (stat == 0) &
            call read_next_token(line, pos, token, coord(2), stat)
         if (stat == 0) &
            call read_next_token(line, pos, token, coord(3), stat)
         if (stat /= 0) then
            call io_error(error, "Cannot read lattice vector", &
               & line, token, filename(unit), lnum, "expected real value")
            return
         end if
         lattice(:, ilat) = coord * aatoau
      end do
      if (.not.cartesian) then
         xyz = matmul(lattice, abc)
      end if
   end if

   if (count(periodic) == 1) then
      call advance_line(unit, line, pos, lnum, stat)
      if (stat == 0) &
         call read_next_token(line, pos, token, coord(1), stat)
      if (stat == 0) &
         call read_next_token(line, pos, token, coord(2), stat)
      if (stat == 0) &
         call read_next_token(line, pos, token, coord(3), stat)
      if (stat /= 0) then
         call io_error(error, "Cannot read lattice vector", &
            & line, token, filename(unit), lnum, "expected real value")
         return
      end if
      if (coord(3) < 1) then
         call io_error(error, "Invalid helical axis rotation order", &
            & line, token, filename(unit), lnum, "expected positive value")
         return
      end if

      ! Store helical axis in *first* lattice vector, however it is not an actual
      ! lattice vector as on would expect but a screw axis
      lattice(:, 1) = [coord(1) * aatoau, coord(2) * pi / 180.0_wp, coord(3)]
   end if

   if (any(periodic)) then
      xyz(:, :) = xyz - spread(origin, 2, natoms)
   end if

   info = structure_info(cartesian=cartesian)
   call new(mol, sym, xyz, lattice=lattice, periodic=periodic, info=info)

contains

subroutine advance_line(unit, line, pos, num, stat)
   integer,intent(in) :: unit
   integer, intent(out) :: pos
   integer, intent(inout) :: num
   character(len=:), allocatable, intent(out) :: line
   integer, intent(out) :: stat
   integer :: ihash

   stat = 0
   do while(stat == 0)
      call next_line(unit, line, pos, num, stat)
      ihash = index(line, '#')
      if (ihash > 0) line = line(:ihash-1)
      if (len_trim(line) > 0) exit
   end do
   line = trim(adjustl(line))
end subroutine advance_line

end subroutine read_genformat