oxiphysics_io/

crystallography_io.rs

// Copyright 2026 COOLJAPAN OU (Team KitaSan)
// SPDX-License-Identifier: Apache-2.0

//! Crystallography file format I/O.
//!
//! Provides parsers and writers for common crystallography formats:
//! - CIF (Crystallographic Information File)
//! - VASP POSCAR/CONTCAR (stub)
//! - XRD pattern simulation and analysis

#![allow(dead_code)]

use crate::{Error, Result};
use std::collections::HashMap;

// ---------------------------------------------------------------------------
// Data structures
// ---------------------------------------------------------------------------

/// A single atom site in a crystal structure.
#[derive(Debug, Clone)]
pub struct CrystalAtom {
    /// Chemical element symbol (e.g. "Fe", "O").
    pub element: String,
    /// Fractional coordinates within the unit cell \[a, b, c\].
    pub fractional_pos: [f64; 3],
    /// Site occupancy (0.0 – 1.0).
    pub occupancy: f64,
    /// Isotropic displacement parameter (B-factor) in Ų.
    pub b_factor: f64,
}

/// A periodic crystal structure described by its lattice and atom list.
#[derive(Debug, Clone)]
pub struct CrystalStructure {
    /// Lattice vectors as row matrix: `lattice[i]` is the i-th basis vector in Å.
    pub lattice: [[f64; 3]; 3],
    /// List of atoms in the asymmetric unit (or full unit cell).
    pub atoms: Vec<CrystalAtom>,
    /// International Tables space group number (1–230).
    pub space_group: u16,
}

/// Parsed X-ray diffraction pattern.
#[derive(Debug, Clone)]
pub struct XrdPattern {
    /// 2θ angles in degrees.
    pub two_theta: Vec<f64>,
    /// Intensity values (arbitrary units).
    pub intensity: Vec<f64>,
}

// ---------------------------------------------------------------------------
// CIF parser
// ---------------------------------------------------------------------------

/// Parser for CIF (Crystallographic Information File) format.
#[derive(Debug, Default)]
pub struct CifParser;

impl CifParser {
    /// Create a new `CifParser`.
    pub fn new() -> Self {
        Self
    }

    /// Parse a CIF file `content` and return a [`CrystalStructure`].
    ///
    /// Handles `_cell_length_a/b/c`, `_cell_angle_alpha/beta/gamma`, and
    /// `_atom_site_*` loop blocks.
    pub fn parse(&self, content: &str) -> Result<CrystalStructure> {
        let mut kv: HashMap<String, String> = HashMap::new();
        let mut in_loop = false;
        let mut loop_headers: Vec<String> = Vec::new();
        let mut loop_rows: Vec<Vec<String>> = Vec::new();
        let mut current_row: Vec<String> = Vec::new();

        for raw_line in content.lines() {
            let line = raw_line.trim();
            // Skip comments and empty lines
            if line.is_empty() || line.starts_with('#') {
                continue;
            }

            if line.eq_ignore_ascii_case("loop_") {
                // Flush any previous loop
                if in_loop && !current_row.is_empty() {
                    loop_rows.push(current_row.clone());
                    current_row.clear();
                }
                in_loop = true;
                loop_headers.clear();
                loop_rows.clear();
                continue;
            }

            if in_loop {
                if line.starts_with('_') {
                    // Still collecting header names
                    let key = line.split_whitespace().next().unwrap_or("").to_lowercase();
                    loop_headers.push(key);
                    continue;
                }
                // Data row
                let tokens = Self::tokenize_line(line);
                for tok in tokens {
                    current_row.push(tok);
                    if current_row.len() == loop_headers.len() {
                        loop_rows.push(current_row.clone());
                        current_row.clear();
                    }
                }
                // Check if the next line starts a new block
                continue;
            }

            // Key-value pair
            if line.starts_with('_') {
                let tokens = Self::tokenize_line(line);
                if tokens.len() >= 2 {
                    kv.insert(tokens[0].to_lowercase(), tokens[1].clone());
                } else if tokens.len() == 1 {
                    // Value on next line — not handled fully; store empty
                    kv.insert(tokens[0].to_lowercase(), String::new());
                }
            }
        }

        // Flush last loop row
        if in_loop && !current_row.is_empty() {
            loop_rows.push(current_row);
        }

        // Extract cell parameters
        let a = Self::parse_float_kv(&kv, "_cell_length_a").unwrap_or(1.0);
        let b = Self::parse_float_kv(&kv, "_cell_length_b").unwrap_or(1.0);
        let c = Self::parse_float_kv(&kv, "_cell_length_c").unwrap_or(1.0);
        let alpha = Self::parse_float_kv(&kv, "_cell_angle_alpha")
            .unwrap_or(90.0)
            .to_radians();
        let beta = Self::parse_float_kv(&kv, "_cell_angle_beta")
            .unwrap_or(90.0)
            .to_radians();
        let gamma = Self::parse_float_kv(&kv, "_cell_angle_gamma")
            .unwrap_or(90.0)
            .to_radians();

        let lattice = Self::build_lattice(a, b, c, alpha, beta, gamma);

        // Space group
        let space_group = kv
            .get("_space_group_it_number")
            .or_else(|| kv.get("_symmetry_int_tables_number"))
            .and_then(|s| s.parse::<u16>().ok())
            .unwrap_or(1);

        // Extract atoms from loop
        let atoms = Self::extract_atoms(&loop_headers, &loop_rows);

        Ok(CrystalStructure {
            lattice,
            atoms,
            space_group,
        })
    }

