readcon-core 0.23.1

use crate::error::ParseError;
use crate::helpers::symbol_to_atomic_number;
use crate::types::{
    AtomDatum, ConFrame, FrameHeader, PreboxHeader, SECTION_CHARGES, SECTION_DISPLACEMENTS,
    SECTION_ENERGIES, SECTION_FORCES, SECTION_MAGMOMS, SECTION_SPINS, SECTION_SPREADS,
    SECTION_VELOCITIES, meta,
};
use serde_json::Value;
use std::iter::Peekable;
use std::sync::Arc;
use std::sync::atomic::{AtomicUsize, Ordering};

/// Line source with peek for section detection (blank-line separators).
///
/// Implemented by the memchr cursor on the hot iterator path and by
/// `Peekable<I>` for tests / ad-hoc callers.
/// Rows of one coordinate block, in file order.
///
/// `atoms` is empty unless the caller asked for per-atom shells. The
/// numeric columns are the bytes a serial `next_coord_row` loop would
/// have stored.
#[doc(hidden)]
pub struct CoordRows {
    pub xyz: Vec<f64>,
    pub ids: Vec<u64>,
    pub flags: Vec<[bool; 3]>,
    pub atoms: Vec<crate::types::AtomDatum>,
    pub consumed: usize,
}

/// Worker count for a coordinate block that is too wide to hand to one task.
///
/// One means the calling thread parses the block. The frame iterator sets
/// this only for a single wide frame, then puts it back.
pub(crate) fn wide_workers() -> usize {
    WIDE_WORKERS.with(|slot| slot.get())
}

pub(crate) fn set_wide_workers(n: usize) {
    WIDE_WORKERS.with(|slot| slot.set(n.max(1)));
}

std::thread_local! {
    static WIDE_WORKERS: std::cell::Cell<usize> = const { std::cell::Cell::new(1) };
}

/// Fill `n` coordinate rows.
///
/// `breaks[i]` is the absolute index of that row's `\n`. Row 0 begins at
/// `row0`. `xyz` receives `n * 3` values. `ids` and `flags` may be null
/// when the caller does not want them. `flags` is the packed column-4 mask.
#[cfg(feature = "parallel")]
pub(crate) fn fill_rows_from_breaks(
    bytes: &[u8],
    breaks: &[usize],
    row0: usize,
    spec_version: u32,
    start_id: u64,
    workers: usize,
    xyz: *mut f64,
    ids: *mut u64,
    flags: *mut u8,
) -> Result<(), ParseError> {
    let n = breaks.len();
    if n == 0 {
        return Ok(());
    }
    let chunk = n.div_ceil(workers.max(1)).max(1);
    let groups = n.div_ceil(chunk);
    let xyz_ptr = xyz as usize;
    let ids_ptr = ids as usize;
    let flags_ptr = flags as usize;
    let fill = || {
        use rayon::prelude::*;
        (0..groups).into_par_iter().try_for_each(|group| {
            let first = group * chunk;
            let last = (first + chunk).min(n);
            for row_index in first..last {
                let start = if row_index == 0 {
                    row0
                } else {
                    breaks[row_index - 1] + 1
                };
                let end = breaks[row_index] + 1;
                if end > bytes.len() || start > end {
                    return Err(ParseError::IncompleteFrame);
                }
                let default_id = start_id + row_index as u64;
                let ((row, fixed, atom_id), _next) =
                    parse_coord_row_at(&bytes[start..end], 0, spec_version, default_id)?;
                unsafe {
                    std::ptr::copy_nonoverlapping(
                        row.as_ptr(),
                        (xyz_ptr as *mut f64).add(row_index * 3),
                        3,
                    );
                    if ids_ptr != 0 {
                        *((ids_ptr as *mut u64).add(row_index)) = atom_id;
                    }
                    if flags_ptr != 0 {
                        let mask = u8::from(fixed[0])
                            | (u8::from(fixed[1]) << 1)
                            | (u8::from(fixed[2]) << 2);
                        *((flags_ptr as *mut u8).add(row_index)) = mask;
                    }
                }
            }
            Ok(())
        })
    };
    on_pool(workers, fill)
}

/// Run `body` on a pool of `workers` threads.
///
/// The pool for a given size is kept for the process. Building one per
/// call was a serial pause on the wide-frame read and the row formatter.
#[cfg(feature = "parallel")]
pub(crate) fn on_pool<T: Send>(workers: usize, body: impl FnOnce() -> T + Send) -> T {
    if workers < 2 {
        return body();
    }
    // The cached pool is leaked for the process. Miri rejects a thread
    // that is still running when the test binary exits, so it builds a
    // pool, runs `body`, and joins.
    #[cfg(miri)]
    {
        return match rayon::ThreadPoolBuilder::new().num_threads(workers).build() {
            Ok(pool) => pool.install(body),
            Err(_) => body(),
        };
    }
    #[cfg(not(miri))]
    {
        let pool = {
            use std::sync::{Mutex, OnceLock};
            static POOLS: OnceLock<Mutex<Vec<(usize, &'static rayon::ThreadPool)>>> =
                OnceLock::new();
            let pools = POOLS.get_or_init(|| Mutex::new(Vec::new()));
            let mut guard = pools.lock().unwrap_or_else(|err| err.into_inner());
            if let Some((_, pool)) = guard.iter().find(|(size, _)| *size == workers) {
                *pool
            } else {
                let next_core = std::sync::Arc::new(std::sync::atomic::AtomicUsize::new(0));
                let built = match rayon::ThreadPoolBuilder::new()
                    .num_threads(workers)
                    .spawn_handler(move |thread| {
                        // One worker per physical core. Leaving the placement to
                        // the kernel put two of these on one core and left
                        // another idle, which is the whole 8-thread gap.
                        let core = next_core.fetch_add(1, std::sync::atomic::Ordering::Relaxed);
                        std::thread::Builder::new().spawn(move || {
                            pin_worker(core);
                            thread.run();
                        })?;
                        Ok(())
                    })
                    .build()
                {
                    Ok(pool) => pool,
                    Err(_) => return body(),
                };
                let pool: &'static rayon::ThreadPool = Box::leak(Box::new(built));
                guard.push((workers, pool));
                pool
            }
        };
        pool.install(body)
    }
}

pub(crate) fn pin_worker(core: usize) {
    // `sched_setaffinity` is Linux. Pinning is a performance hint.
    #[cfg(all(not(miri), target_os = "linux"))]
    unsafe {
        let mut set: libc::cpu_set_t = std::mem::zeroed();
        libc::CPU_ZERO(&mut set);
        libc::CPU_SET(core, &mut set);
        libc::sched_setaffinity(0, std::mem::size_of::<libc::cpu_set_t>(), &set);
    }
    #[cfg(any(miri, not(target_os = "linux")))]
    {
        let _ = core;
    }
}

/// Coordinate rows in `bytes`, which begins at the first row.
///
/// `n` rows are required. The returned `consumed` count is the number of
/// bytes those rows occupy, including their newlines. `symbol` builds the
/// per-atom shells; `None` leaves `atoms` empty. A structural miss returns
/// `Err` and the caller keeps the serial loop.
#[cfg(feature = "parallel")]
pub(crate) fn parallel_coord_rows(
    bytes: &[u8],
    n: usize,
    spec_version: u32,
    start_id: u64,
    workers: usize,
    symbol: Option<&std::sync::Arc<str>>,
) -> Result<CoordRows, ParseError> {
    if n == 0 {
        return Ok(CoordRows {
            xyz: Vec::new(),
            ids: Vec::new(),
            flags: Vec::new(),
            atoms: Vec::new(),
            consumed: 0,
        });
    }
    // Same chunking as the column reader: count newlines, then each worker
    // walks its own run. An offset per line was a serial pass the plain
    // read does not pay.
    let (chunks, consumed) = plan_coord_chunks(bytes, 0, n, workers.max(1))?;
    let mut xyz = Vec::<f64>::with_capacity(n.saturating_mul(3));
    let mut ids = Vec::<u64>::with_capacity(n);
    let mut flags = Vec::<[bool; 3]>::with_capacity(n);
    let mut atoms_raw: Vec<std::mem::MaybeUninit<crate::types::AtomDatum>> = Vec::new();
    if symbol.is_some() {
        atoms_raw = Vec::with_capacity(n);
        unsafe { atoms_raw.set_len(n) };
    }
    let xyz_ptr = xyz.as_mut_ptr() as usize;
    let ids_ptr = ids.as_mut_ptr() as usize;
    let flags_ptr = flags.as_mut_ptr() as usize;
    let atoms_ptr = atoms_raw.as_mut_ptr() as usize;
    let mut row0 = Vec::with_capacity(chunks.len());
    let mut acc = 0usize;
    for (_, nrows) in &chunks {
        row0.push(acc);
        acc += *nrows;
    }
    let fill_chunk = |index: usize| -> Result<(), ParseError> {
        let (start, nrows) = chunks[index];
        let mut row = row0[index];
        let mut pos = start;
        let mut shares = symbol.map(|symbol| symbol_shares(symbol, nrows));
        for _ in 0..nrows {
            let ((xyz, fixed, atom_id), next) =
                parse_coord_row_at(bytes, pos, spec_version, start_id + row as u64)?;
            unsafe {
                std::ptr::copy_nonoverlapping(xyz.as_ptr(), (xyz_ptr as *mut f64).add(row * 3), 3);
                *((ids_ptr as *mut u64).add(row)) = atom_id;
                *((flags_ptr as *mut [bool; 3]).add(row)) = fixed;
                if let Some(shares) = shares.as_mut() {
                    let slot = (xyz_ptr as *mut f64).add(row * 3);
                    let atom =
                        crate::types::AtomDatum::at_slot(shares.take(), slot, fixed, atom_id);
                    (*(atoms_ptr as *mut std::mem::MaybeUninit<crate::types::AtomDatum>).add(row))
                        .write(atom);
                }
            }
            pos = next;
            row += 1;
        }
        Ok(())
    };
    let built = if workers < 2 || chunks.len() < 2 {
        (0..chunks.len()).try_for_each(fill_chunk)
    } else {
        on_pool(workers, || {
            use rayon::prelude::*;
            (0..chunks.len()).into_par_iter().try_for_each(fill_chunk)
        })
    };
    if let Err(err) = built {
        std::mem::forget(atoms_raw);
        return Err(err);
    }
    // SAFETY: every row was written before `built` returned Ok. The
    // length stays 0 on the error path above.
    unsafe {
        xyz.set_len(n * 3);
        ids.set_len(n);
        flags.set_len(n);
    }
    let atoms = if symbol.is_some() {
        let mut raw = std::mem::ManuallyDrop::new(atoms_raw);
        unsafe {
            Vec::from_raw_parts(
                raw.as_mut_ptr() as *mut crate::types::AtomDatum,
                n,
                raw.capacity(),
            )
        }
    } else {
        Vec::new()
    };
    Ok(CoordRows {
        xyz,
        ids,
        flags,
        atoms,
        consumed,
    })
}

/// Walk `nrows` coordinate lines starting at byte `start`.
///
/// `row_base` is the first row's index in `xyz` / `ids` / `flags` / `atoms`.
/// This is the per-member kernel of [`fill_rows_from_chunks`], so a worker
/// can parse a slice as soon as that slice has been inflated.
#[cfg(feature = "parallel")]
pub(crate) fn parse_coord_run(
    bytes: &[u8],
    start: usize,
    nrows: usize,
    spec_version: u32,
    row_base: usize,
    symbol: Option<&std::sync::Arc<str>>,
    xyz: *mut f64,
    ids: *mut u64,
    flags: *mut [bool; 3],
    atoms: *mut std::mem::MaybeUninit<crate::types::AtomDatum>,
) -> Result<usize, ParseError> {
    let mut pos = start;
    let mut shares = symbol.map(|symbol| symbol_shares(symbol, nrows));
    for step in 0..nrows {
        let row = row_base + step;
        let ((xyz_row, fixed, atom_id), next) =
            parse_coord_row_at(bytes, pos, spec_version, row as u64)?;
        unsafe {
            let dest = xyz.add(row * 3);
            std::ptr::copy_nonoverlapping(xyz_row.as_ptr(), dest, 3);
            *ids.add(row) = atom_id;
            *flags.add(row) = fixed;
            if let Some(shares) = shares.as_mut() {
                let atom = crate::types::AtomDatum::at_slot(shares.take(), dest, fixed, atom_id);
                (*atoms.add(row)).write(atom);
            }
        }
        pos = next;
    }
    Ok(pos)
}

#[cfg(feature = "parallel")]
#[allow(dead_code)]
fn coord_line_ranges(
    bytes: &[u8],
    n: usize,
    workers: usize,
) -> Result<Vec<(usize, usize)>, ParseError> {
    let newlines = newline_index(bytes, workers);
    let mut ranges = Vec::with_capacity(n);
    let mut start = 0usize;
    for index in 0..n {
        if start > bytes.len() {
            return Err(ParseError::IncompleteFrame);
        }
        let end = if index < newlines.len() {
            newlines[index] + 1
        } else if index == newlines.len() && start < bytes.len() {
            bytes.len()
        } else {
            return Err(ParseError::IncompleteFrame);
        };
        if end < start {
            return Err(ParseError::IncompleteFrame);
        }
        ranges.push((start, end));
        start = end;
    }
    Ok(ranges)
}

/// Newline positions in `bytes`.
///
/// Chunks are cut on byte boundaries, then each interior cut that lands
/// inside a line is moved to the next newline so a row is not split.
/// Workers count those newlines, a prefix sum places each run, and the
/// same workers write the offsets. The offset array is not zeroed, so
/// the writer is the first touch. `workers < 2` scans on the calling thread.
#[cfg(feature = "parallel")]
pub(crate) fn newline_index(bytes: &[u8], workers: usize) -> Vec<usize> {
    newline_index_at(bytes, 0, workers)
}

/// [`newline_index`] with `base` added to every stored offset.
#[cfg(feature = "parallel")]
pub(crate) fn newline_index_at(bytes: &[u8], base: usize, workers: usize) -> Vec<usize> {
    if bytes.is_empty() || workers < 2 || bytes.len() < (1 << 20) {
        return memchr::memchr_iter(b'\n', bytes)
            .map(|index| base + index)
            .collect();
    }
    let cuts = newline_cuts(bytes, workers);
    newline_index_exact(&cuts, bytes, base)
}

/// One cache line, so adjacent chunk counters do not share a line.
#[cfg(feature = "parallel")]
#[derive(Clone)]
#[repr(C, align(64))]
struct ChunkTally {
    n: usize,
}

#[cfg(feature = "parallel")]
pub(crate) fn newline_cuts(bytes: &[u8], workers: usize) -> Vec<usize> {
    let workers = workers.min(bytes.len()).max(1);
    let mut cuts = Vec::with_capacity(workers + 1);
    cuts.push(0usize);
    for index in 1..workers {
        let mut at = index * bytes.len() / workers;
        // A cut in the middle of a row belongs to the chunk that holds
        // the newline. Move it to the first byte after that newline.
        if at > 0 && at < bytes.len() && bytes[at - 1] != b'\n' {
            if let Some(rel) = memchr::memchr(b'\n', &bytes[at..]) {
                at += rel + 1;
            } else {
                at = bytes.len();
            }
        }
        if cuts.last().copied() != Some(at) {
            cuts.push(at);
        }
    }
    if cuts.last().copied() != Some(bytes.len()) {
        cuts.push(bytes.len());
    }
    cuts
}

/// Count newlines in each `cuts` window. One store per chunk, so the
/// counters do not bounce a cache line during the scan.
#[cfg(feature = "parallel")]
fn count_newline_chunks(bytes: &[u8], cuts: &[usize]) -> Vec<usize> {
    let nchunks = cuts.len().saturating_sub(1);
    let mut tallies = vec![ChunkTally { n: 0 }; nchunks];
    if nchunks == 0 {
        return Vec::new();
    }
    let ptr = tallies.as_mut_ptr() as usize;
    let count = || {
        use rayon::prelude::*;
        (0..nchunks).into_par_iter().for_each(|index| {
            let n = memchr::memchr_iter(b'\n', &bytes[cuts[index]..cuts[index + 1]]).count();
            unsafe {
                (*(ptr as *mut ChunkTally).add(index)).n = n;
            }
        });
    };
    if nchunks < 2 {
        count();
    } else {
        on_pool(nchunks, count);
    }
    tallies.into_iter().map(|tally| tally.n).collect()
}

/// Byte start and row count for each worker that parses `n_rows`
/// coordinate lines beginning at `start`.
///
/// The newline *count* places the rows. Workers then parse their own
/// slice in order, so the read does not materialise an offset per line.
/// The returned end is the first byte after the `n_rows`-th newline.
#[cfg(feature = "parallel")]
pub(crate) fn plan_coord_chunks(
    bytes: &[u8],
    start: usize,
    n_rows: usize,
    workers: usize,
) -> Result<(Vec<(usize, usize)>, usize), ParseError> {
    if n_rows == 0 {
        return Ok((Vec::new(), start));
    }
    if start >= bytes.len() {
        return Err(ParseError::IncompleteFrame);
    }
    let mut window = n_rows.saturating_mul(96).max(1);
    loop {
        let end = start.saturating_add(window).min(bytes.len());
        let region = &bytes[start..end];
        let cuts = newline_cuts(region, workers);
        let counts = count_newline_chunks(region, &cuts);
        let total: usize = counts.iter().sum();
        if total >= n_rows || end == bytes.len() {
            if total < n_rows {
                return Err(ParseError::IncompleteFrame);
            }
            return chunks_from_counts(bytes, start, &cuts, &counts, n_rows);
        }
        window = window.saturating_mul(2);
    }
}

#[cfg(feature = "parallel")]
fn chunks_from_counts(
    bytes: &[u8],
    start: usize,
    cuts: &[usize],
    counts: &[usize],
    n_rows: usize,
) -> Result<(Vec<(usize, usize)>, usize), ParseError> {
    let mut need = n_rows;
    let mut chunks = Vec::new();
    let mut end_abs = start;
    for index in 0..counts.len() {
        if need == 0 {
            break;
        }
        let have = counts[index];
        if have == 0 {
            continue;
        }
        let from = start + cuts[index];
        let to = start + cuts[index + 1];
        let take = have.min(need);
        let chunk_end = if take == have && to > from && bytes[to - 1] == b'\n' {
            to
        } else {
            let mut seen = 0usize;
            let mut found = None;
            for rel in memchr::memchr_iter(b'\n', &bytes[from..to]) {
                seen += 1;
                if seen == take {
                    found = Some(from + rel + 1);
                    break;
                }
            }
            found.ok_or(ParseError::IncompleteFrame)?
        };
        chunks.push((from, take));
        end_abs = chunk_end;
        need -= take;
    }
    if need != 0 {
        return Err(ParseError::IncompleteFrame);
    }
    Ok((chunks, end_abs))
}

/// Parse coordinate rows described by [`plan_coord_chunks`].
///
/// Each chunk is a contiguous run of lines. The worker walks that run
/// and writes `nrows` coordinates at the prefix-summed row.
#[cfg(feature = "parallel")]
pub(crate) fn fill_rows_from_chunks(
    bytes: &[u8],
    chunks: &[(usize, usize)],
    spec_version: u32,
    start_id: u64,
    workers: usize,
    xyz: *mut f64,
    ids: *mut u64,
    flags: *mut u8,
) -> Result<(), ParseError> {
    if chunks.is_empty() {
        return Ok(());
    }
    let mut row0 = Vec::with_capacity(chunks.len());
    let mut acc = 0usize;
    for (_, nrows) in chunks {
        row0.push(acc);
        acc += *nrows;
    }
    let xyz_ptr = xyz as usize;
    let ids_ptr = ids as usize;
    let flags_ptr = flags as usize;
    let fill = || {
        use rayon::prelude::*;
        chunks
            .par_iter()
            .zip(row0.par_iter().copied())
            .try_for_each(|(&(start, nrows), mut row)| {
                let mut pos = start;
                for _ in 0..nrows {
                    let ((xyz, fixed, atom_id), next) =
                        parse_coord_row_at(bytes, pos, spec_version, start_id + row as u64)?;
                    unsafe {
                        std::ptr::copy_nonoverlapping(
                            xyz.as_ptr(),
                            (xyz_ptr as *mut f64).add(row * 3),
                            3,
                        );
                        if ids_ptr != 0 {
                            *((ids_ptr as *mut u64).add(row)) = atom_id;
                        }
                        if flags_ptr != 0 {
                            let mask = u8::from(fixed[0])
                                | (u8::from(fixed[1]) << 1)
                                | (u8::from(fixed[2]) << 2);
                            *((flags_ptr as *mut u8).add(row)) = mask;
                        }
                    }
                    pos = next;
                    row += 1;
                }
                Ok(())
            })
    };
    if workers < 2 {
        let mut row_at = 0usize;
        for &(start, nrows) in chunks {
            let mut pos = start;
            for _ in 0..nrows {
                let ((row_xyz, fixed, atom_id), next) =
                    parse_coord_row_at(bytes, pos, spec_version, start_id + row_at as u64)?;
                unsafe {
                    std::ptr::copy_nonoverlapping(row_xyz.as_ptr(), xyz.add(row_at * 3), 3);
                    if ids_ptr != 0 {
                        *((ids_ptr as *mut u64).add(row_at)) = atom_id;
                    }
                    if flags_ptr != 0 {
                        let mask = u8::from(fixed[0])
                            | (u8::from(fixed[1]) << 1)
                            | (u8::from(fixed[2]) << 2);
                        *((flags_ptr as *mut u8).add(row_at)) = mask;
                    }
                }
                pos = next;
                row_at += 1;
            }
        }
        Ok(())
    } else if chunks.len() < 2 {
        fill()
    } else {
        on_pool(workers, fill)
    }
}

/// Exact index: count, prefix sum, then store into an unzeroed buffer.
#[cfg(feature = "parallel")]
fn newline_index_exact(cuts: &[usize], bytes: &[u8], base: usize) -> Vec<usize> {
    let nchunks = cuts.len() - 1;
    let counts = count_newline_chunks(bytes, cuts);
    let mut slots = vec![0usize; nchunks];
    let mut total = 0usize;
    for index in 0..nchunks {
        slots[index] = total;
        total += counts[index];
    }
    let mut out = Vec::<usize>::with_capacity(total);
    let ptr = out.as_mut_ptr() as usize;
    let store = || {
        use rayon::prelude::*;
        (0..nchunks).into_par_iter().for_each(|index| {
            let from = cuts[index];
            let mut at = slots[index];
            for rel in memchr::memchr_iter(b'\n', &bytes[from..cuts[index + 1]]) {
                unsafe {
                    *((ptr as *mut usize).add(at)) = base + from + rel;
                }
                at += 1;
            }
        });
    };
    if nchunks < 2 {
        store();
    } else {
        on_pool(nchunks.max(1), store);
    }
    // SAFETY: the store wrote every index in `0..total`.
    unsafe { out.set_len(total) };
    out
}

#[cfg(feature = "parallel")]
#[allow(dead_code)]
fn write_coord_group(
    bytes: &[u8],
    group: &[(usize, usize)],
    row0: usize,
    spec_version: u32,
    start_id: u64,
    symbol: Option<&std::sync::Arc<str>>,
    xyz_ptr: usize,
    ids_ptr: usize,
    flags_ptr: usize,
    atoms_ptr: usize,
) -> Result<(), ParseError> {
    let mut shares = symbol.map(|symbol| symbol_shares(symbol, group.len()));
    for (offset, &(start, end)) in group.iter().enumerate() {
        let row_index = row0 + offset;
        let default_id = start_id + row_index as u64;
        let ((row, fixed, atom_id), _next) =
            parse_coord_row_at(&bytes[start..end], 0, spec_version, default_id)?;
        unsafe {
            std::ptr::copy_nonoverlapping(
                row.as_ptr(),
                (xyz_ptr as *mut f64).add(row_index * 3),
                3,
            );
            *((ids_ptr as *mut u64).add(row_index)) = atom_id;
            *((flags_ptr as *mut [bool; 3]).add(row_index)) = fixed;
            if let Some(shares) = shares.as_mut() {
                let slot = (xyz_ptr as *mut f64).add(row_index * 3);
                let atom = crate::types::AtomDatum::at_slot(shares.take(), slot, fixed, atom_id);
                (*(atoms_ptr as *mut std::mem::MaybeUninit<crate::types::AtomDatum>)
                    .add(row_index))
                .write(atom);
            }
        }
    }
    Ok(())
}

pub trait LineStream<'a> {
    fn next_line(&mut self) -> Option<&'a str>;
    fn peek_line(&mut self) -> Option<&'a str>;
    /// Take `n` coordinate rows starting at the cursor, or `None` to keep
    /// the one-row loop. The cursor moves only when this returns `Some`.
    fn try_take_coord_block(
        &mut self,
        _n: usize,
        _spec_version: u32,
        _start_id: u64,
        _symbol: &std::sync::Arc<str>,
        _workers: usize,
    ) -> Option<CoordRows> {
        None
    }
    /// True when every remaining line, including the current peek, is
    /// empty or whitespace. Default is false so a one-line peek cannot
    /// hide a following section.
    fn remainder_is_whitespace(&mut self) -> bool {
        false
    }
    /// Called while a blank line is peeked. True when that blank begins a
    /// section: the next line is a symbol and the line after it is a
    /// component label. A blank that does not is line 1 of the next frame.
    fn blank_line_opens_section(&mut self) -> bool {
        false
    }

    /// With a blank line peeked, true when [`Self::blank_line_opens_section`]
    /// is true. A stream that cannot look ahead keeps the legacy rule and
    /// returns true.
    fn blank_opens_velocity_block(&mut self) -> bool {
        true
    }

    /// One coordinate row, fused with the line scan when the stream can do it.
    ///
    /// The default pulls a `&str` line. [`crate::iterators::MemchrLines`]
    /// parses the row directly from the byte cursor so tokenize and numeric
    /// parse are one pass.
    fn next_coord_row(
        &mut self,
        spec_version: u32,
        default_id: u64,
    ) -> Result<([f64; 3], [bool; 3], u64), ParseError> {
        let line = self.next_line().ok_or(ParseError::IncompleteFrame)?;
        parse_coord_row(line, spec_version, default_id)
    }

    /// One vector-section row (xyz, flag, required id), fused with the line
    /// scan when the stream can do it. End of input is [`ParseError::IncompleteFrame`];
    /// callers that are inside a named section map that to `IncompleteSection`.
    fn next_vector_row(
        &mut self,
        spec_version: u32,
    ) -> Result<([f64; 3], [bool; 3], u64), ParseError> {
        let line = self.next_line().ok_or(ParseError::IncompleteFrame)?;
        parse_vector_section_row(line, spec_version)
    }
}

