sdfrust 0.6.0

A fast, pure-Rust parser for SDF, MOL2, and XYZ chemical structure files
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161

//! RMSD (Root Mean Square Deviation) calculation.

use crate::{Molecule, SdfError};

/// Calculate RMSD between two sets of coordinates.
///
/// The coordinates are provided as slices of (x, y, z) tuples.
/// Both coordinate sets must have the same length.
///
/// # Arguments
///
/// * `coords1` - First set of coordinates
/// * `coords2` - Second set of coordinates
///
/// # Returns
///
/// The RMSD value, or an error if the coordinate sets have different lengths.
///
/// # Example
///
/// ```rust
/// use sdfrust::geometry::rmsd_from_coords;
///
/// let coords1 = vec![(0.0, 0.0, 0.0), (1.0, 0.0, 0.0)];
/// let coords2 = vec![(0.0, 0.0, 0.0), (1.0, 0.0, 0.0)];
/// let rmsd = rmsd_from_coords(&coords1, &coords2).unwrap();
/// assert!(rmsd.abs() < 1e-10);
/// ```
pub fn rmsd_from_coords(
    coords1: &[(f64, f64, f64)],
    coords2: &[(f64, f64, f64)],
) -> Result<f64, SdfError> {
    if coords1.len() != coords2.len() {
        return Err(SdfError::AtomCountMismatch {
            expected: coords1.len(),
            found: coords2.len(),
        });
    }

    if coords1.is_empty() {
        return Ok(0.0);
    }

    let n = coords1.len() as f64;
    let sum_sq: f64 = coords1
        .iter()
        .zip(coords2.iter())
        .map(|((x1, y1, z1), (x2, y2, z2))| {
            let dx = x1 - x2;
            let dy = y1 - y2;
            let dz = z1 - z2;
            dx * dx + dy * dy + dz * dz
        })
        .sum();

    Ok((sum_sq / n).sqrt())
}

/// Calculate RMSD between two molecules.
///
/// Computes the root mean square deviation of atomic positions
/// between two molecules. The molecules must have the same number of atoms.
/// No alignment is performed - atoms are compared directly by index.
///
/// # Arguments
///
/// * `mol1` - First molecule
/// * `mol2` - Second molecule
///
/// # Returns
///
/// The RMSD value in Angstroms, or an error if the molecules have different atom counts.
///
/// # Example
///
/// ```rust
/// use sdfrust::{Molecule, Atom};
/// use sdfrust::geometry::rmsd_to;
///
/// let mut mol1 = Molecule::new("mol1");
/// mol1.atoms.push(Atom::new(0, "C", 0.0, 0.0, 0.0));
/// mol1.atoms.push(Atom::new(1, "C", 1.0, 0.0, 0.0));
///
/// let mut mol2 = Molecule::new("mol2");
/// mol2.atoms.push(Atom::new(0, "C", 0.0, 0.0, 0.0));
/// mol2.atoms.push(Atom::new(1, "C", 1.0, 0.0, 0.0));
///
/// let rmsd = rmsd_to(&mol1, &mol2).unwrap();
/// assert!(rmsd.abs() < 1e-10);
/// ```
pub fn rmsd_to(mol1: &Molecule, mol2: &Molecule) -> Result<f64, SdfError> {
    let coords1: Vec<(f64, f64, f64)> = mol1.atoms.iter().map(|a| (a.x, a.y, a.z)).collect();
    let coords2: Vec<(f64, f64, f64)> = mol2.atoms.iter().map(|a| (a.x, a.y, a.z)).collect();
    rmsd_from_coords(&coords1, &coords2)
}

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

    #[test]
    fn test_rmsd_identical_molecules() {
        let mut mol1 = Molecule::new("mol1");
        mol1.atoms.push(Atom::new(0, "C", 0.0, 0.0, 0.0));
        mol1.atoms.push(Atom::new(1, "C", 1.0, 0.0, 0.0));
        mol1.atoms.push(Atom::new(2, "C", 0.0, 1.0, 0.0));

        let mol2 = mol1.clone();

        let rmsd = rmsd_to(&mol1, &mol2).unwrap();
        assert!(rmsd.abs() < 1e-10);
    }

    #[test]
    fn test_rmsd_translated_molecule() {
        let mut mol1 = Molecule::new("mol1");
        mol1.atoms.push(Atom::new(0, "C", 0.0, 0.0, 0.0));
        mol1.atoms.push(Atom::new(1, "C", 1.0, 0.0, 0.0));

        let mut mol2 = Molecule::new("mol2");
        mol2.atoms.push(Atom::new(0, "C", 1.0, 0.0, 0.0)); // translated by (1, 0, 0)
        mol2.atoms.push(Atom::new(1, "C", 2.0, 0.0, 0.0));

        let rmsd = rmsd_to(&mol1, &mol2).unwrap();
        // Each atom is displaced by 1.0, so RMSD = sqrt(1^2 + 1^2) / sqrt(2) = 1.0
        assert!((rmsd - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_rmsd_different_atom_counts() {
        let mut mol1 = Molecule::new("mol1");
        mol1.atoms.push(Atom::new(0, "C", 0.0, 0.0, 0.0));

        let mut mol2 = Molecule::new("mol2");
        mol2.atoms.push(Atom::new(0, "C", 0.0, 0.0, 0.0));
        mol2.atoms.push(Atom::new(1, "C", 1.0, 0.0, 0.0));

        let result = rmsd_to(&mol1, &mol2);
        assert!(result.is_err());
    }

    #[test]
    fn test_rmsd_empty_molecules() {
        let mol1 = Molecule::new("mol1");
        let mol2 = Molecule::new("mol2");

        let rmsd = rmsd_to(&mol1, &mol2).unwrap();
        assert!(rmsd.abs() < 1e-10);
    }

    #[test]
    fn test_rmsd_from_coords() {
        let coords1 = vec![(0.0, 0.0, 0.0), (1.0, 0.0, 0.0)];
        let coords2 = vec![(0.5, 0.0, 0.0), (1.5, 0.0, 0.0)];

        let rmsd = rmsd_from_coords(&coords1, &coords2).unwrap();
        // Each atom displaced by 0.5
        assert!((rmsd - 0.5).abs() < 1e-10);
    }
}