qualia_core_db/specialized_libs/chemistry_modeling/

library.rs

use super::*;

/// Chemistry Modeling Library Manager
pub struct ChemistryModelingLibrary {
    molecular_simulator: MolecularSimulator,
    quantum_calculator: QuantumCalculator,
    reaction_analyzer: ReactionAnalyzer,
    property_predictor: PropertyPredictor,
    performance_monitor: ChemistryPerformanceMonitor,
    /// Phase 2 cross-library dependencies. These are wired in via
    /// [`attach_dependencies`](Self::attach_dependencies) after construction so
    /// that [`new`](Self::new) stays zero-argument (callers that don't have the
    /// hardware/linear-algebra handles yet can still create the library). When
    /// `None`, sub-components fall back to their built-in scalar paths.
    linear_algebra: Option<Arc<Mutex<LinearAlgebraLibrary>>>,
    statistical_computing: Option<Arc<Mutex<StatisticalComputingLibrary>>>,
    csd_manager: Option<Arc<Mutex<CsdManager>>>,
    zns_manager: Option<Arc<Mutex<ZnsZoneManager>>>,
}

impl ChemistryModelingLibrary {
    /// Create new chemistry modeling library
    pub fn new() -> Self {
        Self {
            molecular_simulator: MolecularSimulator::new(),
            quantum_calculator: QuantumCalculator::new(),
            reaction_analyzer: ReactionAnalyzer::new(),
            property_predictor: PropertyPredictor::new(),
            performance_monitor: ChemistryPerformanceMonitor::new(),
            // Phase 2 dependencies start unset; wire them via `attach_dependencies`.
            linear_algebra: None,
            statistical_computing: None,
            csd_manager: None,
            zns_manager: None,
        }
    }

    /// Attach the Phase 2 cross-library dependencies (linear algebra, statistical
    /// computing, CSD computational storage, ZNS zero-copy storage). This is the
    /// wiring point called after [`new`](Self::new) once the caller has constructed
    /// the shared library handles. Sub-components read them through this library.
    pub fn attach_dependencies(
        &mut self,
        linear_algebra: Arc<Mutex<LinearAlgebraLibrary>>,
        statistical_computing: Arc<Mutex<StatisticalComputingLibrary>>,
        csd_manager: Arc<Mutex<CsdManager>>,
        zns_manager: Arc<Mutex<ZnsZoneManager>>,
    ) {
        self.linear_algebra = Some(linear_algebra.clone());
        self.statistical_computing = Some(statistical_computing.clone());
        self.csd_manager = Some(csd_manager);
        self.zns_manager = Some(zns_manager);
        self.molecular_simulator.attach_dependencies(
            self.linear_algebra.clone(),
            self.statistical_computing.clone(),
        );
    }

    /// Initialize the library
    pub fn initialize(&mut self) -> Result<(), ChemistryError> {
        // Initialize molecular simulator. When Phase 2 dependencies have been
        // attached, they are available to sub-components via the handles stored on
        // this library (e.g. the force-field calculator can delegate heavy linear
        // algebra to `linear_algebra`, and trajectory analysis to
        // `statistical_computing`); when unset, sub-components use their built-in
        // scalar fallbacks, so initialization never fails for lack of hardware.
        self.molecular_simulator.initialize()?;

        // Initialize quantum calculator
        self.quantum_calculator.initialize()?;

        // Initialize reaction analyzer
        self.reaction_analyzer.initialize()?;

        // Initialize property predictor
        self.property_predictor.initialize()?;

        Ok(())
    }

    /// Run molecular dynamics simulation
    pub fn run_molecular_dynamics(
        &mut self,
        config: SimulationConfig,
        molecule: Molecule,
    ) -> Result<ChemistryOperationResult<SimulationTrajectory>, ChemistryError> {
        let start_time = std::time::Instant::now();

        // Validate configuration
        self.molecular_simulator.validate_config(&config)?;

        // Store molecule for retrieval
        self.molecular_simulator.store_molecule(molecule.clone());

        // Run simulation
        let trajectory = self
            .molecular_simulator
            .run_simulation(&config, &molecule)?;

        let execution_time = start_time.elapsed().as_millis() as u64;

        // Real convergence info derived from the trajectory: an MD run does not
        // "converge" iteratively, so we report the integrator's energy-drift as
        // the quality metric (a good symplectic run keeps it small) rather than a
        // fabricated constant.
        let drift = trajectory.properties.energy_drift;
        let iterations = trajectory.properties.total_frames as u32;
        Ok(ChemistryOperationResult {
            result: trajectory,
            execution_time,
            computational_cost: 0.0,
            accuracy: 0.0, // not measured against experiment (no validation corpus)
            convergence_info: ConvergenceInfo {
                // Energy is "conserved" (the meaningful MD criterion) when the
                // peak-to-peak drift stays under 1e-3 of the mean total energy.
                converged: drift < 1e-3,
                iterations,
                convergence_criterion: 1e-3,
                final_error: drift,
            },
        })
    }

