read_xyz_string Subroutine

   pure subroutine read_xyz_string(xyz_string, geom, error)
      !! Parse molecular geometry from XYZ format string
      character(len=*), intent(in) :: xyz_string
      type(geometry_type), intent(out) :: geom
      type(error_t), intent(out) :: error

      character(len=:), allocatable :: lines(:)
      integer :: nlines, iatom, io_stat
      character(len=256) :: element
      real(dp) :: x, y, z

      ! Split into lines
      call split_lines(xyz_string, lines, nlines)

      if (nlines < 2) then
         call error%set(ERROR_PARSE, "XYZ file must have at least 2 lines (natoms + comment)")
         return
      end if

      ! Read number of atoms from first line
      read (lines(1), *, iostat=io_stat) geom%natoms
      if (io_stat /= 0) then
         call error%set(ERROR_PARSE, "Failed to read number of atoms from first line")
         return
      end if

      if (geom%natoms < 0) then
         call error%set(ERROR_PARSE, "Number of atoms must be non-negative")
         return
      end if

      ! Store comment line
      geom%comment = trim(adjustl(lines(2)))

      ! Check we have enough lines
      if (nlines < 2 + geom%natoms) then
         call error%set(ERROR_PARSE, "XYZ file has insufficient lines: expected "// &
                        trim(int_to_string(2 + geom%natoms))//", got "// &
                        trim(int_to_string(nlines)))
         return
      end if

      ! Allocate arrays
      allocate (character(len=MAX_ELEMENT_SYMBOL_LEN) :: geom%elements(geom%natoms))
      allocate (geom%coords(3, geom%natoms))

      ! Read atom data
      do iatom = 1, geom%natoms
         read (lines(2 + iatom), *, iostat=io_stat) element, x, y, z
         if (io_stat /= 0) then
            call error%set(ERROR_PARSE, "Failed to parse atom data on line "// &
                           trim(int_to_string(2 + iatom))//": '"// &
                           trim(lines(2 + iatom))//"'")
            return
         end if

         geom%elements(iatom) = trim(adjustl(element))
         geom%coords(1, iatom) = x
         geom%coords(2, iatom) = y
         geom%coords(3, iatom) = z
      end do

   end subroutine read_xyz_string