    /// Build a lattice matrix from cell parameters (conventional setting).
    fn build_lattice(a: f64, b: f64, c: f64, alpha: f64, beta: f64, gamma: f64) -> [[f64; 3]; 3] {
        let cos_a = alpha.cos();
        let cos_b = beta.cos();
        let cos_g = gamma.cos();
        let sin_g = gamma.sin();

        // Standard triclinic convention
        let cx = cos_b;
        let cy = (cos_a - cos_b * cos_g) / sin_g;
        let cz = (1.0 - cx * cx - cy * cy).max(0.0).sqrt();

        [
            [a, 0.0, 0.0],
            [b * cos_g, b * sin_g, 0.0],
            [c * cx, c * cy, c * cz],
        ]
    }

    /// Tokenize a CIF line, respecting quoted strings.
    fn tokenize_line(line: &str) -> Vec<String> {
        let mut tokens = Vec::new();
        let mut chars = line.char_indices().peekable();
        while let Some((_, c)) = chars.peek().copied() {
            if c.is_whitespace() {
                chars.next();
                continue;
            }
            if c == '\'' || c == '"' {
                chars.next(); // consume opening quote
                let mut tok = String::new();
                for (_, ch) in chars.by_ref() {
                    if ch == c {
                        break;
                    }
                    tok.push(ch);
                }
                tokens.push(tok);
            } else {
                let mut tok = String::new();
                while let Some(&(_, ch)) = chars.peek() {
                    if ch.is_whitespace() {
                        break;
                    }
                    tok.push(ch);
                    chars.next();
                }
                tokens.push(tok);
            }
        }
        tokens
    }

    /// Parse a float from a key-value map entry, ignoring standard uncertainties like "1.234(5)".
    fn parse_float_kv(kv: &HashMap<String, String>, key: &str) -> Option<f64> {
        let s = kv.get(key)?;
        Self::strip_su(s).parse().ok()
    }

    /// Strip standard uncertainty "(n)" from a numeric string.
    fn strip_su(s: &str) -> &str {
        if let Some(pos) = s.find('(') {
            &s[..pos]
        } else {
            s
        }
    }

    /// Extract atom list from loop data.
    fn extract_atoms(headers: &[String], rows: &[Vec<String>]) -> Vec<CrystalAtom> {
        // Find column indices
        let idx = |tag: &str| headers.iter().position(|h| h.contains(tag));

        let i_type_symbol = idx("type_symbol").or_else(|| idx("label"));
        let i_x = idx("fract_x").or_else(|| idx("x_fract"));
        let i_y = idx("fract_y").or_else(|| idx("y_fract"));
        let i_z = idx("fract_z").or_else(|| idx("z_fract"));
        let i_occ = idx("occupancy");
        let i_b = idx("b_iso").or_else(|| idx("u_iso"));

        let mut atoms = Vec::new();
        for row in rows {
            let get = |idx: Option<usize>| -> &str {
                idx.and_then(|i| row.get(i).map(|s| s.as_str()))
                    .unwrap_or("0")
            };

            let element = i_type_symbol
                .and_then(|i| row.get(i))
                .cloned()
                .unwrap_or_default();

            if element.is_empty() {
                continue;
            }

            let x: f64 = Self::strip_su(get(i_x)).parse().unwrap_or(0.0);
            let y: f64 = Self::strip_su(get(i_y)).parse().unwrap_or(0.0);
            let z: f64 = Self::strip_su(get(i_z)).parse().unwrap_or(0.0);
            let occ: f64 = Self::strip_su(get(i_occ)).parse().unwrap_or(1.0);
            let b: f64 = Self::strip_su(get(i_b)).parse().unwrap_or(0.0);

            atoms.push(CrystalAtom {
                element,
                fractional_pos: [x, y, z],
                occupancy: occ,
                b_factor: b,
            });
        }
        atoms
    }
}

// ---------------------------------------------------------------------------
// CrystalStructure methods
// ---------------------------------------------------------------------------

impl CrystalStructure {
    /// Create a simple cubic structure with given lattice parameter.
    pub fn new_cubic(a: f64) -> Self {
        Self {
            lattice: [[a, 0.0, 0.0], [0.0, a, 0.0], [0.0, 0.0, a]],
            atoms: Vec::new(),
            space_group: 1,
        }
    }