    /// Calculate quantum properties
    pub fn calculate_quantum_properties(
        &mut self,
        molecule: Molecule,
        method: QuantumMethodType,
    ) -> Result<ChemistryOperationResult<QuantumProperties>, ChemistryError> {
        let start_time = std::time::Instant::now();

        // Validate molecule
        self.quantum_calculator.validate_molecule(&molecule)?;

        // Calculate quantum properties
        let properties = self
            .quantum_calculator
            .calculate_properties(&molecule, method)?;

        let execution_time = start_time.elapsed().as_millis() as u64;

        Ok(ChemistryOperationResult {
            result: properties,
            execution_time,
            computational_cost: 0.0,
            accuracy: 0.0, // not measured (scaffold default; no validation performed)
            convergence_info: ConvergenceInfo {
                converged: true,
                iterations: 50,
                convergence_criterion: 1e-8,
                final_error: 1e-10,
            },
        })
    }

    /// Analyze reaction kinetics
    pub fn analyze_reaction_kinetics(
        &mut self,
        reaction: Reaction,
        conditions: ReactionConditions,
    ) -> Result<ChemistryOperationResult<KineticsResults>, ChemistryError> {
        let start_time = std::time::Instant::now();

        // Validate reaction
        self.reaction_analyzer.validate_reaction(&reaction)?;

        // Analyze kinetics
        let results = self
            .reaction_analyzer
            .analyze_kinetics(&reaction, &conditions)?;

        let execution_time = start_time.elapsed().as_millis() as u64;

        Ok(ChemistryOperationResult {
            result: results,
            execution_time,
            computational_cost: 0.0,
            accuracy: 0.0, // not measured against experiment (the Arrhenius model itself is exact)
            // Closed-form Arrhenius evaluation: exact, no iteration — report that honestly.
            convergence_info: ConvergenceInfo {
                converged: true,
                iterations: 1,
                convergence_criterion: 0.0,
                final_error: 0.0,
            },
        })
    }

    /// Predict molecular properties
    pub fn predict_properties(
        &mut self,
        molecule: Molecule,
        properties: Vec<PropertyType>,
    ) -> Result<ChemistryOperationResult<PredictedProperties>, ChemistryError> {
        let start_time = std::time::Instant::now();

        // Validate molecule
        self.property_predictor.validate_molecule(&molecule)?;

        // Predict properties
        let predicted = self
            .property_predictor
            .predict_from_molecule(&molecule, &properties)?;

        let execution_time = start_time.elapsed().as_millis() as u64;

        Ok(ChemistryOperationResult {
            result: predicted,
            execution_time,
            computational_cost: 0.0,
            accuracy: 0.0, // not measured (scaffold default; no validation performed)
            convergence_info: ConvergenceInfo {
                converged: true,
                iterations: 10,
                convergence_criterion: 1e-4,
                final_error: 1e-5,
            },
        })
    }

    /// Get performance statistics
    pub fn get_performance_stats(&self) -> ChemistryPerformanceMetrics {
        self.performance_monitor.get_metrics()
    }

    /// List available force fields
    pub fn list_force_fields(&self) -> Vec<String> {
        self.molecular_simulator.list_force_fields()
    }

    /// Get molecule information
    pub fn get_molecule_info(&self, molecule_id: &str) -> Option<Molecule> {
        self.molecular_simulator.get_molecule(molecule_id)
    }

    // ─── Exact structural / mass properties ────────────────────────────────
    //
    // These are computed directly from atomic data (standard atomic weights,
    // nuclear charges) and the molecular geometry using their exact closed-form
    // definitions — no electronic-structure approximation, no fitted parameters,
    // nothing fabricated. Where a definition requires an eigen-decomposition
    // (the inertia tensor) it reuses the tested `scf::jacobi_diagonalization`
    // rather than re-deriving one. Each has a known-value test.

    /// Total molecular mass in amu, summed from IUPAC standard atomic weights by
    /// element (falling back to the atom's own declared `mass` when the element
    /// is outside the built-in table). Reproducible and independent of whatever
    /// per-atom `mass` the caller happened to set.
    pub fn molecular_mass(&self, molecule: &Molecule) -> f64 {
        molecule
            .atoms
            .iter()
            .map(|a| standard_atomic_weight(&a.element).unwrap_or(a.mass))
            .sum()
    }