fn take_vector_row<'a>(
    lines: &mut impl LineStream<'a>,
    spec_version: u32,
    section: &str,
) -> Result<([f64; 3], [bool; 3], u64), ParseError> {
    match lines.next_vector_row(spec_version) {
        Ok(row) => Ok(row),
        Err(ParseError::IncompleteFrame) => Err(ParseError::IncompleteSection(section.to_string())),
        Err(err) => Err(err),
    }
}

impl<'a, I> LineStream<'a> for Peekable<I>
where
    I: Iterator<Item = &'a str> + Clone,
{
    #[inline]
    fn next_line(&mut self) -> Option<&'a str> {
        Iterator::next(self)
    }
    #[inline]
    fn peek_line(&mut self) -> Option<&'a str> {
        Peekable::peek(self).copied()
    }
    fn blank_line_opens_section(&mut self) -> bool {
        let mut probe = self.clone();
        let Some(blank) = probe.next_line() else {
            return false;
        };
        if !blank.trim().is_empty() {
            return false;
        }
        section_follows(&mut probe)
    }

    fn blank_opens_velocity_block(&mut self) -> bool {
        match self.peek_line() {
            Some(line) if line.trim().is_empty() => self.blank_line_opens_section(),
            _ => true,
        }
    }
}

/// True when `line` has the words `of` and `Component` in that order.
pub(crate) fn line_has_component_label(line: &str) -> bool {
    let mut words = line.split_ascii_whitespace();
    while let Some(word) = words.next() {
        if word.eq_ignore_ascii_case("of")
            && words
                .next()
                .is_some_and(|next| next.eq_ignore_ascii_case("component"))
        {
            return true;
        }
    }
    false
}

fn section_follows<'a>(lines: &mut impl LineStream<'a>) -> bool {
    let Some(symbol) = lines.next_line() else {
        return false;
    };
    if symbol.trim().is_empty() || symbol.trim_start().starts_with('{') {
        return false;
    }
    let Some(label) = lines.next_line() else {
        return false;
    };
    line_has_component_label(label)
}

/// Length of the stack buffer of [`parse_line_of_range_f64_stack`]: the most
/// columns an atom line has, which are x, y, z, the flag and the optional id
/// (spec "Coordinates").
pub const NUMERIC_STACK_WIDTH: usize = COORD_ROW_MIN_COLUMNS + 1;

/// Hot-path: parse up to [`NUMERIC_STACK_WIDTH`] whitespace-separated f64s
/// into a stack buffer.
/// Returns count of tokens actually present (before padding).
/// Pads `out[found..max]` from `defaults` when `found < max` and `found >= min`.
///
/// Single-pass over the line bytes: skip ASCII whitespace, then
/// [`fast_float2::parse_partial`] (Eisel–Lemire / SIMD-class decimal kernel)
/// with a token-boundary check. No `SplitWhitespace`, no per-token `&str`,
/// no heap `Vec`. Prefer this over allocating [`parse_line_of_n_f64`] on atom lines.
#[inline]
pub fn parse_line_of_range_f64_stack(
    line: &str,
    min: usize,
    max: usize,
    defaults: &[f64],
    out: &mut [f64; NUMERIC_STACK_WIDTH],
) -> Result<usize, ParseError> {
    debug_assert!(max <= NUMERIC_STACK_WIDTH && min <= max && defaults.len() >= max);
    let bytes = line.as_bytes();
    let n = bytes.len();
    let mut i = 0usize;
    let mut found = 0usize;
    while found < max {
        while i < n && bytes[i].is_ascii_whitespace() {
            i += 1;
        }
        if i >= n {
            break;
        }
        let end = bytes[i..]
            .iter()
            .position(|b| b.is_ascii_whitespace())
            .map(|k| i + k)
            .unwrap_or(n);
        let token = std::str::from_utf8(&bytes[i..end]).unwrap_or("");
        if is_nonfinite_spelling(token) {
            return Err(ParseError::NonFiniteValue(token.to_string()));
        }
        let (val, consumed) = fast_float2::parse_partial::<f64, _>(&bytes[i..])
            .map_err(|_| ParseError::InvalidNumberFormat(format!("invalid float: {token}")))?;
        let next = i + consumed;
        // Reject partial tokens like "1.2abc" (must end at whitespace or EOS).
        if next != end {
            return Err(ParseError::InvalidNumberFormat(format!(
                "invalid float: {token}"
            )));
        }
        if !val.is_finite() {
            return Err(ParseError::NonFiniteValue(token.to_string()));
        }
        out[found] = val;
        found += 1;
        i = next;
    }
    while i < n && bytes[i].is_ascii_whitespace() {
        i += 1;
    }
    // Reject trailing extra tokens (same contract as the old split loop).
    if i < n {
        return Err(ParseError::InvalidVectorLength {
            expected: max,
            found: found + 1,
        });
    }
    if found < min || found > max {
        return Err(ParseError::InvalidVectorLength {
            expected: max,
            found,
        });
    }
    while found < max {
        out[found] = defaults[found];
        found += 1;
    }
    Ok(found)
}

/// Parses a line of whitespace-separated f64 values using fast-float2.
///
/// This is the hot-path parser for coordinate and velocity lines. It uses
/// `fast_float2::parse` instead of `str::parse::<f64>()` for better throughput
/// on the numeric-heavy atom data lines. Fixed-width atom lines use
/// [`parse_line_of_range_f64_stack`] to avoid a heap `Vec` per line.
///
/// # Arguments
///
/// * `line` - A string slice representing a single line of data.
/// * `n` - The exact number of f64 values expected on the line.
pub fn parse_line_of_n_f64(line: &str, n: usize) -> Result<Vec<f64>, ParseError> {
    if n <= NUMERIC_STACK_WIDTH {
        let defaults = [0.0f64; NUMERIC_STACK_WIDTH];
        let mut buf = [0.0f64; NUMERIC_STACK_WIDTH];
        parse_line_of_range_f64_stack(line, n, n, &defaults, &mut buf)?;
        return Ok(buf[..n].to_vec());
    }
    let mut values = Vec::with_capacity(n.min(max_tokens_on_line(line)));
    for token in line.split_ascii_whitespace() {
        values.push(parse_finite_float_token(token)?);
    }
    if values.len() == n {
        Ok(values)
    } else {
        Err(ParseError::InvalidVectorLength {
            expected: n,
            found: values.len(),
        })
    }
}

/// Parses a line of whitespace-separated f64 values, accepting between `min`
/// and `max` values (inclusive). Returns a vector of exactly `max` elements,
/// padding with values from `defaults` when fewer than `max` are present.
///
/// Used for atom lines where column 5 (atom_index) is optional.
/// Prefer [`parse_line_of_range_f64_stack`] on the atom hot path
/// (`max <= NUMERIC_STACK_WIDTH`).
pub fn parse_line_of_range_f64(
    line: &str,
    min: usize,
    max: usize,
    defaults: &[f64],
) -> Result<Vec<f64>, ParseError> {
    if max <= NUMERIC_STACK_WIDTH {
        let mut buf = [0.0f64; NUMERIC_STACK_WIDTH];
        parse_line_of_range_f64_stack(line, min, max, defaults, &mut buf)?;
        return Ok(buf[..max].to_vec());
    }
    let mut values = Vec::with_capacity(max);
    for token in line.split_ascii_whitespace() {
        values.push(parse_finite_float_token(token)?);
    }
    if values.len() < min || values.len() > max {
        return Err(ParseError::InvalidVectorLength {
            expected: max,
            found: values.len(),
        });
    }
    while values.len() < max {
        let idx = values.len();
        values.push(defaults[idx]);
    }
    Ok(values)
}

/// CON-num.1 spellings. A leading sign is optional. Case is ignored.
/// Digit and dot tokens return before any letter compare so the atom
/// hot path does not allocate.
fn is_nonfinite_spelling(token: &str) -> bool {
    let rest = token
        .strip_prefix('+')
        .or_else(|| token.strip_prefix('-'))
        .unwrap_or(token);
    match rest.as_bytes().first() {
        Some(b) if b.is_ascii_digit() || *b == b'.' => return false,
        None => return false,
        _ => {}
    }
    rest.eq_ignore_ascii_case("nan")
        || rest.eq_ignore_ascii_case("inf")
        || rest.eq_ignore_ascii_case("infinity")
}

fn parse_finite_float_token(token: &str) -> Result<f64, ParseError> {
    if is_nonfinite_spelling(token) {
        return Err(ParseError::NonFiniteValue(token.to_string()));
    }
    let (val, consumed) = fast_float2::parse_partial::<f64, _>(token.as_bytes())
        .map_err(|_| ParseError::InvalidNumberFormat(format!("invalid float: {token}")))?;
    if consumed != token.len() {
        return Err(ParseError::InvalidNumberFormat(format!(
            "invalid float: {token}"
        )));
    }
    if !val.is_finite() {
        return Err(ParseError::NonFiniteValue(token.to_string()));
    }
    Ok(val)
}

fn skip_ascii_ws(bytes: &[u8], mut i: usize) -> usize {
    while i < bytes.len() && bytes[i].is_ascii_whitespace() {
        i += 1;
    }
    i
}

fn take_token(line: &str, i: usize) -> Option<(&str, usize)> {
    let bytes = line.as_bytes();
    if i >= bytes.len() {
        return None;
    }
    let mut j = i;
    while j < bytes.len() && !bytes[j].is_ascii_whitespace() {
        j += 1;
    }
    Some((&line[i..j], j))
}

/// Column 4. A sign, a decimal point, or a value outside the spec's set
/// is `InvalidFixedFlag` in every mode (CON-mask.1, CON-mask.5, D8).
fn parse_fixed_flag(token: &str, spec_version: u32) -> Result<[bool; 3], ParseError> {
    let bad = || {
        ParseError::InvalidFixedFlag(format!(
            "column 4 flag {token:?} is not a legal integer for spec {spec_version}"
        ))
    };
    if token.is_empty() || !token.bytes().all(|b| b.is_ascii_digit()) {
        return Err(bad());
    }
    let val: u64 = token.parse::<u64>().map_err(|_| bad())?;
    if spec_version < 2 {
        return match val {
            0 => Ok([false, false, false]),
            1 => Ok([true, true, true]),
            _ => Err(bad()),
        };
    }
    if val > 7 {
        return Err(bad());
    }
    let v = val as u8;
    Ok([v & 1 != 0, v & 2 != 0, v & 4 != 0])
}

/// Column 5. Digits only, parsed as `u64`, never through `f64` (CON-id.1, D9).
/// A sign or a decimal point is `InvalidNumberFormat` at the token level.
fn parse_atom_id_token(token: &str) -> Result<u64, ParseError> {
    if token.is_empty() || !token.bytes().all(|b| b.is_ascii_digit()) {
        return Err(ParseError::InvalidNumberFormat(format!(
            "invalid integer: {token}"
        )));
    }
    token.parse::<u64>().map_err(|_| {
        ParseError::InvalidNumberFormat(format!("integer out of range 0..=2^64-1: {token}"))
    })
}

/// Coordinate row: three finite floats, an integer flag, and an optional id.
/// Four columns take `default_id`, the 0-based position in the frame
/// (CON-id.9). Any other width is `InvalidVectorLength` (CON-num.9).
fn parse_coord_row(
    line: &str,
    spec_version: u32,
    default_id: u64,
) -> Result<([f64; 3], [bool; 3], u64), ParseError> {
    parse_vector_row(line, spec_version, Some(default_id))
}

/// Vector section row. The flag and id columns are required (exact width 5).
fn parse_vector_section_row(
    line: &str,
    spec_version: u32,
) -> Result<([f64; 3], [bool; 3], u64), ParseError> {
    parse_vector_row(line, spec_version, None)
}

fn parse_vector_row(
    line: &str,
    spec_version: u32,
    default_id: Option<u64>,
) -> Result<([f64; 3], [bool; 3], u64), ParseError> {
    let (row, end) = parse_coord_fields(line.as_bytes(), 0, spec_version, default_id)?;
    if end != line.len() {
        let found = 5 + count_tokens_on_line(line.as_bytes(), end);
        return Err(ParseError::InvalidVectorLength { expected: 5, found });
    }
    Ok(row)
}

/// Field separator inside a coordinate row. Newline ends the row.
#[inline(always)]
fn is_field_ws(b: u8) -> bool {
    matches!(b, b' ' | b'\t' | b'\r' | 0x0c)
}

#[inline(always)]
fn skip_field_ws(bytes: &[u8], mut i: usize) -> usize {
    while i < bytes.len() && is_field_ws(bytes[i]) {
        i += 1;
    }
    i
}

fn token_on_line(bytes: &[u8], i: usize) -> String {
    let mut j = i;
    while j < bytes.len() && !is_field_ws(bytes[j]) && bytes[j] != b'\n' {
        j += 1;
    }
    String::from_utf8_lossy(&bytes[i..j]).into_owned()
}

fn count_tokens_on_line(bytes: &[u8], mut i: usize) -> usize {
    let mut n = 0usize;
    loop {
        i = skip_field_ws(bytes, i);
        if i >= bytes.len() || bytes[i] == b'\n' {
            break;
        }
        n += 1;
        while i < bytes.len() && !is_field_ws(bytes[i]) && bytes[i] != b'\n' {
            i += 1;
        }
    }
    n
}

