molcrafts-molrs 0.7.0

//! LAMMPS data file format reader and writer.
//!
//! Implements support for LAMMPS data files as specified in:
//! https://docs.lammps.org/read_data.html
//!
//! # Supported Features
//!
//! - Header section: atom/bond counts, box dimensions, type counts
//! - Atoms section: atomic, full, molecular, charge styles
//! - Bonds section
//! - Type Labels sections: Atom Type Labels, Bond Type Labels
//! - Orthorhombic and triclinic simulation boxes
//!
//! # Unsupported Features (can be added later)
//!
//! - Coeffs sections (Pair Coeffs, Bond Coeffs, etc.)
//! - Angles, Dihedrals, Impropers sections
//! - Velocities section
//! - Special atom styles (ellipsoid, sphere, etc.)
//!
//! # Examples
//!
//! ```no_run
//! use molrs::io::data::lammps_data::{read_lammps_data, write_lammps_data};
//!
//! # fn main() -> std::io::Result<()> {
//! // Read LAMMPS data file
//! let frame = read_lammps_data("system.data")?;
//!
//! // Access atoms
//! let atoms = frame.get("atoms").expect("atoms block");
//! let x = atoms.get_float("x").expect("x coordinates");
//!
//! // Write to new file
//! write_lammps_data("output.data", &frame)?;
//! # Ok(())
//! # }```

use crate::io::reader::{FrameReader, Reader};
use crate::io::writer::FrameWriter;
use molrs::spatial::region::simbox::SimBox;
use molrs::store::block::Block;
use molrs::store::frame::Frame;
use molrs::store::frame_access::FrameAccess;
use molrs::types::{F, I, Pbc3, U};
use ndarray::{Array1, IxDyn, array};
use once_cell::sync::OnceCell;
use std::collections::HashMap;
use std::fs::File;
use std::io::{BufRead, BufReader, Seek, SeekFrom, Write};
use std::path::Path;

// ============================================================================
// Helper Functions
// ============================================================================

/// Parse a line into whitespace-separated tokens
fn tokenize(line: &str) -> Vec<&str> {
    line.split_whitespace().collect()
}

fn err_mapper<E: std::fmt::Display>(e: E) -> std::io::Error {
    std::io::Error::new(std::io::ErrorKind::InvalidData, e.to_string())
}

// ============================================================================
// Data Structures
// ============================================================================

/// Header information from LAMMPS data file
#[derive(Debug, Clone, Default)]
struct LAMMPSHeader {
    num_atoms: usize,
    num_bonds: usize,
    num_angles: usize,
    num_dihedrals: usize,
    num_impropers: usize,
    num_atom_types: usize,
    num_bond_types: usize,
    num_angle_types: usize,
    num_dihedral_types: usize,
    num_improper_types: usize,
    xlo: f64,
    xhi: f64,
    ylo: f64,
    yhi: f64,
    zlo: f64,
    zhi: f64,
    xy: Option<f64>,
    xz: Option<f64>,
    yz: Option<f64>,
}

/// Atom data from Atoms section
#[derive(Debug, Clone)]
struct AtomData {
    id: I,
    molecule_id: Option<I>,
    atom_type: String, // Can be numeric or label
    charge: Option<F>,
    bodyflag: Option<I>, // atom_style="body" only
    mass: Option<F>,     // atom_style="body" only (per-atom mass)
    x: F,
    y: F,
    z: F,
}

/// Bond data from Bonds section
#[derive(Debug, Clone)]
struct BondData {
    _id: I,
    bond_type: String, // Can be numeric or label
    atom_i: I,
    atom_j: I,
}

/// Angle data from Angles section
#[derive(Debug, Clone)]
struct AngleData {
    _id: I,
    angle_type: String, // Can be numeric or label
    atom_i: I,
    atom_j: I,
    atom_k: I,
}

/// Dihedral data from Dihedrals section
#[derive(Debug, Clone)]
struct DihedralData {
    _id: I,
    dihedral_type: String, // Can be numeric or label
    atom_i: I,
    atom_j: I,
    atom_k: I,
    atom_l: I,
}

/// Improper data from Impropers section
#[derive(Debug, Clone)]
struct ImproperData {
    _id: I,
    improper_type: String, // Can be numeric or label
    atom_i: I,
    atom_j: I,
    atom_k: I,
    atom_l: I,
}

// ============================================================================
// Parser Functions
// ============================================================================

/// Parse header section and return the first section line encountered
fn parse_header_with_first_section<R: BufRead>(
    reader: &mut R,
) -> std::io::Result<(LAMMPSHeader, Option<String>)> {
    let mut header = LAMMPSHeader::default();
    let mut line = String::new();

    // Skip first line (comment)
    reader.read_line(&mut line)?;
    line.clear();

    loop {
        let bytes_read = reader.read_line(&mut line)?;
        if bytes_read == 0 {
            return Ok((header, None));
        }

        let trimmed = line.trim();

        // Skip blank lines and comments
        if trimmed.is_empty() || trimmed.starts_with('#') {
            line.clear();
            continue;
        }

        // Check if this looks like a section header (capitalized word followed by optional comment)
        // Section headers: "Atoms", "Bonds", "Atom Type Labels", etc.
        // Header keywords end with lowercase: "atoms", "bonds", "xlo xhi", etc.
        let tokens = tokenize(trimmed);
        if tokens.is_empty() {
            line.clear();
            continue;
        }

        // If first token starts with uppercase and is not followed by a number, it's likely a section
        if let Some(first_char) = tokens[0].chars().next()
            && first_char.is_uppercase()
        {
            // This is a section header, return it along with the header
            return Ok((header, Some(line.clone())));
        }

        // Parse header keywords
        match tokens.last() {
            Some(&"atoms") if tokens.len() >= 2 => {
                header.num_atoms = tokens[0].parse().map_err(err_mapper)?;
            }
            Some(&"bonds") if tokens.len() >= 2 => {
                header.num_bonds = tokens[0].parse().map_err(err_mapper)?;
            }
            Some(&"angles") if tokens.len() >= 2 => {
                header.num_angles = tokens[0].parse().map_err(err_mapper)?;
            }
            Some(&"dihedrals") if tokens.len() >= 2 => {
                header.num_dihedrals = tokens[0].parse().map_err(err_mapper)?;
            }
            Some(&"impropers") if tokens.len() >= 2 => {
                header.num_impropers = tokens[0].parse().map_err(err_mapper)?;
            }
            Some(&"types") if tokens.len() >= 3 => {
                let count: usize = tokens[0].parse().map_err(err_mapper)?;
                match tokens[1] {
                    "atom" => header.num_atom_types = count,
                    "bond" => header.num_bond_types = count,
                    "angle" => header.num_angle_types = count,
                    "dihedral" => header.num_dihedral_types = count,
                    "improper" => header.num_improper_types = count,
                    _ => {}
                }
            }
            Some(&"xhi") if tokens.len() >= 4 && tokens[2] == "xlo" => {
                header.xlo = tokens[0].parse().map_err(err_mapper)?;
                header.xhi = tokens[1].parse().map_err(err_mapper)?;
            }
            Some(&"yhi") if tokens.len() >= 4 && tokens[2] == "ylo" => {
                header.ylo = tokens[0].parse().map_err(err_mapper)?;
                header.yhi = tokens[1].parse().map_err(err_mapper)?;
            }
            Some(&"zhi") if tokens.len() >= 4 && tokens[2] == "zlo" => {
                header.zlo = tokens[0].parse().map_err(err_mapper)?;
                header.zhi = tokens[1].parse().map_err(err_mapper)?;
            }
            Some(&"yz") if tokens.len() >= 6 && tokens[3] == "xy" && tokens[4] == "xz" => {
                header.xy = Some(tokens[0].parse().map_err(err_mapper)?);
                header.xz = Some(tokens[1].parse().map_err(err_mapper)?);
                header.yz = Some(tokens[2].parse().map_err(err_mapper)?);
            }
            _ => {}
        }

        line.clear();
    }
}