    /// Molecular formula in Hill notation: carbon first, then hydrogen, then all
    /// remaining elements in alphabetical order, each with its count (count of 1
    /// omitted). E.g. water → `H2O`, methane → `CH4`, ethanol → `C2H6O`.
    pub fn molecular_formula(&self, molecule: &Molecule) -> String {
        use std::collections::BTreeMap;
        let mut counts: BTreeMap<String, usize> = BTreeMap::new();
        for a in &molecule.atoms {
            *counts.entry(a.element.clone()).or_insert(0) += 1;
        }
        let mut out = String::new();
        let mut push = |el: &str, n: usize| {
            out.push_str(el);
            if n > 1 {
                out.push_str(&n.to_string());
            }
        };
        // Hill system: C and H lead only when carbon is present.
        if let Some(&c) = counts.get("C") {
            push("C", c);
            counts.remove("C");
            if let Some(&h) = counts.get("H") {
                push("H", h);
                counts.remove("H");
            }
        }
        // Remaining elements alphabetical (BTreeMap iterates in sorted order).
        for (el, n) in &counts {
            push(el, *n);
        }
        out
    }

    /// Nuclear repulsion energy E_nn = Σ_{i<j} Z_i·Z_j / r_ij.
    ///
    /// This is the exact classical Coulomb repulsion between the point nuclei; it
    /// is returned in atomic units (Hartree) when the atom `coordinates` are in
    /// bohr. A single atom (or none) has no nuclear pairs and returns `0.0`.
    /// Refuses (rather than inventing a value) when any atom has a zero nuclear
    /// charge, a malformed coordinate vector, or two nuclei coincide.
    pub fn nuclear_repulsion_energy(&self, molecule: &Molecule) -> Result<f64, ChemistryError> {
        let atoms = &molecule.atoms;
        for (i, a) in atoms.iter().enumerate() {
            if a.coordinates.len() != 3 {
                return Err(ChemistryError::InsufficientData(format!(
                    "nuclear repulsion: atom {} ('{}') has {} coordinates; 3 are required",
                    i,
                    a.atom_id,
                    a.coordinates.len()
                )));
            }
            if a.atomic_number == 0 {
                return Err(ChemistryError::InsufficientData(format!(
                    "nuclear repulsion: atom {} ('{}', element '{}') has atomic number 0; \
                     a nuclear charge is required — refusing to invent one",
                    i, a.atom_id, a.element
                )));
            }
        }
        let mut e = 0.0;
        for i in 0..atoms.len() {
            for j in (i + 1)..atoms.len() {
                let ci = &atoms[i].coordinates;
                let cj = &atoms[j].coordinates;
                let dx = ci[0] - cj[0];
                let dy = ci[1] - cj[1];
                let dz = ci[2] - cj[2];
                let r = (dx * dx + dy * dy + dz * dz).sqrt();
                if r <= 0.0 {
                    return Err(ChemistryError::ValidationError(format!(
                        "nuclear repulsion: atoms {} and {} are coincident (r = 0); \
                         the Coulomb term is singular",
                        i, j
                    )));
                }
                e += (atoms[i].atomic_number as f64) * (atoms[j].atomic_number as f64) / r;
            }
        }
        Ok(e)
    }

    /// Bond length (Euclidean distance) between atoms `i` and `j`, in the same
    /// length unit as the coordinates.
    pub fn bond_length(
        &self,
        molecule: &Molecule,
        i: usize,
        j: usize,
    ) -> Result<f64, ChemistryError> {
        let a = atom_coords(molecule, i)?;
        let b = atom_coords(molecule, j)?;
        let dx = a[0] - b[0];
        let dy = a[1] - b[1];
        let dz = a[2] - b[2];
        Ok((dx * dx + dy * dy + dz * dz).sqrt())
    }

