param_get_atoms Subroutine

  subroutine param_get_atoms(lunits)
    !===================================!
    !                                   !
    !!   Fills the atom data block
    !                                   !
    !===================================!

    use w90_constants, only: bohr
    use w90_utility, only: utility_frac_to_cart, utility_cart_to_frac
    use w90_io, only: io_error
    implicit none

    logical, intent(in) :: lunits
    !! Do we expect a first line with the units

    real(kind=dp)     :: atoms_pos_frac_tmp(3, num_atoms)
    real(kind=dp)     :: atoms_pos_cart_tmp(3, num_atoms)
    character(len=20) :: keyword
    integer           :: in, ins, ine, loop, i, line_e, line_s, counter
    integer           :: i_temp, loop2, max_sites, ierr, ic
    logical           :: found_e, found_s, found, frac
    character(len=maxlen) :: dummy, end_st, start_st
    character(len=maxlen) :: ctemp(num_atoms)
    character(len=maxlen) :: atoms_label_tmp(num_atoms)
    logical           :: lconvert

    keyword = "atoms_cart"
    frac = .false.
    call param_get_block_length("atoms_frac", found, i_temp)
    if (found) then
      keyword = "atoms_frac"
      frac = .true.
    end if

    found_s = .false.
    found_e = .false.

    start_st = 'begin '//trim(keyword)
    end_st = 'end '//trim(keyword)

    do loop = 1, num_lines
      ins = index(in_data(loop), trim(keyword))
      if (ins == 0) cycle
      in = index(in_data(loop), 'begin')
      if (in == 0 .or. in > 1) cycle
      line_s = loop
      if (found_s) then
        call io_error('Error: Found '//trim(start_st)//' more than once in input file')
      endif
      found_s = .true.
    end do

    do loop = 1, num_lines
      ine = index(in_data(loop), trim(keyword))
      if (ine == 0) cycle
      in = index(in_data(loop), 'end')
      if (in == 0 .or. in > 1) cycle
      line_e = loop
      if (found_e) then
        call io_error('Error: Found '//trim(end_st)//' more than once in input file')
      endif
      found_e = .true.
    end do

    if (.not. found_e) then
      call io_error('Error: Found '//trim(start_st)//' but no '//trim(end_st)//' in input file')
    end if

    if (line_e <= line_s) then
      call io_error('Error: '//trim(end_st)//' comes before '//trim(start_st)//' in input file')
    end if

    lconvert = .false.
    if (lunits) then
      dummy = in_data(line_s + 1)
      if (index(dummy, 'ang') .ne. 0) then
        lconvert = .false.
      elseif (index(dummy, 'bohr') .ne. 0) then
        lconvert = .true.
      else
        call io_error('Error: Units in block atoms_cart not recognised in param_get_atoms')
      endif
      in_data(line_s) (1:maxlen) = ' '
      line_s = line_s + 1
    endif

    counter = 0
    do loop = line_s + 1, line_e - 1
      dummy = in_data(loop)
      counter = counter + 1
      if (frac) then
        read (dummy, *, err=240, end=240) atoms_label_tmp(counter), (atoms_pos_frac_tmp(i, counter), i=1, 3)
      else
        read (dummy, *, err=240, end=240) atoms_label_tmp(counter), (atoms_pos_cart_tmp(i, counter), i=1, 3)
      end if
    end do

    if (lconvert) atoms_pos_cart_tmp = atoms_pos_cart_tmp*bohr

    in_data(line_s:line_e) (1:maxlen) = ' '

    if (frac) then
      do loop = 1, num_atoms
        call utility_frac_to_cart(atoms_pos_frac_tmp(:, loop), atoms_pos_cart_tmp(:, loop), real_lattice)
      end do
    else
      do loop = 1, num_atoms
        call utility_cart_to_frac(atoms_pos_cart_tmp(:, loop), atoms_pos_frac_tmp(:, loop), recip_lattice)
      end do
    end if

    ! Now we sort the data into the proper structures
    num_species = 1
    ctemp(1) = atoms_label_tmp(1)
    do loop = 2, num_atoms
      do loop2 = 1, loop - 1
        if (trim(atoms_label_tmp(loop)) == trim(atoms_label_tmp(loop2))) exit
        if (loop2 == loop - 1) then
          num_species = num_species + 1
          ctemp(num_species) = atoms_label_tmp(loop)
        end if
      end do
    end do

    allocate (atoms_species_num(num_species), stat=ierr)
    if (ierr /= 0) call io_error('Error allocating atoms_species_num in param_get_atoms')
    allocate (atoms_label(num_species), stat=ierr)
    if (ierr /= 0) call io_error('Error allocating atoms_label in param_get_atoms')
    allocate (atoms_symbol(num_species), stat=ierr)
    if (ierr /= 0) call io_error('Error allocating atoms_symbol in param_get_atoms')
    atoms_species_num(:) = 0

    do loop = 1, num_species
      atoms_label(loop) = ctemp(loop)
      do loop2 = 1, num_atoms
        if (trim(atoms_label(loop)) == trim(atoms_label_tmp(loop2))) then
          atoms_species_num(loop) = atoms_species_num(loop) + 1
        end if
      end do
    end do

    max_sites = maxval(atoms_species_num)
    allocate (atoms_pos_frac(3, max_sites, num_species), stat=ierr)
    if (ierr /= 0) call io_error('Error allocating atoms_pos_frac in param_get_atoms')
    allocate (atoms_pos_cart(3, max_sites, num_species), stat=ierr)
    if (ierr /= 0) call io_error('Error allocating atoms_pos_cart in param_get_atoms')

    do loop = 1, num_species
      counter = 0
      do loop2 = 1, num_atoms
        if (trim(atoms_label(loop)) == trim(atoms_label_tmp(loop2))) then
          counter = counter + 1
          atoms_pos_frac(:, counter, loop) = atoms_pos_frac_tmp(:, loop2)
          atoms_pos_cart(:, counter, loop) = atoms_pos_cart_tmp(:, loop2)
        end if
      end do
    end do

    ! Strip any numeric characters from atoms_label to get atoms_symbol
    do loop = 1, num_species
      atoms_symbol(loop) (1:2) = atoms_label(loop) (1:2)
      ic = ichar(atoms_symbol(loop) (2:2))
      if ((ic .lt. ichar('a')) .or. (ic .gt. ichar('z'))) &
        atoms_symbol(loop) (2:2) = ' '
    end do

    return

240 call io_error('Error: Problem reading block keyword '//trim(keyword))

  end subroutine param_get_atoms