/// Parse Type Labels section and return the next section line if encountered
fn parse_type_labels<R: BufRead>(
    reader: &mut R,
) -> std::io::Result<(HashMap<String, String>, Option<String>)> {
    let mut labels = HashMap::new();
    let mut line = String::new();

    loop {
        let bytes_read = reader.read_line(&mut line)?;
        if bytes_read == 0 {
            return Ok((labels, None));
        }

        let trimmed = line.trim();

        // Stop at blank line or new section
        if trimmed.is_empty() {
            line.clear();
            continue;
        }

        // Check if this is a new section header
        if trimmed.chars().next().is_some_and(|c| c.is_uppercase()) {
            // Check if it looks like a section header (not a type label)
            let tokens = tokenize(trimmed);
            if tokens.len() >= 2 {
                // Could be either "Atoms" or "1 C" (type label)
                // Type labels have numeric first token
                if tokens[0].parse::<i64>().is_err() {
                    // First token is not a number, so it's a section header
                    return Ok((labels, Some(line.clone())));
                }
            } else if tokens.len() == 1 {
                // Single word starting with uppercase is a section header
                return Ok((labels, Some(line.clone())));
            }
        }

        let tokens = tokenize(trimmed);
        if tokens.len() >= 2 {
            let type_id = tokens[0].to_string();
            let label = tokens[1].to_string();
            labels.insert(type_id, label);
        }

        line.clear();
    }
}

/// Parse the `# <style>` hint after an "Atoms" section header line, if any.
fn parse_atoms_style_hint(section_line: &str) -> Option<String> {
    section_line
        .split_once('#')
        .map(|(_, after)| after.split_whitespace().next().unwrap_or("").to_string())
        .filter(|s| !s.is_empty())
}

/// Parse Atoms section
fn parse_atoms<R: BufRead>(
    reader: &mut R,
    num_atoms: usize,
    _atom_type_labels: &HashMap<String, String>,
    style_hint: Option<&str>,
) -> std::io::Result<Vec<AtomData>> {
    let mut atoms = Vec::with_capacity(num_atoms);
    let mut line = String::new();

    while atoms.len() < num_atoms {
        line.clear();
        let bytes_read = reader.read_line(&mut line)?;
        if bytes_read == 0 {
            break;
        }

        let trimmed = line.trim();
        if trimmed.is_empty() || trimmed.starts_with('#') {
            continue;
        }

        let tokens = tokenize(trimmed);

        // atom_style="body" columns: id type bodyflag mass x y z (7), with
        // optional trailing image flags. Body atoms must not be confused with
        // atom_style="full" (id mol type q x y z) which is also 7 columns,
        // so we rely on the section comment hint.
        if matches!(style_hint, Some("body")) {
            if tokens.len() < 7 {
                return Err(err_mapper(format!(
                    "Invalid body-style Atoms line: expected at least 7 columns, got {}",
                    tokens.len()
                )));
            }
            atoms.push(AtomData {
                id: tokens[0].parse::<I>().map_err(err_mapper)?,
                molecule_id: None,
                atom_type: tokens[1].to_string(),
                charge: None,
                bodyflag: Some(tokens[2].parse::<I>().map_err(err_mapper)?),
                mass: Some(tokens[3].parse::<F>().map_err(err_mapper)?),
                x: tokens[4].parse::<F>().map_err(err_mapper)?,
                y: tokens[5].parse::<F>().map_err(err_mapper)?,
                z: tokens[6].parse::<F>().map_err(err_mapper)?,
            });
            line.clear();
            continue;
        }

        // Detect atom style by number of columns
        // atomic: id type x y z (5)
        // charge: id type q x y z (6)
        // molecular: id mol-id type x y z (6)
        // full: id mol-id type q x y z (7)

        let atom = if tokens.len() == 5 {
            // atomic style
            AtomData {
                id: tokens[0].parse::<I>().map_err(err_mapper)?,
                molecule_id: None,
                atom_type: tokens[1].to_string(),
                charge: None,
                bodyflag: None,
                mass: None,
                x: tokens[2].parse::<F>().map_err(err_mapper)?,
                y: tokens[3].parse::<F>().map_err(err_mapper)?,
                z: tokens[4].parse::<F>().map_err(err_mapper)?,
            }
        } else if tokens.len() == 6 {
            // Could be charge or molecular - try to detect
            // If token[2] parses as float, it's charge style
            if tokens[2].parse::<F>().is_ok() {
                // charge style
                AtomData {
                    id: tokens[0].parse::<I>().map_err(err_mapper)?,
                    molecule_id: None,
                    atom_type: tokens[1].to_string(),
                    charge: Some(tokens[2].parse::<F>().map_err(err_mapper)?),
                    bodyflag: None,
                    mass: None,
                    x: tokens[3].parse::<F>().map_err(err_mapper)?,
                    y: tokens[4].parse::<F>().map_err(err_mapper)?,
                    z: tokens[5].parse::<F>().map_err(err_mapper)?,
                }
            } else {
                // molecular style
                AtomData {
                    id: tokens[0].parse::<I>().map_err(err_mapper)?,
                    molecule_id: Some(tokens[1].parse::<I>().map_err(err_mapper)?),
                    atom_type: tokens[2].to_string(),
                    charge: None,
                    bodyflag: None,
                    mass: None,
                    x: tokens[3].parse::<F>().map_err(err_mapper)?,
                    y: tokens[4].parse::<F>().map_err(err_mapper)?,
                    z: tokens[5].parse::<F>().map_err(err_mapper)?,
                }
            }
        } else if tokens.len() >= 7 {
            // full style
            AtomData {
                id: tokens[0].parse::<I>().map_err(err_mapper)?,
                molecule_id: Some(tokens[1].parse::<I>().map_err(err_mapper)?),
                atom_type: tokens[2].to_string(),
                charge: Some(tokens[3].parse::<F>().map_err(err_mapper)?),
                bodyflag: None,
                mass: None,
                x: tokens[4].parse::<F>().map_err(err_mapper)?,
                y: tokens[5].parse::<F>().map_err(err_mapper)?,
                z: tokens[6].parse::<F>().map_err(err_mapper)?,
            }
        } else {
            return Err(err_mapper(format!(
                "Invalid Atoms line: expected 5-7 columns, got {}",
                tokens.len()
            )));
        };

        atoms.push(atom);
        line.clear();
    }

    Ok(atoms)
}