/// Invalid UTF-8 in a numeric row, otherwise a normal number-format error.
fn utf8_or_number(bytes: &[u8], line_start: usize, i: usize) -> ParseError {
    let end = bytes[line_start..]
        .iter()
        .position(|b| *b == b'\n')
        .map(|p| line_start + p)
        .unwrap_or(bytes.len());
    if std::str::from_utf8(&bytes[line_start..end]).is_err() {
        return ParseError::ValidationError("invalid UTF-8".into());
    }
    ParseError::InvalidNumberFormat(format!("invalid float: {}", token_on_line(bytes, i)))
}

#[inline(always)]
fn parse_f64_at(bytes: &[u8], i: usize, line_start: usize) -> Result<(f64, usize), ParseError> {
    if i >= bytes.len() || bytes[i] == b'\n' {
        return Err(ParseError::InvalidVectorLength {
            expected: 5,
            found: 0,
        });
    }
    if bytes[i] >= 0x80 {
        return Err(utf8_or_number(bytes, line_start, i));
    }
    if looks_like_nonfinite(bytes, i) {
        let token = token_on_line(bytes, i);
        if is_nonfinite_spelling(&token) {
            return Err(ParseError::NonFiniteValue(token));
        }
    }
    let (val, consumed) = fast_float2::parse_partial::<f64, _>(&bytes[i..]).map_err(|_| {
        ParseError::InvalidNumberFormat(format!("invalid float: {}", token_on_line(bytes, i)))
    })?;
    let next = i + consumed;
    let boundary = next >= bytes.len() || is_field_ws(bytes[next]) || bytes[next] == b'\n';
    if !boundary || !val.is_finite() {
        let token = token_on_line(bytes, i);
        if is_nonfinite_spelling(&token) || !val.is_finite() {
            return Err(ParseError::NonFiniteValue(token));
        }
        return Err(ParseError::InvalidNumberFormat(format!(
            "invalid float: {token}"
        )));
    }
    Ok((val, next))
}

#[inline(always)]
fn parse_flag_at(
    bytes: &[u8],
    i: usize,
    spec_version: u32,
) -> Result<([bool; 3], usize), ParseError> {
    if i >= bytes.len() || bytes[i] == b'\n' {
        return Err(ParseError::InvalidVectorLength {
            expected: 5,
            found: 3,
        });
    }
    let bad = || {
        ParseError::InvalidFixedFlag(format!(
            "column 4 flag {:?} is not a legal integer for spec {spec_version}",
            token_on_line(bytes, i)
        ))
    };
    if !bytes[i].is_ascii_digit() {
        return Err(bad());
    }
    let mut j = i;
    let mut val: u64 = 0;
    while j < bytes.len() && bytes[j].is_ascii_digit() {
        let d = u64::from(bytes[j] - b'0');
        val = val
            .checked_mul(10)
            .and_then(|v| v.checked_add(d))
            .ok_or_else(bad)?;
        j += 1;
    }
    if j < bytes.len() && !is_field_ws(bytes[j]) && bytes[j] != b'\n' {
        return Err(bad());
    }
    let fixed = if spec_version < 2 {
        match val {
            0 => [false, false, false],
            1 => [true, true, true],
            _ => return Err(bad()),
        }
    } else if val <= 7 {
        let v = val as u8;
        [v & 1 != 0, v & 2 != 0, v & 4 != 0]
    } else {
        return Err(bad());
    };
    Ok((fixed, j))
}

#[inline(always)]
fn parse_id_at(bytes: &[u8], i: usize) -> Result<(u64, usize), ParseError> {
    if i >= bytes.len() || bytes[i] == b'\n' || !bytes[i].is_ascii_digit() {
        return Err(ParseError::InvalidNumberFormat(format!(
            "invalid integer: {}",
            token_on_line(bytes, i)
        )));
    }
    let mut j = i;
    let mut val: u64 = 0;
    while j < bytes.len() && bytes[j].is_ascii_digit() {
        let d = u64::from(bytes[j] - b'0');
        val = val
            .checked_mul(10)
            .and_then(|v| v.checked_add(d))
            .ok_or_else(|| {
                ParseError::InvalidNumberFormat(format!(
                    "integer out of range 0..=2^64-1: {}",
                    token_on_line(bytes, i)
                ))
            })?;
        j += 1;
    }
    if j < bytes.len() && !is_field_ws(bytes[j]) && bytes[j] != b'\n' {
        return Err(ParseError::InvalidNumberFormat(format!(
            "invalid integer: {}",
            token_on_line(bytes, i)
        )));
    }
    Ok((val, j))
}

/// One parsed coordinate row: xyz, fixed flags, and the atom id.
type ParsedCoord = ([f64; 3], [bool; 3], u64);

/// One coordinate or vector row starting at `i`.
///
/// Stops at newline or end of `bytes` and returns the index just past the
/// row (after the newline, when one was present). Floats, the flag, and the
/// id are parsed in that single scan. The success path does not build error
/// strings; those stay on the cold fallback.
#[inline(always)]
fn parse_coord_fields(
    bytes: &[u8],
    i: usize,
    spec_version: u32,
    default_id: Option<u64>,
) -> Result<(ParsedCoord, usize), ParseError> {
    if let Some(parsed) = parse_coord_fields_fast(bytes, i, spec_version, default_id) {
        return Ok(parsed);
    }
    parse_coord_fields_detailed(bytes, i, spec_version, default_id)
}

/// Common coordinate row: three finite decimals, a small flag, an optional id.
#[inline(always)]
fn parse_coord_fields_fast(
    bytes: &[u8],
    i: usize,
    spec_version: u32,
    default_id: Option<u64>,
) -> Option<(ParsedCoord, usize)> {
    let mut i = skip_field_ws(bytes, i);
    let mut xyz = [0.0f64; 3];
    for dest in &mut xyz {
        if i >= bytes.len() || bytes[i] == b'\n' || bytes[i] >= 0x80 {
            return None;
        }
        let (val, next) = crate::fast_dec::parse_f64_hot(bytes, i)?;
        if next < bytes.len() && !is_field_ws(bytes[next]) && bytes[next] != b'\n' {
            return None;
        }
        *dest = val;
        i = skip_field_ws(bytes, next);
    }
    let (fixed, next) = parse_flag_fast(bytes, i, spec_version)?;
    i = skip_field_ws(bytes, next);
    let (atom_id, next) = if i >= bytes.len() || bytes[i] == b'\n' {
        (default_id?, i)
    } else {
        parse_u64_fast(bytes, i)?
    };
    i = skip_field_ws(bytes, next);
    if i < bytes.len() && bytes[i] != b'\n' {
        return None;
    }
    if i < bytes.len() && bytes[i] == b'\n' {
        i += 1;
    }
    Some(((xyz, fixed, atom_id), i))
}

#[inline(always)]
fn parse_u64_fast(bytes: &[u8], i: usize) -> Option<(u64, usize)> {
    if i >= bytes.len() || !bytes[i].is_ascii_digit() {
        return None;
    }
    let mut j = i;
    let mut val: u64 = 0;
    let mut digits = 0u32;
    while j < bytes.len() && bytes[j].is_ascii_digit() {
        digits += 1;
        // 19 digits can overflow u64; the detailed parser reports that.
        if digits > 18 {
            return None;
        }
        val = val * 10 + u64::from(bytes[j] - b'0');
        j += 1;
    }
    if j < bytes.len() && !is_field_ws(bytes[j]) && bytes[j] != b'\n' {
        return None;
    }
    Some((val, j))
}

#[inline(always)]
fn parse_flag_fast(bytes: &[u8], i: usize, spec_version: u32) -> Option<([bool; 3], usize)> {
    let (val, j) = parse_u64_fast(bytes, i)?;
    let fixed = if spec_version < 2 {
        match val {
            0 => [false, false, false],
            1 => [true, true, true],
            _ => return None,
        }
    } else if val <= 7 {
        let v = val as u8;
        [v & 1 != 0, v & 2 != 0, v & 4 != 0]
    } else {
        return None;
    };
    Some((fixed, j))
}

#[cold]
#[inline(never)]
fn parse_coord_fields_detailed(
    bytes: &[u8],
    i: usize,
    spec_version: u32,
    default_id: Option<u64>,
) -> Result<(ParsedCoord, usize), ParseError> {
    let line_start = i;
    let mut i = skip_field_ws(bytes, i);
    if i >= bytes.len() || bytes[i] == b'\n' {
        return Err(ParseError::InvalidVectorLength {
            expected: 5,
            found: 0,
        });
    }
    let mut xyz = [0.0f64; 3];
    for (slot, dest) in xyz.iter_mut().enumerate() {
        if i >= bytes.len() || bytes[i] == b'\n' {
            return Err(ParseError::InvalidVectorLength {
                expected: 5,
                found: slot,
            });
        }
        if bytes[i] >= 0x80 {
            return Err(utf8_or_number(bytes, line_start, i));
        }
        let (val, next) = parse_f64_at(bytes, i, line_start)?;
        *dest = val;
        i = skip_field_ws(bytes, next);
    }
    let (fixed, next) = parse_flag_at(bytes, i, spec_version)?;
    i = skip_field_ws(bytes, next);
    let (atom_id, next) = if i >= bytes.len() || bytes[i] == b'\n' {
        (
            default_id.ok_or(ParseError::InvalidVectorLength {
                expected: 5,
                found: 4,
            })?,
            i,
        )
    } else {
        parse_id_at(bytes, i)?
    };
    i = skip_field_ws(bytes, next);
    if i < bytes.len() && bytes[i] != b'\n' {
        let found = 5 + count_tokens_on_line(bytes, i);
        return Err(ParseError::InvalidVectorLength { expected: 5, found });
    }
    if i < bytes.len() && bytes[i] == b'\n' {
        i += 1;
    }
    Ok(((xyz, fixed, atom_id), i))
}

/// Coordinate row at byte offset `i` in a multi-line buffer.
///
/// The returned index is the start of the next line. ASCII rows are valid
/// UTF-8 without a separate scan; a non-ASCII byte is checked only then.
pub(crate) fn parse_coord_row_at(
    bytes: &[u8],
    i: usize,
    spec_version: u32,
    default_id: u64,
) -> Result<(ParsedCoord, usize), ParseError> {
    if i >= bytes.len() {
        return Err(ParseError::IncompleteFrame);
    }
    parse_coord_fields(bytes, i, spec_version, Some(default_id))
}

/// Vector-section row at `i`. The id column is required.
pub(crate) fn parse_vector_row_at(
    bytes: &[u8],
    i: usize,
    spec_version: u32,
) -> Result<(ParsedCoord, usize), ParseError> {
    if i >= bytes.len() {
        return Err(ParseError::IncompleteFrame);
    }
    parse_coord_fields(bytes, i, spec_version, None)
}

/// Scalar section row: one finite float, a flag, and an id. Exact width 3.
fn parse_scalar_section_row(
    line: &str,
    spec_version: u32,
) -> Result<(f64, [bool; 3], u64), ParseError> {
    let bytes = line.as_bytes();
    let i = skip_ascii_ws(bytes, 0);
    let (value, i) = parse_finite_at(line, i)?;
    let i = skip_ascii_ws(bytes, i);
    let (flag_tok, i) = take_token(line, i).ok_or(ParseError::InvalidVectorLength {
        expected: 3,
        found: 1,
    })?;
    let fixed = parse_fixed_flag(flag_tok, spec_version)?;
    let i = skip_ascii_ws(bytes, i);
    let (id_tok, i) = take_token(line, i).ok_or(ParseError::InvalidVectorLength {
        expected: 3,
        found: 2,
    })?;
    let atom_id = parse_atom_id_token(id_tok)?;
    let i = skip_ascii_ws(bytes, i);
    if i < bytes.len() {
        let found = 3 + line[i..].split_ascii_whitespace().count();
        return Err(ParseError::InvalidVectorLength { expected: 3, found });
    }
    Ok((value, fixed, atom_id))
}

fn parse_finite_at(line: &str, i: usize) -> Result<(f64, usize), ParseError> {
    let bytes = line.as_bytes();
    if i >= bytes.len() {
        return Err(ParseError::InvalidVectorLength {
            expected: 3,
            found: 0,
        });
    }
    // "infinity" parses as the prefix "inf" in fast_float2. Name that
    // spelling before the numeric kernel so the leftover is not a format error.
    if looks_like_nonfinite(bytes, i) {
        let token = take_token(line, i).map(|(tok, _)| tok).unwrap_or("");
        if is_nonfinite_spelling(token) {
            return Err(ParseError::NonFiniteValue(token.to_string()));
        }
    }
    let (val, consumed) = fast_float2::parse_partial::<f64, _>(&bytes[i..]).map_err(|_| {
        let token = take_token(line, i).map(|(tok, _)| tok).unwrap_or("");
        ParseError::InvalidNumberFormat(format!("invalid float: {token}"))
    })?;
    let next = i + consumed;
    let boundary = next == bytes.len() || bytes[next].is_ascii_whitespace();
    if !boundary || !val.is_finite() {
        let token = take_token(line, i).map(|(tok, _)| tok).unwrap_or("");
        if is_nonfinite_spelling(token) || !val.is_finite() {
            return Err(ParseError::NonFiniteValue(token.to_string()));
        }
        return Err(ParseError::InvalidNumberFormat(format!(
            "invalid float: {token}"
        )));
    }
    Ok((val, next))
}

/// True when the token may be a non-finite spelling. Digit and dot
/// tokens, the common atom-line case, are false after one byte.
fn looks_like_nonfinite(bytes: &[u8], i: usize) -> bool {
    let mut j = i;
    if bytes[j] == b'+' || bytes[j] == b'-' {
        j += 1;
        if j >= bytes.len() {
            return false;
        }
    }
    let b = bytes[j];
    !b.is_ascii_digit() && b != b'.'
}

/// Case-insensitive, after collapsing each run of spaces (CON-file.5, D17).
fn labels_equal(found: &str, expected: &str) -> bool {
    let mut found_words = found.split_ascii_whitespace();
    let mut expected_words = expected.split_ascii_whitespace();
    loop {
        match (found_words.next(), expected_words.next()) {
            (Some(a), Some(b)) => {
                if !a.eq_ignore_ascii_case(b) {
                    return false;
                }
            }
            (None, None) => return true,
            _ => return false,
        }
    }
}

/// True when `line`, after the same folding, starts with every word of `prefix`.
fn label_has_prefix(line: &str, prefix: &str) -> bool {
    let mut found_words = line.split_ascii_whitespace();
    for word in prefix.split_ascii_whitespace() {
        match found_words.next() {
            Some(got) if got.eq_ignore_ascii_case(word) => {}
            _ => return false,
        }
    }
    true
}

fn token_count(line: &str) -> usize {
    line.split_ascii_whitespace().count()
}

/// Most whitespace-separated tokens `line` can hold.
///
/// `k` tokens take at least one byte each and one separator between each
/// pair, `2k - 1` bytes in all.
fn max_tokens_on_line(line: &str) -> usize {
    line.len().div_ceil(2)
}

fn json_u64(value: &Value) -> Option<u64> {
    value.as_u64()
}

fn validate_section_arity_value(value: &Value) -> Result<(), ParseError> {
    let obj = value.as_object().ok_or_else(|| {
        metadata_json_error("section_arity must be an object mapping section names to integers")
    })?;
    for (name, entry) in obj {
        let n = json_u64(entry).ok_or_else(|| {
            metadata_json_error(format!("section_arity.{name} must be an integer >= 1"))
        })?;
        if n < 1 {
            return Err(ParseError::InvalidCount(format!(
                "section_arity.{name} must be >= 1"
            )));
        }
        // Reserved names have a fixed value-field count (CON-section.13).
        if let Some(expect) = crate::types::reserved_section_value_arity(name)
            && n != expect
        {
            return Err(ParseError::InvalidVectorLength {
                expected: usize::try_from(expect).unwrap_or(usize::MAX),
                found: usize::try_from(n).unwrap_or(usize::MAX),
            });
        }
    }
    Ok(())
}

/// Integer counts. A sign, a decimal point, or any other non-digit token
/// is `InvalidCount` (CON-error.3). Zero is a legal token here. The caller
/// rejects a per-type count of 0 when the type count is 1 or more (D4).
pub(crate) fn parse_count_line(line: &str, n: usize) -> Result<Vec<usize>, ParseError> {
    let mut values = Vec::with_capacity(n.min(max_tokens_on_line(line)));
    for token in line.split_ascii_whitespace() {
        if token.is_empty() || !token.bytes().all(|b| b.is_ascii_digit()) {
            return Err(ParseError::InvalidCount(format!(
                "count field is not an unsigned integer: {token}"
            )));
        }
        let val: u64 = token.parse::<u64>().map_err(|_| {
            ParseError::InvalidCount(format!("count field is out of range: {token}"))
        })?;
        let val = usize::try_from(val).map_err(|_| {
            ParseError::InvalidCount(format!("count field is out of range: {token}"))
        })?;
        values.push(val);
    }
    if values.len() == n {
        Ok(values)
    } else {
        Err(ParseError::InvalidVectorLength {
            expected: n,
            found: values.len(),
        })
    }
}

/// True when `line` has no tokens. Spaces and tabs count as empty.
pub(crate) fn line_is_whitespace(line: &str) -> bool {
    line.split_ascii_whitespace().next().is_none()
}

/// Lines 8 and 9 of a type count of 0.
///
/// A token on either line is `InvalidCount`. At end of file, a missing
/// line is empty (CON-count.4), and that frame is not `IncompleteHeader`
/// (CON-frame.3). Mid-file, both lines are required.
pub(crate) fn read_zero_atom_tail<'a>(
    mut next_line: impl FnMut() -> Option<&'a str>,
) -> Result<(), ParseError> {
    let bad = || {
        ParseError::InvalidCount("a zero-atom frame requires whitespace-only lines 8 and 9".into())
    };
    match next_line() {
        None => Ok(()),
        Some(line) if line_is_whitespace(line) => match next_line() {
            None => Ok(()),
            Some(line) if line_is_whitespace(line) => Ok(()),
            Some(_) => Err(bad()),
        },
        Some(_) => Err(bad()),
    }
}

fn metadata_json_error(message: impl Into<String>) -> ParseError {
    ParseError::InvalidMetadataJson(message.into())
}

fn validate_metadata_number(key: &str, value: &Value) -> Result<(), ParseError> {
    if value.as_f64().is_some() {
        Ok(())
    } else {
        Err(metadata_json_error(format!(
            "{key} must be a finite number"
        )))
    }
}

fn validate_metadata_integer(key: &str, value: &Value) -> Result<(), ParseError> {
    if value.as_u64().is_some() {
        Ok(())
    } else {
        Err(metadata_json_error(format!(
            "{key} must be a non-negative integer"
        )))
    }
}

/// Validates a parsed metadata JSON object against the spec v2 schema.
///
/// Type-checks the `validate` and `sections` keys, then runs per-key
/// schema checks for the recommended metadata keys. Used by the
/// builder/Python `set_metadata_json` paths to fail fast on malformed
/// input, and by the parser when the file requested strict validation
/// via `"validate": true`.
pub fn validate_metadata_schema(
    json_obj: &serde_json::Map<String, Value>,
) -> Result<(), ParseError> {
    let strict_requested = match json_obj.get(meta::VALIDATE) {
        Some(Value::Bool(b)) => *b,
        Some(_) => return Err(metadata_json_error("validate must be a boolean")),
        None => false,
    };

    match json_obj.get(meta::SECTIONS) {
        Some(Value::Array(values)) => {
            if values.iter().any(|entry| !entry.is_string()) {
                return Err(metadata_json_error("sections must be an array of strings"));
            }
        }
        Some(_) => return Err(metadata_json_error("sections must be an array of strings")),
        None if strict_requested => {
            return Err(metadata_json_error(
                "validate=true requires a sections array, even when empty",
            ));
        }
        None => {}
    }

    for (key, value) in json_obj {
        match key.as_str() {
            meta::CON_SPEC_VERSION | meta::SECTIONS | meta::VALIDATE => {}
            meta::CON_SPEC_MINOR => {
                if json_u64(value).is_none() {
                    return Err(metadata_json_error(
                        "con_spec_minor must be a JSON integer >= 0",
                    ));
                }
            }
            meta::SECTION_ARITY => validate_section_arity_value(value)?,
            meta::ENERGY
            | meta::TIME
            | meta::TIMESTEP
            | meta::CONVERGENCE_FMAX
            | meta::CONVERGENCE_ENERGY
            | meta::FMAX => validate_metadata_number(key, value)?,
            meta::FRAME_INDEX | meta::NEB_BEAD | meta::NEB_BAND => {
                validate_metadata_integer(key, value)?
            }
            meta::GENERATOR if !value.is_string() => {
                return Err(metadata_json_error("generator must be a string"));
            }
            meta::GENERATOR => {}
            meta::UNITS | meta::POTENTIAL if !value.is_object() => {
                return Err(metadata_json_error(format!("{key} must be an object")));
            }
            meta::POTENTIAL => {
                if let Some(potential_type) = value.get("type")
                    && !potential_type.is_string()
                {
                    return Err(metadata_json_error("potential.type must be a string"));
                }
            }
            meta::UNITS => {
                // Full v3 checks (required keys) run when con_spec_version >= 3.
            }
            meta::PBC => validate_pbc_metadata(value)?,
            meta::BONDS => validate_bonds_metadata(value)?,
            meta::LATTICE_VECTORS => validate_lattice_vectors_metadata(value)?,
            meta::CONVERGED if !value.is_boolean() => {
                return Err(metadata_json_error("converged must be a boolean"));
            }
            meta::CONVERGED => {}
            _ => {}
        }
    }

    Ok(())
}