    /// Convert all atom fractional coordinates to Cartesian coordinates (Å).
    pub fn to_cartesian(&self) -> Vec<[f64; 3]> {
        self.atoms
            .iter()
            .map(|atom| {
                let [u, v, w] = atom.fractional_pos;
                let x = u * self.lattice[0][0] + v * self.lattice[1][0] + w * self.lattice[2][0];
                let y = u * self.lattice[0][1] + v * self.lattice[1][1] + w * self.lattice[2][1];
                let z = u * self.lattice[0][2] + v * self.lattice[1][2] + w * self.lattice[2][2];
                [x, y, z]
            })
            .collect()
    }

    /// Volume of the unit cell in ų (scalar triple product of basis vectors).
    pub fn volume(&self) -> f64 {
        let a = self.lattice[0];
        let b = self.lattice[1];
        let c = self.lattice[2];
        // a · (b × c)
        let bxc = [
            b[1] * c[2] - b[2] * c[1],
            b[2] * c[0] - b[0] * c[2],
            b[0] * c[1] - b[1] * c[0],
        ];
        (a[0] * bxc[0] + a[1] * bxc[1] + a[2] * bxc[2]).abs()
    }

    /// Crystal density in g/cm³.
    ///
    /// `molar_mass` is the molar mass of one formula unit in g/mol.
    /// The number of formula units Z is taken as the atom count.
    pub fn density(&self, molar_mass: f64) -> f64 {
        const NA: f64 = 6.022_140_76e23; // Avogadro
        const ANG3_TO_CM3: f64 = 1e-24;
        let z = self.atoms.len() as f64;
        let vol_cm3 = self.volume() * ANG3_TO_CM3;
        (z * molar_mass) / (NA * vol_cm3)
    }

    /// Compute the reciprocal lattice vectors (2π convention).
    ///
    /// Returns `[[b1x,b1y,b1z\],[b2x,b2y,b2z],[b3x,b3y,b3z]]` in Å⁻¹.
    pub fn reciprocal_lattice(&self) -> [[f64; 3]; 3] {
        let a1 = self.lattice[0];
        let a2 = self.lattice[1];
        let a3 = self.lattice[2];
        let vol = self.volume();

        let cross = |u: [f64; 3], v: [f64; 3]| -> [f64; 3] {
            [
                u[1] * v[2] - u[2] * v[1],
                u[2] * v[0] - u[0] * v[2],
                u[0] * v[1] - u[1] * v[0],
            ]
        };
        let scale = |s: f64, u: [f64; 3]| -> [f64; 3] { [s * u[0], s * u[1], s * u[2]] };

        let b1 = scale(2.0 * std::f64::consts::PI / vol, cross(a2, a3));
        let b2 = scale(2.0 * std::f64::consts::PI / vol, cross(a3, a1));
        let b3 = scale(2.0 * std::f64::consts::PI / vol, cross(a1, a2));

        [b1, b2, b3]
    }

    /// Interplanar d-spacing for Miller indices (h,k,l) in Å.
    ///
    /// Uses the general triclinic formula via the reciprocal lattice.
    pub fn d_spacing(&self, h: i32, k: i32, l: i32) -> f64 {
        let rec = self.reciprocal_lattice();
        // G = h*b1 + k*b2 + l*b3
        let hf = h as f64;
        let kf = k as f64;
        let lf = l as f64;
        let gx = hf * rec[0][0] + kf * rec[1][0] + lf * rec[2][0];
        let gy = hf * rec[0][1] + kf * rec[1][1] + lf * rec[2][1];
        let gz = hf * rec[0][2] + kf * rec[1][2] + lf * rec[2][2];
        let g_mag = (gx * gx + gy * gy + gz * gz).sqrt();
        if g_mag < 1e-12 {
            f64::INFINITY
        } else {
            2.0 * std::f64::consts::PI / g_mag
        }
    }

    /// Compute the kinematic structure factor F(h,k,l) = (Re, Im).
    ///
    /// Uses atomic scattering factors approximated by atomic number (simplified).
    pub fn structure_factor(&self, h: i32, k: i32, l: i32) -> (f64, f64) {
        let mut re = 0.0_f64;
        let mut im = 0.0_f64;

        for atom in &self.atoms {
            let [u, v, w] = atom.fractional_pos;
            let phase = 2.0 * std::f64::consts::PI * (h as f64 * u + k as f64 * v + l as f64 * w);
            // Simplified scattering factor: f ≈ atomic_number(element) * exp(-B*sin²θ/λ²)
            let f = Self::approx_scattering_factor(&atom.element) * atom.occupancy;
            re += f * phase.cos();
            im += f * phase.sin();
        }
        (re, im)
    }