/// Parse Bonds section
fn parse_bonds<R: BufRead>(
    reader: &mut R,
    num_bonds: usize,
    _bond_type_labels: &HashMap<String, String>,
) -> std::io::Result<Vec<BondData>> {
    let mut bonds = Vec::with_capacity(num_bonds);
    let mut line = String::new();

    while bonds.len() < num_bonds {
        line.clear();
        let bytes_read = reader.read_line(&mut line)?;
        if bytes_read == 0 {
            break;
        }

        let trimmed = line.trim();
        if trimmed.is_empty() || trimmed.starts_with('#') {
            continue;
        }

        let tokens = tokenize(trimmed);
        if tokens.len() < 4 {
            return Err(err_mapper(format!(
                "Invalid Bonds line: expected 4 columns, got {}",
                tokens.len()
            )));
        }

        let bond = BondData {
            _id: tokens[0].parse::<I>().map_err(err_mapper)?,
            bond_type: tokens[1].to_string(),
            atom_i: tokens[2].parse::<I>().map_err(err_mapper)?,
            atom_j: tokens[3].parse::<I>().map_err(err_mapper)?,
        };

        bonds.push(bond);
        line.clear();
    }

    Ok(bonds)
}

/// Parse Angles section
fn parse_angles<R: BufRead>(
    reader: &mut R,
    num_angles: usize,
    _angle_type_labels: &HashMap<String, String>,
) -> std::io::Result<Vec<AngleData>> {
    let mut angles = Vec::with_capacity(num_angles);
    let mut line = String::new();

    while angles.len() < num_angles {
        line.clear();
        let bytes_read = reader.read_line(&mut line)?;
        if bytes_read == 0 {
            break;
        }

        let trimmed = line.trim();
        if trimmed.is_empty() || trimmed.starts_with('#') {
            continue;
        }

        let tokens = tokenize(trimmed);
        if tokens.len() < 5 {
            return Err(err_mapper(format!(
                "Invalid Angles line: expected 5 columns, got {}",
                tokens.len()
            )));
        }

        let angle = AngleData {
            _id: tokens[0].parse::<I>().map_err(err_mapper)?,
            angle_type: tokens[1].to_string(),
            atom_i: tokens[2].parse::<I>().map_err(err_mapper)?,
            atom_j: tokens[3].parse::<I>().map_err(err_mapper)?,
            atom_k: tokens[4].parse::<I>().map_err(err_mapper)?,
        };

        angles.push(angle);
    }

    Ok(angles)
}

/// Parse Dihedrals section
fn parse_dihedrals<R: BufRead>(
    reader: &mut R,
    num_dihedrals: usize,
    _dihedral_type_labels: &HashMap<String, String>,
) -> std::io::Result<Vec<DihedralData>> {
    let mut dihedrals = Vec::with_capacity(num_dihedrals);
    let mut line = String::new();

    while dihedrals.len() < num_dihedrals {
        line.clear();
        let bytes_read = reader.read_line(&mut line)?;
        if bytes_read == 0 {
            break;
        }

        let trimmed = line.trim();
        if trimmed.is_empty() || trimmed.starts_with('#') {
            continue;
        }

        let tokens = tokenize(trimmed);
        if tokens.len() < 6 {
            return Err(err_mapper(format!(
                "Invalid Dihedrals line: expected 6 columns, got {}",
                tokens.len()
            )));
        }

        let dihedral = DihedralData {
            _id: tokens[0].parse::<I>().map_err(err_mapper)?,
            dihedral_type: tokens[1].to_string(),
            atom_i: tokens[2].parse::<I>().map_err(err_mapper)?,
            atom_j: tokens[3].parse::<I>().map_err(err_mapper)?,
            atom_k: tokens[4].parse::<I>().map_err(err_mapper)?,
            atom_l: tokens[5].parse::<I>().map_err(err_mapper)?,
        };

        dihedrals.push(dihedral);
    }

    Ok(dihedrals)
}

/// Parse Impropers section. Column layout matches dihedrals:
/// `id type i j k l` with `j` the central atom by LAMMPS convention.
fn parse_impropers<R: BufRead>(
    reader: &mut R,
    num_impropers: usize,
    _improper_type_labels: &HashMap<String, String>,
) -> std::io::Result<Vec<ImproperData>> {
    let mut impropers = Vec::with_capacity(num_impropers);
    let mut line = String::new();

    while impropers.len() < num_impropers {
        line.clear();
        let bytes_read = reader.read_line(&mut line)?;
        if bytes_read == 0 {
            break;
        }

        let trimmed = line.trim();
        if trimmed.is_empty() || trimmed.starts_with('#') {
            continue;
        }

        let tokens = tokenize(trimmed);
        if tokens.len() < 6 {
            return Err(err_mapper(format!(
                "Invalid Impropers line: expected 6 columns, got {}",
                tokens.len()
            )));
        }

        impropers.push(ImproperData {
            _id: tokens[0].parse::<I>().map_err(err_mapper)?,
            improper_type: tokens[1].to_string(),
            atom_i: tokens[2].parse::<I>().map_err(err_mapper)?,
            atom_j: tokens[3].parse::<I>().map_err(err_mapper)?,
            atom_k: tokens[4].parse::<I>().map_err(err_mapper)?,
            atom_l: tokens[5].parse::<I>().map_err(err_mapper)?,
        });
    }

    Ok(impropers)
}