fn validate_pbc_metadata(value: &Value) -> Result<(), ParseError> {
    let Some(values) = value.as_array() else {
        return Err(metadata_json_error("pbc must be a length-3 boolean array"));
    };
    if values.len() != 3 || values.iter().any(|entry| !entry.is_boolean()) {
        return Err(metadata_json_error("pbc must be a length-3 boolean array"));
    }
    Ok(())
}

/// Validate optional `bonds` frame topology metadata.
///
/// Each element is either a length-2 non-negative integer pair `[i, j]` or an
/// object `{"i": i, "j": j, "order"?: integer}`. Indices are 0-based into
/// `atom_data` order (parser does not yet know atom count at metadata time;
/// bounds are enforced when projecting to chemfiles / selection).
fn validate_bonds_metadata(value: &Value) -> Result<(), ParseError> {
    let Some(items) = value.as_array() else {
        return Err(metadata_json_error("bonds must be an array"));
    };
    for (idx, item) in items.iter().enumerate() {
        if let Some(pair) = item.as_array() {
            if pair.len() != 2 {
                return Err(metadata_json_error(format!(
                    "bonds[{idx}] pair must have exactly two indices"
                )));
            }
            for (k, entry) in pair.iter().enumerate() {
                let Some(n) = entry.as_u64() else {
                    return Err(metadata_json_error(format!(
                        "bonds[{idx}][{k}] must be a non-negative integer"
                    )));
                };
                if n > u32::MAX as u64 {
                    return Err(metadata_json_error(format!(
                        "bonds[{idx}][{k}] index exceeds u32"
                    )));
                }
            }
            continue;
        }
        if let Some(obj) = item.as_object() {
            for key in ["i", "j"] {
                let Some(entry) = obj.get(key) else {
                    return Err(metadata_json_error(format!(
                        "bonds[{idx}] object must include \"{key}\""
                    )));
                };
                let Some(n) = entry.as_u64() else {
                    return Err(metadata_json_error(format!(
                        "bonds[{idx}].{key} must be a non-negative integer"
                    )));
                };
                if n > u32::MAX as u64 {
                    return Err(metadata_json_error(format!(
                        "bonds[{idx}].{key} index exceeds u32"
                    )));
                }
            }
            if let Some(order) = obj.get("order")
                && order.as_i64().is_none()
            {
                return Err(metadata_json_error(format!(
                    "bonds[{idx}].order must be an integer when present"
                )));
            }
            continue;
        }
        return Err(metadata_json_error(format!(
            "bonds[{idx}] must be [i, j] or {{\"i\": i, \"j\": j, \"order\"?: ...}}"
        )));
    }
    Ok(())
}

fn validate_lattice_vectors_metadata(value: &Value) -> Result<(), ParseError> {
    let Some(rows) = value.as_array() else {
        return Err(metadata_json_error(
            "lattice_vectors must be a 3x3 numeric array",
        ));
    };
    if rows.len() != 3 {
        return Err(metadata_json_error(
            "lattice_vectors must be a 3x3 numeric array",
        ));
    }
    for row in rows {
        let Some(entries) = row.as_array() else {
            return Err(metadata_json_error(
                "lattice_vectors must be a 3x3 numeric array",
            ));
        };
        if entries.len() != 3 || entries.iter().any(|entry| entry.as_f64().is_none()) {
            return Err(metadata_json_error(
                "lattice_vectors must be a 3x3 numeric array",
            ));
        }
    }
    Ok(())
}

/// Parses a line of whitespace-separated values into a vector of a specific type.
///
/// This generic helper function takes a string slice, splits it by whitespace,
/// and attempts to parse each substring into the target type `T`. The type `T`
/// must implement `std::str::FromStr`.
///
/// # Arguments
///
/// * `line` - A string slice representing a single line of data.
/// * `n` - The exact number of values expected on the line.
///
/// # Errors
///
/// * `ParseError::InvalidVectorLength` if the number of parsed values is not equal to `n`.
/// * Propagates any error from the `parse()` method of the target type `T`.
///
/// # Example
///
/// ```
/// use readcon_core::parser::parse_line_of_n;
/// let line = "10.5 20.0 30.5";
/// let values: Vec<f64> = parse_line_of_n(line, 3).unwrap();
/// assert_eq!(values, vec![10.5, 20.0, 30.5]);
///
/// let result = parse_line_of_n::<i32>(line, 2);
/// assert!(result.is_err());
/// ```
pub fn parse_line_of_n<T: std::str::FromStr>(line: &str, n: usize) -> Result<Vec<T>, ParseError>
where
    ParseError: From<<T as std::str::FromStr>::Err>,
{
    let values: Vec<T> = line
        .split_whitespace()
        .map(|s| s.parse::<T>())
        .collect::<Result<_, _>>()?;

    if values.len() == n {
        Ok(values)
    } else {
        Err(ParseError::InvalidVectorLength {
            expected: n,
            found: values.len(),
        })
    }
}

/// Parses the 9-line header of a `.con` file frame from an iterator.
///
/// This function consumes the next 9 lines from the given line iterator to
/// construct a `FrameHeader`. The iterator is advanced by 9 lines on success.
///
/// # Arguments
///
/// * `lines` - A mutable reference to an iterator that yields string slices.
///
/// # Errors
///
/// * `ParseError::IncompleteHeader` if the iterator has fewer than 9 lines remaining.
/// * Propagates any errors from `parse_line_of_n` if the numeric data within
///   the header is malformed.
///
/// # Panics
///
/// This function will panic if the intermediate vectors for box dimensions or angles,
/// after being successfully parsed, cannot be converted into fixed-size arrays.
/// This should not happen if `parse_line_of_n` is used correctly with `n=3`.
pub fn parse_frame_header<'a>(
    lines: &mut impl Iterator<Item = &'a str>,
) -> Result<FrameHeader, ParseError> {
    let prebox1 = lines
        .next()
        .ok_or(ParseError::IncompleteHeader)?
        .to_string();
    // CON-file.2: a leading BOM is rejected. It would sit on line 1.
    if prebox1.starts_with('\u{feff}') {
        return Err(ParseError::ValidationError(
            "file begins with a byte-order mark".into(),
        ));
    }
    let prebox2_raw = lines.next().ok_or(ParseError::IncompleteHeader)?;

    // Line 2: if it starts with '{', parse as JSON metadata (spec v2+).
    // Otherwise treat as a legacy (pre-v2) file with spec_version = 1.
    // One JSON object and nothing else (CON-line2.1). serde_json::from_str
    // rejects trailing junk after the value.
    let trimmed = prebox2_raw.trim();
    let (spec_version, sections, validate, sections_declared) =
        if trimmed == "{\"con_spec_version\":2}" {
            // The line written by every stock trajectory. Same object the
            // JSON parser would build: major 2, no extra keys.
            (2, Vec::new(), false, false)
        } else if let Some(ver) = simple_metadata_version(trimmed) {
            // A header whose only keys are the spec version, an optional
            // minor, and a units object. Every byte is still read. Any
            // other key falls through to the general parser.
            (ver, Vec::new(), false, false)
        } else if trimmed.starts_with('{') {
            let json_val: serde_json::Value = serde_json::from_str(trimmed)
                .map_err(|e| ParseError::InvalidMetadataJson(e.to_string()))?;
            let json_obj = json_val.as_object().ok_or_else(|| {
                ParseError::InvalidMetadataJson("expected a JSON object".to_string())
            })?;
            let ver_value = json_obj
                .get(meta::CON_SPEC_VERSION)
                .ok_or(ParseError::MissingSpecVersion)?;
            let ver = json_u64(ver_value).ok_or_else(|| {
                ParseError::InvalidMetadataJson(
                    "con_spec_version must be a JSON integer, the major version".into(),
                )
            })?;
            if ver > u64::from(u32::MAX) {
                return Err(ParseError::UnsupportedVersion(u32::MAX));
            }
            let ver = ver as u32;
            // JSON majors 0 and 1 are UnsupportedVersion (CON-line2.5).
            // An unknown major is UnsupportedVersion (CON-line2.6).
            // Unknown minors of a known major are accepted (CON-line2.7, CON-status.1).
            if ver == 0 || ver == 1 || ver > crate::CON_SPEC_VERSION {
                return Err(ParseError::UnsupportedVersion(ver));
            }
            if ver >= 3 {
                match json_obj.get(meta::UNITS) {
                    Some(u) => crate::units::validate_v3_units_metadata(u)
                        .map_err(|e| ParseError::InvalidUnits(format!("v3 units: {e}")))?,
                    None => {
                        return Err(ParseError::InvalidUnits(
                        "con_spec_version >= 3 requires metadata \"units\" with length and energy"
                            .into(),
                    ));
                    }
                }
            }

            // Single pass over the JSON object: collect sections, capture the
            // validate flag, copy the rest into metadata. Folds the previous
            // pre-extract get(validate) + re-iterate pattern into one walk.
            let mut sections: Vec<String> = Vec::new();
            let mut sections_declared = false;
            let mut validate = false;
            for (k, v) in json_obj {
                match k.as_str() {
                    meta::CON_SPEC_VERSION => {}
                    meta::CON_SPEC_MINOR => {
                        if json_u64(v).is_none() {
                            return Err(metadata_json_error(
                                "con_spec_minor must be a JSON integer >= 0",
                            ));
                        }
                    }
                    meta::SECTION_ARITY => {
                        validate_section_arity_value(v)?;
                    }
                    meta::SECTIONS => {
                        sections_declared = true;
                        let arr = v.as_array().ok_or_else(|| {
                            metadata_json_error("sections must be an array of strings")
                        })?;
                        sections.reserve(arr.len());
                        for entry in arr {
                            let s = entry.as_str().ok_or_else(|| {
                                metadata_json_error("sections must be an array of strings")
                            })?;
                            sections.push(s.to_string());
                        }
                    }
                    meta::VALIDATE => {
                        validate = match v {
                            Value::Bool(b) => *b,
                            _ => return Err(metadata_json_error("validate must be a boolean")),
                        };
                    }
                    _ => {}
                }
            }

            // Strict-mode schema check fires only when the file requested
            // it. Hot-path parses (validate=false) skip the per-key match.
            if validate {
                validate_metadata_schema(json_obj)?;
            }

            (ver, sections, validate, sections_declared)
        } else {
            // Legacy file: no JSON metadata line.
            (1_u32, Vec::new(), false, false)
        };
    let prebox2 = prebox2_raw.to_string();

    let boxl_vec = parse_line_of_n_f64(lines.next().ok_or(ParseError::IncompleteHeader)?, 3)?;
    let angles_vec = parse_line_of_n_f64(lines.next().ok_or(ParseError::IncompleteHeader)?, 3)?;
    let postbox1 = lines
        .next()
        .ok_or(ParseError::IncompleteHeader)?
        .to_string();
    let postbox2 = lines
        .next()
        .ok_or(ParseError::IncompleteHeader)?
        .to_string();
    let natm_types_line = lines.next().ok_or(ParseError::IncompleteHeader)?;
    let natm_types = parse_count_line(natm_types_line, 1)?[0];
    // A type count of 0 is a zero-atom frame (CON-count.1, CON-count.4).
    // A per-type count of 0 is InvalidCount only when the type count is >= 1.
    let (natms_per_type, masses_per_type) = if natm_types == 0 {
        read_zero_atom_tail(|| lines.next())?;
        (Vec::new(), Vec::new())
    } else {
        let natms_per_type = parse_count_line(
            lines.next().ok_or(ParseError::IncompleteHeader)?,
            natm_types,
        )?;
        if natms_per_type.contains(&0) {
            return Err(ParseError::InvalidCount(
                "each per-type count must be >= 1".into(),
            ));
        }
        let masses_per_type = parse_line_of_n_f64(
            lines.next().ok_or(ParseError::IncompleteHeader)?,
            natm_types,
        )?;
        (natms_per_type, masses_per_type)
    };
    if validate {
        validate_header_geometry(&boxl_vec, &angles_vec)?;
        validate_masses(&masses_per_type)?;
    }
    let shared_line: std::sync::Arc<str> = std::sync::Arc::from(prebox2);
    let metadata = crate::types::FrameMeta::from_shared(std::sync::Arc::clone(&shared_line));
    Ok(FrameHeader {
        prebox_header: PreboxHeader {
            user: prebox1,
            metadata_line: shared_line,
        },
        boxl: boxl_vec.try_into().unwrap(),
        angles: angles_vec.try_into().unwrap(),
        postbox_header: [postbox1, postbox2],
        natm_types,
        natms_per_type,
        masses_per_type,
        spec_version,
        metadata,
        sections,
        strict_validation: validate,
        sections_declared,
    })
}

/// Spec line limited to `con_spec_version`, optional `con_spec_minor`, and `units`.
///
/// Returns the major version when the object is exactly that shape and, for
/// version 3 or newer, `length` and `energy` are valid units. Anything else
/// is `None` so the general JSON parser can accept or reject it.
fn simple_metadata_version(line: &str) -> Option<u32> {
    let bytes = line.as_bytes();
    let mut i = 0usize;
    let skip = |i: &mut usize| {
        while *i < bytes.len() && bytes[*i].is_ascii_whitespace() {
            *i += 1;
        }
    };
    let eat = |i: &mut usize, byte: u8| -> bool {
        skip(i);
        if *i < bytes.len() && bytes[*i] == byte {
            *i += 1;
            true
        } else {
            false
        }
    };
    let quoted = |i: &mut usize| -> Option<&str> {
        skip(i);
        if *i >= bytes.len() || bytes[*i] != b'"' {
            return None;
        }
        *i += 1;
        let start = *i;
        while *i < bytes.len() && bytes[*i] != b'"' {
            if bytes[*i] == b'\\' {
                return None;
            }
            *i += 1;
        }
        if *i >= bytes.len() {
            return None;
        }
        let text = std::str::from_utf8(&bytes[start..*i]).ok()?;
        *i += 1;
        Some(text)
    };
    let integer = |i: &mut usize| -> Option<u64> {
        skip(i);
        let start = *i;
        if *i < bytes.len() && bytes[*i] == b'-' {
            return None;
        }
        while *i < bytes.len() && bytes[*i].is_ascii_digit() {
            *i += 1;
        }
        if start == *i || (bytes[start] == b'0' && *i > start + 1) {
            return None;
        }
        std::str::from_utf8(&bytes[start..*i]).ok()?.parse().ok()
    };
    if !eat(&mut i, b'{') {
        return None;
    }
    let mut version = None;
    let mut saw_minor = false;
    let mut length = None;
    let mut energy = None;
    let mut extra_unit = false;
    loop {
        skip(&mut i);
        if i < bytes.len() && bytes[i] == b'}' {
            i += 1;
            break;
        }
        if version.is_some() || saw_minor || length.is_some() || energy.is_some() || extra_unit {
            if !eat(&mut i, b',') {
                return None;
            }
        }
        let key = quoted(&mut i)?;
        if !eat(&mut i, b':') {
            return None;
        }
        match key {
            "con_spec_version" if version.is_none() => version = Some(integer(&mut i)?),
            "con_spec_minor" if !saw_minor => {
                integer(&mut i)?;
                saw_minor = true;
            }
            "units" if length.is_none() && energy.is_none() && !extra_unit => {
                if !eat(&mut i, b'{') {
                    return None;
                }
                loop {
                    skip(&mut i);
                    if i < bytes.len() && bytes[i] == b'}' {
                        i += 1;
                        break;
                    }
                    if length.is_some() || energy.is_some() || extra_unit {
                        if !eat(&mut i, b',') {
                            return None;
                        }
                    }
                    let ukey = quoted(&mut i)?;
                    if !eat(&mut i, b':') {
                        return None;
                    }
                    let value = quoted(&mut i)?;
                    match ukey {
                        "length" if length.is_none() => length = Some(value),
                        "energy" if energy.is_none() => energy = Some(value),
                        "mass" | "time" | "velocity" | "force" => {
                            if crate::units::validate_unit_for_quantity(ukey, value).is_err() {
                                return None;
                            }
                            extra_unit = true;
                        }
                        _ => return None,
                    }
                }
            }
            _ => return None,
        }
    }
    skip(&mut i);
    if i != bytes.len() {
        return None;
    }
    let ver = u32::try_from(version?).ok()?;
    if ver == 0 || ver == 1 || ver > crate::CON_SPEC_VERSION {
        return None;
    }
    if ver >= 3 {
        let length = length?;
        let energy = energy?;
        if length.trim().is_empty() || energy.trim().is_empty() {
            return None;
        }
        crate::units::validate_unit_for_quantity("length", length).ok()?;
        crate::units::validate_unit_for_quantity("energy", energy).ok()?;
    }
    Some(ver)
}

/// Parses a complete frame from a `.con` file, including its header and atomic data.
///
/// This function first parses the complete frame header and then uses the information within it
/// (specifically the number of atom types and atoms per type) to parse the subsequent
/// atom coordinate blocks.
///
/// Each atom takes one line, so the atom vectors reserve no more rows than the
/// upper bound of `lines.size_hint()`. With no upper bound, as for `str::lines`,
/// they grow as the rows arrive.
///
/// # Arguments
///
/// * `lines` - A mutable reference to an iterator that yields string slices for the frame.
///
/// # Errors
///
/// * `ParseError::IncompleteFrame` if the iterator ends before all expected
///   atomic data has been read.
/// * Propagates any errors from the underlying calls to `parse_frame_header` and
///   `parse_line_of_n`.
///
/// # Example
///
/// ```
/// use readcon_core::parser::parse_single_frame;
///
/// let frame_text = r#"
///Generated by test
///{"con_spec_version":2}
///10.0 10.0 10.0
///90.0 90.0 90.0
///POSTBOX LINE 1
///POSTBOX LINE 2
///2
///1 1
///12.011 1.008
///C
///Coordinates of Component 1
///1.0 1.0 1.0 0 1
///H
///Coordinates of Component 2
///2.0 2.0 2.0 0 2
/// "#;
///
/// let mut lines = frame_text.trim().lines();
/// let con_frame = parse_single_frame(&mut lines).unwrap();
///
/// assert_eq!(con_frame.header.natm_types, 2);
/// assert_eq!(con_frame.atom_data.len(), 2);
/// assert_eq!(&*con_frame.atom_data[0].symbol, "C");
/// assert_eq!(con_frame.atom_data[1].atom_id, 2);
/// ```
pub fn parse_single_frame<'a>(
    lines: &mut impl Iterator<Item = &'a str>,
) -> Result<ConFrame, ParseError> {
    // The caller owns the iterator and may not be able to clone it. Section
    // detection, which needs a lookahead clone, runs later on a `LineStream`
    // that can probe (`MemchrLines`, or a `Peekable` whose source is `Clone`).
    let max_rows = lines.size_hint().1;
    parse_single_frame_ls(&mut IterStream(lines), max_rows)
}

