fill_element_basis Subroutine

   pure subroutine fill_element_basis(basis_string, element_name, atom_basis, error)
      !! Fill in the shell data for a specific element from a GAMESS formatted basis string
      character(len=*), intent(in) :: basis_string
      character(len=*), intent(in) :: element_name
      type(atomic_basis_type), intent(inout) :: atom_basis
      type(error_t), intent(out) :: error

      integer :: line_start, line_end, line_type
      character(len=256) :: line
      logical :: in_data_block, in_target_element
      character(len=1) :: ang_mom
      integer :: nfunc, func_num, ishell, ifunc
      real(dp) :: exponent, coeff_s, coeff_p
      logical :: has_p

      ! L shell handling: we split into two shells, need to track both
      logical :: reading_l_shell
      integer :: l_shell_s_idx, l_shell_p_idx
      integer :: stat

      in_data_block = .false.
      in_target_element = .false.
      ishell = 0
      reading_l_shell = .false.

      line_start = 1
      do while (line_start <= len(basis_string))
         call get_next_line(basis_string, line_start, line, line_end)
         if (line_end == 0) exit

         line = adjustl(line)
         line_type = classify_line(line)

         select case (line_type)
            ! case (LINE_UNKNOWN)
            !   if (index(line, '$DATA') > 0) then
            !     in_data_block = .true.
            !   else if (index(line, '$END') > 0) then
            !     exit
            !   end if

         case (LINE_ATOM)
            if (strings_equal(line, element_name)) then
               in_target_element = .true.
            else
               if (in_target_element) exit
               in_target_element = .false.
            end if

         case (LINE_SHELL)
            if (in_target_element) then
               ! Parse shell header
               call parse_shell_header(line, ang_mom, nfunc, stat)
               if (stat /= 0) then
                  call error%set(ERROR_PARSE, "Failed to parse shell header: "//trim(line))
                  return
               end if

               if (ang_mom == 'L') then
                  ! L shell: create two shells (S and P)
                  reading_l_shell = .true.

                  ishell = ishell + 1
                  l_shell_s_idx = ishell
                  atom_basis%shells(ishell)%ang_mom = 0  ! S
                  call atom_basis%shells(ishell)%allocate_arrays(nfunc)

                  ishell = ishell + 1
                  l_shell_p_idx = ishell
                  atom_basis%shells(ishell)%ang_mom = 1  ! P
                  call atom_basis%shells(ishell)%allocate_arrays(nfunc)

                  ifunc = 0  ! Reset function counter
               else
                  ! Regular shell
                  reading_l_shell = .false.
                  ishell = ishell + 1

                  ! Set angular momentum (S=0, P=1, D=2, F=3, G=4, H=5, I=6)
                  atom_basis%shells(ishell)%ang_mom = ang_mom_char_to_int(ang_mom)

                  call atom_basis%shells(ishell)%allocate_arrays(nfunc)
                  ifunc = 0
               end if
            end if

         case (LINE_FUNCTION)
            if (in_target_element) then
               call parse_function_line(line, func_num, exponent, coeff_s, coeff_p, has_p, stat)
               if (stat /= 0) then
                  call error%set(ERROR_PARSE, "Failed to parse function line: "//trim(line))
                  return
               end if

               ifunc = ifunc + 1

               if (reading_l_shell) then
               if (.not. has_p) then
                  call error%set(ERROR_PARSE, "L shell requires both S and P coefficients")
                  return
               end if
               ! Store in both S and P shells
               atom_basis%shells(l_shell_s_idx)%exponents(ifunc) = exponent
               atom_basis%shells(l_shell_s_idx)%coefficients(ifunc) = coeff_s

               atom_basis%shells(l_shell_p_idx)%exponents(ifunc) = exponent
               atom_basis%shells(l_shell_p_idx)%coefficients(ifunc) = coeff_p
               else
               ! Store in current shell
               atom_basis%shells(ishell)%exponents(ifunc) = exponent
               atom_basis%shells(ishell)%coefficients(ifunc) = coeff_s
               end if
            end if

         case default
            ! Skip unknown line types (e.g., LINE_UNKNOWN, blank lines, comments)
            continue
         end select

         line_start = line_end
      end do

   end subroutine fill_element_basis