oxiphysics-io 0.1.2

//! Auto-generated module
//!
//! 🤖 Generated with [SplitRS](https://github.com/cool-japan/splitrs)

use super::super::types::{
    AmberAngle, AmberBond, AmberDihedral, AmberEnergyFrame, AmberTopology, FfAngleParam,
    FfBondParam, LjParameter, MdcrdFrame,
};

/// Extract the raw text body of a `%FLAG `name` section.
///
/// Returns the lines after `%FORMAT(...)` up to the next `%FLAG` or EOF.
pub fn parse_flag_section(content: &str, flag_name: &str) -> Option<String> {
    let target = format!("%FLAG {}", flag_name.to_uppercase());
    let lines: Vec<&str> = content.lines().collect();
    let flag_pos = lines.iter().position(|l| l.trim() == target.trim())?;
    let data_start = flag_pos + 2;
    let mut body = String::new();
    for line in lines.iter().skip(data_start) {
        let trimmed = line.trim();
        if trimmed.starts_with("%FLAG") || trimmed.starts_with("%VERSION") {
            break;
        }
        body.push_str(line);
        body.push('\n');
    }
    Some(body)
}
/// Parse Fortran-formatted real numbers (fixed-width or space-separated).
pub fn parse_fortran_reals(s: &str) -> Vec<f64> {
    s.split_whitespace()
        .filter_map(|tok| tok.replace('D', "E").replace('d', "e").parse::<f64>().ok())
        .collect()
}
/// Parse Fortran-formatted integers (space-separated).
pub fn parse_fortran_ints(s: &str) -> Vec<i64> {
    s.split_whitespace()
        .filter_map(|tok| tok.parse::<i64>().ok())
        .collect()
}
/// Parse Fortran-formatted string tokens (fixed-width 4-char columns or whitespace separated).
pub fn parse_fortran_strings(s: &str) -> Vec<String> {
    let mut result = Vec::new();
    for line in s.lines() {
        if line.trim().is_empty() {
            continue;
        }
        let chars: Vec<char> = line.chars().collect();
        if chars.len() >= 4 {
            let mut i = 0;
            while i + 4 <= chars.len() {
                let tok: String = chars[i..i + 4].iter().collect();
                let trimmed = tok.trim().to_string();
                if !trimmed.is_empty() {
                    result.push(trimmed);
                }
                i += 4;
            }
        } else {
            for tok in line.split_whitespace() {
                result.push(tok.to_string());
            }
        }
    }
    result
}
/// Build `AmberBond` entries from a flat integer array (i*3, j*3, type_idx triplets).
pub fn build_bonds(ints: &[i64], force_k: &[f64], equil_r: &[f64]) -> Vec<AmberBond> {
    let mut bonds = Vec::new();
    let mut idx = 0;
    while idx + 2 < ints.len() {
        let i = (ints[idx] / 3) as usize;
        let j = (ints[idx + 1] / 3) as usize;
        let type_idx = (ints[idx + 2] - 1).max(0) as usize;
        let k = force_k.get(type_idx).copied().unwrap_or(0.0);
        let r0 = equil_r.get(type_idx).copied().unwrap_or(0.0);
        bonds.push(AmberBond { i, j, k, r0 });
        idx += 3;
    }
    bonds
}
/// Build `AmberAngle` entries from a flat integer array (i*3, j*3, k*3, type_idx quads).
pub fn build_angles(ints: &[i64], force_k: &[f64], equil_t: &[f64]) -> Vec<AmberAngle> {
    let mut angles = Vec::new();
    let mut idx = 0;
    while idx + 3 < ints.len() {
        let i = (ints[idx] / 3) as usize;
        let j = (ints[idx + 1] / 3) as usize;
        let k_idx = (ints[idx + 3] - 1).max(0) as usize;
        let k = force_k.get(k_idx).copied().unwrap_or(0.0);
        let theta0 = equil_t.get(k_idx).copied().unwrap_or(0.0);
        angles.push(AmberAngle {
            i,
            j,
            k_idx,
            k,
            theta0,
        });
        idx += 4;
    }
    angles
}
/// Parse Lennard-Jones A/B coefficient pairs from a raw prmtop string.
///
/// Looks for the `%FLAG LENNARD_JONES_ACOEF` and `%FLAG LENNARD_JONES_BCOEF`
/// sections and returns one [`LjParameter`] per pair.
///
/// This is a best-effort parser; it handles well-formed AMBER 7+ files only.
pub fn extract_lj_parameters(prmtop: &str) -> Vec<LjParameter> {
    let a_vals = extract_flag_floats(prmtop, "LENNARD_JONES_ACOEF");
    let b_vals = extract_flag_floats(prmtop, "LENNARD_JONES_BCOEF");
    let n = a_vals.len().min(b_vals.len());
    (0..n)
        .map(|i| LjParameter {
            type_a: i,
            type_b: i,
            a_coeff: a_vals[i],
            b_coeff: b_vals[i],
        })
        .collect()
}
/// Parse dihedral entries from `%FLAG DIHEDRALS_INC_HYDROGEN` and
/// `%FLAG DIHEDRALS_WITHOUT_HYDROGEN` sections.
///
/// Returns a flat list of [`AmberDihedral`]s in the order they appear.
pub fn extract_dihedrals(prmtop: &str) -> Vec<AmberDihedral> {
    let mut out = Vec::new();
    for flag in &["DIHEDRALS_INC_HYDROGEN", "DIHEDRALS_WITHOUT_HYDROGEN"] {
        let raw = extract_flag_integers(prmtop, flag);
        let chunks = raw.chunks_exact(5);
        for chunk in chunks {
            let is_improper = chunk[2] < 0;
            out.push(AmberDihedral {
                i: (chunk[0].abs() / 3) as usize,
                j: (chunk[1].abs() / 3) as usize,
                k: (chunk[2].abs() / 3) as usize,
                l: (chunk[3].abs() / 3) as usize,
                v_n: 0.0,
                gamma: 0.0,
                n: chunk[4].abs(),
                is_improper,
            });
        }
    }
    out
}
/// Parse a sequence of energy frames from an AMBER mdout-style string.
///
/// Only lines containing `NSTEP =` and `Etot` are considered; the rest are
/// skipped.  Returns one frame per `NSTEP` block.
pub fn parse_energy_frames(mdout: &str) -> Vec<AmberEnergyFrame> {
    let mut frames = Vec::new();
    let mut current_step: Option<u64> = None;
    let mut current_time: f64 = 0.0;
    let mut current_temp: f64 = 0.0;
    let mut e_tot = 0.0f64;
    let mut e_kin = 0.0f64;
    let mut e_pot = 0.0f64;
    let mut have_energy = false;
    for line in mdout.lines() {
        let t = line.trim();
        if t.contains("NSTEP =") {
            if let Some(step) = current_step.take() {
                if have_energy {
                    frames.push(AmberEnergyFrame {
                        step,
                        time_ps: current_time,
                        temp_k: current_temp,
                        e_tot,
                        e_kin,
                        e_pot,
                    });
                }
                have_energy = false;
            }
            current_step = parse_field(t, "NSTEP =").and_then(|s| s.parse().ok());
            current_time = parse_field(t, "TIME(PS) =")
                .and_then(|s| s.parse().ok())
                .unwrap_or(0.0);
            current_temp = parse_field(t, "TEMP(K) =")
                .and_then(|s| s.parse().ok())
                .unwrap_or(0.0);
        } else if t.starts_with("Etot") {
            e_tot = parse_field(t, "Etot   =")
                .and_then(|s| s.parse().ok())
                .unwrap_or(0.0);
            e_kin = parse_field(t, "EKtot   =")
                .and_then(|s| s.parse().ok())
                .unwrap_or(0.0);
            e_pot = parse_field(t, "EPtot   =")
                .and_then(|s| s.parse().ok())
                .unwrap_or(0.0);
            have_energy = true;
        }
    }
    if let Some(step) = current_step
        && have_energy
    {
        frames.push(AmberEnergyFrame {
            step,
            time_ps: current_time,
            temp_k: current_temp,
            e_tot,
            e_kin,
            e_pot,
        });
    }
    frames
}
/// Write a minimal AMBER restart (.inpcrd) string from positions and optional
/// velocities.
///
/// Each coordinate triplet is formatted to 12.7f in groups of 2 per line
/// (6 values per line), which matches the standard AMBER inpcrd layout.
pub fn write_inpcrd(
    title: &str,
    positions: &[[f64; 3]],
    velocities: Option<&[[f64; 3]]>,
) -> String {
    let mut out = String::new();
    out.push_str(title);
    out.push('\n');
    out.push_str(&format!("{:6}\n", positions.len()));
    let flat_pos: Vec<f64> = positions.iter().flat_map(|p| p.iter().copied()).collect();
    for chunk in flat_pos.chunks(6) {
        let line: String = chunk.iter().map(|v| format!("{:12.7}", v)).collect();
        out.push_str(&line);
        out.push('\n');
    }
    if let Some(vels) = velocities {
        let flat_vel: Vec<f64> = vels.iter().flat_map(|v| v.iter().copied()).collect();
        for chunk in flat_vel.chunks(6) {
            let line: String = chunk.iter().map(|v| format!("{:12.7}", v)).collect();
            out.push_str(&line);
            out.push('\n');
        }
    }
    out
}
/// Minimal built-in AMBER99SB bond parameter table.
///
/// Only the most common backbone bonds are included for testing/demo purposes.
pub fn amber99sb_bonds() -> Vec<FfBondParam> {
    vec![
        FfBondParam {
            type_a: "N".into(),
            type_b: "H".into(),
            k: 434.0,
            r0: 1.010,
        },
        FfBondParam {
            type_a: "N".into(),
            type_b: "CT".into(),
            k: 337.0,
            r0: 1.449,
        },
        FfBondParam {
            type_a: "CT".into(),
            type_b: "CT".into(),
            k: 310.0,
            r0: 1.526,
        },
        FfBondParam {
            type_a: "CT".into(),
            type_b: "HC".into(),
            k: 340.0,
            r0: 1.090,
        },
        FfBondParam {
            type_a: "CT".into(),
            type_b: "C".into(),
            k: 317.0,
            r0: 1.522,
        },
        FfBondParam {
            type_a: "C".into(),
            type_b: "O".into(),
            k: 570.0,
            r0: 1.229,
        },
        FfBondParam {
            type_a: "C".into(),
            type_b: "N".into(),
            k: 490.0,
            r0: 1.335,
        },
        FfBondParam {
            type_a: "OH".into(),
            type_b: "HO".into(),
            k: 553.0,
            r0: 0.960,
        },
    ]
}
/// Minimal built-in AMBER99SB angle parameter table.
pub fn amber99sb_angles() -> Vec<FfAngleParam> {
    vec![
        FfAngleParam {
            type_i: "H".into(),
            type_j: "N".into(),
            type_k: "H".into(),
            k: 35.0,
            theta0_deg: 120.0,
        },
        FfAngleParam {
            type_i: "H".into(),
            type_j: "N".into(),
            type_k: "CT".into(),
            k: 50.0,
            theta0_deg: 118.0,
        },
        FfAngleParam {
            type_i: "N".into(),
            type_j: "CT".into(),
            type_k: "C".into(),
            k: 63.0,
            theta0_deg: 111.2,
        },
        FfAngleParam {
            type_i: "N".into(),
            type_j: "CT".into(),
            type_k: "CT".into(),
            k: 80.0,
            theta0_deg: 111.2,
        },
        FfAngleParam {
            type_i: "CT".into(),
            type_j: "C".into(),
            type_k: "O".into(),
            k: 80.0,
            theta0_deg: 120.4,
        },
        FfAngleParam {
            type_i: "CT".into(),
            type_j: "C".into(),
            type_k: "N".into(),
            k: 70.0,
            theta0_deg: 116.6,
        },
        FfAngleParam {
            type_i: "CT".into(),
            type_j: "CT".into(),
            type_k: "HC".into(),
            k: 50.0,
            theta0_deg: 109.5,
        },
        FfAngleParam {
            type_i: "HC".into(),
            type_j: "CT".into(),
            type_k: "HC".into(),
            k: 35.0,
            theta0_deg: 109.5,
        },
    ]
}
/// Compute summary statistics over a slice of energy frames.
///
/// Returns `(mean_etot, min_etot, max_etot)`.
pub fn energy_summary(frames: &[AmberEnergyFrame]) -> (f64, f64, f64) {
    if frames.is_empty() {
        return (0.0, 0.0, 0.0);
    }
    let vals: Vec<f64> = frames.iter().map(|f| f.e_tot).collect();
    let mean = vals.iter().sum::<f64>() / vals.len() as f64;
    let min = vals.iter().cloned().fold(f64::INFINITY, f64::min);
    let max = vals.iter().cloned().fold(f64::NEG_INFINITY, f64::max);
    (mean, min, max)
}
pub(super) fn parse_field<'a>(line: &'a str, key: &str) -> Option<&'a str> {
    let pos = line.find(key)?;
    let rest = line[pos + key.len()..].trim_start();
    let end = rest.find([' ', '\t', '\n']).unwrap_or(rest.len());
    Some(&rest[..end])
}
pub(super) fn extract_flag_floats(prmtop: &str, flag: &str) -> Vec<f64> {
    let marker = format!("%FLAG {flag}");
    let mut in_section = false;
    let mut out = Vec::new();
    for line in prmtop.lines() {
        let t = line.trim();
        if t == marker {
            in_section = true;
            continue;
        }
        if in_section {
            if t.starts_with("%FORMAT") {
                continue;
            }
            if t.starts_with('%') {
                break;
            }
            for tok in t.split_whitespace() {
                if let Ok(v) = tok.parse::<f64>() {
                    out.push(v);
                }
            }
        }
    }
    out
}
pub(super) fn extract_flag_integers(prmtop: &str, flag: &str) -> Vec<i32> {
    let marker = format!("%FLAG {flag}");
    let mut in_section = false;
    let mut out = Vec::new();
    for line in prmtop.lines() {
        let t = line.trim();
        if t == marker {
            in_section = true;
            continue;
        }
        if in_section {
            if t.starts_with("%FORMAT") {
                continue;
            }
            if t.starts_with('%') {
                break;
            }
            for tok in t.split_whitespace() {
                if let Ok(v) = tok.parse::<i32>() {
                    out.push(v);
                }
            }
        }
    }
    out
}
/// Strip AMBER-style inline comment (everything from `;` onward).
pub fn strip_amber_comment(s: &str) -> &str {
    if let Some(pos) = s.find(';') {
        s[..pos].trim()
    } else {
        s.trim()
    }
}
/// Collect all floating-point numbers from a string, ignoring non-numeric tokens.
pub fn collect_numbers(s: &str) -> Vec<f64> {
    s.split_whitespace()
        .filter_map(|tok| tok.parse::<f64>().ok())
        .collect()
}
/// Parse a residue selection like "1", "1-10", "1,3,5", "1-5,8,10-12".
pub fn parse_residue_selection(s: &str) -> Result<Vec<usize>, String> {
    let mut result = Vec::new();
    for token in s.split(',') {
        let t = token.trim();
        if t.is_empty() {
            continue;
        }
        if let Some(dash_pos) = t.find('-') {
            let start: usize = t[..dash_pos]
                .trim()
                .parse()
                .map_err(|_| format!("Invalid range start: '{}'", &t[..dash_pos]))?;
            let end: usize = t[dash_pos + 1..]
                .trim()
                .parse()
                .map_err(|_| format!("Invalid range end: '{}'", &t[dash_pos + 1..]))?;
            if end < start {
                return Err(format!("Range end {} < start {}", end, start));
            }
            result.extend(start..=end);
        } else {
            let n: usize = t
                .parse()
                .map_err(|_| format!("Invalid residue number: '{}'", t))?;
            result.push(n);
        }
    }
    Ok(result)
}
/// Parse AMBER MDCRD trajectory data from a string.
///
/// MDCRD format: title line, then coordinates in groups of 10 per line
/// (8.3f format), with optional box at the end of each frame.
pub fn parse_mdcrd(s: &str, n_atoms: usize, has_box: bool) -> Vec<MdcrdFrame> {
    let mut frames = Vec::new();
    let mut lines = s.lines();
    let _ = lines.next();
    let n_coords = n_atoms * 3;
    let coords_per_line = 10_usize;
    let coords_lines = n_coords.div_ceil(coords_per_line);
    let mut remaining_lines: Vec<&str> = lines.collect();
    let box_lines = if has_box { 1 } else { 0 };
    let frame_lines = coords_lines + box_lines;
    let mut i = 0;
    while i + coords_lines <= remaining_lines.len() {
        let mut coords: Vec<f64> = Vec::with_capacity(n_coords);
        for line in &remaining_lines[i..i + coords_lines] {
            for tok in line.split_whitespace() {
                if let Ok(v) = tok.parse::<f64>() {
                    coords.push(v);
                    if coords.len() == n_coords {
                        break;
                    }
                }
            }
        }
        let box_dims = if has_box && i + frame_lines <= remaining_lines.len() {
            let box_line = remaining_lines[i + coords_lines];
            let vals: Vec<f64> = box_line
                .split_whitespace()
                .filter_map(|t| t.parse().ok())
                .collect();
            if vals.len() >= 3 {
                Some([vals[0], vals[1], vals[2]])
            } else {
                None
            }
        } else {
            None
        };
        if coords.len() >= n_coords {
            frames.push(MdcrdFrame { coords, box_dims });
        }
        i += frame_lines;
        if frame_lines == 0 {
            break;
        }
        remaining_lines = remaining_lines[i..].to_vec();
        i = 0;
    }
    frames
}
/// Serialises a minimal AMBER prmtop-format string from an `AmberTopology`.
///
/// The output contains `%VERSION`, `%FLAG TITLE`, `%FLAG POINTERS`,
/// `%FLAG ATOM_NAME`, `%FLAG CHARGE`, and `%FLAG MASS` sections so that it
/// can be round-tripped through `AmberTopology::from_prmtop_str`.
pub fn write_prmtop(topo: &AmberTopology) -> String {
    let mut out = String::new();
    out.push_str("%VERSION  VERSION_STAMP = V0001.000  DATE = 01/01/26  00:00:00\n");
    out.push_str("%FLAG TITLE\n");
    out.push_str("%FORMAT(20a4)\n");
    out.push_str(&format!("{}\n", topo.title));
    out.push_str("%FLAG POINTERS\n");
    out.push_str("%FORMAT(10I8)\n");
    let n_atoms = topo.atoms.len() as i64;
    let n_types = 1i64;
    let n_bonds = topo.bonds.len() as i64;
    let n_bonds_nh = n_bonds;
    let n_angles = topo.angles.len() as i64;
    let pointers: Vec<i64> = vec![
        n_atoms, n_types, n_bonds, n_bonds_nh, n_angles, n_angles, 0, 0, 0, 0,
    ];
    for (i, p) in pointers.iter().enumerate() {
        out.push_str(&format!("{:8}", p));
        if (i + 1) % 10 == 0 {
            out.push('\n');
        }
    }
    if !pointers.len().is_multiple_of(10) {
        out.push('\n');
    }
    out.push_str("%FLAG ATOM_NAME\n");
    out.push_str("%FORMAT(20a4)\n");
    for (i, atom) in topo.atoms.iter().enumerate() {
        out.push_str(&format!("{:<4}", &atom.name[..atom.name.len().min(4)]));
        if (i + 1) % 20 == 0 {
            out.push('\n');
        }
    }
    if !topo.atoms.is_empty() && !topo.atoms.len().is_multiple_of(20) {
        out.push('\n');
    }
    pub(super) const AMBER_CHARGE_SCALE: f64 = 18.2223;
    out.push_str("%FLAG CHARGE\n");
    out.push_str("%FORMAT(5E16.8)\n");
    for (i, atom) in topo.atoms.iter().enumerate() {
        out.push_str(&format!("{:16.8E}", atom.charge * AMBER_CHARGE_SCALE));
        if (i + 1) % 5 == 0 {
            out.push('\n');
        }
    }
    if !topo.atoms.is_empty() && !topo.atoms.len().is_multiple_of(5) {
        out.push('\n');
    }
    out.push_str("%FLAG MASS\n");
    out.push_str("%FORMAT(5E16.8)\n");
    for (i, atom) in topo.atoms.iter().enumerate() {
        out.push_str(&format!("{:16.8E}", atom.mass));
        if (i + 1) % 5 == 0 {
            out.push('\n');
        }
    }
    if !topo.atoms.is_empty() && !topo.atoms.len().is_multiple_of(5) {
        out.push('\n');
    }
    out
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::amber::types::*;
    use std::str::FromStr;
    fn minimal_prmtop() -> String {
        let n_charge = -0.4157 * 18.2223;
        let h_charge = 0.2719 * 18.2223;
        format!(
            "%VERSION  VERSION_STAMP = V0001.000\n\
             %FLAG TITLE\n\
             %FORMAT(20a4)\n\
             test molecule\n\
             %FLAG POINTERS\n\
             %FORMAT(10I8)\n\
               2    0    1    0    1    0    0    0    0    0\n\
               0    0    0    0    0    0    0    0    0    0\n\
               0    0    0    0    0    0    0    0    0    0\n\
               0\n\
             %FLAG ATOM_NAME\n\
             %FORMAT(20a4)\n\
             N   H   \n\
             %FLAG CHARGE\n\
             %FORMAT(5E16.8)\n\
               {:.8E}  {:.8E}\n\
             %FLAG MASS\n\
             %FORMAT(5E16.8)\n\
                1.40100000E+01   1.00800000E+00\n\
             %FLAG AMBER_ATOM_TYPE\n\
             %FORMAT(20a4)\n\
             N   H   \n\
             %FLAG RESIDUE_LABEL\n\
             %FORMAT(20a4)\n\
             ALA \n\
             %FLAG RESIDUE_POINTER\n\
             %FORMAT(10I8)\n\
               1\n\
             %FLAG BONDS_INC_HYDROGEN\n\
             %FORMAT(10I8)\n\
               0    3    1\n\
             %FLAG BONDS_WITHOUT_HYDROGEN\n\
             %FORMAT(10I8)\n\
             %FLAG BOND_FORCE_CONSTANT\n\
             %FORMAT(5E16.8)\n\
                4.34000000E+02\n\
             %FLAG BOND_EQUIL_VALUE\n\
             %FORMAT(5E16.8)\n\
                1.01000000E+00\n\
             %FLAG ANGLES_INC_HYDROGEN\n\
             %FORMAT(10I8)\n\
             %FLAG ANGLES_WITHOUT_HYDROGEN\n\
             %FORMAT(10I8)\n\
             %FLAG ANGLE_FORCE_CONSTANT\n\
             %FORMAT(5E16.8)\n\
             %FLAG ANGLE_EQUIL_VALUE\n\
             %FORMAT(5E16.8)\n",
            n_charge, h_charge
        )
    }
    fn parse_minimal() -> AmberTopology {
        AmberTopology::from_prmtop_str(&minimal_prmtop()).unwrap()
    }
    #[test]
    fn test_parse_flag_section_found() {
        let s = "%FLAG TITLE\n%FORMAT(20a4)\nhello\n%FLAG OTHER\n%FORMAT\ndata\n";
        let result = parse_flag_section(s, "TITLE");
        assert!(result.is_some());
        assert!(result.unwrap().contains("hello"));
    }
    #[test]
    fn test_parse_flag_section_not_found() {
        let s = "%FLAG TITLE\n%FORMAT(20a4)\nhello\n";
        let result = parse_flag_section(s, "MISSING");
        assert!(result.is_none());
    }
    #[test]
    fn test_parse_fortran_reals_basic() {
        let s = "   1.23000000E+00  -4.56000000E+00\n";
        let vals = parse_fortran_reals(s);
        assert_eq!(vals.len(), 2);
        assert!((vals[0] - 1.23).abs() < 1e-6);
        assert!((vals[1] + 4.56).abs() < 1e-6);
    }
    #[test]
    fn test_parse_fortran_reals_d_exponent() {
        let s = "1.0D+00  2.5D-01\n";
        let vals = parse_fortran_reals(s);
        assert_eq!(vals.len(), 2);
        assert!((vals[1] - 0.25).abs() < 1e-9);
    }
    #[test]
    fn test_parse_fortran_ints_basic() {
        let s = "   2    0    1    0\n";
        let vals = parse_fortran_ints(s);
        assert_eq!(vals, vec![2, 0, 1, 0]);
    }
    #[test]
    fn test_parse_title() {
        let top = parse_minimal();
        assert!(top.title.contains("test molecule"), "title='{}'", top.title);
    }
    #[test]
    fn test_n_atoms() {
        let top = parse_minimal();
        assert_eq!(top.n_atoms, 2);
    }
    #[test]
    fn test_atom_names() {
        let top = parse_minimal();
        assert_eq!(top.atoms[0].name, "N");
        assert_eq!(top.atoms[1].name, "H");
    }
    #[test]
    fn test_atom_masses() {
        let top = parse_minimal();
        assert!(
            (top.atoms[0].mass - 14.01).abs() < 0.01,
            "N mass={}",
            top.atoms[0].mass
        );
        assert!(
            (top.atoms[1].mass - 1.008).abs() < 0.01,
            "H mass={}",
            top.atoms[1].mass
        );
    }
    #[test]
    fn test_atom_charges_scaled() {
        let top = parse_minimal();
        assert!(
            (top.atoms[0].charge - (-0.4157)).abs() < 0.001,
            "N charge={:.4}",
            top.atoms[0].charge
        );
    }
    #[test]
    fn test_atom_types() {
        let top = parse_minimal();
        assert_eq!(top.atoms[0].atom_type, "N");
        assert_eq!(top.atoms[1].atom_type, "H");
    }
    #[test]
    fn test_residue_names() {
        let top = parse_minimal();
        let res = top.residue_names();
        assert!(res.contains(&"ALA".to_string()), "residues={res:?}");
    }
    #[test]
    fn test_bonds_parsed() {
        let top = parse_minimal();
        assert_eq!(top.bonds.len(), 1, "bonds={}", top.bonds.len());
        assert_eq!(top.bonds[0].i, 0);
        assert_eq!(top.bonds[0].j, 1);
    }
    #[test]
    fn test_bond_force_constant() {
        let top = parse_minimal();
        assert!((top.bonds[0].k - 434.0).abs() < 1.0, "k={}", top.bonds[0].k);
    }
    #[test]
    fn test_bond_equil_value() {
        let top = parse_minimal();
        assert!(
            (top.bonds[0].r0 - 1.01).abs() < 0.01,
            "r0={}",
            top.bonds[0].r0
        );
    }
    #[test]
    fn test_total_charge() {
        let top = parse_minimal();
        let q = top.total_charge();
        assert!((q - (-0.1438)).abs() < 0.01, "total_charge={q:.4}");
    }
    #[test]
    fn test_total_mass() {
        let top = parse_minimal();
        let m = top.total_mass();
        assert!((m - 15.018).abs() < 0.1, "total_mass={m:.3}");
    }
    #[test]
    fn test_write_summary_contains_atoms() {
        let top = parse_minimal();
        let summary = top.write_summary();
        assert!(summary.contains("Atoms"), "summary={summary}");
        assert!(summary.contains('2'), "summary={summary}");
    }
    #[test]
    fn test_from_prmtop_empty_string() {
        let top = AmberTopology::from_prmtop_str("").unwrap();
        assert_eq!(top.n_atoms, 0);
        assert!(top.atoms.is_empty());
    }
    fn minimal_inpcrd() -> String {
        "test coords\n    2\n   1.0000000   2.0000000   3.0000000   4.0000000   5.0000000   6.0000000\n"
            .to_string()
    }
    #[test]
    fn test_amber_coords_parse_basic() {
        let coords = AmberCoordinates::from_str(&minimal_inpcrd()).unwrap();
        assert_eq!(coords.n_atoms, 2);
        assert_eq!(coords.coords.len(), 6);
    }
    #[test]
    fn test_amber_coords_position() {
        let coords = AmberCoordinates::from_str(&minimal_inpcrd()).unwrap();
        let p0 = coords.position(0);
        assert!((p0[0] - 1.0).abs() < 1e-6);
        assert!((p0[1] - 2.0).abs() < 1e-6);
        assert!((p0[2] - 3.0).abs() < 1e-6);
    }
    #[test]
    fn test_amber_coords_position_atom1() {
        let coords = AmberCoordinates::from_str(&minimal_inpcrd()).unwrap();
        let p1 = coords.position(1);
        assert!((p1[0] - 4.0).abs() < 1e-6);
    }
    #[test]
    fn test_amber_coords_no_velocities() {
        let coords = AmberCoordinates::from_str(&minimal_inpcrd()).unwrap();
        assert!(coords.velocities.is_none());
    }
    #[test]
    fn test_amber_coords_restart_roundtrip() {
        let coords = AmberCoordinates::from_str(&minimal_inpcrd()).unwrap();
        let restart = coords.to_restart_string();
        assert!(restart.contains("test coords"));
        assert!(restart.contains("2"));
    }
    #[test]
    fn test_amber_coords_empty() {
        let result = AmberCoordinates::from_str("");
        assert!(result.is_err());
    }
    #[test]
    fn test_amber_coords_with_velocities() {
        let s = "test\n    1\n   1.0000000   2.0000000   3.0000000   0.1000000   0.2000000   0.3000000\n";
        let coords = AmberCoordinates::from_str(s).unwrap();
        assert_eq!(coords.n_atoms, 1);
        assert!(coords.velocities.is_some());
        let vels = coords.velocities.unwrap();
        assert!((vels[0] - 0.1).abs() < 1e-6);
    }
    #[test]
    fn test_extract_lj_empty() {
        let params = extract_lj_parameters("");
        assert!(params.is_empty());
    }
    #[test]
    fn test_extract_dihedrals_empty() {
        let dihedrals = extract_dihedrals("");
        assert!(dihedrals.is_empty());
    }
    #[test]
    fn test_atom_type_names() {
        let top = parse_minimal();
        let types = top.atom_type_names();
        assert!(types.contains(&"N".to_string()));
        assert!(types.contains(&"H".to_string()));
    }
    #[test]
    fn test_n_atom_types() {
        let top = parse_minimal();
        assert_eq!(top.n_atom_types(), 2);
    }
    #[test]
    fn test_atoms_in_residue() {
        let top = parse_minimal();
        let ala_atoms = top.atoms_in_residue("ALA");
        assert_eq!(ala_atoms.len(), 2);
    }
    #[test]
    fn test_atoms_in_residue_missing() {
        let top = parse_minimal();
        let gly_atoms = top.atoms_in_residue("GLY");
        assert!(gly_atoms.is_empty());
    }
    #[test]
    fn test_extract_lj_basic() {
        let prmtop = "%FLAG LENNARD_JONES_ACOEF\n%FORMAT(5E16.8)\n1.0  2.0  3.0\n%FLAG LENNARD_JONES_BCOEF\n%FORMAT(5E16.8)\n0.5  1.0  1.5\n";
        let params = extract_lj_parameters(prmtop);
        assert_eq!(params.len(), 3);
        assert!((params[0].a_coeff - 1.0).abs() < 1e-10);
        assert!((params[0].b_coeff - 0.5).abs() < 1e-10);
    }
    #[test]
    fn test_lj_r_min() {
        let lj = LjParameter {
            type_a: 0,
            type_b: 0,
            a_coeff: 2.0,
            b_coeff: 2.0,
        };
        let r_min = lj.r_min();
        assert!((r_min - 2.0f64.powf(1.0 / 6.0)).abs() < 1e-10);
    }
    #[test]
    fn test_lj_epsilon() {
        let lj = LjParameter {
            type_a: 0,
            type_b: 0,
            a_coeff: 4.0,
            b_coeff: 4.0,
        };
        let eps = lj.epsilon();
        assert!((eps - 1.0).abs() < 1e-10);
    }
    #[test]
    fn test_lj_zero_b_coeff_r_min() {
        let lj = LjParameter {
            type_a: 0,
            type_b: 0,
            a_coeff: 1.0,
            b_coeff: 0.0,
        };
        assert_eq!(lj.r_min(), 0.0);
    }
    #[test]
    fn test_extract_dihedrals_basic() {
        let prmtop = "%FLAG DIHEDRALS_INC_HYDROGEN\n%FORMAT(10I8)\n0  3  6  9  1\n";
        let dihedrals = extract_dihedrals(prmtop);
        assert_eq!(dihedrals.len(), 1);
        assert_eq!(dihedrals[0].n, 1);
        assert!(!dihedrals[0].is_improper);
    }
    #[test]
    fn test_extract_dihedrals_improper_flag() {
        let prmtop = "%FLAG DIHEDRALS_WITHOUT_HYDROGEN\n%FORMAT(10I8)\n0  3  -6  9  2\n";
        let dihedrals = extract_dihedrals(prmtop);
        assert_eq!(dihedrals.len(), 1);
        assert!(dihedrals[0].is_improper);
    }
    fn sample_mdout() -> &'static str {
        " NSTEP =       10  TIME(PS) =  0.010000  TEMP(K) =  298.5  PRESS =    0.0\n\
         Etot   =    -100.00  EKtot   =     50.00  EPtot   =    -150.00\n"
    }
    #[test]
    fn test_parse_energy_frames_single() {
        let frames = parse_energy_frames(sample_mdout());
        assert_eq!(frames.len(), 1);
        assert_eq!(frames[0].step, 10);
        assert!((frames[0].time_ps - 0.01).abs() < 1e-6);
        assert!((frames[0].temp_k - 298.5).abs() < 1e-3);
        assert!((frames[0].e_tot - (-100.0)).abs() < 1e-6);
        assert!((frames[0].e_kin - 50.0).abs() < 1e-6);
        assert!((frames[0].e_pot - (-150.0)).abs() < 1e-6);
    }
    #[test]
    fn test_parse_energy_frames_empty() {
        let frames = parse_energy_frames("");
        assert!(frames.is_empty());
    }
    #[test]
    fn test_parse_energy_frames_two_frames() {
        let mdout = " NSTEP =        1  TIME(PS) =  0.001000  TEMP(K) =  300.0  PRESS =   0.0\n\
                     Etot   =   -200.00  EKtot   =   100.00  EPtot   =   -300.00\n\
                     \n\
                     NSTEP =        2  TIME(PS) =  0.002000  TEMP(K) =  301.0  PRESS =   0.0\n\
                     Etot   =   -190.00  EKtot   =   110.00  EPtot   =   -300.00\n";
        let frames = parse_energy_frames(mdout);
        assert_eq!(frames.len(), 2);
        assert_eq!(frames[1].step, 2);
    }
    #[test]
    fn test_energy_summary_basic() {
        let frames = vec![
            AmberEnergyFrame {
                step: 1,
                time_ps: 0.0,
                temp_k: 300.0,
                e_tot: -100.0,
                e_kin: 50.0,
                e_pot: -150.0,
            },
            AmberEnergyFrame {
                step: 2,
                time_ps: 0.1,
                temp_k: 300.0,
                e_tot: -200.0,
                e_kin: 50.0,
                e_pot: -250.0,
            },
        ];
        let (mean, min, max) = energy_summary(&frames);
        assert!((mean - (-150.0)).abs() < 1e-10);
        assert!((min - (-200.0)).abs() < 1e-10);
        assert!((max - (-100.0)).abs() < 1e-10);
    }
    #[test]
    fn test_energy_summary_empty() {
        let (mean, min, max) = energy_summary(&[]);
        assert_eq!(mean, 0.0);
        assert_eq!(min, 0.0);
        assert_eq!(max, 0.0);
    }
    #[test]
    fn test_write_inpcrd_basic() {
        let positions = vec![[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]];
        let s = write_inpcrd("test restart", &positions, None);
        assert!(s.starts_with("test restart"));
        assert!(s.contains("2"));
        assert!(s.contains("1.0000000"));
    }
    #[test]
    fn test_write_inpcrd_with_velocities() {
        let positions = vec![[0.0, 0.0, 0.0]];
        let velocities = vec![[0.1, 0.2, 0.3]];
        let s = write_inpcrd("vel test", &positions, Some(&velocities));
        assert!(s.contains("0.1000000"));
    }
    #[test]
    fn test_amber99sb_bonds_not_empty() {
        let bonds = amber99sb_bonds();
        assert!(!bonds.is_empty());
        let ct_ct = bonds.iter().find(|b| b.type_a == "CT" && b.type_b == "CT");
        assert!(ct_ct.is_some());
    }
    #[test]
    fn test_amber99sb_angles_not_empty() {
        let angles = amber99sb_angles();
        assert!(!angles.is_empty());
        for a in &angles {
            assert!(
                a.k > 0.0,
                "k should be positive for {}-{}-{}",
                a.type_i,
                a.type_j,
                a.type_k
            );
            assert!(a.theta0_deg > 0.0);
        }
    }
}