/// Iterator adapter so [`parse_single_frame`] can satisfy [`LineStream`]
/// without requiring [`Clone`].
struct IterStream<I>(I);

impl<'a, I> Iterator for IterStream<I>
where
    I: Iterator<Item = &'a str>,
{
    type Item = &'a str;
    fn next(&mut self) -> Option<&'a str> {
        self.0.next()
    }
}

impl<'a, I> LineStream<'a> for IterStream<I>
where
    I: Iterator<Item = &'a str>,
{
    fn next_line(&mut self) -> Option<&'a str> {
        Iterator::next(self)
    }
    fn peek_line(&mut self) -> Option<&'a str> {
        None
    }
}

struct CoordShell {
    atoms: Vec<crate::types::AtomDatum>,
    coords: Vec<f64>,
    ids: Vec<u64>,
}

thread_local! {
    static COORD_SHELLS: std::cell::RefCell<Vec<CoordShell>> =
        const { std::cell::RefCell::new(Vec::new()) };
}

/// Allocate `n_frames` coordinate shells on this thread before the parse.
///
/// One tight reservation, then each frame takes a shell. Doing the
/// reservations inside the parse mixed the allocator with the float
/// parse and made eight threads slower than the bytes they write.
pub(crate) fn prestash_coord_shells(n_frames: usize, n_atoms: usize) {
    // Small frames (the 32-atom strong-scaling files) are cheaper to
    // allocate inside the parse. A burst of shells ahead of those frames
    // made eight threads slower. Wide atoms, as in the weak-scaling
    // fixture, are the case the burst is for.
    if n_frames == 0 || n_atoms < 128 {
        return;
    }
    COORD_SHELLS.with(|pool| {
        let mut pool = pool.borrow_mut();
        pool.reserve(n_frames);
        for _ in 0..n_frames {
            pool.push(CoordShell {
                atoms: Vec::with_capacity(n_atoms),
                coords: Vec::with_capacity(n_atoms.saturating_mul(3)),
                ids: Vec::with_capacity(n_atoms),
            });
        }
    });
}

pub(crate) fn clear_coord_shells() {
    COORD_SHELLS.with(|pool| pool.borrow_mut().clear());
}

fn take_coord_shell(n_atoms: usize) -> Option<CoordShell> {
    COORD_SHELLS.with(|pool| {
        let shell = pool.borrow_mut().pop()?;
        if shell.atoms.capacity() >= n_atoms
            && shell.coords.capacity() >= n_atoms.saturating_mul(3)
            && shell.ids.capacity() >= n_atoms
        {
            Some(shell)
        } else {
            None
        }
    })
}

/// Columns every atom line has: x, y, z and the flag (spec "Coordinates",
/// CON-num.9). The id column is optional (CON-id.9).
const COORD_ROW_MIN_COLUMNS: usize = 4;

/// Fewest bytes an atom line takes, as in `0 0 0 0`: one byte for each of
/// the [`COORD_ROW_MIN_COLUMNS`] columns and one separator between each
/// pair. Newlines and the symbol and label lines only add bytes, so `n`
/// unread bytes hold at most `n / MIN_COORD_ROW_BYTES` atom lines.
pub(crate) const MIN_COORD_ROW_BYTES: usize = 2 * COORD_ROW_MIN_COLUMNS - 1;

/// Frame parse over a [`LineStream`].
///
/// Coordinate rows go through [`LineStream::next_coord_row`] so the memchr
/// cursor can fuse the newline scan with numeric parse. Positions, ids, and
/// flags land in flat buffers; [`AtomDatum`] values are built from those
/// buffers afterwards (one `Vec` allocation, no per-atom heap traffic).
///
/// `max_rows` bounds the atom lines the rest of the input can hold, and is
/// `None` when that is unknown. The buffers reserve no more atoms than it
/// allows, so a header that declares more atoms than the input holds
/// cannot force a large allocation. The declared counts still decide how
/// many rows are read, and so which error is reported.
pub(crate) fn parse_single_frame_ls<'a, L>(
    lines: &mut L,
    max_rows: Option<usize>,
) -> Result<ConFrame, ParseError>
where
    L: LineStream<'a> + Iterator<Item = &'a str>,
{
    let header = parse_frame_header(lines)?;
    // A count the unread input can hold stays in this function, on the same
    // reservation and `Arc` bump as before the bound existed. A count that
    // cannot fit is a separate cold function so it does not shape this loop.
    let Some(total_atoms) = fitting_atom_count(&header.natms_per_type, max_rows) else {
        return parse_oversized_frame(lines, header, max_rows);
    };
    let validate = header.strict_validation;
    use crate::storage_dtype::{ElementKind, FloatArray2};
    let dt = header.metadata.dtypes();
    let f64_positions = dt.positions == ElementKind::Float64;
    let shell = take_coord_shell(total_atoms);
    let (mut pos_flat, mut ids, mut atom_data) = match shell {
        Some(shell) => (shell.coords, shell.ids, shell.atoms),
        None => (
            Vec::with_capacity(total_atoms.saturating_mul(3)),
            Vec::with_capacity(total_atoms),
            Vec::with_capacity(total_atoms),
        ),
    };
    let mut positions_other = if f64_positions {
        None
    } else {
        Some(FloatArray2::zeros(dt.positions, total_atoms, 3))
    };
    let mut global_atom_idx: u64 = 0;
    for (type_idx, num_atoms) in header.natms_per_type.iter().enumerate() {
        let symbol_line = lines.next().ok_or(ParseError::IncompleteFrame)?;
        let symbol: Arc<str> = Arc::from(symbol_line.trim());
        let coord_label = lines.next().ok_or(ParseError::IncompleteFrame)?;
        if validate {
            validate_coordinate_component(type_idx, symbol.as_ref(), coord_label)?;
        }
        global_atom_idx = ingest_coord_rows(
            lines,
            *num_atoms,
            header.spec_version,
            global_atom_idx,
            &symbol,
            CoordSink {
                pos_flat: &mut pos_flat,
                ids: &mut ids,
                positions_other: &mut positions_other,
                atom_data: &mut atom_data,
            },
        )?;
    }
    let positions = if f64_positions {
        FloatArray2::from_f64_row_major(total_atoms, 3, pos_flat)
    } else {
        positions_other.expect("non-f64 positions allocated")
    };
    Ok(crate::types::con_frame_coords_only_with_ids(
        header, atom_data, positions, ids,
    ))
}

/// Atom total when it can use [`Vec::with_capacity`] safely.
///
/// `None` means the counts wrap, exceed `isize::MAX / 4`, or exceed the
/// unread atom lines. A `u32::MAX` count with no size hint is `None`: that
/// reservation is tens of GiB.
#[inline]
fn fitting_atom_count(counts: &[usize], max_rows: Option<usize>) -> Option<usize> {
    let mut total = 0usize;
    for &count in counts {
        total = total.checked_add(count)?;
    }
    if total > (isize::MAX as usize) / 4 {
        return None;
    }
    match max_rows {
        Some(rows) if total > rows => None,
        None if total > (1 << 20) => None,
        _ => Some(total),
    }
}

/// A header declared more atoms than the rest of the input can hold.
///
/// Reservations stay within `max_rows`. Rows are still read up to the
/// declared count, so a short file is `IncompleteFrame`.
#[cold]
#[inline(never)]
fn parse_oversized_frame<'a, L>(
    lines: &mut L,
    header: FrameHeader,
    max_rows: Option<usize>,
) -> Result<ConFrame, ParseError>
where
    L: LineStream<'a> + Iterator<Item = &'a str>,
{
    let validate = header.strict_validation;
    let total_atoms = header
        .natms_per_type
        .iter()
        .fold(0usize, |sum, &n| sum.saturating_add(n));
    let reserve_atoms = max_rows.map_or(0, |rows| total_atoms.min(rows));
    use crate::storage_dtype::FloatArray2;
    let dt = header.metadata.dtypes();
    let mut pos_flat = Vec::new();
    let _ = pos_flat.try_reserve_exact(reserve_atoms.saturating_mul(3));
    let mut ids = Vec::new();
    let _ = ids.try_reserve_exact(reserve_atoms);
    let mut atom_data: Vec<AtomDatum> = Vec::new();
    let _ = atom_data.try_reserve_exact(reserve_atoms);

    let mut global_atom_idx: u64 = 0;
    for (type_idx, num_atoms) in header.natms_per_type.iter().enumerate() {
        let symbol_line = lines.next().ok_or(ParseError::IncompleteFrame)?;
        let symbol: Arc<str> = Arc::from(symbol_line.trim());
        let coord_label = lines.next().ok_or(ParseError::IncompleteFrame)?;
        if validate {
            validate_coordinate_component(type_idx, symbol.as_ref(), coord_label)?;
        }
        global_atom_idx = ingest_coord_rows_cloned(
            lines,
            *num_atoms,
            header.spec_version,
            global_atom_idx,
            &symbol,
            &mut pos_flat,
            &mut ids,
            &mut atom_data,
        )?;
    }
    let mut positions = FloatArray2::from_f64_row_major(total_atoms, 3, pos_flat);
    positions.project_to(dt.positions);
    Ok(crate::types::con_frame_coords_only_with_ids(
        header, atom_data, positions, ids,
    ))
}

struct CoordSink<'a> {
    pos_flat: &'a mut Vec<f64>,
    ids: &'a mut Vec<u64>,
    positions_other: &'a mut Option<crate::storage_dtype::FloatArray2>,
    atom_data: &'a mut Vec<AtomDatum>,
}

/// Parse `n` coordinate rows into flat SoA buffers.
fn ingest_coord_rows<'a, L>(
    lines: &mut L,
    n: usize,
    spec_version: u32,
    global_atom_idx: u64,
    symbol: &Arc<str>,
    sink: CoordSink<'_>,
) -> Result<u64, ParseError>
where
    L: LineStream<'a>,
{
    let CoordSink {
        pos_flat,
        ids,
        positions_other,
        atom_data,
    } = sink;
    if positions_other.is_none() {
        return ingest_coord_rows_f64(
            lines,
            n,
            spec_version,
            global_atom_idx,
            symbol,
            pos_flat,
            ids,
            atom_data,
        );
    }
    let mut global_atom_idx = global_atom_idx;
    for _ in 0..n {
        let (xyz, fixed, atom_id) = lines.next_coord_row(spec_version, global_atom_idx)?;
        if let Some(pos) = positions_other.as_mut() {
            pos.set_f64_row(ids.len(), xyz);
        }
        let index = pos_flat.len() / 3;
        pos_flat.extend_from_slice(&xyz);
        let slot = unsafe { pos_flat.as_mut_ptr().add(index * 3) };
        ids.push(atom_id);
        atom_data.push(AtomDatum::at_slot(Arc::clone(symbol), slot, fixed, atom_id));
        global_atom_idx += 1;
    }
    Ok(global_atom_idx)
}

/// Coordinate rows for the f64 storage path.
///
/// The symbol `Arc` is bumped once for the whole block. Cloning it per atom
/// was a visible share of the coordinate parse.
#[inline(never)]
#[allow(clippy::too_many_arguments)]
fn ingest_coord_rows_f64<'a, L>(
    lines: &mut L,
    n: usize,
    spec_version: u32,
    mut global_atom_idx: u64,
    symbol: &Arc<str>,
    pos_flat: &mut Vec<f64>,
    ids: &mut Vec<u64>,
    atom_data: &mut Vec<AtomDatum>,
) -> Result<u64, ParseError>
where
    L: LineStream<'a>,
{
    let workers = wide_workers();
    if workers >= 2 && n >= 4096 {
        if let Some(mut block) =
            lines.try_take_coord_block(n, spec_version, global_atom_idx, symbol, workers)
        {
            let base = pos_flat.len();
            // Keep a reservation that covers the whole frame. Replacing it
            // with a shorter block makes the next species group reallocate
            // under live pointers.
            if pos_flat.is_empty() && block.xyz.capacity() >= pos_flat.capacity() {
                *pos_flat = block.xyz;
            } else {
                pos_flat.extend_from_slice(&block.xyz);
            }
            let ptr = pos_flat.as_mut_ptr();
            for (i, atom) in block.atoms.iter_mut().enumerate() {
                atom.retarget(unsafe { ptr.add(base + i * 3) });
            }
            if ids.is_empty() {
                *ids = block.ids;
            } else {
                ids.extend_from_slice(&block.ids);
            }
            if atom_data.is_empty() {
                *atom_data = block.atoms;
            } else {
                atom_data.extend(block.atoms);
            }
            return Ok(global_atom_idx + n as u64);
        }
    }
    let mut shares = symbol_shares(symbol, n);
    for _ in 0..n {
        let (xyz, fixed, atom_id) = lines.next_coord_row(spec_version, global_atom_idx)?;
        let sym = shares.take();
        let index = pos_flat.len() / 3;
        pos_flat.extend_from_slice(&xyz);
        let slot = unsafe { pos_flat.as_mut_ptr().add(index * 3) };
        ids.push(atom_id);
        atom_data.push(AtomDatum::at_slot(sym, slot, fixed, atom_id));
        global_atom_idx += 1;
    }
    Ok(global_atom_idx)
}

/// Coordinate rows when the declared count must not bump an `Arc` by `n`.
fn ingest_coord_rows_cloned<'a, L>(
    lines: &mut L,
    n: usize,
    spec_version: u32,
    mut global_atom_idx: u64,
    symbol: &Arc<str>,
    pos_flat: &mut Vec<f64>,
    ids: &mut Vec<u64>,
    atom_data: &mut Vec<AtomDatum>,
) -> Result<u64, ParseError>
where
    L: LineStream<'a>,
{
    for _ in 0..n {
        let (xyz, fixed, atom_id) = lines.next_coord_row(spec_version, global_atom_idx)?;
        pos_flat.extend_from_slice(&xyz);
        ids.push(atom_id);
        atom_data.push(AtomDatum::coords(
            Arc::clone(symbol),
            xyz[0],
            xyz[1],
            xyz[2],
            fixed,
            atom_id,
        ));
        global_atom_idx += 1;
    }
    Ok(global_atom_idx)
}

/// `n` extra strong references to `symbol`, taken one at a time.
///
/// `Arc<str>` stores the header (strong count, weak count) immediately
/// before the bytes. One `fetch_add` covers the block; each [`SymbolShares::take`]
/// turns one of those counts into an `Arc`. Leftovers drop with the guard.
struct SymbolShares {
    data: *const u8,
    len: usize,
    left: usize,
}

impl SymbolShares {
    fn raw(&self) -> *const str {
        std::ptr::slice_from_raw_parts(self.data, self.len) as *const str
    }

    fn take(&mut self) -> Arc<str> {
        debug_assert!(self.left > 0);
        self.left -= 1;
        unsafe { Arc::from_raw(self.raw()) }
    }
}

impl Drop for SymbolShares {
    fn drop(&mut self) {
        for _ in 0..self.left {
            drop(unsafe { Arc::from_raw(self.raw()) });
        }
        self.left = 0;
    }
}

fn symbol_shares(symbol: &Arc<str>, n: usize) -> SymbolShares {
    if n == 0 {
        return SymbolShares {
            data: std::ptr::null(),
            len: 0,
            left: 0,
        };
    }
    let data = Arc::as_ptr(symbol);
    let offset = arc_str_header(symbol.as_ref());
    // SAFETY: `Arc::as_ptr` points at the `str` inside the allocation, and
    // the two atomic counts sit in front of it. This is the same subtraction
    // `Arc::increment_strong_count` performs, done once for `n` clones.
    let bytes = data as *const u8;
    let strong = unsafe { bytes.sub(offset) as *const AtomicUsize };
    unsafe { (*strong).fetch_add(n, Ordering::Relaxed) };
    SymbolShares {
        data: bytes,
        len: symbol.len(),
        left: n,
    }
}

fn arc_str_header(text: &str) -> usize {
    let header = 2 * std::mem::size_of::<usize>();
    let align = std::mem::align_of_val(text);
    let rem = header % align;
    if rem == 0 {
        header
    } else {
        header + (align - rem)
    }
}

#[cfg(test)]
mod symbol_share_tests {
    use super::symbol_shares;
    use std::sync::Arc;

    #[test]
    fn one_bump_matches_n_clones() {
        for text in ["", "Cu", "Gold", "a-much-longer-component-symbol"] {
            let symbol: Arc<str> = Arc::from(text);
            let mut shares = symbol_shares(&symbol, 3);
            assert_eq!(Arc::strong_count(&symbol), 4);
            let owned = [shares.take(), shares.take(), shares.take()];
            assert!(owned.iter().all(|arc| arc.as_ref() == text));
            drop(owned);
            assert_eq!(Arc::strong_count(&symbol), 1);
            let mut partial = symbol_shares(&symbol, 2);
            let one = partial.take();
            drop(partial);
            drop(one);
            assert_eq!(Arc::strong_count(&symbol), 1);
        }
    }
}

fn validate_header_geometry(boxl: &[f64], angles: &[f64]) -> Result<(), ParseError> {
    if boxl
        .iter()
        .any(|length| !length.is_finite() || *length <= 0.0)
        || angles
            .iter()
            .any(|angle| !angle.is_finite() || *angle <= 0.0 || *angle >= 180.0)
    {
        return Err(ParseError::ValidationError(
            "cell geometry must have positive lengths and angles between 0 and 180 degrees"
                .to_string(),
        ));
    }
    Ok(())
}

fn validate_masses(masses_per_type: &[f64]) -> Result<(), ParseError> {
    if masses_per_type
        .iter()
        .any(|mass| !mass.is_finite() || *mass <= 0.0)
    {
        return Err(ParseError::ValidationError(
            "component masses must be positive".to_string(),
        ));
    }
    Ok(())
}

fn validate_coordinate_component(
    type_idx: usize,
    symbol: &str,
    label: &str,
) -> Result<(), ParseError> {
    let expected_label = format!("Coordinates of Component {}", type_idx + 1);
    if !labels_equal(label, &expected_label) {
        return Err(ParseError::ValidationError(format!(
            "expected coordinate label {expected_label:?}, found {label:?}"
        )));
    }
    if symbol != "X" && symbol_to_atomic_number(symbol) == 0 {
        return Err(ParseError::ValidationError(format!(
            "unknown component symbol {symbol}"
        )));
    }
    Ok(())
}

fn validate_section_component(
    section: &str,
    type_idx: usize,
    atom_idx: usize,
    symbol: &str,
    label: &str,
    header: &FrameHeader,
    atom_data: &[AtomDatum],
) -> Result<(), ParseError> {
    let expected_label = format!("{section} of Component {}", type_idx + 1);
    if !labels_equal(label, &expected_label) {
        return Err(ParseError::ValidationError(format!(
            "expected section label {expected_label:?}, found {label:?}"
        )));
    }

    if header.natms_per_type[type_idx] == 0 {
        return Ok(());
    }

    let expected_symbol = atom_data
        .get(atom_idx)
        .map(|atom| atom.symbol.as_ref())
        .ok_or_else(|| {
            ParseError::ValidationError(format!(
                "{section} component {} has no coordinate atom to validate against",
                type_idx + 1
            ))
        })?;
    if symbol != expected_symbol {
        return Err(ParseError::ValidationError(format!(
            "{section} component {} symbol mismatch: expected {expected_symbol}, found {symbol}",
            type_idx + 1
        )));
    }

    Ok(())
}

fn validate_section_atom_identity(
    section: &str,
    atom_idx: usize,
    fixed: [bool; 3],
    atom_id: u64,
    atom_data: &[AtomDatum],
) -> Result<(), ParseError> {
    let atom = atom_data.get(atom_idx).ok_or_else(|| {
        ParseError::ValidationError(format!(
            "{section} row {atom_idx} has no coordinate atom to validate against"
        ))
    })?;

    if atom.fixed() != fixed {
        return Err(ParseError::ValidationError(format!(
            "{section} row {atom_idx} fixed mask mismatch for atom_id {}",
            atom.atom_id
        )));
    }
    if atom.atom_id != atom_id {
        return Err(ParseError::ValidationError(format!(
            "{section} row {atom_idx} atom_id mismatch: expected {}, found {atom_id}",
            atom.atom_id
        )));
    }

    Ok(())
}