/// Build Frame from parsed data
#[allow(clippy::too_many_arguments)]
fn build_frame(
    header: &LAMMPSHeader,
    atoms: Vec<AtomData>,
    bonds: Vec<BondData>,
    angles: Vec<AngleData>,
    dihedrals: Vec<DihedralData>,
    impropers: Vec<ImproperData>,
    atom_type_labels: HashMap<String, String>,
    bond_type_labels: HashMap<String, String>,
    angle_type_labels: HashMap<String, String>,
    dihedral_type_labels: HashMap<String, String>,
    improper_type_labels: HashMap<String, String>,
) -> std::io::Result<Frame> {
    let mut frame = Frame::new();

    // Build atoms block
    if !atoms.is_empty() {
        let n = atoms.len();
        let mut atom_block = Block::new();

        let mut x_vec = Vec::with_capacity(n);
        let mut y_vec = Vec::with_capacity(n);
        let mut z_vec = Vec::with_capacity(n);
        let mut id_vec = Vec::with_capacity(n);
        let mut type_vec = Vec::with_capacity(n);
        let mut charge_vec = Vec::with_capacity(n);
        let mut mol_id_vec = Vec::with_capacity(n);
        let mut bodyflag_vec = Vec::with_capacity(n);
        let mut mass_vec = Vec::with_capacity(n);

        let mut has_charge = false;
        let mut has_mol_id = false;
        let mut has_bodyflag = false;
        let mut has_mass = false;

        // Build atom ID to index map for bonds
        let mut atom_id_map = HashMap::new();

        for (idx, atom) in atoms.iter().enumerate() {
            x_vec.push(atom.x);
            y_vec.push(atom.y);
            z_vec.push(atom.z);
            id_vec.push(atom.id);

            // Convert type label to numeric if it's a label
            let type_num = if let Ok(num) = atom.atom_type.parse::<I>() {
                num
            } else {
                // It's a label, find the numeric type
                atom_type_labels
                    .iter()
                    .find(|(_, label)| *label == &atom.atom_type)
                    .and_then(|(id, _)| id.parse::<I>().ok())
                    .unwrap_or(1)
            };
            type_vec.push(type_num);

            if let Some(q) = atom.charge {
                charge_vec.push(q);
                has_charge = true;
            } else {
                charge_vec.push(0.0 as F);
            }

            if let Some(mol_id) = atom.molecule_id {
                mol_id_vec.push(mol_id);
                has_mol_id = true;
            } else {
                mol_id_vec.push(0);
            }

            if let Some(bf) = atom.bodyflag {
                bodyflag_vec.push(bf);
                has_bodyflag = true;
            } else {
                bodyflag_vec.push(0);
            }

            if let Some(m) = atom.mass {
                mass_vec.push(m);
                has_mass = true;
            } else {
                mass_vec.push(0.0 as F);
            }

            atom_id_map.insert(atom.id, idx as U);
        }

        // Insert coordinate arrays
        let to_arr_float = |v: Vec<F>| -> std::io::Result<ndarray::ArrayD<F>> {
            Ok(Array1::from_vec(v)
                .into_shape_with_order(IxDyn(&[n]))
                .map_err(err_mapper)?
                .into_dyn())
        };
        let to_arr_int = |v: Vec<I>| -> std::io::Result<ndarray::ArrayD<I>> {
            Ok(Array1::from_vec(v)
                .into_shape_with_order(IxDyn(&[n]))
                .map_err(err_mapper)?
                .into_dyn())
        };

        atom_block
            .insert("x", to_arr_float(x_vec)?)
            .map_err(err_mapper)?;
        atom_block
            .insert("y", to_arr_float(y_vec)?)
            .map_err(err_mapper)?;
        atom_block
            .insert("z", to_arr_float(z_vec)?)
            .map_err(err_mapper)?;
        atom_block
            .insert("id", to_arr_int(id_vec)?)
            .map_err(err_mapper)?;
        atom_block
            .insert("type", to_arr_int(type_vec)?)
            .map_err(err_mapper)?;

        if has_charge {
            atom_block
                .insert("charge", to_arr_float(charge_vec)?)
                .map_err(err_mapper)?;
        }
        if has_mol_id {
            atom_block
                .insert("molecule_id", to_arr_int(mol_id_vec)?)
                .map_err(err_mapper)?;
        }
        if has_bodyflag {
            atom_block
                .insert("bodyflag", to_arr_int(bodyflag_vec)?)
                .map_err(err_mapper)?;
        }
        if has_mass {
            atom_block
                .insert("mass", to_arr_float(mass_vec)?)
                .map_err(err_mapper)?;
        }

        frame.insert("atoms", atom_block);

        // Build bonds block
        if !bonds.is_empty() {
            let bn = bonds.len();
            let mut bond_block = Block::new();

            let mut atom_i_vec = Vec::with_capacity(bn);
            let mut atom_j_vec = Vec::with_capacity(bn);
            let mut bond_type_vec = Vec::with_capacity(bn);

            for bond in &bonds {
                // Map atom IDs to indices
                let idx_i = *atom_id_map.get(&bond.atom_i).ok_or_else(|| {
                    err_mapper(format!("Bond references unknown atom ID: {}", bond.atom_i))
                })?;
                let idx_j = *atom_id_map.get(&bond.atom_j).ok_or_else(|| {
                    err_mapper(format!("Bond references unknown atom ID: {}", bond.atom_j))
                })?;

                atom_i_vec.push(idx_i);
                atom_j_vec.push(idx_j);

                // Convert bond type label to numeric if it's a label
                let type_num = if let Ok(num) = bond.bond_type.parse::<I>() {
                    num
                } else {
                    bond_type_labels
                        .iter()
                        .find(|(_, label)| *label == &bond.bond_type)
                        .and_then(|(id, _)| id.parse::<I>().ok())
                        .unwrap_or(1)
                };
                bond_type_vec.push(type_num);
            }

            let mk_barr_u = |v: Vec<U>| -> std::io::Result<ndarray::ArrayD<U>> {
                Array1::from_vec(v)
                    .into_shape_with_order(IxDyn(&[bn]))
                    .map_err(err_mapper)
                    .map(|a| a.into_dyn())
            };
            let mk_barr_i = |v: Vec<I>| -> std::io::Result<ndarray::ArrayD<I>> {
                Array1::from_vec(v)
                    .into_shape_with_order(IxDyn(&[bn]))
                    .map_err(err_mapper)
                    .map(|a| a.into_dyn())
            };

            bond_block
                .insert("atomi", mk_barr_u(atom_i_vec)?)
                .map_err(err_mapper)?;
            bond_block
                .insert("atomj", mk_barr_u(atom_j_vec)?)
                .map_err(err_mapper)?;
            bond_block
                .insert("type", mk_barr_i(bond_type_vec)?)
                .map_err(err_mapper)?;

            frame.insert("bonds", bond_block);
        }

        // Build angles block
        if !angles.is_empty() {
            let an = angles.len();
            let mut angle_block = Block::new();

            let mut atom_i_vec = Vec::with_capacity(an);
            let mut atom_j_vec = Vec::with_capacity(an);
            let mut atom_k_vec = Vec::with_capacity(an);
            let mut angle_type_vec = Vec::with_capacity(an);

            for angle in &angles {
                let idx_i = *atom_id_map.get(&angle.atom_i).ok_or_else(|| {
                    err_mapper(format!(
                        "Angle references unknown atom ID: {}",
                        angle.atom_i
                    ))
                })?;
                let idx_j = *atom_id_map.get(&angle.atom_j).ok_or_else(|| {
                    err_mapper(format!(
                        "Angle references unknown atom ID: {}",
                        angle.atom_j
                    ))
                })?;
                let idx_k = *atom_id_map.get(&angle.atom_k).ok_or_else(|| {
                    err_mapper(format!(
                        "Angle references unknown atom ID: {}",
                        angle.atom_k
                    ))
                })?;

                atom_i_vec.push(idx_i);
                atom_j_vec.push(idx_j);
                atom_k_vec.push(idx_k);

                let type_num = if let Ok(num) = angle.angle_type.parse::<I>() {
                    num
                } else {
                    angle_type_labels
                        .iter()
                        .find(|(_, label)| *label == &angle.angle_type)
                        .and_then(|(id, _)| id.parse::<I>().ok())
                        .unwrap_or(1)
                };
                angle_type_vec.push(type_num);
            }

            let mk_aarr_u = |v: Vec<U>| -> std::io::Result<ndarray::ArrayD<U>> {
                Array1::from_vec(v)
                    .into_shape_with_order(IxDyn(&[an]))
                    .map_err(err_mapper)
                    .map(|a| a.into_dyn())
            };
            let mk_aarr_i = |v: Vec<I>| -> std::io::Result<ndarray::ArrayD<I>> {
                Array1::from_vec(v)
                    .into_shape_with_order(IxDyn(&[an]))
                    .map_err(err_mapper)
                    .map(|a| a.into_dyn())
            };

            angle_block
                .insert("atomi", mk_aarr_u(atom_i_vec)?)
                .map_err(err_mapper)?;
            angle_block
                .insert("atomj", mk_aarr_u(atom_j_vec)?)
                .map_err(err_mapper)?;
            angle_block
                .insert("atomk", mk_aarr_u(atom_k_vec)?)
                .map_err(err_mapper)?;
            angle_block
                .insert("type", mk_aarr_i(angle_type_vec)?)
                .map_err(err_mapper)?;

            frame.insert("angles", angle_block);
        }

        // Build dihedrals block
        if !dihedrals.is_empty() {
            let dn = dihedrals.len();
            let mut dihedral_block = Block::new();

            let mut atom_i_vec = Vec::with_capacity(dn);
            let mut atom_j_vec = Vec::with_capacity(dn);
            let mut atom_k_vec = Vec::with_capacity(dn);
            let mut atom_l_vec = Vec::with_capacity(dn);
            let mut dihedral_type_vec = Vec::with_capacity(dn);

            for dihedral in &dihedrals {
                let idx_i = *atom_id_map.get(&dihedral.atom_i).ok_or_else(|| {
                    err_mapper(format!(
                        "Dihedral references unknown atom ID: {}",
                        dihedral.atom_i
                    ))
                })?;
                let idx_j = *atom_id_map.get(&dihedral.atom_j).ok_or_else(|| {
                    err_mapper(format!(
                        "Dihedral references unknown atom ID: {}",
                        dihedral.atom_j
                    ))
                })?;
                let idx_k = *atom_id_map.get(&dihedral.atom_k).ok_or_else(|| {
                    err_mapper(format!(
                        "Dihedral references unknown atom ID: {}",
                        dihedral.atom_k
                    ))
                })?;
                let idx_l = *atom_id_map.get(&dihedral.atom_l).ok_or_else(|| {
                    err_mapper(format!(
                        "Dihedral references unknown atom ID: {}",
                        dihedral.atom_l
                    ))
                })?;

                atom_i_vec.push(idx_i);
                atom_j_vec.push(idx_j);
                atom_k_vec.push(idx_k);
                atom_l_vec.push(idx_l);

                let type_num = if let Ok(num) = dihedral.dihedral_type.parse::<I>() {
                    num
                } else {
                    dihedral_type_labels
                        .iter()
                        .find(|(_, label)| *label == &dihedral.dihedral_type)
                        .and_then(|(id, _)| id.parse::<I>().ok())
                        .unwrap_or(1)
                };
                dihedral_type_vec.push(type_num);
            }

            let mk_darr_u = |v: Vec<U>| -> std::io::Result<ndarray::ArrayD<U>> {
                Array1::from_vec(v)
                    .into_shape_with_order(IxDyn(&[dn]))
                    .map_err(err_mapper)
                    .map(|a| a.into_dyn())
            };
            let mk_darr_i = |v: Vec<I>| -> std::io::Result<ndarray::ArrayD<I>> {
                Array1::from_vec(v)
                    .into_shape_with_order(IxDyn(&[dn]))
                    .map_err(err_mapper)
                    .map(|a| a.into_dyn())
            };

            dihedral_block
                .insert("atomi", mk_darr_u(atom_i_vec)?)
                .map_err(err_mapper)?;
            dihedral_block
                .insert("atomj", mk_darr_u(atom_j_vec)?)
                .map_err(err_mapper)?;
            dihedral_block
                .insert("atomk", mk_darr_u(atom_k_vec)?)
                .map_err(err_mapper)?;
            dihedral_block
                .insert("atoml", mk_darr_u(atom_l_vec)?)
                .map_err(err_mapper)?;
            dihedral_block
                .insert("type", mk_darr_i(dihedral_type_vec)?)
                .map_err(err_mapper)?;

            frame.insert("dihedrals", dihedral_block);
        }

        // Build impropers block
        if !impropers.is_empty() {
            let in_ = impropers.len();
            let mut improper_block = Block::new();

            let mut atom_i_vec = Vec::with_capacity(in_);
            let mut atom_j_vec = Vec::with_capacity(in_);
            let mut atom_k_vec = Vec::with_capacity(in_);
            let mut atom_l_vec = Vec::with_capacity(in_);
            let mut improper_type_vec = Vec::with_capacity(in_);

            for improper in &impropers {
                let idx_i = *atom_id_map.get(&improper.atom_i).ok_or_else(|| {
                    err_mapper(format!(
                        "Improper references unknown atom ID: {}",
                        improper.atom_i
                    ))
                })?;
                let idx_j = *atom_id_map.get(&improper.atom_j).ok_or_else(|| {
                    err_mapper(format!(
                        "Improper references unknown atom ID: {}",
                        improper.atom_j
                    ))
                })?;
                let idx_k = *atom_id_map.get(&improper.atom_k).ok_or_else(|| {
                    err_mapper(format!(
                        "Improper references unknown atom ID: {}",
                        improper.atom_k
                    ))
                })?;
                let idx_l = *atom_id_map.get(&improper.atom_l).ok_or_else(|| {
                    err_mapper(format!(
                        "Improper references unknown atom ID: {}",
                        improper.atom_l
                    ))
                })?;

                atom_i_vec.push(idx_i);
                atom_j_vec.push(idx_j);
                atom_k_vec.push(idx_k);
                atom_l_vec.push(idx_l);

                let type_num = if let Ok(num) = improper.improper_type.parse::<I>() {
                    num
                } else {
                    improper_type_labels
                        .iter()
                        .find(|(_, label)| *label == &improper.improper_type)
                        .and_then(|(id, _)| id.parse::<I>().ok())
                        .unwrap_or(1)
                };
                improper_type_vec.push(type_num);
            }

            let mk_iarr_u = |v: Vec<U>| -> std::io::Result<ndarray::ArrayD<U>> {
                Array1::from_vec(v)
                    .into_shape_with_order(IxDyn(&[in_]))
                    .map_err(err_mapper)
                    .map(|a| a.into_dyn())
            };
            let mk_iarr_i = |v: Vec<I>| -> std::io::Result<ndarray::ArrayD<I>> {
                Array1::from_vec(v)
                    .into_shape_with_order(IxDyn(&[in_]))
                    .map_err(err_mapper)
                    .map(|a| a.into_dyn())
            };

            improper_block
                .insert("atomi", mk_iarr_u(atom_i_vec)?)
                .map_err(err_mapper)?;
            improper_block
                .insert("atomj", mk_iarr_u(atom_j_vec)?)
                .map_err(err_mapper)?;
            improper_block
                .insert("atomk", mk_iarr_u(atom_k_vec)?)
                .map_err(err_mapper)?;
            improper_block
                .insert("atoml", mk_iarr_u(atom_l_vec)?)
                .map_err(err_mapper)?;
            improper_block
                .insert("type", mk_iarr_i(improper_type_vec)?)
                .map_err(err_mapper)?;

            frame.insert("impropers", improper_block);
        }
    }

    // Build SimBox from the header when the file actually declares a box.
    // Origin = [xlo, ylo, zlo]; box vectors derive from hi-lo lengths plus
    // optional triclinic tilt factors (xy, xz, yz). LAMMPS data files carry no
    // boundary flags — assume fully periodic.
    let lx = header.xhi - header.xlo;
    let ly = header.yhi - header.ylo;
    let lz = header.zhi - header.zlo;
    if lx > 0.0 && ly > 0.0 && lz > 0.0 {
        let origin = array![header.xlo, header.ylo, header.zlo];
        let pbc: Pbc3 = [true, true, true];

        let simbox = if let (Some(xy), Some(xz), Some(yz)) = (header.xy, header.xz, header.yz) {
            let h = array![[lx, xy, xz], [0.0, ly, yz], [0.0, 0.0, lz]];
            SimBox::new(h, origin, pbc).map_err(|e| err_mapper(format!("{:?}", e)))?
        } else {
            SimBox::ortho(array![lx, ly, lz], origin, pbc)
                .map_err(|e| err_mapper(format!("{:?}", e)))?
        };
        frame.simbox = Some(simbox);
    }

    // Add type labels to metadata
    if !atom_type_labels.is_empty() {
        let labels_str = atom_type_labels
            .iter()
            .map(|(id, label)| format!("{}:{}", id, label))
            .collect::<Vec<_>>()
            .join(",");
        frame
            .meta
            .insert("atom_type_labels".to_string(), labels_str);
    }

    if !bond_type_labels.is_empty() {
        let labels_str = bond_type_labels
            .iter()
            .map(|(id, label)| format!("{}:{}", id, label))
            .collect::<Vec<_>>()
            .join(",");
        frame
            .meta
            .insert("bond_type_labels".to_string(), labels_str);
    }

    if !angle_type_labels.is_empty() {
        let labels_str = angle_type_labels
            .iter()
            .map(|(id, label)| format!("{}:{}", id, label))
            .collect::<Vec<_>>()
            .join(",");
        frame
            .meta
            .insert("angle_type_labels".to_string(), labels_str);
    }

    if !dihedral_type_labels.is_empty() {
        let labels_str = dihedral_type_labels
            .iter()
            .map(|(id, label)| format!("{}:{}", id, label))
            .collect::<Vec<_>>()
            .join(",");
        frame
            .meta
            .insert("dihedral_type_labels".to_string(), labels_str);
    }

    if !improper_type_labels.is_empty() {
        let labels_str = improper_type_labels
            .iter()
            .map(|(id, label)| format!("{}:{}", id, label))
            .collect::<Vec<_>>()
            .join(",");
        frame
            .meta
            .insert("improper_type_labels".to_string(), labels_str);
    }

    Ok(frame)
}

