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
237    #[cfg(feature = "parallel")]
238    {
239        Ok(surface::sasa::compute_sasa_parallel(
240            elements,
241            &positions,
242            probe_radius,
243            None,
244        ))
245    }
246
247    #[cfg(not(feature = "parallel"))]
248    {
249        Ok(surface::sasa::compute_sasa(
250            elements,
251            &positions,
252            probe_radius,
253            None,
254        ))
255    }
256}
257
258/// Compute Mulliken & Löwdin population analysis from atomic elements and positions.
259pub fn compute_population(
260    elements: &[u8],
261    positions: &[[f64; 3]],
262) -> Result<population::PopulationResult, String> {
263    let eht_result = eht::solve_eht(elements, positions, None)?;
264    Ok(population::compute_population(
265        elements,
266        positions,
267        &eht_result.coefficients,
268        eht_result.n_electrons,
269    ))
270}
271
272/// Compute molecular dipole moment (Debye) from atomic elements and positions.
273pub fn compute_dipole(
274    elements: &[u8],
275    positions: &[[f64; 3]],
276) -> Result<dipole::DipoleResult, String> {
277    let eht_result = eht::solve_eht(elements, positions, None)?;
278    Ok(dipole::compute_dipole_from_eht(
279        elements,
280        positions,
281        &eht_result.coefficients,
282        eht_result.n_electrons,
283    ))
284}
285
286/// Compute ESP grid from atomic elements, positions and Mulliken charges.
287pub fn compute_esp(
288    elements: &[u8],
289    positions: &[[f64; 3]],
290    spacing: f64,
291    padding: f64,
292) -> Result<esp::EspGrid, String> {
293    let pop = compute_population(elements, positions)?;
294    #[cfg(feature = "parallel")]
295    {
296        Ok(esp::compute_esp_grid_parallel(
297            elements,
298            positions,
299            &pop.mulliken_charges,
300            spacing,
301            padding,
302        ))
303    }
304
305    #[cfg(not(feature = "parallel"))]
306    {
307        Ok(esp::compute_esp_grid(
308            elements,
309            positions,
310            &pop.mulliken_charges,
311            spacing,
312            padding,
313        ))
314    }
315}
316
317/// Compute DOS/PDOS from EHT orbital energies.
318pub fn compute_dos(
319    elements: &[u8],
320    positions: &[[f64; 3]],
321    sigma: f64,
322    e_min: f64,
323    e_max: f64,
324    n_points: usize,
325) -> Result<dos::DosResult, String> {
326    let eht_result = eht::solve_eht(elements, positions, None)?;
327    let flat: Vec<f64> = positions.iter().flat_map(|p| p.iter().copied()).collect();
328
329    #[cfg(feature = "parallel")]
330    {
331        Ok(dos::compute_pdos_parallel(
332            elements,
333            &flat,
334            &eht_result.energies,
335            &eht_result.coefficients,
336            eht_result.n_electrons,
337            sigma,
338            e_min,
339            e_max,
340            n_points,
341        ))
342    }
343
344    #[cfg(not(feature = "parallel"))]
345    {
346        Ok(dos::compute_pdos(
347            elements,
348            &flat,
349            &eht_result.energies,
350            &eht_result.coefficients,
351            eht_result.n_electrons,
352            sigma,
353            e_min,
354            e_max,
355            n_points,
356        ))
357    }
358}
359
360/// Compute RMSD between two coordinate sets after Kabsch alignment.
361pub fn compute_rmsd(coords: &[f64], reference: &[f64]) -> f64 {
362    alignment::compute_rmsd(coords, reference)
363}
364
365/// Compute UFF force field energy for a molecule.
366pub fn compute_uff_energy(smiles: &str, coords: &[f64]) -> Result<f64, String> {
367    let mol = graph::Molecule::from_smiles(smiles)?;
368    let n = mol.graph.node_count();
369    if coords.len() != n * 3 {
370        return Err(format!("coords length {} != 3 * atoms {}", coords.len(), n));
371    }
372    let ff = forcefield::builder::build_uff_force_field(&mol);
373    let mut gradient = vec![0.0f64; n * 3];
374    let energy = ff.compute_system_energy_and_gradients(coords, &mut gradient);
375    Ok(energy)
376}
377
378/// Create a periodic unit cell from lattice parameters (a, b, c in Å; α, β, γ in degrees).
379pub fn create_unit_cell(
380    a: f64,
381    b: f64,
382    c: f64,
383    alpha: f64,
384    beta: f64,
385    gamma: f64,
386) -> materials::UnitCell {
387    materials::UnitCell::from_parameters(&materials::CellParameters {
388        a,
389        b,
390        c,
391        alpha,
392        beta,
393        gamma,
394    })
395}
396
397/// Assemble a framework crystal structure from node/linker SBUs on a topology.
398///
399/// Returns the assembled crystal structure as a JSON-serializable `CrystalStructure`.
400pub fn assemble_framework(
401    node: &materials::Sbu,
402    linker: &materials::Sbu,
403    topology: &materials::Topology,
404    cell: &materials::UnitCell,
405) -> materials::CrystalStructure {
406    materials::assemble_framework(node, linker, topology, cell)
407}
sci_form/lib.rs

sci_form/
lib.rs