/// Consume the single blank line that opens a section.
///
/// A blank line followed only by whitespace to EOF does not open a section
/// (CON-frame.5). The lines stay in the stream.
fn frame_minor(header: &FrameHeader) -> u64 {
    header.metadata.spec_minor()
}

fn section_row_arity(header: &FrameHeader, name: &str) -> Option<usize> {
    let obj = header.metadata.get(meta::SECTION_ARITY)?.as_object()?;
    let n = obj.get(name)?.as_u64()?;
    if n >= 1 { Some(n as usize) } else { None }
}

/// Skip one unknown section of a newer minor.
///
/// `arity` is the value-field count (CON-section.12). Each data row holds
/// that many values, then the flag, and an optional id, so the row has
/// `arity + 1` or `arity + 2` tokens (CON-section.14). The label is the
/// section name under CON-file.5, then `of Component N` (CON-section.15).
fn skip_unknown_section<'a>(
    lines: &mut impl LineStream<'a>,
    header: &FrameHeader,
    name: &str,
    arity: usize,
) -> Result<(), ParseError> {
    if !open_section(lines)? {
        return Err(ParseError::IncompleteSection(name.to_string()));
    }
    let with_flag = arity.saturating_add(1);
    let with_id = arity.saturating_add(2);
    for (type_idx, &num_atoms) in header.natms_per_type.iter().enumerate() {
        if lines.next_line().is_none() {
            return Err(ParseError::IncompleteSection(name.to_string()));
        }
        let label = lines
            .next_line()
            .ok_or_else(|| ParseError::IncompleteSection(name.to_string()))?;
        let expected = format!("{name} of Component {}", type_idx + 1);
        if !labels_equal(label, &expected) {
            return Err(ParseError::IncompleteSection(name.to_string()));
        }
        for _ in 0..num_atoms {
            let row = lines
                .next_line()
                .ok_or_else(|| ParseError::IncompleteSection(name.to_string()))?;
            let found = token_count(row);
            if found != with_flag && found != with_id {
                let expected = if found < with_flag {
                    with_flag
                } else {
                    with_id
                };
                return Err(ParseError::InvalidVectorLength { expected, found });
            }
        }
    }
    Ok(())
}

fn open_section<'a>(lines: &mut impl LineStream<'a>) -> Result<bool, ParseError> {
    match lines.peek_line() {
        Some(line) if line.trim().is_empty() => {
            if lines.remainder_is_whitespace() {
                return Ok(false);
            }
            lines.next_line();
            Ok(true)
        }
        _ => Ok(false),
    }
}

/// Attempts to parse an optional velocity section following coordinate blocks.
///
/// In `.convel` files, after all coordinate blocks there is a blank separator line
/// followed by per-component velocity blocks with the same structure as coordinate
/// blocks (symbol line, "Velocities of Component N" line, then atom lines with
/// `vx vy vz fixed atomID`).
///
/// This function peeks at the next line. If it is blank (or contains only whitespace),
/// it consumes the blank line and writes each velocity into `column` (row-major
/// `n * 3` `f64`s, one allocation). If the next line is not blank (or is absent),
/// no velocities are parsed and `column` is left untouched.
///
/// Returns `Ok(true)` if velocities were found and parsed, `Ok(false)` otherwise.
pub fn parse_velocity_section<'a>(
    lines: &mut impl LineStream<'a>,
    header: &FrameHeader,
    atom_data: &mut [AtomDatum],
    column: &mut Vec<f64>,
) -> Result<bool, ParseError> {
    let validate = header.strict_validation;
    // A blank line followed only by whitespace to EOF is not a velocity
    // block (CON-frame.5, D2). Leave those lines for the iterator.
    match lines.peek_line() {
        Some(line) if line.trim().is_empty() => {
            if lines.remainder_is_whitespace() {
                return Ok(false);
            }
            // CON-section.2 / the legacy fallback: the blank counts only
            // when it begins a complete velocity block. Otherwise it is
            // line 1 of the next frame (CON-frame.7, CON-header.4).
            if !lines.blank_line_opens_section() {
                return Ok(false);
            }
            lines.next_line();
        }
        _ => return Ok(false),
    }

    let n = atom_data.len();
    column.clear();
    column.resize(n.saturating_mul(3), 0.0);
    let mut atom_idx: usize = 0;
    for (type_idx, &num_atoms) in header.natms_per_type.iter().enumerate() {
        // Symbol line
        let symbol = lines
            .next_line()
            .ok_or(ParseError::IncompleteSection(SECTION_VELOCITIES.into()))?
            .trim();

        // "Velocities of Component N" line
        let comp_line = lines
            .next_line()
            .ok_or(ParseError::IncompleteSection(SECTION_VELOCITIES.into()))?;
        if !label_has_prefix(comp_line, "Velocities of Component") {
            return Err(ParseError::IncompleteSection(SECTION_VELOCITIES.into()));
        }
        if validate {
            validate_section_component(
                "Velocities",
                type_idx,
                atom_idx,
                symbol,
                comp_line,
                header,
                atom_data,
            )?;
        }

        for _ in 0..num_atoms {
            let (vals, fixed, atom_id) =
                take_vector_row(lines, header.spec_version, SECTION_VELOCITIES)?;
            if validate {
                validate_section_atom_identity("velocities", atom_idx, fixed, atom_id, atom_data)?;
            }
            if atom_idx < n {
                let base = atom_idx * 3;
                column[base] = vals[0];
                column[base + 1] = vals[1];
                column[base + 2] = vals[2];
            }
            atom_idx += 1;
        }
    }

    Ok(true)
}

/// Attempts to parse a force section following coordinate (and optional velocity) blocks.
///
/// Force sections mirror velocity sections: a blank separator line followed by per-component
/// force blocks (symbol line, "Forces of Component N" line, then atom lines with
/// `fx fy fz fixed_flag atom_id`).
///
/// Returns `Ok(true)` if forces were found and parsed, `Ok(false)` otherwise.
pub fn parse_force_section<'a>(
    lines: &mut impl LineStream<'a>,
    header: &FrameHeader,
    atom_data: &mut [AtomDatum],
    column: &mut Vec<f64>,
) -> Result<bool, ParseError> {
    let validate = header.strict_validation;
    if !open_section(lines)? {
        return Ok(false);
    }

    let n = atom_data.len();
    column.clear();
    column.resize(n.saturating_mul(3), 0.0);
    let mut atom_idx: usize = 0;
    for (type_idx, &num_atoms) in header.natms_per_type.iter().enumerate() {
        let symbol = lines
            .next_line()
            .ok_or(ParseError::IncompleteSection(SECTION_FORCES.into()))?
            .trim();

        let comp_line = lines
            .next_line()
            .ok_or(ParseError::IncompleteSection(SECTION_FORCES.into()))?;
        if !label_has_prefix(comp_line, "Forces of Component") {
            return Err(ParseError::IncompleteSection(SECTION_FORCES.into()));
        }
        if validate {
            validate_section_component(
                "Forces", type_idx, atom_idx, symbol, comp_line, header, atom_data,
            )?;
        }

        for _ in 0..num_atoms {
            let (vals, fixed, atom_id) =
                take_vector_row(lines, header.spec_version, SECTION_FORCES)?;
            if validate {
                validate_section_atom_identity("forces", atom_idx, fixed, atom_id, atom_data)?;
            }
            if atom_idx < n {
                let base = atom_idx * 3;
                column[base] = vals[0];
                column[base + 1] = vals[1];
                column[base + 2] = vals[2];
            }
            atom_idx += 1;
        }
    }

    Ok(true)
}

/// Attempts to parse an energies section following coordinate (and optional
/// velocity / force) blocks.
///
/// Energy sections mirror force sections but with one scalar per atom:
/// blank separator, then per-component blocks of (symbol, "Energies of
/// Component N", and atom lines `e fixed_flag atom_id`). The two
/// trailing identity columns are optional and used only for strict
/// validation; in non-strict mode any whitespace after the energy is
/// ignored.
///
/// Returns `Ok(true)` if energies were found and parsed, `Ok(false)`
/// otherwise.
pub fn parse_energy_section<'a>(
    lines: &mut impl LineStream<'a>,
    header: &FrameHeader,
    atom_data: &mut [AtomDatum],
) -> Result<bool, ParseError> {
    let validate = header.strict_validation;
    if !open_section(lines)? {
        return Ok(false);
    }

    let mut atom_idx: usize = 0;
    for (type_idx, &num_atoms) in header.natms_per_type.iter().enumerate() {
        let symbol = lines
            .next_line()
            .ok_or(ParseError::IncompleteSection(SECTION_ENERGIES.into()))?
            .trim();

        let comp_line = lines
            .next_line()
            .ok_or(ParseError::IncompleteSection(SECTION_ENERGIES.into()))?;
        if !label_has_prefix(comp_line, "Energies of Component") {
            return Err(ParseError::IncompleteSection(SECTION_ENERGIES.into()));
        }
        if validate {
            validate_section_component(
                "Energies", type_idx, atom_idx, symbol, comp_line, header, atom_data,
            )?;
        }

        for _ in 0..num_atoms {
            let energy_line = lines
                .next_line()
                .ok_or(ParseError::IncompleteSection(SECTION_ENERGIES.into()))?;
            let (value, fixed, atom_id) =
                parse_scalar_section_row(energy_line, header.spec_version)?;
            if validate {
                validate_section_atom_identity("energies", atom_idx, fixed, atom_id, atom_data)?;
            }
            if atom_idx < atom_data.len() {
                atom_data[atom_idx].set_energy(Some(value));
            }
            atom_idx += 1;
        }
    }

    Ok(true)
}

/// Parses declared sections from a frame's header metadata.
///
/// If `header.sections` is non-empty (v2 file with `"sections"` key in JSON),
/// parses each declared section in order. Otherwise falls back to legacy
/// blank-separator velocity detection.
pub fn parse_declared_sections<'a>(
    lines: &mut impl LineStream<'a>,
    frame: &mut ConFrame,
) -> Result<usize, ParseError> {
    let mut applied = 0usize;
    // A zero-atom frame has no coordinate block and no section block.
    // The legacy velocity fallback must not run (CON-count.1).
    if frame.header.natm_types == 0 {
        if !frame.header.sections.is_empty() {
            return Err(ParseError::InvalidCount(
                "a zero-atom frame cannot declare sections".into(),
            ));
        }
        return Ok(0);
    }
    if !frame.header.sections_declared && frame.header.sections.is_empty() {
        // Legacy: try velocity detection via blank separator. A blank
        // line can also be the empty line 1 of the next frame, so look for
        // the velocity label before committing (CON-frame.5).
        let mut column = Vec::new();
        let found = lines.blank_opens_velocity_block()
            && parse_velocity_section(lines, &frame.header, &mut frame.atom_data, &mut column)?;
        if found {
            frame.adopt_velocity_column(column);
            frame.header.sections.push(SECTION_VELOCITIES.into());
            applied = 1;
        }
    } else {
        let sections = std::mem::take(&mut frame.header.sections);
        for section in &sections {
            match section.as_str() {
                SECTION_VELOCITIES => {
                    let mut column = Vec::new();
                    let found = parse_velocity_section(
                        lines,
                        &frame.header,
                        &mut frame.atom_data,
                        &mut column,
                    )?;
                    if !found {
                        return Err(ParseError::IncompleteSection(SECTION_VELOCITIES.into()));
                    }
                    frame.adopt_velocity_column(column);
                    applied += 1;
                }
                SECTION_FORCES => {
                    let mut column = Vec::new();
                    let found = parse_force_section(
                        lines,
                        &frame.header,
                        &mut frame.atom_data,
                        &mut column,
                    )?;
                    if !found {
                        return Err(ParseError::IncompleteSection(SECTION_FORCES.into()));
                    }
                    frame.adopt_force_column(column);
                    applied += 1;
                }
                SECTION_ENERGIES => {
                    let found = parse_energy_section(lines, &frame.header, &mut frame.atom_data)?;
                    if !found {
                        return Err(ParseError::IncompleteSection(SECTION_ENERGIES.into()));
                    }
                    applied += 1;
                }
                SECTION_CHARGES => {
                    let found = parse_charge_section(lines, &frame.header, &mut frame.atom_data)?;
                    if !found {
                        return Err(ParseError::IncompleteSection(SECTION_CHARGES.into()));
                    }
                    applied += 1;
                }
                SECTION_SPINS => {
                    let found = parse_spin_section(lines, &frame.header, &mut frame.atom_data)?;
                    if !found {
                        return Err(ParseError::IncompleteSection(SECTION_SPINS.into()));
                    }
                    applied += 1;
                }
                SECTION_MAGMOMS => {
                    let found = parse_magmom_section(lines, &frame.header, &mut frame.atom_data)?;
                    if !found {
                        return Err(ParseError::IncompleteSection(SECTION_MAGMOMS.into()));
                    }
                    applied += 1;
                }
                SECTION_DISPLACEMENTS => {
                    let found =
                        parse_displacement_section(lines, &frame.header, &mut frame.atom_data)?;
                    if !found {
                        return Err(ParseError::IncompleteSection(SECTION_DISPLACEMENTS.into()));
                    }
                    applied += 1;
                }
                SECTION_SPREADS => {
                    let found = parse_spread_section(lines, &frame.header, &mut frame.atom_data)?;
                    if !found {
                        return Err(ParseError::IncompleteSection(SECTION_SPREADS.into()));
                    }
                    applied += 1;
                }
                other => {
                    if frame_minor(&frame.header) > u64::from(crate::CON_SPEC_MINOR)
                        && let Some(arity) = section_row_arity(&frame.header, other)
                    {
                        skip_unknown_section(lines, &frame.header, other, arity)?;
                        continue;
                    }
                    return Err(ParseError::UnknownSection(other.to_string()));
                }
            }
        }
        frame.header.sections = sections;
    }
    reject_velocities_without_time(&frame.header)?;
    Ok(applied)
}

/// CON-units.11. Versions 1, 2, and 3.0 do not require `units.time`
/// (CON-units.8, CON-units.9). From minor 1 the key is required when a
/// velocities section is present.
fn reject_velocities_without_time(header: &FrameHeader) -> Result<(), ParseError> {
    if header.spec_version < 3 || frame_minor(header) < 1 {
        return Ok(());
    }
    if !header
        .sections
        .iter()
        .any(|name| name == SECTION_VELOCITIES)
    {
        return Ok(());
    }
    let time_ok = header
        .metadata
        .get(meta::UNITS)
        .and_then(|units| units.get("time"))
        .and_then(|value| value.as_str())
        .is_some_and(|unit| !unit.trim().is_empty());
    if time_ok {
        Ok(())
    } else {
        Err(ParseError::InvalidUnits(
            "minor >= 1 requires units.time when a velocities section is present".into(),
        ))
    }
}

/// Scalar per-atom section: blank separator + per-type blocks
/// (`symbol`, `"X of Component N"`, lines `value fixed atom_id`).
fn parse_scalar_atom_section<'a>(
    lines: &mut impl LineStream<'a>,
    header: &FrameHeader,
    atom_data: &mut [AtomDatum],
    section_name: &str,
    component_label: &str,
    mut set_value: impl FnMut(&mut AtomDatum, f64),
) -> Result<bool, ParseError> {
    let validate = header.strict_validation;
    if !open_section(lines)? {
        return Ok(false);
    }

    let mut atom_idx: usize = 0;
    for (type_idx, &num_atoms) in header.natms_per_type.iter().enumerate() {
        let symbol = lines
            .next_line()
            .ok_or_else(|| ParseError::IncompleteSection(section_name.into()))?
            .trim();

        let comp_line = lines
            .next_line()
            .ok_or_else(|| ParseError::IncompleteSection(section_name.into()))?;
        if !label_has_prefix(comp_line, component_label) {
            return Err(ParseError::IncompleteSection(section_name.into()));
        }
        if validate {
            validate_section_component(
                component_label.trim_end_matches(" of Component").trim(),
                type_idx,
                atom_idx,
                symbol,
                comp_line,
                header,
                atom_data,
            )?;
        }

        for _ in 0..num_atoms {
            let data_line = lines
                .next_line()
                .ok_or_else(|| ParseError::IncompleteSection(section_name.into()))?;
            let (value, fixed, atom_id) = parse_scalar_section_row(data_line, header.spec_version)?;
            if validate {
                validate_section_atom_identity(section_name, atom_idx, fixed, atom_id, atom_data)?;
            }
            if atom_idx < atom_data.len() {
                set_value(&mut atom_data[atom_idx], value);
            }
            atom_idx += 1;
        }
    }
    Ok(true)
}

pub fn parse_charge_section<'a>(
    lines: &mut impl LineStream<'a>,
    header: &FrameHeader,
    atom_data: &mut [AtomDatum],
) -> Result<bool, ParseError> {
    parse_scalar_atom_section(
        lines,
        header,
        atom_data,
        SECTION_CHARGES,
        "Charges of Component",
        |a, v| a.set_charge(Some(v)),
    )
}

pub fn parse_spin_section<'a>(
    lines: &mut impl LineStream<'a>,
    header: &FrameHeader,
    atom_data: &mut [AtomDatum],
) -> Result<bool, ParseError> {
    parse_scalar_atom_section(
        lines,
        header,
        atom_data,
        SECTION_SPINS,
        "Spins of Component",
        |a, v| a.set_spin(Some(v)),
    )
}

fn parse_vector_named_section<'a>(
    lines: &mut impl LineStream<'a>,
    header: &FrameHeader,
    atom_data: &mut [AtomDatum],
    names: (&str, &str, &str),
    mut set_value: impl FnMut(&mut AtomDatum, [f64; 3]),
    mut check_row: impl FnMut(usize, [f64; 3]) -> Result<(), ParseError>,
) -> Result<bool, ParseError> {
    let (section_name, component_label, label_kind) = names;
    let validate = header.strict_validation;
    if !open_section(lines)? {
        return Ok(false);
    }
    let mut atom_idx: usize = 0;
    for (type_idx, &num_atoms) in header.natms_per_type.iter().enumerate() {
        let symbol = lines
            .next_line()
            .ok_or_else(|| ParseError::IncompleteSection(section_name.into()))?
            .trim();
        let comp_line = lines
            .next_line()
            .ok_or_else(|| ParseError::IncompleteSection(section_name.into()))?;
        if !label_has_prefix(comp_line, component_label) {
            return Err(ParseError::IncompleteSection(section_name.into()));
        }
        if validate {
            validate_section_component(
                label_kind, type_idx, atom_idx, symbol, comp_line, header, atom_data,
            )?;
        }
        for _ in 0..num_atoms {
            let (vals, fixed, atom_id) = take_vector_row(lines, header.spec_version, section_name)?;
            check_row(atom_idx, vals)?;
            if validate {
                validate_section_atom_identity(section_name, atom_idx, fixed, atom_id, atom_data)?;
            }
            if atom_idx < atom_data.len() {
                set_value(&mut atom_data[atom_idx], vals);
            }
            atom_idx += 1;
        }
    }
    Ok(true)
}

/// Magmoms: 3-vector per atom, same layout as velocities/forces.
pub fn parse_magmom_section<'a>(
    lines: &mut impl LineStream<'a>,
    header: &FrameHeader,
    atom_data: &mut [AtomDatum],
) -> Result<bool, ParseError> {
    parse_vector_named_section(
        lines,
        header,
        atom_data,
        (SECTION_MAGMOMS, "Magmoms of Component", "Magmoms"),
        |a, v| a.set_magmom(Some(v)),
        |_, _| Ok(()),
    )
}

/// Displacements: 3-vector per atom (Angstrom), same layout as velocities/forces.
pub fn parse_displacement_section<'a>(
    lines: &mut impl LineStream<'a>,
    header: &FrameHeader,
    atom_data: &mut [AtomDatum],
) -> Result<bool, ParseError> {
    parse_vector_named_section(
        lines,
        header,
        atom_data,
        (
            SECTION_DISPLACEMENTS,
            "Displacements of Component",
            "Displacements",
        ),
        |a, v| a.set_displacement(Some(v)),
        |_, _| Ok(()),
    )
}