    /// Approximate atomic scattering factor by element symbol (f(0)).
    fn approx_scattering_factor(element: &str) -> f64 {
        // Very simplified: use atomic number as proxy for f(0)
        match element.trim_end_matches(|c: char| c.is_ascii_digit()) {
            "H" => 1.0,
            "He" => 2.0,
            "Li" => 3.0,
            "Be" => 4.0,
            "B" => 5.0,
            "C" => 6.0,
            "N" => 7.0,
            "O" => 8.0,
            "F" => 9.0,
            "Na" => 11.0,
            "Mg" => 12.0,
            "Al" => 13.0,
            "Si" => 14.0,
            "P" => 15.0,
            "S" => 16.0,
            "Cl" => 17.0,
            "Ar" => 18.0,
            "K" => 19.0,
            "Ca" => 20.0,
            "Ti" => 22.0,
            "V" => 23.0,
            "Cr" => 24.0,
            "Mn" => 25.0,
            "Fe" => 26.0,
            "Co" => 27.0,
            "Ni" => 28.0,
            "Cu" => 29.0,
            "Zn" => 30.0,
            "Ga" => 31.0,
            "Ge" => 32.0,
            "As" => 33.0,
            "Se" => 34.0,
            "Br" => 35.0,
            "Sr" => 38.0,
            "Mo" => 42.0,
            "Ag" => 47.0,
            "Sn" => 50.0,
            "Ba" => 56.0,
            "La" => 57.0,
            "W" => 74.0,
            "Pt" => 78.0,
            "Au" => 79.0,
            "Hg" => 80.0,
            "Pb" => 82.0,
            "U" => 92.0,
            _ => 10.0, // fallback
        }
    }
}

// ---------------------------------------------------------------------------
// CIF writer
// ---------------------------------------------------------------------------

/// Write a [`CrystalStructure`] to a CIF-formatted string.
pub fn write_cif(crystal: &CrystalStructure) -> String {
    let mut out = String::new();

    out.push_str("# CIF written by oxiphysics-io\n");
    out.push_str("data_oxiphysics\n\n");

    // Cell parameters
    let a = crystal.lattice[0];
    let b = crystal.lattice[1];
    let c = crystal.lattice[2];

    let a_len = (a[0] * a[0] + a[1] * a[1] + a[2] * a[2]).sqrt();
    let b_len = (b[0] * b[0] + b[1] * b[1] + b[2] * b[2]).sqrt();
    let c_len = (c[0] * c[0] + c[1] * c[1] + c[2] * c[2]).sqrt();

    let dot = |u: [f64; 3], v: [f64; 3]| u[0] * v[0] + u[1] * v[1] + u[2] * v[2];
    let alpha = if b_len > 0.0 && c_len > 0.0 {
        (dot(b, c) / (b_len * c_len))
            .clamp(-1.0, 1.0)
            .acos()
            .to_degrees()
    } else {
        90.0
    };
    let beta = if a_len > 0.0 && c_len > 0.0 {
        (dot(a, c) / (a_len * c_len))
            .clamp(-1.0, 1.0)
            .acos()
            .to_degrees()
    } else {
        90.0
    };
    let gamma = if a_len > 0.0 && b_len > 0.0 {
        (dot(a, b) / (a_len * b_len))
            .clamp(-1.0, 1.0)
            .acos()
            .to_degrees()
    } else {
        90.0
    };

    out.push_str(&format!("_cell_length_a    {a_len:.6}\n"));
    out.push_str(&format!("_cell_length_b    {b_len:.6}\n"));
    out.push_str(&format!("_cell_length_c    {c_len:.6}\n"));
    out.push_str(&format!("_cell_angle_alpha {alpha:.4}\n"));
    out.push_str(&format!("_cell_angle_beta  {beta:.4}\n"));
    out.push_str(&format!("_cell_angle_gamma {gamma:.4}\n"));
    out.push_str(&format!(
        "_symmetry_int_tables_number {}\n\n",
        crystal.space_group
    ));

    // Atom loop
    out.push_str("loop_\n");
    out.push_str("  _atom_site_type_symbol\n");
    out.push_str("  _atom_site_fract_x\n");
    out.push_str("  _atom_site_fract_y\n");
    out.push_str("  _atom_site_fract_z\n");
    out.push_str("  _atom_site_occupancy\n");
    out.push_str("  _atom_site_b_iso\n");

    for atom in &crystal.atoms {
        out.push_str(&format!(
            "  {:4} {:10.6} {:10.6} {:10.6} {:6.4} {:8.4}\n",
            atom.element,
            atom.fractional_pos[0],
            atom.fractional_pos[1],
            atom.fractional_pos[2],
            atom.occupancy,
            atom.b_factor,
        ));
    }

    out
}

// ---------------------------------------------------------------------------
// VASP POSCAR reader (stub)
// ---------------------------------------------------------------------------

/// Reader for VASP POSCAR/CONTCAR format.
#[derive(Debug, Default)]
pub struct VaspReader;

impl VaspReader {
    /// Create a new `VaspReader`.
    pub fn new() -> Self {
        Self
    }