// ============================================================================
// Reader
// ============================================================================

/// LAMMPS data file reader implementing the unified `Reader` interface.
///
/// LAMMPS data files are single-frame files; use `read_frame()` to load the frame.
/// and `read(n)` for n > 0 returns `Ok(None)`.
pub struct LAMMPSDataReader<R: BufRead + Seek> {
    reader: R,
    frame: OnceCell<Option<Frame>>,
    returned: bool,
}

impl<R: BufRead + Seek> LAMMPSDataReader<R> {
    /// Create a new LAMMPS data file reader
    pub fn new(reader: R) -> Self {
        Self {
            reader,
            frame: OnceCell::new(),
            returned: false,
        }
    }

    /// Parse the entire file and return the frame
    fn parse_file(&mut self) -> std::io::Result<Option<Frame>> {
        self.reader.seek(SeekFrom::Start(0))?;

        // Parse header returns when it hits a section header line
        // We need to handle that first line
        let (header, first_section_line) = parse_header_with_first_section(&mut self.reader)?;

        let mut atom_type_labels = HashMap::new();
        let mut bond_type_labels = HashMap::new();
        let mut angle_type_labels = HashMap::new();
        let mut dihedral_type_labels = HashMap::new();
        let mut atoms = Vec::new();
        let mut bonds = Vec::new();
        let mut angles = Vec::new();
        let mut dihedrals = Vec::new();

        // Process the first section line that parse_header found
        let mut next_section = first_section_line;

        let mut impropers: Vec<ImproperData> = Vec::new();
        let mut improper_type_labels: HashMap<String, String> = HashMap::new();

        while let Some(line) = next_section.take() {
            let trimmed = line.trim();
            if trimmed.starts_with("Atom Type Labels") {
                let (labels, next) = parse_type_labels(&mut self.reader)?;
                atom_type_labels = labels;
                next_section = next;
            } else if trimmed.starts_with("Bond Type Labels") {
                let (labels, next) = parse_type_labels(&mut self.reader)?;
                bond_type_labels = labels;
                next_section = next;
            } else if trimmed.starts_with("Angle Type Labels") {
                let (labels, next) = parse_type_labels(&mut self.reader)?;
                angle_type_labels = labels;
                next_section = next;
            } else if trimmed.starts_with("Dihedral Type Labels") {
                let (labels, next) = parse_type_labels(&mut self.reader)?;
                dihedral_type_labels = labels;
                next_section = next;
            } else if trimmed.starts_with("Improper Type Labels") {
                let (labels, next) = parse_type_labels(&mut self.reader)?;
                improper_type_labels = labels;
                next_section = next;
            } else if trimmed.starts_with("Atoms") {
                let hint = parse_atoms_style_hint(trimmed);
                atoms = parse_atoms(
                    &mut self.reader,
                    header.num_atoms,
                    &atom_type_labels,
                    hint.as_deref(),
                )?;
                break; // Atoms section doesn't return next section
            } else if trimmed.starts_with("Bonds") {
                bonds = parse_bonds(&mut self.reader, header.num_bonds, &bond_type_labels)?;
                break; // Bonds section doesn't return next section
            } else if trimmed.starts_with("Angles") {
                angles = parse_angles(&mut self.reader, header.num_angles, &angle_type_labels)?;
                break; // Angles section doesn't return next section
            } else if trimmed.starts_with("Dihedrals") {
                dihedrals = parse_dihedrals(
                    &mut self.reader,
                    header.num_dihedrals,
                    &dihedral_type_labels,
                )?;
                break; // Dihedrals section doesn't return next section
            } else if trimmed.starts_with("Impropers") {
                impropers = parse_impropers(
                    &mut self.reader,
                    header.num_impropers,
                    &improper_type_labels,
                )?;
                break; // Impropers section doesn't return next section
            } else {
                break; // Unknown section
            }
        }

        // Continue reading remaining sections
        let mut line = String::new();
        loop {
            line.clear();
            let bytes_read = self.reader.read_line(&mut line)?;
            if bytes_read == 0 {
                break;
            }

            let trimmed = line.trim();

            // Skip blank lines and comments
            if trimmed.is_empty() || trimmed.starts_with('#') {
                continue;
            }

            // Check for section headers
            if trimmed.starts_with("Atom Type Labels") {
                let (labels, _) = parse_type_labels(&mut self.reader)?;
                atom_type_labels = labels;
            } else if trimmed.starts_with("Bond Type Labels") {
                let (labels, _) = parse_type_labels(&mut self.reader)?;
                bond_type_labels = labels;
            } else if trimmed.starts_with("Angle Type Labels") {
                let (labels, _) = parse_type_labels(&mut self.reader)?;
                angle_type_labels = labels;
            } else if trimmed.starts_with("Dihedral Type Labels") {
                let (labels, _) = parse_type_labels(&mut self.reader)?;
                dihedral_type_labels = labels;
            } else if trimmed.starts_with("Improper Type Labels") {
                let (labels, _) = parse_type_labels(&mut self.reader)?;
                improper_type_labels = labels;
            } else if trimmed.starts_with("Atoms") {
                let hint = parse_atoms_style_hint(trimmed);
                atoms = parse_atoms(
                    &mut self.reader,
                    header.num_atoms,
                    &atom_type_labels,
                    hint.as_deref(),
                )?;
            } else if trimmed.starts_with("Bonds") {
                bonds = parse_bonds(&mut self.reader, header.num_bonds, &bond_type_labels)?;
            } else if trimmed.starts_with("Angles") {
                angles = parse_angles(&mut self.reader, header.num_angles, &angle_type_labels)?;
            } else if trimmed.starts_with("Dihedrals") {
                dihedrals = parse_dihedrals(
                    &mut self.reader,
                    header.num_dihedrals,
                    &dihedral_type_labels,
                )?;
            } else if trimmed.starts_with("Impropers") {
                impropers = parse_impropers(
                    &mut self.reader,
                    header.num_impropers,
                    &improper_type_labels,
                )?;
            }
            // Ignore other sections (Masses, Velocities, Bodies, etc.)
        }

        if atoms.is_empty() && header.num_atoms > 0 {
            return Err(err_mapper("No atoms found in file"));
        }

        let frame = build_frame(
            &header,
            atoms,
            bonds,
            angles,
            dihedrals,
            impropers,
            atom_type_labels,
            bond_type_labels,
            angle_type_labels,
            dihedral_type_labels,
            improper_type_labels,
        )?;
        Ok(Some(frame))
    }
}

