sci_form/
lib.rs

1// Scientific/numerical code patterns that are idiomatic in this domain
2#![allow(clippy::too_many_arguments)]
3#![allow(clippy::needless_range_loop)]
4
5pub mod alignment;
6pub mod charges;
7pub mod conformer;
8pub mod dipole;
9pub mod distgeom;
10pub mod dos;
11pub mod eht;
12pub mod esp;
13pub mod etkdg;
14pub mod forcefield;
15pub mod graph;
16pub mod materials;
17pub mod optimization;
18pub mod population;
19pub mod smarts;
20pub mod smiles;
21pub mod surface;
22pub mod transport;
23
24use serde::{Deserialize, Serialize};
25
26// ─── Public API Types ────────────────────────────────────────────────────────
27
28/// Result of a 3D conformer generation for a single molecule.
29#[derive(Debug, Clone, Serialize, Deserialize)]
30pub struct ConformerResult {
31    /// Input SMILES string.
32    pub smiles: String,
33    /// Number of atoms in the molecule.
34    pub num_atoms: usize,
35    /// Flat xyz coordinates: [x0, y0, z0, x1, y1, z1, ...].
36    /// Empty on failure.
37    pub coords: Vec<f64>,
38    /// Atom elements (atomic numbers) in the same order as coords.
39    pub elements: Vec<u8>,
40    /// Bond list as (start_atom, end_atom, order_string).
41    pub bonds: Vec<(usize, usize, String)>,
42    /// Error message if generation failed.
43    pub error: Option<String>,
44    /// Generation time in milliseconds.
45    pub time_ms: f64,
46}
47
48/// Configuration for conformer generation.
49#[derive(Debug, Clone, Serialize, Deserialize)]
50pub struct ConformerConfig {
51    /// RNG seed (same seed = reproducible output).
52    pub seed: u64,
53    /// Number of threads for batch processing (0 = auto-detect).
54    pub num_threads: usize,
55}
56
57impl Default for ConformerConfig {
58    fn default() -> Self {
59        Self {
60            seed: 42,
61            num_threads: 0,
62        }
63    }
64}
65
66// ─── Public API Functions ────────────────────────────────────────────────────
67
68/// Library version string.
69pub fn version() -> String {
70    format!("sci-form {}", env!("CARGO_PKG_VERSION"))
71}
72
73/// Generate a 3D conformer from a SMILES string.
74pub fn embed(smiles: &str, seed: u64) -> ConformerResult {
75    #[cfg(not(target_arch = "wasm32"))]
76    let start = std::time::Instant::now();
77
78    let mol = match graph::Molecule::from_smiles(smiles) {
79        Ok(m) => m,
80        Err(e) => {
81            return ConformerResult {
82                smiles: smiles.to_string(),
83                num_atoms: 0,
84                coords: vec![],
85                elements: vec![],
86                bonds: vec![],
87                error: Some(e),
88                #[cfg(not(target_arch = "wasm32"))]
89                time_ms: start.elapsed().as_secs_f64() * 1000.0,
90                #[cfg(target_arch = "wasm32")]
91                time_ms: 0.0,
92            };
93        }
94    };
95
96    let n = mol.graph.node_count();
97    let elements: Vec<u8> = (0..n)
98        .map(|i| mol.graph[petgraph::graph::NodeIndex::new(i)].element)
99        .collect();
100    let bonds: Vec<(usize, usize, String)> = mol
101        .graph
102        .edge_indices()
103        .map(|e| {
104            let (a, b) = mol.graph.edge_endpoints(e).unwrap();
105            let order = match mol.graph[e].order {
106                graph::BondOrder::Single => "SINGLE",
107                graph::BondOrder::Double => "DOUBLE",
108                graph::BondOrder::Triple => "TRIPLE",
109                graph::BondOrder::Aromatic => "AROMATIC",
110                graph::BondOrder::Unknown => "UNKNOWN",
111            };
112            (a.index(), b.index(), order.to_string())
113        })
114        .collect();
115
116    match conformer::generate_3d_conformer(&mol, seed) {
117        Ok(coords) => {
118            let mut flat = Vec::with_capacity(n * 3);
119            for i in 0..n {
120                flat.push(coords[(i, 0)] as f64);
121                flat.push(coords[(i, 1)] as f64);
122                flat.push(coords[(i, 2)] as f64);
123            }
124            ConformerResult {
125                smiles: smiles.to_string(),
126                num_atoms: n,
127                coords: flat,
128                elements,
129                bonds,
130                error: None,
131                #[cfg(not(target_arch = "wasm32"))]
132                time_ms: start.elapsed().as_secs_f64() * 1000.0,
133                #[cfg(target_arch = "wasm32")]
134                time_ms: 0.0,
135            }
136        }
137        Err(e) => ConformerResult {
138            smiles: smiles.to_string(),
139            num_atoms: n,
140            coords: vec![],
141            elements,
142            bonds,
143            error: Some(e),
144            #[cfg(not(target_arch = "wasm32"))]
145            time_ms: start.elapsed().as_secs_f64() * 1000.0,
146            #[cfg(target_arch = "wasm32")]
147            time_ms: 0.0,
148        },
149    }
150}
151
152/// Batch-embed multiple SMILES in parallel.
153///
154/// Uses rayon thread pool for parallel processing.
155/// `config.num_threads = 0` means auto-detect CPU count.
156#[cfg(feature = "parallel")]
157pub fn embed_batch(smiles_list: &[&str], config: &ConformerConfig) -> Vec<ConformerResult> {
158    use rayon::prelude::*;
159
160    if config.num_threads > 0 {
161        let pool = rayon::ThreadPoolBuilder::new()
162            .num_threads(config.num_threads)
163            .build()
164            .unwrap();
165        pool.install(|| {
166            smiles_list
167                .par_iter()
168                .map(|smi| embed(smi, config.seed))
169                .collect()
170        })
171    } else {
172        smiles_list
173            .par_iter()
174            .map(|smi| embed(smi, config.seed))
175            .collect()
176    }
177}
178
179/// Batch-embed multiple SMILES sequentially (no rayon dependency).
180#[cfg(not(feature = "parallel"))]
181pub fn embed_batch(smiles_list: &[&str], config: &ConformerConfig) -> Vec<ConformerResult> {
182    smiles_list
183        .iter()
184        .map(|smi| embed(smi, config.seed))
185        .collect()
186}
187
188/// Parse a SMILES string and return molecular structure (no 3D generation).
189pub fn parse(smiles: &str) -> Result<graph::Molecule, String> {
190    graph::Molecule::from_smiles(smiles)
191}
192
193/// Compute Gasteiger-Marsili partial charges from a SMILES string.
194///
195/// Parses the molecule, extracts bonds and elements, then runs 6 iterations
196/// of electronegativity equalization.
197pub fn compute_charges(smiles: &str) -> Result<charges::gasteiger::ChargeResult, String> {
198    let mol = graph::Molecule::from_smiles(smiles)?;
199    let n = mol.graph.node_count();
200    let elements: Vec<u8> = (0..n)
201        .map(|i| mol.graph[petgraph::graph::NodeIndex::new(i)].element)
202        .collect();
203    let formal_charges: Vec<i8> = (0..n)
204        .map(|i| mol.graph[petgraph::graph::NodeIndex::new(i)].formal_charge)
205        .collect();
206    let bonds: Vec<(usize, usize)> = mol
207        .graph
208        .edge_indices()
209        .map(|e| {
210            let (a, b) = mol.graph.edge_endpoints(e).unwrap();
211            (a.index(), b.index())
212        })
213        .collect();
214    charges::gasteiger::gasteiger_marsili_charges(&elements, &bonds, &formal_charges, 6)
215}
216
217/// Compute solvent-accessible surface area from SMILES + 3D coordinates.
218///
219/// The `coords_flat` parameter is a flat [x0,y0,z0, x1,y1,z1,...] array.
220pub fn compute_sasa(
221    elements: &[u8],
222    coords_flat: &[f64],
223    probe_radius: Option<f64>,
224) -> Result<surface::sasa::SasaResult, String> {
225    if coords_flat.len() != elements.len() * 3 {
226        return Err(format!(
227            "coords length {} != 3 * elements {}",
228            coords_flat.len(),
229            elements.len()
230        ));
231    }
232    let positions: Vec<[f64; 3]> = coords_flat
233        .chunks_exact(3)
234        .map(|c| [c[0], c[1], c[2]])
235        .collect();
236    Ok(surface::sasa::compute_sasa(
237        elements,
238        &positions,
239        probe_radius,
240        None,
241    ))
242}
243
244/// Compute Mulliken & Löwdin population analysis from atomic elements and positions.
245pub fn compute_population(
246    elements: &[u8],
247    positions: &[[f64; 3]],
248) -> Result<population::PopulationResult, String> {
249    let eht_result = eht::solve_eht(elements, positions, None)?;
250    Ok(population::compute_population(
251        elements,
252        positions,
253        &eht_result.coefficients,
254        eht_result.n_electrons,
255    ))
256}
257
258/// Compute molecular dipole moment (Debye) from atomic elements and positions.
259pub fn compute_dipole(
260    elements: &[u8],
261    positions: &[[f64; 3]],
262) -> Result<dipole::DipoleResult, String> {
263    let eht_result = eht::solve_eht(elements, positions, None)?;
264    Ok(dipole::compute_dipole_from_eht(
265        elements,
266        positions,
267        &eht_result.coefficients,
268        eht_result.n_electrons,
269    ))
270}
271
272/// Compute ESP grid from atomic elements, positions and Mulliken charges.
273pub fn compute_esp(
274    elements: &[u8],
275    positions: &[[f64; 3]],
276    spacing: f64,
277    padding: f64,
278) -> Result<esp::EspGrid, String> {
279    let pop = compute_population(elements, positions)?;
280    Ok(esp::compute_esp_grid(
281        elements,
282        positions,
283        &pop.mulliken_charges,
284        spacing,
285        padding,
286    ))
287}
288
289/// Compute DOS/PDOS from EHT orbital energies.
290pub fn compute_dos(
291    elements: &[u8],
292    positions: &[[f64; 3]],
293    sigma: f64,
294    e_min: f64,
295    e_max: f64,
296    n_points: usize,
297) -> Result<dos::DosResult, String> {
298    let eht_result = eht::solve_eht(elements, positions, None)?;
299    let flat: Vec<f64> = positions.iter().flat_map(|p| p.iter().copied()).collect();
300    Ok(dos::compute_pdos(
301        elements,
302        &flat,
303        &eht_result.energies,
304        &eht_result.coefficients,
305        eht_result.n_electrons,
306        sigma,
307        e_min,
308        e_max,
309        n_points,
310    ))
311}
312
313/// Compute RMSD between two coordinate sets after Kabsch alignment.
314pub fn compute_rmsd(coords: &[f64], reference: &[f64]) -> f64 {
315    alignment::compute_rmsd(coords, reference)
316}
317
318/// Compute UFF force field energy for a molecule.
319pub fn compute_uff_energy(smiles: &str, coords: &[f64]) -> Result<f64, String> {
320    let mol = graph::Molecule::from_smiles(smiles)?;
321    let n = mol.graph.node_count();
322    if coords.len() != n * 3 {
323        return Err(format!("coords length {} != 3 * atoms {}", coords.len(), n));
324    }
325    let ff = forcefield::builder::build_uff_force_field(&mol);
326    let mut gradient = vec![0.0f64; n * 3];
327    let energy = ff.compute_system_energy_and_gradients(coords, &mut gradient);
328    Ok(energy)
329}
330
331/// Create a periodic unit cell from lattice parameters (a, b, c in Å; α, β, γ in degrees).
332pub fn create_unit_cell(
333    a: f64,
334    b: f64,
335    c: f64,
336    alpha: f64,
337    beta: f64,
338    gamma: f64,
339) -> materials::UnitCell {
340    materials::UnitCell::from_parameters(&materials::CellParameters {
341        a,
342        b,
343        c,
344        alpha,
345        beta,
346        gamma,
347    })
348}
349
350/// Assemble a framework crystal structure from node/linker SBUs on a topology.
351///
352/// Returns the assembled crystal structure as a JSON-serializable `CrystalStructure`.
353pub fn assemble_framework(
354    node: &materials::Sbu,
355    linker: &materials::Sbu,
356    topology: &materials::Topology,
357    cell: &materials::UnitCell,
358) -> materials::CrystalStructure {
359    materials::assemble_framework(node, linker, topology, cell)
360}
sci_form/lib.rs

sci_form/
lib.rs