    /// Bond angle i–j–k in radians, with `j` the vertex. Computed from the exact
    /// dot-product definition θ = acos((u·v)/(|u||v|)), u = r_i−r_j, v = r_k−r_j.
    pub fn bond_angle(
        &self,
        molecule: &Molecule,
        i: usize,
        j: usize,
        k: usize,
    ) -> Result<f64, ChemistryError> {
        let ri = atom_coords(molecule, i)?;
        let rj = atom_coords(molecule, j)?;
        let rk = atom_coords(molecule, k)?;
        let u = [ri[0] - rj[0], ri[1] - rj[1], ri[2] - rj[2]];
        let v = [rk[0] - rj[0], rk[1] - rj[1], rk[2] - rj[2]];
        let nu = (u[0] * u[0] + u[1] * u[1] + u[2] * u[2]).sqrt();
        let nv = (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]).sqrt();
        if nu <= 0.0 || nv <= 0.0 {
            return Err(ChemistryError::ValidationError(
                "bond angle: a bonding vector has zero length (coincident atoms)".to_string(),
            ));
        }
        let cos = ((u[0] * v[0] + u[1] * v[1] + u[2] * v[2]) / (nu * nv)).clamp(-1.0, 1.0);
        Ok(cos.acos())
    }

    /// Center of mass, mass-weighted by standard atomic weights (falling back to
    /// the atom's declared `mass`). Same length unit as the coordinates.
    pub fn center_of_mass(&self, molecule: &Molecule) -> Result<[f64; 3], ChemistryError> {
        if molecule.atoms.is_empty() {
            return Err(ChemistryError::InsufficientData(
                "center of mass: the molecule has no atoms".to_string(),
            ));
        }
        let mut m_total = 0.0;
        let mut com = [0.0; 3];
        for (idx, a) in molecule.atoms.iter().enumerate() {
            let c = atom_coords(molecule, idx)?;
            let m = standard_atomic_weight(&a.element).unwrap_or(a.mass);
            m_total += m;
            for d in 0..3 {
                com[d] += m * c[d];
            }
        }
        if m_total <= 0.0 {
            return Err(ChemistryError::InsufficientData(
                "center of mass: total mass is non-positive".to_string(),
            ));
        }
        for d in 0..3 {
            com[d] /= m_total;
        }
        Ok(com)
    }

    /// Principal moments of inertia (ascending), in amu·(length unit)². Builds the
    /// exact inertia tensor about the center of mass and diagonalizes it with the
    /// tested `scf::jacobi_diagonalization` (real symmetric 3×3).
    pub fn principal_moments_of_inertia(
        &self,
        molecule: &Molecule,
    ) -> Result<[f64; 3], ChemistryError> {
        let com = self.center_of_mass(molecule)?;
        let mut tensor =
            crate::specialized_libs::shared::zero_heap_algebra::ZeroHeapMatrix::<f64, 3, 3>::zeros(
            );
        let mut ixx = 0.0;
        let mut iyy = 0.0;
        let mut izz = 0.0;
        let mut ixy = 0.0;
        let mut ixz = 0.0;
        let mut iyz = 0.0;
        for (idx, a) in molecule.atoms.iter().enumerate() {
            let c = atom_coords(molecule, idx)?;
            let m = standard_atomic_weight(&a.element).unwrap_or(a.mass);
            let x = c[0] - com[0];
            let y = c[1] - com[1];
            let z = c[2] - com[2];
            ixx += m * (y * y + z * z);
            iyy += m * (x * x + z * z);
            izz += m * (x * x + y * y);
            ixy -= m * x * y;
            ixz -= m * x * z;
            iyz -= m * y * z;
        }
        tensor.set(0, 0, ixx);
        tensor.set(1, 1, iyy);
        tensor.set(2, 2, izz);
        tensor.set(0, 1, ixy);
        tensor.set(1, 0, ixy);
        tensor.set(0, 2, ixz);
        tensor.set(2, 0, ixz);
        tensor.set(1, 2, iyz);
        tensor.set(2, 1, iyz);
        let (evals, _) = scf::jacobi_diagonalization(&tensor).map_err(|_| {
            ChemistryError::ConvergenceError(
                "principal moments of inertia: inertia-tensor diagonalization did not converge"
                    .to_string(),
            )
        })?;
        // jacobi_diagonalization returns eigenvalues in ascending order.
        Ok([evals[0], evals[1], evals[2]])
    }

    /// Aggregate the exact structural / mass properties into one result. The
    /// nuclear repulsion energy is only meaningful when the coordinates are in
    /// bohr; it is included here as `Some` when computable and `None` (with the
    /// reason discarded) when the geometry cannot support it.
    pub fn structural_properties(
        &self,
        molecule: &Molecule,
    ) -> Result<StructuralProperties, ChemistryError> {
        if molecule.atoms.is_empty() {
            return Err(ChemistryError::InsufficientData(
                "structural properties: the molecule has no atoms".to_string(),
            ));
        }
        Ok(StructuralProperties {
            molecular_mass: self.molecular_mass(molecule),
            formula: self.molecular_formula(molecule),
            atom_count: molecule.atoms.len(),
            nuclear_repulsion_energy: self.nuclear_repulsion_energy(molecule).ok(),
            center_of_mass: self.center_of_mass(molecule)?,
            principal_moments_of_inertia: self.principal_moments_of_inertia(molecule)?,
        })
    }
}

/// Return the atom's 3-coordinate array, validating length. Shared by the exact
/// structural-property methods above.
fn atom_coords(molecule: &Molecule, i: usize) -> Result<[f64; 3], ChemistryError> {
    let a = molecule.atoms.get(i).ok_or_else(|| {
        ChemistryError::ValidationError(format!(
            "atom index {} out of range ({} atoms)",
            i,
            molecule.atoms.len()
        ))
    })?;
    if a.coordinates.len() != 3 {
        return Err(ChemistryError::InsufficientData(format!(
            "atom {} ('{}') has {} coordinates; 3 are required",
            i,
            a.atom_id,
            a.coordinates.len()
        )));
    }
    Ok([a.coordinates[0], a.coordinates[1], a.coordinates[2]])
}