/// Spreads: per-atom root-mean-square spread `[sx, sy, sz]` (Angstrom, each
/// non-negative), same layout as velocities/forces.
pub fn parse_spread_section<'a>(
    lines: &mut impl LineStream<'a>,
    header: &FrameHeader,
    atom_data: &mut [AtomDatum],
) -> Result<bool, ParseError> {
    parse_vector_named_section(
        lines,
        header,
        atom_data,
        (SECTION_SPREADS, "Spreads of Component", "Spreads"),
        |a, v| a.set_spread(Some(v)),
        |atom_idx, vals| {
            if !vals.iter().all(|v| v.is_finite()) {
                return Err(ParseError::NonFiniteValue(format!(
                    "spreads atom {atom_idx}"
                )));
            }
            if vals.iter().any(|v| *v < 0.0) {
                return Err(crate::types::invalid_spread(atom_idx));
            }
            Ok(())
        },
    )
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::iterators::ConFrameIterator;

    #[test]
    fn test_parse_line_of_n_success() {
        let line = "1.0 2.5 -3.0";
        let values = parse_line_of_n::<f64>(line, 3).unwrap();
        assert_eq!(values, vec![1.0, 2.5, -3.0]);
    }

    #[test]
    fn test_parse_line_of_n_too_short() {
        let line = "1.0 2.5";
        let result = parse_line_of_n::<f64>(line, 3);
        assert!(result.is_err());
        assert!(matches!(
            result.unwrap_err(),
            ParseError::InvalidVectorLength {
                expected: 3,
                found: 2
            }
        ));
    }

    #[test]
    fn test_parse_line_of_n_too_long() {
        let line = "1.0 2.5 -3.0 4.0";
        let result = parse_line_of_n::<f64>(line, 3);
        assert!(result.is_err());
        assert!(matches!(
            result.unwrap_err(),
            ParseError::InvalidVectorLength {
                expected: 3,
                found: 4
            }
        ));
    }

    #[test]
    fn test_parse_line_of_n_invalid_float() {
        let line = "1.0 abc -3.0";
        let result = parse_line_of_n::<f64>(line, 3);
        assert!(result.is_err());
        assert!(matches!(
            result.unwrap_err(),
            ParseError::InvalidNumberFormat(_)
        ));
    }

    #[test]
    fn test_parse_frame_header_success() {
        let lines = [
            "PREBOX1",
            "{\"con_spec_version\":2}",
            "10.0 20.0 30.0",
            "90.0 90.0 90.0",
            "POSTBOX1",
            "POSTBOX2",
            "2",
            "1 1",
            "12.011 1.008",
        ];
        let mut line_it = lines.iter().copied();
        match parse_frame_header(&mut line_it) {
            Ok(header) => {
                assert_eq!(header.prebox_header.user, "PREBOX1");
                assert_eq!(header.spec_version, 2);
                assert_eq!(header.boxl, [10.0, 20.0, 30.0]);
                assert_eq!(header.angles, [90.0, 90.0, 90.0]);
                assert_eq!(header.postbox_header, ["POSTBOX1", "POSTBOX2"]);
                assert_eq!(header.natm_types, 2);
                assert_eq!(header.natms_per_type, vec![1, 1]);
                assert_eq!(header.masses_per_type, vec![12.011, 1.008]);
            }
            Err(e) => {
                panic!("Parsing failed when it should have succeeded. Error: {e:?}");
            }
        }
    }

    #[test]
    fn test_parse_frame_header_missing_line() {
        let lines = [
            "PREBOX1",
            "{\"con_spec_version\":2}",
            "10.0 20.0 30.0",
            "90.0 90.0 90.0",
            "POSTBOX1",
            "POSTBOX2",
            "2",
            "1 1",
            // Missing masses_per_type
        ];
        let mut line_it = lines.iter().copied();
        let result = parse_frame_header(&mut line_it);
        assert!(result.is_err());
        assert!(matches!(result.unwrap_err(), ParseError::IncompleteHeader));
    }

    #[test]
    fn test_parse_frame_header_missing_spec_version() {
        let lines = [
            "PREBOX1",
            "{}",
            "10.0 20.0 30.0",
            "90.0 90.0 90.0",
            "POSTBOX1",
            "POSTBOX2",
            "2",
            "1 1",
            "12.011 1.008",
        ];
        let mut line_it = lines.iter().copied();
        let result = parse_frame_header(&mut line_it);
        assert!(result.is_err());
        assert!(matches!(
            result.unwrap_err(),
            ParseError::MissingSpecVersion
        ));
    }

    #[test]
    fn test_parse_frame_header_legacy_no_json() {
        // A non-JSON line 2 is treated as a legacy (v1) file.
        let lines = [
            "PREBOX1",
            "0.0000 TIME",
            "10.0 20.0 30.0",
            "90.0 90.0 90.0",
            "POSTBOX1",
            "POSTBOX2",
            "2",
            "1 1",
            "12.011 1.008",
        ];
        let mut line_it = lines.iter().copied();
        let header = parse_frame_header(&mut line_it).unwrap();
        assert_eq!(header.spec_version, 1);
        assert!(header.metadata.is_empty());
    }

    #[test]
    fn test_parse_frame_header_malformed_json() {
        // Line 2 starts with '{' but is not valid JSON -- this IS an error.
        let lines = [
            "PREBOX1",
            "{broken json",
            "10.0 20.0 30.0",
            "90.0 90.0 90.0",
            "POSTBOX1",
            "POSTBOX2",
            "2",
            "1 1",
            "12.011 1.008",
        ];
        let mut line_it = lines.iter().copied();
        let result = parse_frame_header(&mut line_it);
        assert!(result.is_err());
        assert!(matches!(
            result.unwrap_err(),
            ParseError::InvalidMetadataJson(_)
        ));
    }

    #[test]
    fn test_parse_frame_header_unsupported_version() {
        let lines = [
            "PREBOX1",
            "{\"con_spec_version\":999}",
            "10.0 20.0 30.0",
            "90.0 90.0 90.0",
            "POSTBOX1",
            "POSTBOX2",
            "2",
            "1 1",
            "12.011 1.008",
        ];
        let mut line_it = lines.iter().copied();
        let result = parse_frame_header(&mut line_it);
        assert!(result.is_err());
        assert!(matches!(
            result.unwrap_err(),
            ParseError::UnsupportedVersion(999)
        ));
    }

    #[test]
    fn test_v3_missing_units_rejected() {
        let lines = [
            "PREBOX1",
            "{\"con_spec_version\":3}",
            "10.0 20.0 30.0",
            "90.0 90.0 90.0",
            "POSTBOX1",
            "POSTBOX2",
            "1",
            "1",
            "1.0",
        ];
        let mut line_it = lines.iter().copied();
        let err = parse_frame_header(&mut line_it).unwrap_err();
        let msg = err.to_string();
        assert!(
            msg.contains("units") || msg.contains("v3"),
            "expected units error, got {msg}"
        );
    }

    #[test]
    fn test_v3_invalid_units_rejected() {
        let lines = [
            "PREBOX1",
            r#"{"con_spec_version":3,"units":{"length":"eV","energy":"angstrom"}}"#,
            "10.0 20.0 30.0",
            "90.0 90.0 90.0",
            "POSTBOX1",
            "POSTBOX2",
            "1",
            "1",
            "1.0",
        ];
        let mut line_it = lines.iter().copied();
        assert!(parse_frame_header(&mut line_it).is_err());
    }

    #[test]
    fn test_v3_valid_units_exposes_length_energy() {
        let lines = [
            "PREBOX1",
            r#"{"con_spec_version":3,"units":{"length":"angstrom","energy":"eV","mass":"amu","time":"fs"}}"#,
            "10.0 20.0 30.0",
            "90.0 90.0 90.0",
            "POSTBOX1",
            "POSTBOX2",
            "1",
            "1",
            "1.0",
        ];
        let mut line_it = lines.iter().copied();
        let header = parse_frame_header(&mut line_it).unwrap();
        assert_eq!(header.spec_version, 3);
        assert_eq!(header.length_unit(), Some("angstrom"));
        assert_eq!(header.energy_unit(), Some("eV"));
        let f = header.conversion_factor_to("length", "nm").unwrap();
        assert!((f - 0.1).abs() < 1e-12);
    }

    #[test]
    fn test_parse_frame_header_extra_metadata_preserved() {
        let lines = [
            "PREBOX1",
            "{\"con_spec_version\":2,\"generator\":\"test\"}",
            "10.0 20.0 30.0",
            "90.0 90.0 90.0",
            "POSTBOX1",
            "POSTBOX2",
            "2",
            "1 1",
            "12.011 1.008",
        ];
        let mut line_it = lines.iter().copied();
        let header = parse_frame_header(&mut line_it).unwrap();
        assert_eq!(header.spec_version, 2);
        assert_eq!(
            header.metadata.get("generator"),
            Some(&serde_json::Value::String("test".to_string()))
        );
    }

    #[test]
    fn test_parse_frame_header_invalid_natms_per_type() {
        let lines = [
            "PREBOX1",
            "{\"con_spec_version\":2}",
            "10.0 20.0 30.0",
            "90.0 90.0 90.0",
            "POSTBOX1",
            "POSTBOX2",
            "2",
            "1 1 1", // 3 values, but natm_types is 2
            "12.011 1.008",
        ];
        let mut line_it = lines.iter().copied();
        let result = parse_frame_header(&mut line_it);
        assert!(result.is_err());
        assert!(matches!(
            result.unwrap_err(),
            ParseError::InvalidVectorLength {
                expected: 2,
                found: 3
            }
        ));
    }

    #[test]
    fn test_parse_single_frame_success() {
        let lines = vec![
            "PREBOX1",
            "{\"con_spec_version\":2}",
            "10.0 20.0 30.0",
            "90.0 90.0 90.0",
            "POSTBOX1",
            "POSTBOX2",
            "2",
            "3 3",
            "12.011 1.008",
            "1",
            "Coordinates of Component 1",
            "0.0 0.0 0.0 0 1",
            "1.0940 0.0 0.0 0 2",
            "-0.5470 0.9499 0.0 0 3",
            "2",
            "Coordinates of Component 2",
            "5.0 5.0 5.0 0 4",
            "6.0940 5.0 5.0 0 5",
            "5.5470 5.9499 5.0 0 6",
        ];
        let mut line_it = lines.iter().copied();
        let frame = parse_single_frame(&mut line_it).unwrap();

        assert_eq!(frame.header.natm_types, 2);
        assert_eq!(frame.header.natms_per_type, vec![3, 3]);
        assert_eq!(frame.header.masses_per_type, vec![12.011, 1.008]);
        assert_eq!(frame.atom_data.len(), 6);
        assert_eq!(&*frame.atom_data[0].symbol, "1");
        assert_eq!(frame.atom_data[0].atom_id, 1);
        assert_eq!(&*frame.atom_data[5].symbol, "2");
        assert_eq!(frame.atom_data[5].atom_id, 6);
    }

    #[test]
    fn test_column_4_value_1_is_x_only_on_spec_2() {
        let lines = [
            "Generated by probe",
            "{\"con_spec_version\":2}",
            "20.0 20.0 20.0",
            "90.0 90.0 90.0",
            "",
            "",
            "1",
            "1",
            "63.5",
            "Cu",
            "Coordinates of component   1",
            "1.0 2.0 3.0 1 1",
        ];
        let mut line_it = lines.iter().copied();
        let frame = parse_single_frame(&mut line_it).unwrap();
        assert_eq!(frame.header.spec_version, 2);
        assert_eq!(frame.atom_data[0].fixed(), [true, false, false]);
    }

    #[test]
    fn test_column_4_value_1_is_all_fixed_on_spec_1() {
        let lines = [
            "Generated by probe",
            "0.0000 TIME",
            "20.0 20.0 20.0",
            "90.0 90.0 90.0",
            "",
            "",
            "1",
            "1",
            "63.5",
            "Cu",
            "Coordinates of component   1",
            "1.0 2.0 3.0 1 1",
        ];
        let mut line_it = lines.iter().copied();
        let frame = parse_single_frame(&mut line_it).unwrap();
        assert_eq!(frame.header.spec_version, 1);
        assert_eq!(frame.atom_data[0].fixed(), [true, true, true]);
    }

    #[test]
    fn test_parse_single_frame_missing_line() {
        // With a valid header but truncated atom data, we get IncompleteFrame.
        let lines = vec![
            "PREBOX1",
            "{\"con_spec_version\":2}",
            "10.0 20.0 30.0",
            "90.0 90.0 90.0",
            "POSTBOX1",
            "POSTBOX2",
            "2",
            "3 3",
            "12.011 1.008",
            "1",
            "Coordinates of Component 1",
            "0.0 0.0 0.0 0 1",
            "1.0940 0.0 0.0 0 2",
            "-0.5470 0.9499 0.0 0 3",
            // Missing Component 2 entirely
        ];
        let mut line_it = lines.iter().copied();
        let result = parse_single_frame(&mut line_it);
        assert!(result.is_err());
        assert!(matches!(result.unwrap_err(), ParseError::IncompleteFrame));
    }

    #[test]
    fn test_parse_single_frame_missing_atom_index_defaults_sequential() {
        // Column 5 (atom_index) is optional; when absent, defaults to sequential.
        let lines = vec![
            "PREBOX1",
            "{\"con_spec_version\":2}",
            "10.0 20.0 30.0",
            "90.0 90.0 90.0",
            "POSTBOX1",
            "POSTBOX2",
            "2",
            "3 3",
            "12.011 1.008",
            "1",
            "Coordinates of Component 1",
            "0.0 0.0 0.0 0 1",
            "1.0940 0.0 0.0 0 2",
            "-0.5470 0.9499 0.0 0 3",
            "2",
            "Coordinates of Component 2",
            "5.0 5.0 5.0 0",       // No atom_index: defaults to 3
            "6.0940 5.0 5.0 0 10", // Explicit atom_index: 10
            "5.5470 5.9499 5.0 0", // No atom_index: defaults to 5
        ];
        let mut line_it = lines.iter().copied();
        let frame = parse_single_frame(&mut line_it).unwrap();
        assert_eq!(frame.atom_data.len(), 6);
        // First type: explicit atom_index values
        assert_eq!(frame.atom_data[0].atom_id, 1);
        assert_eq!(frame.atom_data[1].atom_id, 2);
        assert_eq!(frame.atom_data[2].atom_id, 3);
        // Second type: mixed explicit and defaulted
        assert_eq!(frame.atom_data[3].atom_id, 3); // defaulted (global idx 3)
        assert_eq!(frame.atom_data[4].atom_id, 10); // explicit
        assert_eq!(frame.atom_data[5].atom_id, 5); // defaulted (global idx 5)
    }

    #[test]
    fn test_parse_single_frame_too_few_columns_fails() {
        // Only 3 columns (missing fixed_flag too) should still fail.
        let lines = [
            "PREBOX1",
            "{\"con_spec_version\":2}",
            "10.0 20.0 30.0",
            "90.0 90.0 90.0",
            "POSTBOX1",
            "POSTBOX2",
            "1",
            "1",
            "12.011",
            "C",
            "Coordinates of Component 1",
            "0.0 0.0 0.0", // Only 3 values
        ];
        let mut line_it = lines.iter().copied();
        let result = parse_single_frame(&mut line_it);
        assert!(result.is_err());
        assert!(matches!(
            result.unwrap_err(),
            ParseError::InvalidVectorLength {
                expected: 5,
                found: 3
            }
        ));
    }

    #[test]
    fn test_parse_velocity_section_present() {
        let lines = vec![
            "PREBOX1",
            "{\"con_spec_version\":2}",
            "10.0 20.0 30.0",
            "90.0 90.0 90.0",
            "POSTBOX1",
            "POSTBOX2",
            "2",
            "1 1",
            "63.546 1.008",
            "Cu",
            "Coordinates of Component 1",
            "0.0 0.0 0.0 1 0",
            "H",
            "Coordinates of Component 2",
            "1.0 2.0 3.0 0 1",
            "",
            "Cu",
            "Velocities of Component 1",
            "0.1 0.2 0.3 1 0",
            "H",
            "Velocities of Component 2",
            "0.4 0.5 0.6 0 1",
        ];
        let mut line_it = lines.iter().copied().peekable();
        // Parse the frame first (consuming 15 lines)
        let mut frame =
            parse_single_frame(&mut line_it).expect("coordinate parsing should succeed");
        assert!(!frame.has_velocities());

        // Now parse the velocity section
        let mut column = Vec::new();
        let has_vel = parse_velocity_section(
            &mut line_it,
            &frame.header,
            &mut frame.atom_data,
            &mut column,
        )
        .expect("velocity parsing should succeed");
        assert!(has_vel);
        frame.adopt_velocity_column(column);
        assert_eq!(frame.atom_data[0].velocity(), Some([0.1, 0.2, 0.3]));
        assert_eq!(frame.atom_data[1].velocity(), Some([0.4, 0.5, 0.6]));
    }

    #[test]
    fn test_validate_true_accepts_matching_section_identity() {
        let text = r#"
PREBOX1
{"con_spec_version":2,"sections":["velocities"],"validate":true}
10.0 20.0 30.0
90.0 90.0 90.0
POSTBOX1
POSTBOX2
2
1 1
63.546 1.008
Cu
Coordinates of Component 1
0.0 0.0 0.0 5 0
H
Coordinates of Component 2
1.0 2.0 3.0 0 1

Cu
Velocities of Component 1
0.1 0.2 0.3 5 0
H
Velocities of Component 2
0.4 0.5 0.6 0 1
"#;
        let mut iter = ConFrameIterator::new(text.trim());
        let frame = iter.next().unwrap().unwrap();

        assert!(frame.has_velocities());
        assert_eq!(
            frame
                .header
                .metadata
                .get("validate")
                .and_then(|v| v.as_bool()),
            Some(true)
        );
    }

    #[test]
    fn test_validate_true_rejects_section_atom_id_mismatch() {
        let text = r#"
PREBOX1
{"con_spec_version":2,"sections":["velocities"],"validate":true}
10.0 20.0 30.0
90.0 90.0 90.0
POSTBOX1
POSTBOX2
1
1
63.546
Cu
Coordinates of Component 1
0.0 0.0 0.0 5 0

Cu
Velocities of Component 1
0.1 0.2 0.3 5 99
"#;
        let mut iter = ConFrameIterator::new(text.trim());
        let err = iter.next().unwrap().unwrap_err();

        assert!(matches!(err, ParseError::ValidationError(_)));
        assert!(err.to_string().contains("atom_id mismatch"));
    }

    #[test]
    fn test_validate_true_rejects_section_symbol_mismatch() {
        let text = r#"
PREBOX1
{"con_spec_version":2,"sections":["forces"],"validate":true}
10.0 20.0 30.0
90.0 90.0 90.0
POSTBOX1
POSTBOX2
1
1
63.546
Cu
Coordinates of Component 1
0.0 0.0 0.0 5 0

H
Forces of Component 1
0.1 0.2 0.3 5 0
"#;
        let mut iter = ConFrameIterator::new(text.trim());
        let err = iter.next().unwrap().unwrap_err();

        assert!(matches!(err, ParseError::ValidationError(_)));
        assert!(err.to_string().contains("symbol mismatch"));
    }

    #[test]
    fn test_validate_absent_allows_legacy_duplicate_identity_mismatch() {
        let text = r#"
PREBOX1
{"con_spec_version":2,"sections":["velocities"]}
10.0 20.0 30.0
90.0 90.0 90.0
POSTBOX1
POSTBOX2
1
1
63.546
Cu
Coordinates of Component 1
0.0 0.0 0.0 5 0

Cu
Velocities of Component 1
0.1 0.2 0.3 0 99
"#;
        let mut iter = ConFrameIterator::new(text.trim());
        let frame = iter.next().unwrap().unwrap();

        assert!(frame.has_velocities());
        assert_eq!(frame.atom_data[0].atom_id, 0);
        assert_eq!(frame.atom_data[0].fixed(), [true, false, true]);
    }

    #[test]
    fn test_validate_must_be_boolean_when_present() {
        let lines = [
            "PREBOX1",
            "{\"con_spec_version\":2,\"validate\":\"yes\",\"sections\":[]}",
            "10.0 20.0 30.0",
            "90.0 90.0 90.0",
            "POSTBOX1",
            "POSTBOX2",
            "1",
            "1",
            "12.011",
        ];
        let mut line_it = lines.iter().copied();
        let err = parse_frame_header(&mut line_it).unwrap_err();

        assert!(matches!(err, ParseError::InvalidMetadataJson(_)));
        assert!(err.to_string().contains("validate"));
    }

    #[test]
    fn test_validate_true_requires_sections_key() {
        let lines = [
            "PREBOX1",
            "{\"con_spec_version\":2,\"validate\":true}",
            "10.0 20.0 30.0",
            "90.0 90.0 90.0",
            "POSTBOX1",
            "POSTBOX2",
            "1",
            "1",
            "12.011",
        ];
        let mut line_it = lines.iter().copied();
        let err = parse_frame_header(&mut line_it).unwrap_err();

        assert!(matches!(err, ParseError::InvalidMetadataJson(_)));
        assert!(err.to_string().contains("sections"));
    }

    #[test]
    fn test_sections_must_be_string_array_when_present() {
        let lines = [
            "PREBOX1",
            "{\"con_spec_version\":2,\"sections\":[\"velocities\",7]}",
            "10.0 20.0 30.0",
            "90.0 90.0 90.0",
            "POSTBOX1",
            "POSTBOX2",
            "1",
            "1",
            "12.011",
        ];
        let mut line_it = lines.iter().copied();
        let err = parse_frame_header(&mut line_it).unwrap_err();

        assert!(matches!(err, ParseError::InvalidMetadataJson(_)));
        assert!(err.to_string().contains("sections"));
    }

    #[test]
    fn test_validate_true_rejects_non_integer_coordinate_identity_columns() {
        let text = r#"
PREBOX1
{"con_spec_version":2,"sections":[],"validate":true}
10.0 20.0 30.0
90.0 90.0 90.0
POSTBOX1
POSTBOX2
1
1
63.546
Cu
Coordinates of Component 1
0.0 0.0 0.0 5.0 0
"#;
        let mut iter = ConFrameIterator::new(text.trim());
        let err = iter.next().unwrap().unwrap_err();

        assert!(matches!(err, ParseError::InvalidFixedFlag(_)));
        assert!(err.to_string().contains("5.0"));
    }

    #[test]
    fn test_validate_true_rejects_non_exact_coordinate_label() {
        let text = r#"
PREBOX1
{"con_spec_version":2,"sections":[],"validate":true}
10.0 20.0 30.0
90.0 90.0 90.0
POSTBOX1
POSTBOX2
1
1
63.546
Cu
Coordinates Component 1
0.0 0.0 0.0 5 0
"#;
        let mut iter = ConFrameIterator::new(text.trim());
        let err = iter.next().unwrap().unwrap_err();

        assert!(matches!(err, ParseError::ValidationError(_)));
        assert!(err.to_string().contains("Coordinates of Component 1"));
    }

    #[test]
    fn test_validate_true_rejects_unknown_component_symbol() {
        let text = r#"
PREBOX1
{"con_spec_version":2,"sections":[],"validate":true}
10.0 20.0 30.0
90.0 90.0 90.0
POSTBOX1
POSTBOX2
1
1
63.546
Qq
Coordinates of Component 1
0.0 0.0 0.0 0 0
"#;
        let mut iter = ConFrameIterator::new(text.trim());
        let err = iter.next().unwrap().unwrap_err();

        assert!(matches!(err, ParseError::ValidationError(_)));
        assert!(err.to_string().contains("symbol"));
    }

    #[test]
    fn test_declared_section_must_be_present() {
        let text = r#"
PREBOX1
{"con_spec_version":2,"sections":["velocities"]}
10.0 20.0 30.0
90.0 90.0 90.0
POSTBOX1
POSTBOX2
1
1
63.546
Cu
Coordinates of Component 1
0.0 0.0 0.0 0 0
"#;
        let mut iter = ConFrameIterator::new(text.trim());
        let err = iter.next().unwrap().unwrap_err();

        assert!(matches!(
            err,
            ParseError::IncompleteSection(ref name) if name == "velocities"
        ));
    }

    #[test]
    fn test_non_finite_cell_geometry_rejected_in_strict_mode() {
        let lines = [
            "PREBOX1",
            "{\"con_spec_version\":2,\"sections\":[],\"validate\":true}",
            "10.0 NaN 30.0",
            "90.0 90.0 90.0",
            "POSTBOX1",
            "POSTBOX2",
            "1",
            "1",
            "12.011",
        ];
        let mut line_it = lines.iter().copied();
        let err = parse_frame_header(&mut line_it).unwrap_err();

        assert!(matches!(err, ParseError::NonFiniteValue(_)));
        assert!(err.to_string().contains("NaN"));
    }

    #[test]
    fn test_validate_true_rejects_non_physical_cell_geometry() {
        let lines = [
            "PREBOX1",
            "{\"con_spec_version\":2,\"sections\":[],\"validate\":true}",
            "0.0 20.0 30.0",
            "90.0 180.0 90.0",
            "POSTBOX1",
            "POSTBOX2",
            "1",
            "1",
            "12.011",
        ];
        let mut line_it = lines.iter().copied();
        let err = parse_frame_header(&mut line_it).unwrap_err();

        assert!(matches!(err, ParseError::ValidationError(_)));
        assert!(err.to_string().contains("cell"));
    }

    #[test]
    fn test_validate_true_rejects_reserved_metadata_type_mismatch() {
        let lines = [
            "PREBOX1",
            "{\"con_spec_version\":2,\"sections\":[],\"validate\":true,\"energy\":\"low\"}",
            "10.0 20.0 30.0",
            "90.0 90.0 90.0",
            "POSTBOX1",
            "POSTBOX2",
            "1",
            "1",
            "12.011",
        ];
        let mut line_it = lines.iter().copied();
        let err = parse_frame_header(&mut line_it).unwrap_err();

        assert!(matches!(err, ParseError::InvalidMetadataJson(_)));
        assert!(err.to_string().contains("energy"));
    }

    #[test]
    fn test_validate_true_rejects_malformed_bonds() {
        let lines = [
            "PREBOX1",
            "{\"con_spec_version\":2,\"sections\":[],\"validate\":true,\"bonds\":[[0]]}",
            "10.0 20.0 30.0",
            "90.0 90.0 90.0",
            "POSTBOX1",
            "POSTBOX2",
            "1",
            "1",
            "12.011",
        ];
        let mut line_it = lines.iter().copied();
        let err = parse_frame_header(&mut line_it).unwrap_err();
        assert!(matches!(err, ParseError::InvalidMetadataJson(_)));
        assert!(err.to_string().contains("bonds"));
    }

    #[test]
    fn test_bonds_metadata_round_trip_in_header() {
        use crate::types::{Bond, meta};
        let lines = [
            "PREBOX1",
            "{\"con_spec_version\":2,\"bonds\":[[0,1],{\"i\":0,\"j\":2,\"order\":1}]}",
            "10.0 20.0 30.0",
            "90.0 90.0 90.0",
            "POSTBOX1",
            "POSTBOX2",
            "1",
            "1",
            "12.011",
        ];
        let mut line_it = lines.iter().copied();
        let header = parse_frame_header(&mut line_it).expect("header");
        let bonds = header.bonds();
        assert_eq!(bonds.len(), 2);
        assert_eq!(bonds[0], Bond::new(0, 1));
        assert_eq!(bonds[1].i, 0);
        assert_eq!(bonds[1].j, 2);
        assert_eq!(bonds[1].order, Some(1));
        assert!(header.metadata.contains_key(meta::BONDS));
    }

    #[test]
    fn test_parse_velocity_section_absent() {
        let lines = [
            "PREBOX1",
            "{\"con_spec_version\":2}",
            "10.0 20.0 30.0",
            "90.0 90.0 90.0",
            "POSTBOX1",
            "POSTBOX2",
            "1",
            "1",
            "12.011",
            "C",
            "Coordinates of Component 1",
            "0.0 0.0 0.0 0 1",
        ];
        let mut line_it = lines.iter().copied().peekable();
        let mut frame = parse_single_frame(&mut line_it).expect("parse should succeed");
        let mut column = Vec::new();
        let has_vel = parse_velocity_section(
            &mut line_it,
            &frame.header,
            &mut frame.atom_data,
            &mut column,
        )
        .expect("should succeed with no velocities");
        assert!(!has_vel);
        assert!(column.is_empty());
        assert_eq!(frame.atom_data[0].velocity(), None);
    }

    /// The velocity column is the storage `atom.velocity` reads. Writing the
    /// atom updates that buffer in place, and the bytes match the file.
    #[test]
    fn velocity_column_aliases_soa_bytes() {
        let text = "\
bulk
{\"con_spec_version\":2}
4.0 5.0 6.0
90.0 90.0 90.0
0 0
0 0
1
2
63.546
Cu
Coordinates of Component 1
0.0 0.0 0.0 7 0
1.0 0.0 0.0 0 1

Cu
Velocities of Component 1
0.125 -0.25 0.5 7 0
-0.125 0.25 -0.5 0 1
";
        let mut frame = ConFrameIterator::new(text)
            .next()
            .expect("frame")
            .expect("parse");
        let soa = frame.velocities.f64_slice().expect("f64 column");
        let mut packed = Vec::with_capacity(soa.len());
        for atom in &frame.atom_data {
            packed.extend_from_slice(&atom.velocity().expect("velocity"));
        }
        assert_eq!(soa, packed.as_slice());
        assert_eq!(soa, &[0.125, -0.25, 0.5, -0.125, 0.25, -0.5][..]);
        assert_eq!(frame.atom_data[0].velocity_slot(), soa.as_ptr());
        frame.atom_data[0].set_velocity(Some([9.0, 8.0, 7.0]));
        assert_eq!(frame.velocities.as_f64_row(0), [9.0, 8.0, 7.0]);
        assert_eq!(frame.velocities.as_f64_row(1), [-0.125, 0.25, -0.5]);
    }

    #[test]
    fn force_column_aliases_soa_bytes() {
        let text = "\
bulk
{\"con_spec_version\":2,\"sections\":[\"forces\"]}
4.0 5.0 6.0
90.0 90.0 90.0
0 0
0 0
1
2
63.546
Cu
Coordinates of Component 1
0.0 0.0 0.0 0 0
1.0 1.0 1.0 7 1

Cu
Forces of Component 1
0.1 0.2 0.3 0 0
-0.4 0.5 -0.6 7 1
";
        let mut frame = ConFrameIterator::new(text)
            .next()
            .expect("frame")
            .expect("parse");
        let soa = frame.forces.f64_slice().expect("f64 column");
        let mut packed = Vec::with_capacity(soa.len());
        for atom in &frame.atom_data {
            packed.extend_from_slice(&atom.force().expect("force"));
        }
        assert_eq!(soa, packed.as_slice());
        assert_eq!(soa, &[0.1, 0.2, 0.3, -0.4, 0.5, -0.6][..]);
        assert_eq!(frame.atom_data[0].force_slot(), soa.as_ptr());
        frame.atom_data[1].set_force(Some([1.0, 2.0, 3.0]));
        assert_eq!(frame.forces.as_f64_row(1), [1.0, 2.0, 3.0]);
    }

    /// Million-atom `.convel` parse. Release builds must stay under the
    /// per-atom box path (about 0.7 s on this host).
    #[test]
    #[ignore = "times a million-atom convel parse"]
    fn microbench_convel_column_parse() {
        let n = 1_000_000usize;
        let mut text = String::with_capacity(n * 64);
        text.push_str("bench\n{\"con_spec_version\":2}\n4 5 6\n90 90 90\n0 0\n0 0\n1\n");
        text.push_str(&format!("{n}\n63.546\nCu\nCoordinates of Component 1\n"));
        for i in 0..n {
            text.push_str(&format!("0 0 0 0 {i}\n"));
        }
        text.push_str("\nCu\nVelocities of Component 1\n");
        for i in 0..n {
            let v = (i % 1000) as f64 * 0.001;
            text.push_str(&format!("{v} 0 0 0 {i}\n"));
        }
        let started = std::time::Instant::now();
        let frame = ConFrameIterator::new(&text)
            .next()
            .expect("frame")
            .expect("parse");
        let elapsed = started.elapsed().as_secs_f64();
        eprintln!("convel column parse {n} atoms {elapsed:.4} s");
        let soa = frame.velocities.f64_slice().expect("f64");
        assert_eq!(frame.atom_data[0].velocity_slot(), soa.as_ptr());
        assert_eq!(frame.atom_data[n - 1].velocity().unwrap()[0], 0.999);
        assert_eq!(soa[(n - 1) * 3], 0.999);
        if !cfg!(debug_assertions) {
            assert!(elapsed < 0.35, "convel parse {elapsed:.4} s");
        }
        let real = std::path::Path::new("bench/h2h/work/thread-reads/con-h2h.convel");
        if let Ok(bytes) = std::fs::read(real) {
            let started = std::time::Instant::now();
            let frame = ConFrameIterator::from_bytes(&bytes)
                .next()
                .expect("frame")
                .expect("parse");
            let elapsed = started.elapsed().as_secs_f64();
            eprintln!("real convel column parse {elapsed:.4} s");
            let soa = frame.velocities.f64_slice().expect("f64");
            assert_eq!(frame.atom_data[0].velocity_slot(), soa.as_ptr());
            assert_eq!(
                frame.atom_data[0].velocity().unwrap(),
                frame.velocities.as_f64_row(0)
            );
        }
    }

    #[test]
    fn test_parse_line_of_range_f64_exact() {
        let vals = parse_line_of_range_f64("1.0 2.0 3.0 0.0 42", 4, 5, &[0.0; 5]).unwrap();
        assert_eq!(vals, vec![1.0, 2.0, 3.0, 0.0, 42.0]);
    }

    #[test]
    fn test_byte_scan_stack_parses_realistic_coord_line() {
        // Typical CON atom line: x y z fixed_mask [atom_id]
        let line = "   0.63939999999999997    0.90449999999999997   -0.00009999999999977 1    0";
        let defaults = [0.0, 0.0, 0.0, 0.0, 99.0];
        let mut buf = [0.0f64; NUMERIC_STACK_WIDTH];
        let n = parse_line_of_range_f64_stack(line, 4, 5, &defaults, &mut buf).unwrap();
        assert_eq!(n, 5);
        assert!((buf[0] - 0.6394).abs() < 1e-4);
        assert!((buf[1] - 0.9045).abs() < 1e-4);
        assert_eq!(buf[3] as u8, 1);
        assert_eq!(buf[4] as u64, 0);
        // trailing junk must error
        let bad = "1.0 2.0 3.0 1 0 EXTRA";
        assert!(parse_line_of_range_f64_stack(bad, 4, 5, &defaults, &mut buf).is_err());
        // partial non-boundary token must error
        let glued = "1.0 2.0 3.0abc 1 0";
        assert!(parse_line_of_range_f64_stack(glued, 4, 5, &defaults, &mut buf).is_err());
    }

    #[test]
    fn test_parse_line_of_range_f64_stack_matches_vec_api() {
        let defaults = [0.0, 0.0, 0.0, 0.0, 7.0];
        let line = "1.0 2.0 3.0 0.0";
        let mut buf = [0.0f64; NUMERIC_STACK_WIDTH];
        parse_line_of_range_f64_stack(line, 4, 5, &defaults, &mut buf).unwrap();
        let via_vec = parse_line_of_range_f64(line, 4, 5, &defaults).unwrap();
        assert_eq!(&buf[..5], via_vec.as_slice());
        assert_eq!(buf[4], 7.0);
    }

    #[test]
    fn test_parse_line_of_range_f64_padded() {
        let defaults = [0.0, 0.0, 0.0, 0.0, 99.0];
        let vals = parse_line_of_range_f64("1.0 2.0 3.0 0.0", 4, 5, &defaults).unwrap();
        assert_eq!(vals, vec![1.0, 2.0, 3.0, 0.0, 99.0]);
    }

    #[test]
    fn test_parse_line_of_range_f64_too_few() {
        let result = parse_line_of_range_f64("1.0 2.0 3.0", 4, 5, &[0.0; 5]);
        assert!(result.is_err());
    }

    #[test]
    fn test_parse_line_of_range_f64_too_many() {
        let result = parse_line_of_range_f64("1.0 2.0 3.0 0.0 5.0 6.0", 4, 5, &[0.0; 5]);
        assert!(result.is_err());
    }

    #[test]
    fn test_parse_all_four_column_lines() {
        // All atom lines have only 4 columns; atom_index defaults to sequential.
        let lines = [
            "PREBOX1",
            "{\"con_spec_version\":2}",
            "10.0 10.0 10.0",
            "90.0 90.0 90.0",
            "POSTBOX1",
            "POSTBOX2",
            "1",
            "3",
            "12.011",
            "C",
            "Coordinates of Component 1",
            "0.0 0.0 0.0 0",
            "1.0 0.0 0.0 0",
            "2.0 0.0 0.0 1",
        ];
        let mut line_it = lines.iter().copied();
        let frame = parse_single_frame(&mut line_it).unwrap();
        assert_eq!(frame.atom_data[0].atom_id, 0);
        assert_eq!(frame.atom_data[1].atom_id, 1);
        assert_eq!(frame.atom_data[2].atom_id, 2);
        assert!(frame.atom_data[2].is_fixed());
    }

    #[test]
    fn test_positions_take_the_storage_dtype_of_the_metadata() {
        use crate::storage_dtype::{ElementKind, FloatArray2};
        let rows = [[0.5, -1.25, 3.0], [1.0e-3, 2.5e4, -7.75]];
        let atom_lines: Vec<String> = rows
            .iter()
            .enumerate()
            .map(|(id, [x, y, z])| format!("{x} {y} {z} 0 {id}"))
            .collect();
        for &kind in ElementKind::all_hosted() {
            let metadata = format!(
                "{{\"con_spec_version\":3,\"units\":{{\"length\":\"angstrom\",\"energy\":\"eV\"}},\
                 \"storage_dtypes\":{{\"positions\":\"{}\"}}}}",
                kind.as_str()
            );
            let mut lines = vec![
                "PREBOX1",
                metadata.as_str(),
                "10.0 10.0 10.0",
                "90.0 90.0 90.0",
                "POSTBOX1",
                "POSTBOX2",
                "1",
                "2",
                "1.008",
                "H",
                "Coordinates of Component 1",
            ];
            lines.extend(atom_lines.iter().map(String::as_str));
            let frame = parse_single_frame(&mut lines.iter().copied()).unwrap();
            assert_eq!(frame.positions.kind(), kind);
            assert_eq!(frame.positions.nrows(), rows.len());
            for (i, &xyz) in rows.iter().enumerate() {
                let mut expected = FloatArray2::zeros(kind, 1, 3);
                expected.set_f64_row(0, xyz);
                assert_eq!(
                    frame.positions.as_f64_row(i),
                    expected.as_f64_row(0),
                    "{kind:?} row {i}"
                );
            }
        }
    }
}

#[cfg(all(test, feature = "parallel"))]
mod newline_index_bench {
    use super::newline_index_at;

    #[test]
    fn index_matches_memchr() {
        let mut bytes = vec![b'x'; 80_000];
        for index in (20..bytes.len()).step_by(23) {
            bytes[index] = b'\n';
        }
        let serial: Vec<usize> = memchr::memchr_iter(b'\n', &bytes).collect();
        assert_eq!(newline_index_at(&bytes, 0, 8), serial);
        // Above the 1 MiB serial cutoff, including a buffer of bare newlines.
        let mut long = Vec::new();
        while long.len() < (1 << 20) + 64 {
            long.extend_from_slice(b"01234567 1.25 -2.5 3.5 0 7\n");
        }
        let long_serial: Vec<usize> = memchr::memchr_iter(b'\n', &long).collect();
        assert_eq!(newline_index_at(&long, 40, 8), {
            let mut shifted = long_serial.clone();
            for offset in &mut shifted {
                *offset += 40;
            }
            shifted
        });
        let bare = vec![b'\n'; (1 << 20) + 64];
        let bare_serial: Vec<usize> = memchr::memchr_iter(b'\n', &bare).collect();
        assert_eq!(newline_index_at(&bare, 0, 8), bare_serial);
    }
}