    /// Parse a POSCAR/CONTCAR file and return a [`CrystalStructure`].
    ///
    /// Currently a stub — handles the header, scale, lattice vectors, and
    /// Cartesian/Direct coordinate blocks.
    pub fn parse_poscar(&self, content: &str) -> Result<CrystalStructure> {
        let mut lines = content.lines().peekable();

        // Line 1: comment
        let _comment = lines.next().unwrap_or("").trim().to_string();
        // Line 2: scale factor
        let scale: f64 = lines
            .next()
            .unwrap_or("1.0")
            .trim()
            .parse()
            .map_err(|_| Error::Parse("POSCAR: invalid scale factor".into()))?;

        // Lines 3-5: lattice vectors
        let mut lattice = [[0.0_f64; 3]; 3];
        for row in lattice.iter_mut() {
            let line = lines
                .next()
                .ok_or_else(|| Error::Parse("POSCAR: missing lattice vector".into()))?;
            let vals: Vec<f64> = line
                .split_whitespace()
                .filter_map(|s| s.parse().ok())
                .collect();
            if vals.len() < 3 {
                return Err(Error::Parse("POSCAR: bad lattice vector".into()));
            }
            row[0] = vals[0] * scale;
            row[1] = vals[1] * scale;
            row[2] = vals[2] * scale;
        }

        // Line 6: element symbols (VASP5+) or direct species count
        let species_line = lines.next().unwrap_or("").trim().to_string();
        let (elements, count_line) = if species_line
            .chars()
            .next()
            .map(|c| c.is_alphabetic())
            .unwrap_or(false)
        {
            let elems: Vec<String> = species_line
                .split_whitespace()
                .map(|s| s.to_string())
                .collect();
            let cl = lines.next().unwrap_or("").trim().to_string();
            (elems, cl)
        } else {
            // VASP4: no element names
            (vec!["X".to_string()], species_line.clone())
        };

        let counts: Vec<usize> = count_line
            .split_whitespace()
            .filter_map(|s| s.parse().ok())
            .collect();

        // Line 7: Selective dynamics or coordinate type
        let coord_type_line = lines.next().unwrap_or("").trim().to_string();
        let coord_type_line = if coord_type_line.to_lowercase().starts_with('s') {
            lines.next().unwrap_or("").trim().to_string()
        } else {
            coord_type_line
        };
        let is_cartesian = coord_type_line.to_lowercase().starts_with('c')
            || coord_type_line.to_lowercase().starts_with('k');

        // Atoms
        let mut atoms = Vec::new();
        for (ei, &count) in counts.iter().enumerate() {
            let element = elements.get(ei).cloned().unwrap_or_else(|| "X".to_string());
            for _ in 0..count {
                let line = lines.next().unwrap_or("0 0 0");
                let vals: Vec<f64> = line
                    .split_whitespace()
                    .take(3)
                    .filter_map(|s| s.parse().ok())
                    .collect();
                let mut pos = [
                    vals.first().copied().unwrap_or(0.0),
                    vals.get(1).copied().unwrap_or(0.0),
                    vals.get(2).copied().unwrap_or(0.0),
                ];
                if is_cartesian {
                    // Convert Cartesian to fractional (simplified — assumes orthogonal)
                    pos[0] /= lattice[0][0].abs().max(1e-12);
                    pos[1] /= lattice[1][1].abs().max(1e-12);
                    pos[2] /= lattice[2][2].abs().max(1e-12);
                }
                atoms.push(CrystalAtom {
                    element: element.clone(),
                    fractional_pos: pos,
                    occupancy: 1.0,
                    b_factor: 0.0,
                });
            }
        }

        Ok(CrystalStructure {
            lattice,
            atoms,
            space_group: 1,
        })
    }
}

// ---------------------------------------------------------------------------
// XRD pattern analysis
// ---------------------------------------------------------------------------

impl XrdPattern {
    /// Create an XRD pattern from parallel 2θ and intensity arrays.
    pub fn new(two_theta: Vec<f64>, intensity: Vec<f64>) -> Self {
        Self {
            two_theta,
            intensity,
        }
    }

    /// Find peak positions (local maxima) using a simple derivative method.
    ///
    /// Returns indices into `two_theta`/`intensity` where peaks occur.
    pub fn find_peaks(&self) -> Vec<usize> {
        let n = self.intensity.len();
        if n < 3 {
            return vec![];
        }
        let mut peaks = Vec::new();
        for i in 1..n - 1 {
            if self.intensity[i] > self.intensity[i - 1]
                && self.intensity[i] > self.intensity[i + 1]
            {
                peaks.push(i);
            }
        }
        peaks
    }

    /// Convert 2θ peak positions to d-spacings using Bragg's law: d = λ/(2 sin θ).
    ///
    /// Uses Cu Kα radiation (λ = 1.5406 Å) by default.
    pub fn d_spacings(&self) -> Vec<f64> {
        const LAMBDA: f64 = 1.5406; // Å, Cu Kα
        self.find_peaks()
            .iter()
            .map(|&i| {
                let theta = self.two_theta[i].to_radians() / 2.0;
                let sin_t = theta.sin();
                if sin_t > 1e-12 {
                    LAMBDA / (2.0 * sin_t)
                } else {
                    f64::INFINITY
                }
            })
            .collect()
    }
}

// ---------------------------------------------------------------------------
// XRD simulation
// ---------------------------------------------------------------------------

