omgkit_match/
props.rs

1//! 把一个净化过的分子预计算成**逐原子逐键的查询性质**。
2//!
3//! # 为什么要预计算
4//!
5//! 匹配是回溯搜索:同一个原子会被反复拿来与不同的查询原子比对,一个 20 原子
6//! 的模式配到 40 原子的分子上,单个目标原子被求值几十上百次很常见。
7//!
8//! 环成员数、最小环大小这些量要遍历环集才能算出来。放在求值里现算,等于把
9//! "遍历环集"塞进了搜索的最内层循环 —— 无论搜索本身多快都救不回来。
10//!
11//! # 环的三个量来自两处
12//!
13//! `R`(所属环个数)与 `x`(环键数)要数**具体的环**,只能从环集来;
14//! `r`(最小环大小)环感知那一步已经算好了,直接取。
15//!
16//! 这两处必须用**同一个环集**,否则会出现"`R0` 但 `r6`"这种自相矛盾的性质。
17//! 所以本模块自己跑一遍环感知和环集,不接受调用方分别传进来。
18//!
19//! # 氢数要把**图里的氢原子**也算上
20//!
21//! SMARTS 的 `H<n>` 数的是"这个原子上一共有几个氢",不管氢是记成计数还是
22//! 画成图里的一个节点。`[H]N(C)C` 的氮:显式氢 0、隐式氢 0,但 `[N;H1]`
23//! 命中它 —— 那个氢是它的**邻居**。
24//!
25//! 漏掉这一项的失效方式很安静:绝大多数分子不写显式氢,测试全绿,直到遇上
26//! 一条 `[H]N(...)` 才莫名其妙地不命中。
27
28use omgkit_chem::{perceive_rings, ring_set};
29use omgkit_core::{AtomFlags, BondOrder, MolBuilder};
30use omgkit_io::smarts::{AtomProps, BondProps};
31
32/// 一个分子的全部查询性质,按原子/键下标索引。
33#[derive(Debug, Clone)]
34pub struct MolProps {
35    /// 逐原子
36    pub atoms: Vec<AtomProps>,
37    /// 逐键。配位键的朝向留给匹配时按方向决定,这里 `dative_forward` 恒为真。
38    pub bonds: Vec<BondProps>,
39}
40
41impl MolProps {
42    /// 从一个**已净化**的分子预计算。
43    ///
44    /// 分子必须已经跑过价键计算与芳香性感知 —— 隐式氢、芳香标志都直接取
45    /// 分子上的字段。没净化的分子算出来的性质是错的,而且错得很安静:
46    /// 隐式氢全是 0,芳香标志是"作者声称"而非感知结果。
47    #[must_use]
48    pub fn compute(mol: &MolBuilder) -> Self {
49        // 环感知与环集必须来自同一次计算,见模块文档
50        let mut scratch = mol.clone();
51        let rings = perceive_rings(&mut scratch);
52        let cycles = ring_set(&scratch);
53
54        let n = mol.num_atoms();
55        let mut ring_count = vec![0u32; n];
56        let mut ring_bond_count = vec![0u32; n];
57        for ring in &cycles {
58            for &a in &ring.atoms {
59                ring_count[a as usize] += 1;
60            }
61        }
62        // 环键数要按**键**去重:同一条键属于多个环时,对端点只算一次
63        let mut bond_is_ring = vec![false; mol.num_bonds()];
64        for ring in &cycles {
65            for &b in &ring.bonds {
66                bond_is_ring[b as usize] = true;
67            }
68        }
69        for (bi, &in_ring) in bond_is_ring.iter().enumerate() {
70            if in_ring {
71                let b = mol.bonds()[bi];
72                ring_bond_count[b.begin as usize] += 1;
73                ring_bond_count[b.end as usize] += 1;
74            }
75        }
76
77        let atoms = (0..n)
78            .map(|i| {
79                let a = mol.atoms()[i];
80                AtomProps {
81                    atomic_num: a.atomic_num,
82                    aromatic: a.flags.contains(AtomFlags::AROMATIC),
83                    charge: i32::from(a.formal_charge),
84                    isotope: a.isotope,
85                    degree: mol.degree(i as u32) as u32,
86                    total_hs: u32::from(a.num_explicit_hs)
87                        + u32::from(a.num_implicit_hs)
88                        + neighbour_hydrogens(mol, i as u32),
89                    implicit_hs: u32::from(a.num_implicit_hs),
90                    // **总价只有一处实现。** 先前这里自己写了一份"键级和四舍五入
91                    // 加氢数",少了 core 那份的芳香价回落:那一步只对**两根**芳香键
92                    // 成立,稠合位有三根,4.5 进位成 5,于是萘的两个稠合碳被判成
93                    // 5 价 —— 任何用 `[v4]` 挑碳的 SMARTS 在稠环上都漏掉稠合位。
94                    valence: omgkit_core::valence::total_valence_nonstrict(mol, i as u32)
95                        .max(0)
96                        .unsigned_abs(),
97                    ring_count: ring_count[i],
98                    min_ring_size: u32::from(rings.atom_min_ring_size[i]),
99                    ring_bonds: ring_bond_count[i],
100                    chiral_tag: a.chiral_tag,
101                    atom_map: a.atom_map,
102                }
103            })
104            .collect();
105
106        let bonds = (0..mol.num_bonds())
107            .map(|i| {
108                let b = mol.bonds()[i];
109                BondProps {
110                    order: b.order,
111                    in_ring: rings.bond_in_ring[i],
112                    direction: b.direction,
113                    dative_forward: true,
114                }
115            })
116            .collect();
117
118        Self { atoms, bonds }
119    }
120}
121
122/// 邻居里画成独立节点的氢原子数。
123///
124/// 只数**真的氢**:同位素(氘、氚)也算,通配原子不算。
125fn neighbour_hydrogens(mol: &MolBuilder, atom: u32) -> u32 {
126    mol.neighbors(atom)
127        .filter(|&(other, _)| mol.atoms()[other as usize].atomic_num == 1)
128        .count() as u32
129}
130
131/// 该原子上的键在 `BondOrder::Dative` 时,给体是不是 `from` 端。
132#[must_use]
133pub fn dative_points_from(mol: &MolBuilder, bond: u32, from: u32) -> bool {
134    let b = mol.bonds()[bond as usize];
135    b.order != BondOrder::Dative || b.begin == from
136}
omgkit_match/props.rs

omgkit_match/
props.rs