impl<R: BufRead + Seek> Reader for LAMMPSDataReader<R> {
    type R = R;
    type Frame = Frame;

    fn new(reader: R) -> Self {
        Self::new(reader)
    }
}

impl<R: BufRead + Seek> FrameReader for LAMMPSDataReader<R> {
    fn read_frame(&mut self) -> std::io::Result<Option<Self::Frame>> {
        if self.returned {
            return Ok(None);
        }
        // Parse on first access and cache
        if self.frame.get().is_none() {
            let frame = self.parse_file()?;
            let _ = self.frame.set(frame);
        }

        self.returned = true;
        Ok(self.frame.get().unwrap().clone())
    }
}

// ============================================================================
// Writer
// ============================================================================

/// LAMMPS data file writer implementing the `FrameWriter` interface.
pub struct LAMMPSDataWriter<W: Write> {
    writer: W,
}

impl<W: Write> crate::io::writer::Writer for LAMMPSDataWriter<W> {
    type W = W;
    type FrameLike = Frame;

    fn new(writer: W) -> Self {
        Self { writer }
    }
}

impl<W: Write> LAMMPSDataWriter<W> {
    /// Create a new LAMMPS data file writer
    pub fn new(writer: W) -> Self {
        Self { writer }
    }
}

impl<W: Write> FrameWriter for LAMMPSDataWriter<W> {
    fn write_frame(&mut self, frame: &Frame) -> std::io::Result<()> {
        write_lammps_data_frame(&mut self.writer, frame)
    }
}