/// Simulate a powder XRD pattern for the given crystal.
///
/// `wavelength` in Å, `theta_range` in degrees (2θ_min, 2θ_max), `n_pts` sample points.
pub fn simulate_xrd(
    crystal: &CrystalStructure,
    wavelength: f64,
    theta_range: (f64, f64),
    n_pts: usize,
) -> XrdPattern {
    let (two_theta_min, two_theta_max) = theta_range;
    let step = if n_pts > 1 {
        (two_theta_max - two_theta_min) / (n_pts - 1) as f64
    } else {
        0.0
    };

    let two_theta: Vec<f64> = (0..n_pts)
        .map(|i| two_theta_min + i as f64 * step)
        .collect();
    let mut intensity = vec![0.0_f64; n_pts];

    // Consider reflections up to h,k,l = ±max_hkl
    let max_hkl = 5_i32;
    let vol = crystal.volume();
    if vol < 1e-12 {
        return XrdPattern::new(two_theta, intensity);
    }

    for h in -max_hkl..=max_hkl {
        for k in -max_hkl..=max_hkl {
            for l in -max_hkl..=max_hkl {
                if h == 0 && k == 0 && l == 0 {
                    continue;
                }
                let d = crystal.d_spacing(h, k, l);
                if d.is_infinite() || d < 0.5 {
                    continue;
                }

                // Bragg angle
                let sin_theta = wavelength / (2.0 * d);
                if sin_theta > 1.0 {
                    continue;
                }
                let two_theta_hkl = 2.0 * sin_theta.asin().to_degrees();

                // Find nearest grid point
                if two_theta_hkl < two_theta_min || two_theta_hkl > two_theta_max {
                    continue;
                }
                let idx = ((two_theta_hkl - two_theta_min) / step.max(1e-12)) as usize;
                if idx >= n_pts {
                    continue;
                }

                // Structure factor intensity
                let (fre, fim) = crystal.structure_factor(h, k, l);
                let f2 = fre * fre + fim * fim;

                // Lorentz-polarization factor
                let theta = two_theta_hkl.to_radians() / 2.0;
                let lp = (1.0 + (2.0 * theta).cos().powi(2))
                    / (theta.sin().powi(2) * (2.0 * theta).cos()).abs().max(1e-12);

                // Spread over a few points with a Gaussian (FWHM ~ 0.1°)
                let sigma_pts = 2.0;
                for di in -(3 * sigma_pts as i64)..=(3 * sigma_pts as i64) {
                    let j = idx as i64 + di;
                    if j < 0 || j >= n_pts as i64 {
                        continue;
                    }
                    let dt = di as f64 / sigma_pts;
                    intensity[j as usize] += f2 * lp * (-0.5 * dt * dt).exp();
                }
            }
        }
    }

    XrdPattern::new(two_theta, intensity)
}

