gchemol_core/
data.rs

1// imports
2
3//===============================================================================#
4//   DESCRIPTION: provide basic chemical data for atoms and molecules
5//
6//       OPTIONS:  ---
7//  REQUIREMENTS:  ---
8//         NOTES:  ---
9//        AUTHOR:  Wenping Guo <ybyygu@gmail.com>
10//       LICENCE:  GPL version 3
11//       CREATED:  <2018-04-12 Thu 14:40>
12//===============================================================================#
13
14use crate::{Atom, AtomKind, Molecule};
15
16// radii
17// Element radii data taking from: https://mendeleev.readthedocs.io/en/stable/data.html
18
19// Data in columns:
20// covalent_radii_single covalent_radii_double, covalent_radii_triple, vdw_radii
21
22
23const RADII_DATA: [[f64; 4]; 118] = [
24    [0.32, 0.32, 0.32, 1.1],
25    [0.46, 0.46, 0.46, 1.4],
26    [1.33, 1.24, 1.24, 1.82],
27    [1.02, 0.9, 0.85, 1.53],
28    [0.85, 0.78, 0.73, 1.92],
29    [0.75, 0.67, 0.6, 1.7],
30    [0.71, 0.6, 0.54, 1.55],
31    [0.63, 0.57, 0.53, 1.52],
32    [0.64, 0.59, 0.53, 1.47],
33    [0.67, 0.96, 0.96, 1.54],
34    [1.55, 1.6, 1.6, 2.27],
35    [1.39, 1.32, 1.27, 1.73],
36    [1.26, 1.13, 1.11, 1.84],
37    [1.16, 1.07, 1.02, 2.1],
38    [1.11, 1.02, 0.94, 1.8],
39    [1.03, 0.94, 0.95, 1.8],
40    [0.99, 0.95, 0.93, 1.75],
41    [0.96, 1.07, 0.96, 1.88],
42    [1.96, 1.93, 1.93, 2.75],
43    [1.71, 1.47, 1.33, 2.31],
44    [1.48, 1.16, 1.14, 2.15],
45    [1.36, 1.17, 1.08, 2.11],
46    [1.34, 1.12, 1.06, 2.07],
47    [1.22, 1.11, 1.03, 2.06],
48    [1.19, 1.05, 1.03, 2.05],
49    [1.16, 1.09, 1.02, 2.04],
50    [1.11, 1.03, 0.96, 2.0],
51    [1.1, 1.01, 1.01, 1.97],
52    [1.12, 1.15, 1.2, 1.96],
53    [1.18, 1.2, 1.2, 2.01],
54    [1.24, 1.17, 1.21, 1.87],
55    [1.21, 1.11, 1.14, 2.11],
56    [1.21, 1.14, 1.06, 1.85],
57    [1.16, 1.07, 1.07, 1.9],
58    [1.14, 1.09, 1.1, 1.85],
59    [1.17, 1.21, 1.08, 2.02],
60    [2.1, 2.02, 2.02, 3.03],
61    [1.85, 1.57, 1.39, 2.49],
62    [1.63, 1.3, 1.24, 2.32],
63    [1.54, 1.27, 1.21, 2.23],
64    [1.47, 1.25, 1.16, 2.18],
65    [1.38, 1.21, 1.13, 2.17],
66    [1.28, 1.2, 1.1, 2.16],
67    [1.25, 1.14, 1.03, 2.13],
68    [1.25, 1.1, 1.06, 2.1],
69    [1.2, 1.17, 1.12, 2.1],
70    [1.28, 1.39, 1.37, 2.11],
71    [1.36, 1.44, 1.44, 2.18],
72    [1.42, 1.36, 1.46, 1.93],
73    [1.4, 1.3, 1.32, 2.17],
74    [1.4, 1.33, 1.27, 2.06],
75    [1.36, 1.28, 1.21, 2.06],
76    [1.33, 1.29, 1.25, 1.98],
77    [1.31, 1.35, 1.22, 2.16],
78    [2.32, 2.09, 2.09, 3.43],
79    [1.96, 1.61, 1.49, 2.68],
80    [1.8, 1.39, 1.39, 2.43],
81    [1.63, 1.37, 1.31, 2.42],
82    [1.76, 1.38, 1.28, 2.4],
83    [1.74, 1.37, 1.37, 2.39],
84    [1.73, 1.35, 1.35, 2.38],
85    [1.72, 1.34, 1.34, 2.36],
86    [1.68, 1.34, 1.34, 2.35],
87    [1.69, 1.35, 1.32, 2.34],
88    [1.68, 1.35, 1.35, 2.33],
89    [1.67, 1.33, 1.33, 2.31],
90    [1.66, 1.33, 1.33, 2.3],
91    [1.65, 1.33, 1.33, 2.29],
92    [1.64, 1.31, 1.31, 2.27],
93    [1.7, 1.29, 1.29, 2.26],
94    [1.62, 1.31, 1.31, 2.24],
95    [1.52, 1.28, 1.22, 2.23],
96    [1.46, 1.26, 1.19, 2.22],
97    [1.37, 1.2, 1.15, 2.18],
98    [1.31, 1.19, 1.1, 2.16],
99    [1.29, 1.16, 1.09, 2.16],
100    [1.22, 1.15, 1.07, 2.13],
101    [1.23, 1.12, 1.1, 2.13],
102    [1.24, 1.21, 1.23, 2.14],
103    [1.33, 1.42, 1.42, 2.23],
104    [1.44, 1.42, 1.5, 1.96],
105    [1.44, 1.35, 1.37, 2.02],
106    [1.51, 1.41, 1.35, 2.07],
107    [1.45, 1.35, 1.29, 1.97],
108    [1.47, 1.38, 1.38, 2.02],
109    [1.42, 1.45, 1.33, 2.2],
110    [2.23, 2.18, 2.18, 3.48],
111    [2.01, 1.73, 1.59, 2.83],
112    [1.86, 1.53, 1.4, 2.47],
113    [1.75, 1.43, 1.36, 2.45],
114    [1.69, 1.38, 1.29, 2.43],
115    [1.7, 1.34, 1.18, 2.41],
116    [1.71, 1.36, 1.16, 2.39],
117    [1.72, 1.35, 1.35, 2.43],
118    [1.66, 1.35, 1.35, 2.44],
119    [1.66, 1.36, 1.36, 2.45],
120    [1.68, 1.39, 1.39, 2.44],
121    [1.68, 1.4, 1.4, 2.45],
122    [1.65, 1.4, 1.4, 2.45],
123    [1.67, 1.67, 1.67, 2.45],
124    [1.73, 1.39, 1.39, 2.46],
125    [1.76, 1.76, 1.76, 2.46],
126    [1.61, 1.41, 1.41, 2.46],
127    [1.57, 1.4, 1.31, 2.46],
128    [1.49, 1.36, 1.26, 2.46],
129    [1.43, 1.28, 1.21, 2.46],
130    [1.41, 1.28, 1.19, 2.46],
131    [1.34, 1.25, 1.18, 2.46],
132    [1.29, 1.25, 1.13, 2.46],
133    [1.28, 1.16, 1.18, 2.46],
134    [1.21, 1.16, 1.18, 2.46],
135    [1.22, 1.37, 1.3, 2.46],
136    [1.36, 1.36, 1.36, 2.46],
137    [1.43, 1.43, 1.43, 2.46],
138    [1.62, 1.62, 1.62, 2.46],
139    [1.75, 1.75, 1.75, 2.46],
140    [1.65, 1.65, 1.65, 2.46],
141    [1.57, 1.57, 1.57, 2.46],
142];
143
144// Following data are taken from jmol for auto bonding, which works better for
145// detecting chemical bonds, especially for metals.
146const BONDING_RADII: [f64; 109] = [
147    0.230, 0.930, 0.680, 0.350, 0.830, 0.680, 0.680, 0.680, 0.640, 1.120, 0.970, 1.100, 1.350, 1.200, 0.750, 1.020, 0.990, 1.570, 1.330,
148    0.990, 1.440, 1.470, 1.330, 1.350, 1.350, 1.340, 1.330, 1.500, 1.520, 1.450, 1.220, 1.170, 1.210, 1.220, 1.210, 1.910, 1.470, 1.120,
149    1.780, 1.560, 1.480, 1.470, 1.350, 1.400, 1.450, 1.500, 1.590, 1.690, 1.630, 1.460, 1.460, 1.470, 1.400, 1.980, 1.670, 1.340, 1.870,
150    1.830, 1.820, 1.810, 1.800, 1.800, 1.990, 1.790, 1.760, 1.750, 1.740, 1.730, 1.720, 1.940, 1.720, 1.570, 1.430, 1.370, 1.350, 1.370,
151    1.320, 1.500, 1.500, 1.700, 1.550, 1.540, 1.540, 1.680, 1.700, 2.400, 2.000, 1.900, 1.880, 1.790, 1.610, 1.580, 1.550, 1.530, 1.510,
152    1.500, 1.500, 1.500, 1.500, 1.500, 1.500, 1.500, 1.500, 1.600, 1.600, 1.600, 1.600, 1.600, 1.600,
153];
154
155/// Return covalent radius for single, double, or triple bonds
156fn get_cov_radius(element_number: usize, bond_order: usize) -> Option<f64> {
157    if element_number > 0 && element_number <= RADII_DATA.len() {
158        // only for single, double, or triple bond
159        if bond_order > 0 && bond_order <= 3 {
160            let irow = element_number - 1;
161            let icol = bond_order - 1;
162
163            let r = RADII_DATA[irow][icol];
164            return Some(r);
165        }
166    }
167
168    None
169}
170
171/// Return Van der Waals radius if any
172fn get_vdw_radius(element_number: usize) -> Option<f64> {
173    // column index to vdw radii
174    let icol = 3;
175    if element_number > 0 && element_number <= RADII_DATA.len() {
176        let irow = element_number - 1;
177        let r = RADII_DATA[irow][icol];
178        return Some(r);
179    }
180    None
181}
182
183// Return a radius for auto bond
184fn get_bonding_radius(element_number: usize) -> Option<f64> {
185    BONDING_RADII.get(element_number - 1).copied()
186}
187
188// masses
189// # get from python mendeleev package
190
191const MASSES_DATA: [f64; 118] = [
192    1.008000, 4.002602, 6.940000, 9.012183, 10.810000, 12.011000, 14.007000, 15.999000, 18.998403, 20.179700,
193    22.989769, 24.305000, 26.981538, 28.085000, 30.973762, 32.060000, 35.450000, 39.948000, 39.098300, 40.078000,
194    44.955908, 47.867000, 50.941500, 51.996100, 54.938044, 55.845000, 58.933194, 58.693400, 63.546000, 65.380000,
195    69.723000, 72.630000, 74.921595, 78.971000, 79.904000, 83.798000, 85.467800, 87.620000, 88.905840, 91.224000,
196    92.906370, 95.950000, 97.907210, 101.070000, 102.905500, 106.420000, 107.868200, 112.414000, 114.818000,
197    118.710000, 121.760000, 127.600000, 126.904470, 131.293000, 132.905452, 137.327000, 138.905470, 140.116000,
198    140.907660, 144.242000, 144.912760, 150.360000, 151.964000, 157.250000, 158.925350, 162.500000, 164.930330,
199    167.259000, 168.934220, 173.045000, 174.966800, 178.490000, 180.947880, 183.840000, 186.207000, 190.230000,
200    192.217000, 195.084000, 196.966569, 200.592000, 204.380000, 207.200000, 208.980400, 209.000000, 210.000000,
201    222.000000, 223.000000, 226.000000, 227.000000, 232.037700, 231.035880, 238.028910, 237.000000, 244.000000,
202    243.000000, 247.000000, 247.000000, 251.000000, 252.000000, 257.000000, 258.000000, 259.000000, 262.000000,
203    267.000000, 268.000000, 271.000000, 274.000000, 269.000000, 276.000000, 281.000000, 281.000000, 285.000000,
204    286.000000, 289.000000, 288.000000, 293.000000, 294.000000, 294.000000,
205];
206
207fn get_atom_mass(atom: &Atom) -> Option<f64> {
208    match atom.kind() {
209        AtomKind::Element(n) => {
210            assert!(*n > 0, "invalid element number");
211            Some(MASSES_DATA[n - 1])
212        }
213        AtomKind::Dummy(_) => None,
214    }
215}
216
217// atom
218// #+name: a8beb3dd
219
220/// Core data for `Atom`
221impl Atom {
222    /// Access covalent radius of atom. Return None if no data available or atom is dummy.
223    pub fn get_cov_radius(&self) -> Option<f64> {
224        get_cov_radius(self.number(), 1)
225    }
226
227    /// Access Van der Waals radius of atom.
228    /// Return None if no data available
229    pub fn get_vdw_radius(&self) -> Option<f64> {
230        get_vdw_radius(self.number())
231    }
232
233    /// Return bonding radius of this element if available
234    pub fn get_bonding_radius(&self) -> Option<f64> {
235        get_bonding_radius(self.number())
236    }
237
238    #[deprecated(note = "use get_cov_radius instead")]
239    /// Access covalent radius of atom. Return None if no data available or atom is dummy.
240    pub fn cov_radius(&self) -> Option<f64> {
241        self.get_cov_radius()
242    }
243
244    /// Get mass in atomic mass unit. Return None if atom is dummy.
245    pub fn get_mass(&self) -> Option<f64> {
246        self.mass.or(get_atom_mass(self))
247    }
248
249    #[deprecated(note = "Use get_vdw_radius instead")]
250    /// Access Van der Waals radius of atom.
251    /// Return None if no data available
252    pub fn vdw_radius(&self) -> Option<f64> {
253        self.get_vdw_radius()
254    }
255}
256
257// molecule
258// #+name: 190e0574
259
260impl Molecule {
261    /// Returns `Molecule` created from the internal database (mainly for tests).
262    ///
263    /// Available names: `CH4`, `H2O`, `HCN`
264    pub fn from_database(name: &str) -> Self {
265        let ch4 = "
266C  -0.0000   -0.0000    0.0000
267H   1.0900   -0.0000    0.0000
268H  -0.3633    1.0277    0.0000
269H  -0.3633   -0.5138    0.8900
270H  -0.3633   -0.5138   -0.8900 ";
271
272        let h2o = "
273O    -1.4689     2.1375     0.0000
274H    -0.5089     2.1375     0.0000
275H    -1.7894     3.0424     0.0000
276";
277
278        let hcn = "
279H		-2.5671751	1.2900188	0.0000000
280C		-1.4971751	1.2900188	0.0000000
281N		-0.3505751	1.2900188	0.0000000
282";
283
284        match name {
285            "HCN" => {
286                let atoms = hcn.trim().lines().filter_map(|line| line.parse::<Atom>().ok());
287                Molecule::from_atoms(atoms)
288            }
289            "H2O" => {
290                let atoms = h2o.trim().lines().filter_map(|line| line.parse::<Atom>().ok());
291                Molecule::from_atoms(atoms)
292            }
293            "CH4" => {
294                let atoms = ch4.trim().lines().filter_map(|line| line.parse::<Atom>().ok());
295                Molecule::from_atoms(atoms)
296            }
297            _ => unimplemented!(),
298        }
299    }
300}
301
302// test
303
304#[test]
305fn test_atom_data() {
306    // carbon atom
307    let atom1 = Atom::new("C", [-0.90203687, 0.62555259, 0.0081889]);
308    // dummy atom
309    let atom2 = Atom::new("X", [-0.54538244, -0.38325741, 0.0081889]);
310
311    assert!(atom1.get_cov_radius().is_some());
312    assert!(atom1.get_vdw_radius().is_some());
313    assert!(atom2.get_cov_radius().is_none());
314
315    // atom mass
316    assert_eq!(atom1.get_mass(), Some(12.011));
317    assert_eq!(atom2.get_mass(), None);
318}
gchemol_core/data.rs

gchemol_core/
data.rs