/// Write a single frame in LAMMPS data file format.
///
/// Accepts any type implementing [`FrameAccess`], including both [`Frame`] and
/// [`FrameView`](crate::io::frame_view::FrameView).
fn write_lammps_data_frame<W: Write>(
    writer: &mut W,
    frame: &impl FrameAccess,
) -> std::io::Result<()> {
    // Write comment line
    writeln!(writer, "# LAMMPS data file generated by molrs")?;
    writeln!(writer)?;

    // Get atom count and type info via FrameAccess
    let atom_types = frame
        .get_int("atoms", "type")
        .ok_or_else(|| err_mapper("Atoms block must contain 'type' column"))?;
    let num_atoms = atom_types.shape().first().copied().unwrap_or(0);
    let num_atom_types = atom_types.iter().max().copied().unwrap_or(1) as usize;

    // Get bonds info
    let bonds_type_view = frame.get_int("bonds", "type");
    let num_bonds = frame
        .visit_block("bonds", |b| b.nrows().unwrap_or(0))
        .unwrap_or(0);
    let num_bond_types = bonds_type_view
        .as_ref()
        .map(|bt| bt.iter().max().copied().unwrap_or(1) as usize)
        .unwrap_or(0);

    // Write header
    writeln!(writer, "{} atoms", num_atoms)?;
    if num_bonds > 0 {
        writeln!(writer, "{} bonds", num_bonds)?;
    }
    writeln!(writer, "{} atom types", num_atom_types)?;
    if num_bond_types > 0 {
        writeln!(writer, "{} bond types", num_bond_types)?;
    }
    writeln!(writer)?;

    // Write box dimensions from the canonical SimBox. Fall back to a unit box
    // if the frame has none — LAMMPS data files require box bounds in the
    // header, and a default is nicer than failing the whole write.
    let simbox_ref = frame.simbox_ref();
    let (box_origin, box_lengths, tilts) = if let Some(sb) = simbox_ref {
        let origin_view = sb.origin_view();
        let lengths = sb.lengths();
        let tilts = sb.tilts();
        (
            [origin_view[0], origin_view[1], origin_view[2]],
            [lengths[0], lengths[1], lengths[2]],
            [tilts[0], tilts[1], tilts[2]],
        )
    } else {
        ([0.0; 3], [1.0; 3], [0.0; 3])
    };

    writeln!(
        writer,
        "{} {} xlo xhi",
        box_origin[0],
        box_origin[0] + box_lengths[0]
    )?;
    writeln!(
        writer,
        "{} {} ylo yhi",
        box_origin[1],
        box_origin[1] + box_lengths[1]
    )?;
    writeln!(
        writer,
        "{} {} zlo zhi",
        box_origin[2],
        box_origin[2] + box_lengths[2]
    )?;

    // Write tilt factors if the box is triclinic (any nonzero tilt).
    if tilts.iter().any(|&t| t != 0.0) {
        writeln!(writer, "{} {} {} xy xz yz", tilts[0], tilts[1], tilts[2])?;
    }
    writeln!(writer)?;

    let meta = frame.meta_ref();

    // Write type labels if present
    if let Some(atom_labels_str) = meta.get("atom_type_labels") {
        writeln!(writer, "Atom Type Labels")?;
        writeln!(writer)?;
        for pair in atom_labels_str.split(',') {
            let parts: Vec<&str> = pair.split(':').collect();
            if parts.len() == 2 {
                writeln!(writer, "{} {}", parts[0], parts[1])?;
            }
        }
        writeln!(writer)?;
    }

    if let Some(bond_labels_str) = meta.get("bond_type_labels") {
        writeln!(writer, "Bond Type Labels")?;
        writeln!(writer)?;
        for pair in bond_labels_str.split(',') {
            let parts: Vec<&str> = pair.split(':').collect();
            if parts.len() == 2 {
                writeln!(writer, "{} {}", parts[0], parts[1])?;
            }
        }
        writeln!(writer)?;
    }

    // Write Atoms section
    let x = frame
        .get_float("atoms", "x")
        .ok_or_else(|| err_mapper("Missing 'x' column"))?;
    let y = frame
        .get_float("atoms", "y")
        .ok_or_else(|| err_mapper("Missing 'y' column"))?;
    let z = frame
        .get_float("atoms", "z")
        .ok_or_else(|| err_mapper("Missing 'z' column"))?;
    let ids = frame
        .get_int("atoms", "id")
        .ok_or_else(|| err_mapper("Missing 'id' column"))?;

    let charges = frame.get_float("atoms", "charge");
    let mol_ids = frame.get_int("atoms", "molecule_id");
    let has_charge = charges.is_some();
    let has_mol_id = mol_ids.is_some();

    // Determine atom style
    let atom_style = if has_mol_id && has_charge {
        "full"
    } else if has_mol_id {
        "molecular"
    } else if has_charge {
        "charge"
    } else {
        "atomic"
    };

    writeln!(writer, "Atoms # {}", atom_style)?;
    writeln!(writer)?;

    for i in 0..num_atoms {
        write!(writer, "{}", ids[i])?;

        if let Some(ref mol_ids) = mol_ids {
            write!(writer, " {}", mol_ids[i])?;
        }

        write!(writer, " {}", atom_types[i])?;

        if let Some(ref charges) = charges {
            write!(writer, " {}", charges[i])?;
        }

        writeln!(writer, " {} {} {}", x[i], y[i], z[i])?;
    }
    writeln!(writer)?;

    // Write Bonds section if present
    if num_bonds > 0 {
        writeln!(writer, "Bonds")?;
        writeln!(writer)?;

        let atom_i = frame
            .get_uint("bonds", "atomi")
            .ok_or_else(|| err_mapper("Missing 'atomi' column"))?;
        let atom_j = frame
            .get_uint("bonds", "atomj")
            .ok_or_else(|| err_mapper("Missing 'atomj' column"))?;
        let bond_types = bonds_type_view.ok_or_else(|| err_mapper("Missing 'type' column"))?;

        // Convert indices back to IDs
        for i in 0..num_bonds {
            let id_i = ids[atom_i[i] as usize];
            let id_j = ids[atom_j[i] as usize];
            writeln!(writer, "{} {} {} {}", i + 1, bond_types[i], id_i, id_j)?;
        }
        writeln!(writer)?;
    }

    Ok(())
}