// ---------------------------------------------------------------------------
// Tests
// ---------------------------------------------------------------------------

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

    // --- CifParser ---

    #[test]
    fn test_cif_parser_minimal_cubic() {
        let cif = r#"
data_test
_cell_length_a 5.0
_cell_length_b 5.0
_cell_length_c 5.0
_cell_angle_alpha 90.0
_cell_angle_beta  90.0
_cell_angle_gamma 90.0
_symmetry_int_tables_number 225
"#;
        let parser = CifParser::new();
        let crystal = parser.parse(cif).unwrap();
        assert!((crystal.volume() - 125.0).abs() < 0.01);
        assert_eq!(crystal.space_group, 225);
    }

    #[test]
    fn test_cif_parser_space_group_fallback() {
        let cif = "data_x\n_cell_length_a 3.0\n_cell_length_b 3.0\n_cell_length_c 3.0\n";
        let parser = CifParser::new();
        let crystal = parser.parse(cif).unwrap();
        assert_eq!(crystal.space_group, 1);
    }

    #[test]
    fn test_cif_parser_with_standard_uncertainty() {
        let cif =
            "data_x\n_cell_length_a 4.123(5)\n_cell_length_b 4.123(5)\n_cell_length_c 4.123(5)\n";
        let parser = CifParser::new();
        let crystal = parser.parse(cif).unwrap();
        assert!((crystal.lattice[0][0] - 4.123).abs() < 0.001);
    }

    #[test]
    fn test_volume_cubic() {
        let crystal = CrystalStructure::new_cubic(3.0);
        assert!((crystal.volume() - 27.0).abs() < 1e-10);
    }

    #[test]
    fn test_volume_orthorhombic() {
        let crystal = CrystalStructure {
            lattice: [[2.0, 0.0, 0.0], [0.0, 3.0, 0.0], [0.0, 0.0, 4.0]],
            atoms: vec![],
            space_group: 1,
        };
        assert!((crystal.volume() - 24.0).abs() < 1e-10);
    }

    #[test]
    fn test_to_cartesian_fractional_origin() {
        let mut crystal = CrystalStructure::new_cubic(4.0);
        crystal.atoms.push(CrystalAtom {
            element: "Fe".into(),
            fractional_pos: [0.0, 0.0, 0.0],
            occupancy: 1.0,
            b_factor: 0.0,
        });
        let cart = crystal.to_cartesian();
        assert!((cart[0][0]).abs() < 1e-12);
        assert!((cart[0][1]).abs() < 1e-12);
        assert!((cart[0][2]).abs() < 1e-12);
    }

    #[test]
    fn test_to_cartesian_fractional_half() {
        let mut crystal = CrystalStructure::new_cubic(4.0);
        crystal.atoms.push(CrystalAtom {
            element: "O".into(),
            fractional_pos: [0.5, 0.5, 0.5],
            occupancy: 1.0,
            b_factor: 0.0,
        });
        let cart = crystal.to_cartesian();
        assert!((cart[0][0] - 2.0).abs() < 1e-10);
        assert!((cart[0][1] - 2.0).abs() < 1e-10);
        assert!((cart[0][2] - 2.0).abs() < 1e-10);
    }

    #[test]
    fn test_d_spacing_cubic_100() {
        // For cubic with a=3: d(1,0,0) = a = 3 Å
        let crystal = CrystalStructure::new_cubic(3.0);
        let d = crystal.d_spacing(1, 0, 0);
        assert!((d - 3.0).abs() < 0.001, "d_100 = {d}");
    }

    #[test]
    fn test_d_spacing_cubic_110() {
        // d(1,1,0) = a/√2
        let a = 4.0;
        let crystal = CrystalStructure::new_cubic(a);
        let d = crystal.d_spacing(1, 1, 0);
        let expected = a / 2.0_f64.sqrt();
        assert!(
            (d - expected).abs() < 0.001,
            "d_110 = {d}, expected {expected}"
        );
    }

    #[test]
    fn test_d_spacing_cubic_111() {
        // d(1,1,1) = a/√3
        let a = 3.0;
        let crystal = CrystalStructure::new_cubic(a);
        let d = crystal.d_spacing(1, 1, 1);
        let expected = a / 3.0_f64.sqrt();
        assert!(
            (d - expected).abs() < 0.001,
            "d_111 = {d}, expected {expected}"
        );
    }

    #[test]
    fn test_reciprocal_lattice_cubic() {
        let a = 2.0 * std::f64::consts::PI;
        let crystal = CrystalStructure::new_cubic(a);
        let rec = crystal.reciprocal_lattice();
        // For cubic: b1 = 2π/a * x̂ → magnitude = 1.0
        assert!((rec[0][0] - 1.0).abs() < 1e-10);
        assert!((rec[1][1] - 1.0).abs() < 1e-10);
        assert!((rec[2][2] - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_structure_factor_single_atom_at_origin() {
        let mut crystal = CrystalStructure::new_cubic(4.0);
        crystal.atoms.push(CrystalAtom {
            element: "Fe".into(),
            fractional_pos: [0.0, 0.0, 0.0],
            occupancy: 1.0,
            b_factor: 0.0,
        });
        // Any hkl: phase = 0, so F = f(Fe) + 0i
        let (re, im) = crystal.structure_factor(1, 0, 0);
        assert!(re > 0.0, "real part should be positive");
        assert!(im.abs() < 1e-10, "imaginary part should be 0");
    }

    #[test]
    fn test_density_positive() {
        let mut crystal = CrystalStructure::new_cubic(3.615); // copper
        for _ in 0..4 {
            crystal.atoms.push(CrystalAtom {
                element: "Cu".into(),
                fractional_pos: [0.0, 0.0, 0.0],
                occupancy: 1.0,
                b_factor: 0.0,
            });
        }
        let rho = crystal.density(63.546); // Cu molar mass
        assert!(rho > 0.0, "density must be positive");
        // copper density ~8.9 g/cm³
        assert!(
            rho > 5.0 && rho < 15.0,
            "density {rho} out of expected range"
        );
    }

    #[test]
    fn test_write_cif_roundtrip_cell() {
        let crystal = CrystalStructure::new_cubic(5.0);
        let cif = write_cif(&crystal);
        assert!(cif.contains("_cell_length_a"));
        assert!(cif.contains("5.000000"));
    }

    #[test]
    fn test_write_cif_contains_atom_loop() {
        let mut crystal = CrystalStructure::new_cubic(4.0);
        crystal.atoms.push(CrystalAtom {
            element: "Na".into(),
            fractional_pos: [0.0, 0.0, 0.0],
            occupancy: 1.0,
            b_factor: 0.5,
        });
        let cif = write_cif(&crystal);
        assert!(cif.contains("loop_"));
        assert!(cif.contains("Na"));
    }

    #[test]
    fn test_vasp_reader_simple_poscar() {
        let poscar = "Fe BCC\n1.0\n2.87 0.0 0.0\n0.0 2.87 0.0\n0.0 0.0 2.87\nFe\n2\nDirect\n0.0 0.0 0.0\n0.5 0.5 0.5\n";
        let reader = VaspReader::new();
        let crystal = reader.parse_poscar(poscar).unwrap();
        assert_eq!(crystal.atoms.len(), 2);
        assert_eq!(crystal.atoms[0].element, "Fe");
    }

    #[test]
    fn test_vasp_reader_lattice_scale() {
        let poscar =
            "Test\n2.0\n1.0 0.0 0.0\n0.0 1.0 0.0\n0.0 0.0 1.0\nH\n1\nDirect\n0.0 0.0 0.0\n";
        let reader = VaspReader::new();
        let crystal = reader.parse_poscar(poscar).unwrap();
        // Scale=2 applied to lattice
        assert!((crystal.lattice[0][0] - 2.0).abs() < 1e-10);
    }

    #[test]
    fn test_xrd_pattern_find_peaks() {
        // Manually create a pattern with clear peaks
        let two_theta: Vec<f64> = (0..20).map(|i| i as f64).collect();
        let mut intensity = vec![0.0; 20];
        intensity[5] = 100.0;
        intensity[10] = 200.0;
        let pattern = XrdPattern::new(two_theta, intensity);
        let peaks = pattern.find_peaks();
        assert!(peaks.contains(&5), "peak at index 5");
        assert!(peaks.contains(&10), "peak at index 10");
    }

    #[test]
    fn test_xrd_d_spacings_bragg() {
        // 2θ = 44.5° for Fe(110) with Cu Kα → d ~ 2.03 Å
        let two_theta = vec![40.0, 44.5, 50.0];
        let intensity = vec![5.0, 100.0, 5.0];
        let pattern = XrdPattern::new(two_theta, intensity);
        let ds = pattern.d_spacings();
        assert!(!ds.is_empty());
        // d for 2θ=44.5° ≈ 2.03 Å
        assert!(ds[0] > 1.5 && ds[0] < 3.0, "d = {}", ds[0]);
    }

    #[test]
    fn test_simulate_xrd_returns_correct_size() {
        let crystal = CrystalStructure::new_cubic(3.0);
        let pattern = simulate_xrd(&crystal, 1.5406, (20.0, 80.0), 100);
        assert_eq!(pattern.two_theta.len(), 100);
        assert_eq!(pattern.intensity.len(), 100);
    }

    #[test]
    fn test_simulate_xrd_non_negative_intensity() {
        let mut crystal = CrystalStructure::new_cubic(3.615);
        crystal.atoms.push(CrystalAtom {
            element: "Cu".into(),
            fractional_pos: [0.0, 0.0, 0.0],
            occupancy: 1.0,
            b_factor: 0.0,
        });
        let pattern = simulate_xrd(&crystal, 1.5406, (20.0, 80.0), 200);
        for &i in &pattern.intensity {
            assert!(i >= 0.0, "negative intensity: {i}");
        }
    }

    #[test]
    fn test_simulate_xrd_has_peaks() {
        let mut crystal = CrystalStructure::new_cubic(3.615);
        crystal.atoms.push(CrystalAtom {
            element: "Cu".into(),
            fractional_pos: [0.0, 0.0, 0.0],
            occupancy: 1.0,
            b_factor: 0.0,
        });
        let pattern = simulate_xrd(&crystal, 1.5406, (20.0, 100.0), 500);
        let peaks = pattern.find_peaks();
        assert!(!peaks.is_empty(), "simulated XRD should have peaks");
    }

    #[test]
    fn test_approx_scattering_factor_known() {
        assert!((CrystalStructure::approx_scattering_factor("Fe") - 26.0).abs() < 0.5);
        assert!((CrystalStructure::approx_scattering_factor("O") - 8.0).abs() < 0.5);
        assert!((CrystalStructure::approx_scattering_factor("Cu") - 29.0).abs() < 0.5);
    }

    #[test]
    fn test_build_lattice_cubic_is_diagonal() {
        let a = 4.0;
        let pi2 = std::f64::consts::PI / 2.0;
        let lat = CifParser::build_lattice(a, a, a, pi2, pi2, pi2);
        assert!((lat[0][0] - a).abs() < 1e-10);
        assert!(lat[0][1].abs() < 1e-10);
        assert!(lat[0][2].abs() < 1e-10);
        assert!((lat[1][1] - a).abs() < 1e-8, "b[1] = {}", lat[1][1]);
    }

    #[test]
    fn test_cif_tokenize_quoted() {
        let tokens = CifParser::tokenize_line("'hello world' 3.14");
        assert_eq!(tokens[0], "hello world");
        assert_eq!(tokens[1], "3.14");
    }

    #[test]
    fn test_strip_su() {
        assert_eq!(CifParser::strip_su("1.234(5)"), "1.234");
        assert_eq!(CifParser::strip_su("5.000"), "5.000");
    }

    #[test]
    fn test_d_spacing_000_is_infinity() {
        let crystal = CrystalStructure::new_cubic(3.0);
        let d = crystal.d_spacing(0, 0, 0);
        assert!(d.is_infinite());
    }

    #[test]
    fn test_crystal_structure_no_atoms_cartesian() {
        let crystal = CrystalStructure::new_cubic(5.0);
        assert!(crystal.to_cartesian().is_empty());
    }

    #[test]
    fn test_cif_parser_empty_returns_default() {
        let parser = CifParser::new();
        let crystal = parser.parse("").unwrap();
        assert!((crystal.volume() - 1.0).abs() < 1e-10); // 1×1×1
    }
}