// ============================================================================
// Helper Functions
// ============================================================================

/// Read a LAMMPS data file from a path
pub fn read_lammps_data<P: AsRef<Path>>(path: P) -> std::io::Result<Frame> {
    let file = File::open(path)?;
    let mut reader = LAMMPSDataReader::new(BufReader::new(file));
    reader
        .read_frame()?
        .ok_or_else(|| err_mapper("No frame found in LAMMPS data file"))
}

/// Write a Frame to a LAMMPS data file.
///
/// Accepts any type implementing [`FrameAccess`], including both [`Frame`] and
/// [`FrameView`](crate::io::frame_view::FrameView).
pub fn write_lammps_data<P: AsRef<Path>>(path: P, frame: &impl FrameAccess) -> std::io::Result<()> {
    let file = File::create(path)?;
    let mut writer = std::io::BufWriter::new(file);
    write_lammps_data_frame(&mut writer, frame)
}

// ============================================================================
// Streaming
// ============================================================================

use crate::io::streaming::{FrameIndexBuilder, FrameIndexEntry};
use std::io::Cursor;

/// Parse a LAMMPS data frame from a tightly-bounded byte slice. The slice
/// must contain the entire LAMMPS data file (data files are single-frame).
pub fn parse_frame_bytes(bytes: &[u8]) -> std::io::Result<Frame> {
    let cursor = Cursor::new(bytes);
    let mut reader = LAMMPSDataReader::new(cursor);
    reader
        .read_frame()?
        .ok_or_else(|| err_mapper("No frame found in LAMMPS data slice"))
}

/// Streaming frame indexer for LAMMPS data files.
///
/// LAMMPS data files contain exactly one frame; this indexer simply tracks
/// total bytes seen and emits a single `FrameIndexEntry { 0, bytes_seen }`
/// from `finish`.
///
/// `byte_len` is `u32` so files larger than 4 GiB are rejected — that's a
/// future-spec lift.
pub struct LammpsDataIndexBuilder {
    bytes_seen: u64,
}

impl Default for LammpsDataIndexBuilder {
    fn default() -> Self {
        Self::new()
    }
}

impl LammpsDataIndexBuilder {
    pub fn new() -> Self {
        Self { bytes_seen: 0 }
    }
}

impl FrameIndexBuilder for LammpsDataIndexBuilder {
    fn feed(&mut self, chunk: &[u8], global_offset: u64) {
        self.bytes_seen = global_offset.saturating_add(chunk.len() as u64);
    }

    fn drain(&mut self) -> Vec<FrameIndexEntry> {
        Vec::new()
    }

    fn finish(self: Box<Self>) -> std::io::Result<Vec<FrameIndexEntry>> {
        if self.bytes_seen == 0 {
            return Ok(Vec::new());
        }
        if self.bytes_seen > u32::MAX as u64 {
            return Err(std::io::Error::new(
                std::io::ErrorKind::InvalidData,
                "LAMMPS data file exceeds 4 GiB",
            ));
        }
        Ok(vec![FrameIndexEntry {
            byte_offset: 0,
            byte_len: self.bytes_seen as u32,
        }])
    }

    fn bytes_seen(&self) -> u64 {
        self.bytes_seen
    }
}

#[cfg(test)]
mod streaming_tests {
    use super::*;

    const TINY_DATA: &str = concat!(
        "LAMMPS data file\n",
        "\n",
        "2 atoms\n",
        "1 atom types\n",
        "\n",
        "0.0 10.0 xlo xhi\n",
        "0.0 10.0 ylo yhi\n",
        "0.0 10.0 zlo zhi\n",
        "\n",
        "Atoms\n",
        "\n",
        "1 1 0.0 0.0 0.0\n",
        "2 1 1.0 0.0 0.0\n",
    );

    fn build_chunked(bytes: &[u8], cs: usize) -> Vec<FrameIndexEntry> {
        let mut b = Box::new(LammpsDataIndexBuilder::new());
        let mut off: u64 = 0;
        let mut out: Vec<FrameIndexEntry> = Vec::new();
        for piece in bytes.chunks(cs.max(1)) {
            b.feed(piece, off);
            off += piece.len() as u64;
            out.extend(b.drain());
        }
        out.extend(b.finish().expect("finish"));
        out
    }

    #[test]
    fn lammps_data_streaming_single_frame() {
        let bytes = TINY_DATA.as_bytes();
        let one = build_chunked(bytes, bytes.len());
        for cs in [1usize, 7, 31, 64, bytes.len()] {
            let chunked = build_chunked(bytes, cs);
            assert_eq!(one, chunked, "chunk size {}", cs);
        }
        assert_eq!(one.len(), 1);
        assert_eq!(one[0].byte_offset, 0);
        assert_eq!(one[0].byte_len as usize, bytes.len());
        let frame = parse_frame_bytes(bytes).expect("parse");
        assert_eq!(frame.get("atoms").unwrap().nrows().unwrap(), 2);
    }
}