use crate::params::{self, Source};
use crate::smooth::Bounds;
use omgkit_core::MolBuilder;
#[must_use]
pub fn vdw_frac(topo_dist: usize) -> f64 {
match topo_dist {
0..=4 => 0.70,
5 => 0.85,
_ => 1.00,
}
}
pub const MAX_UPPER: f64 = 1000.0;
pub const DIST12_TOL: f64 = 0.01;
pub const DIST13_TOL: f64 = 0.04;
pub const DIST14_TOL: f64 = 0.06;
pub const DIST15_TOL: f64 = 0.08;
pub const BOND_REL: f64 = 0.012;
pub const ANGLE_TOL: f64 = 2.5;
#[must_use]
pub fn ring_internal_torsion(size: usize, aromatic: bool, sp3: bool) -> (f64, f64) {
type Row = (f64, f64, f64);
static T: std::sync::OnceLock<std::collections::HashMap<(usize, bool, bool), Row>> =
std::sync::OnceLock::new();
let t = T.get_or_init(|| {
let mut m = std::collections::HashMap::new();
for line in include_str!("../data/mmff.ringtorsion.tsv").lines() {
if line.starts_with('#') {
continue;
}
let f: Vec<&str> = line.split('\t').collect();
if f.len() < 7 {
continue;
}
let (Ok(sz), Ok(ar), Ok(sp), Ok(p05), Ok(p95)) = (
f[0].parse::<usize>(),
f[1].parse::<u8>(),
f[2].parse::<u8>(),
f[5].parse::<f64>(),
f[6].parse::<f64>(),
) else {
continue;
};
let Ok(med) = f[4].parse::<f64>() else {
continue;
};
m.insert((sz, ar == 1, sp == 1), (med, p05, p95));
}
m
});
t.get(&(size, aromatic, sp3))
.map_or((0.0, 180.0), |&(med, p05, p95)| {
if aromatic {
(med, med)
} else {
(p05, p95)
}
})
}
fn stereo_path_torsion(mol: &MolBuilder, bidx: usize, i: u32, k: u32, j: u32) -> Option<f64> {
use omgkit_core::BondStereo;
let bd = mol.bonds().get(bidx)?;
let same_side = match bd.stereo {
BondStereo::Z | BondStereo::Cis => true,
BondStereo::E | BondStereo::Trans => false,
BondStereo::None => return None,
};
let (ra, rb) = (bd.stereo_atoms[0], bd.stereo_atoms[1]);
if ra == omgkit_core::BondData::NO_STEREO_ATOM || rb == omgkit_core::BondData::NO_STEREO_ATOM {
return None;
}
let (ref_i, ref_j) = if k == bd.begin { (ra, rb) } else { (rb, ra) };
let flips = usize::from(i != ref_i) + usize::from(j != ref_j);
let cis = if flips % 2 == 0 {
same_side
} else {
!same_side
};
Some(if cis { 0.0 } else { 180.0 })
}
struct Rings<'a> {
sets: &'a [Vec<u32>],
aromatic: &'a [bool],
sp3: &'a [bool],
}
fn ring_path_torsion(
mol: &MolBuilder,
rings: &Rings<'_>,
i: u32,
k: u32,
l: u32,
j: u32,
) -> Option<(f64, f64)> {
let (ring_sets, ring_aromatic, ring_sp3) = (rings.sets, rings.aromatic, rings.sp3);
let mut best: Option<(usize, usize)> = None; for (r, set) in ring_sets.iter().enumerate() {
let has = |a: u32| set.binary_search(&a).is_ok();
if has(k) && has(l) && best.map_or(true, |(sz, _)| set.len() < sz) {
best = Some((set.len(), r));
}
}
let (size, r) = best?;
let sp2 = |a: u32| mol.atoms()[a as usize].hybridization == omgkit_core::Hybridization::Sp2;
let conj_planar = size <= 8 && sp2(k) && sp2(l);
if !ring_aromatic[r] && !ring_sp3[r] && !conj_planar {
return None;
}
let set = &ring_sets[r];
let has = |a: u32| set.binary_search(&a).is_ok();
let sp2_or_in_ring = |end: u32, attach: u32| has(end) || sp2(attach);
if !sp2_or_in_ring(i, k) || !sp2_or_in_ring(j, l) {
return None;
}
let (t_lo, t_hi) = if conj_planar && !ring_aromatic[r] {
(0.0, 0.0)
} else {
ring_internal_torsion(size, ring_aromatic[r], ring_sp3[r])
};
Some(if has(i) == has(j) {
(t_lo, t_hi)
} else {
(180.0 - t_hi, 180.0 - t_lo)
})
}
#[must_use]
pub fn torsion_envelope(shared_ring: usize, aromatic: bool) -> (f64, f64) {
if shared_ring == 0 {
return (0.0, 180.0);
}
if aromatic {
return (0.0, 0.0);
}
match shared_ring {
3 | 4 => (0.0, 0.0),
5 => (0.0, 40.0),
6 => (0.0, 60.0),
_ => (0.0, 90.0),
}
}
#[derive(Debug, Clone, Copy, Default, PartialEq, Eq)]
pub struct Stats {
pub atoms: usize,
pub n12: usize,
pub n13: usize,
pub n14: usize,
pub n13_conflict: usize,
pub n14_degenerate: usize,
pub bond_table: usize,
pub bond_relaxed: usize,
pub bond_model: usize,
pub angle_table: usize,
pub angle_relaxed: usize,
pub angle_model: usize,
pub angle_envelope: usize,
}
#[must_use]
pub fn coord_angle_envelope(degree: usize) -> (f64, f64) {
match degree {
0 | 1 => (180.0, 180.0),
2 => (180.0, 180.0),
3 => (120.0, 120.0),
4 => (109.4712, 109.4712),
5 | 6 => (90.0, 180.0),
7 => (72.0, 180.0),
_ => (70.0, 180.0),
}
}
fn third_side(a: f64, b: f64, theta_deg: f64) -> f64 {
let t = theta_deg.to_radians();
(a * a + b * b - 2.0 * a * b * t.cos()).max(0.0).sqrt()
}
fn bond_range(mol: &MolBuilder, i: u32, j: u32, min_ring: usize, st: &mut Stats) -> (f64, f64) {
let ord = mol
.neighbors(i)
.find(|(y, _)| *y == j)
.map_or(omgkit_core::BondOrder::Single, |(_, bi)| {
mol.bonds()[bi as usize].order
});
let p = params::bond_length(
mol.atoms()[i as usize].atomic_num,
mol.atoms()[j as usize].atomic_num,
ord,
min_ring,
);
match p.source {
Source::Table => st.bond_table += 1,
Source::RingRelaxed => st.bond_relaxed += 1,
Source::Model => st.bond_model += 1,
}
if p.source == Source::Model {
(p.value * (1.0 - BOND_REL), p.value * (1.0 + BOND_REL))
} else {
(p.value - DIST12_TOL, p.value + DIST12_TOL)
}
}
fn angle_range(
mol: &MolBuilder,
c: u32,
ring_self: usize,
ring_shared: usize,
st: &mut Stats,
) -> (f64, f64) {
let deg = mol.neighbors(c).count();
if deg >= 5 {
st.angle_envelope += 1;
return coord_angle_envelope(deg);
}
let p = params::angle(
mol.atoms()[c as usize].atomic_num,
deg,
mol.atoms()[c as usize]
.flags
.contains(omgkit_core::AtomFlags::AROMATIC),
ring_self,
ring_shared,
);
match p.source {
Source::Table => st.angle_table += 1,
Source::RingRelaxed => st.angle_relaxed += 1,
Source::Model => st.angle_model += 1,
}
if p.source == Source::Model {
(p.lo.to_degrees(), p.hi.to_degrees())
} else {
let v = p.value.to_degrees();
((v - ANGLE_TOL).max(1.0), (v + ANGLE_TOL).min(180.0))
}
}
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum Kind {
None,
B12,
B13,
B14,
}
fn tighten(b: &mut Bounds, i: usize, j: usize, lo: f64, hi: f64) -> bool {
let (clo, chi) = (b.lower(i, j), b.upper(i, j));
let (nlo, nhi) = (clo.max(lo), chi.min(hi));
if nlo <= nhi {
b.set_lower(i, j, nlo);
b.set_upper(i, j, nhi);
return false;
}
b.set_lower(i, j, clo.min(lo));
b.set_upper(i, j, chi.max(hi));
true
}
fn ring_info(mol: &MolBuilder) -> (Vec<usize>, Vec<Vec<u32>>) {
let n = mol.num_atoms();
let rings = omgkit_chem::sssr::ring_set(mol);
let mut min_ring = vec![0usize; n];
let mut sets: Vec<Vec<u32>> = Vec::with_capacity(rings.len());
for r in &rings {
let sz = r.atoms.len();
for a in &r.atoms {
let a = *a as usize;
if a < n && (min_ring[a] == 0 || sz < min_ring[a]) {
min_ring[a] = sz;
}
}
let mut v = r.atoms.clone();
v.sort_unstable();
sets.push(v);
}
(min_ring, sets)
}
fn shared_ring(sets: &[Vec<u32>], atoms: &[u32]) -> usize {
let mut best = 0usize;
for s in sets {
if atoms.iter().all(|a| s.binary_search(a).is_ok()) && (best == 0 || s.len() < best) {
best = s.len();
}
}
best
}
#[must_use]
pub fn build(mol: &MolBuilder) -> (Bounds, Stats) {
let n = mol.num_atoms();
let mut st = Stats {
atoms: n,
..Stats::default()
};
let mut b = Bounds::new(n, 0.0, MAX_UPPER);
if n == 0 {
return (b, st);
}
let (min_ring, ring_sets) = ring_info(mol);
let aromatic_ring: Vec<bool> = ring_sets
.iter()
.map(|set| {
set.iter().all(|a| {
mol.atoms()[*a as usize]
.flags
.contains(omgkit_core::AtomFlags::AROMATIC)
})
})
.collect();
let sp3_ring: Vec<bool> = ring_sets
.iter()
.map(|set| {
set.iter()
.all(|a| mol.atoms()[*a as usize].hybridization == omgkit_core::Hybridization::Sp3)
})
.collect();
let mut kind = vec![Kind::None; n * n];
for bd in mol.bonds() {
let (i, j) = (bd.begin as usize, bd.end as usize);
let r = shared_ring(&ring_sets, &[bd.begin, bd.end]);
let (lo, hi) = bond_range(mol, bd.begin, bd.end, r, &mut st);
tighten(&mut b, i, j, lo, hi);
kind[i * n + j] = Kind::B12;
kind[j * n + i] = Kind::B12;
st.n12 += 1;
}
for (c, &self_ring) in min_ring.iter().enumerate().take(n) {
let Ok(cu) = u32::try_from(c) else { continue };
let nb: Vec<u32> = mol.neighbors(cu).map(|(y, _)| y).collect();
if nb.len() < 2 {
continue;
}
for x in 0..nb.len() {
for y in (x + 1)..nb.len() {
let (i, j) = (nb[x] as usize, nb[y] as usize);
if i == j {
continue;
}
let shared = shared_ring(&ring_sets, &[nb[x], cu, nb[y]]);
let (ang_lo, ang_hi) = angle_range(mol, cu, self_ring, shared, &mut st);
let b1 = ((b.lower(i, c)) + (b.upper(i, c))) / 2.0;
let b2 = ((b.lower(c, j)) + (b.upper(c, j))) / 2.0;
let (lo, hi) = if ang_hi - ang_lo > 2.0 * ANGLE_TOL + 1e-9 {
(
third_side(b1, b2, ang_lo) - DIST13_TOL,
third_side(b1, b2, ang_hi) + DIST13_TOL,
)
} else {
let d = third_side(b1, b2, ((ang_lo) + (ang_hi)) / 2.0);
(d - DIST13_TOL, d + DIST13_TOL)
};
if kind[i * n + j] == Kind::B12 {
continue;
}
if tighten(&mut b, i, j, lo, hi) {
st.n13_conflict += 1;
}
if kind[i * n + j] == Kind::None {
kind[i * n + j] = Kind::B13;
kind[j * n + i] = Kind::B13;
}
st.n13 += 1;
}
}
}
for (bidx, bd) in mol.bonds().iter().enumerate() {
let (k, l) = (bd.begin, bd.end);
let nk: Vec<u32> = mol
.neighbors(k)
.map(|(y, _)| y)
.filter(|y| *y != l)
.collect();
let nl: Vec<u32> = mol
.neighbors(l)
.map(|(y, _)| y)
.filter(|y| *y != k)
.collect();
for &i in &nk {
for &j in &nl {
if i == j {
continue; }
let (iu, ju) = (i as usize, j as usize);
let (ku, lu) = (k as usize, l as usize);
if kind[iu * n + ju] != Kind::None {
continue;
}
let mid = |x: usize, y: usize| (b.lower(x, y) + b.upper(x, y)) / 2.0;
let Some((cis, trans)) = torsion_span(
mid(iu, ku),
mid(ku, lu),
mid(lu, ju),
mid(iu, lu),
mid(ku, ju),
) else {
st.n14_degenerate += 1;
continue;
};
let arom = mol.bonds()[bidx].order == omgkit_core::BondOrder::Aromatic;
let (t_lo, t_hi) = stereo_path_torsion(mol, bidx, i, k, j)
.map(|t| (t, t))
.or_else(|| {
ring_path_torsion(
mol,
&Rings {
sets: &ring_sets,
aromatic: &aromatic_ring,
sp3: &sp3_ring,
},
i,
k,
l,
j,
)
})
.unwrap_or_else(|| {
torsion_envelope(shared_ring(&ring_sets, &[i, k, l, j]), arom)
});
let f = |t: f64| {
let c = t.to_radians().cos();
let (a2, b2) = (cis * cis, trans * trans);
(a2 + (b2 - a2) * (1.0 - c) / 2.0).max(0.0).sqrt()
};
tighten(&mut b, iu, ju, f(t_lo) - DIST14_TOL, f(t_hi) + DIST14_TOL);
kind[iu * n + ju] = Kind::B14;
kind[ju * n + iu] = Kind::B14;
st.n14 += 1;
}
}
}
let mut topo = vec![6u8; n * n];
for start in 0..n {
let mut d = vec![u8::MAX; n];
d[start] = 0;
let mut q = std::collections::VecDeque::from([start]);
while let Some(x) = q.pop_front() {
if d[x] >= 6 {
continue;
}
let Ok(xu) = u32::try_from(x) else { continue };
for (y, _) in mol.neighbors(xu) {
let y = y as usize;
if y < n && d[y] == u8::MAX {
d[y] = d[x] + 1;
q.push_back(y);
}
}
}
for j in 0..n {
topo[start * n + j] = d[j].min(6);
}
}
for i in 0..n {
for j in (i + 1)..n {
if kind[i * n + j] != Kind::None {
continue;
}
let d0 = vdw_frac(topo[i * n + j] as usize)
* (params::vdw_radius(mol.atoms()[i].atomic_num)
+ params::vdw_radius(mol.atoms()[j].atomic_num));
if d0 > b.lower(i, j) {
b.set_lower(i, j, d0);
}
}
}
(b, st)
}
fn torsion_span(d_ik: f64, d_kl: f64, d_lj: f64, d_il: f64, d_kj: f64) -> Option<(f64, f64)> {
if d_kl <= 1e-9 {
return None;
}
let zi = (d_ik * d_ik + d_kl * d_kl - d_il * d_il) / (2.0 * d_kl);
let ri2 = d_ik * d_ik - zi * zi;
let zj = (d_lj * d_lj + d_kl * d_kl - d_kj * d_kj) / (2.0 * d_kl);
let rj2 = d_lj * d_lj - zj * zj;
let (ri, rj) = (ri2.max(0.0).sqrt(), rj2.max(0.0).sqrt());
let dz = d_kl - zi - zj;
let cis = (dz * dz + (ri - rj) * (ri - rj)).max(0.0).sqrt();
let trans = (dz * dz + (ri + rj) * (ri + rj)).max(0.0).sqrt();
(cis.is_finite() && trans.is_finite()).then_some((cis, trans))
}
#[cfg(test)]
mod tests {
use super::*;
fn prep(smi: &str) -> MolBuilder {
let mut m = omgkit_io::smiles::parse(smi).expect("SMILES 该解析得了");
omgkit_chem::pipeline::sanitize(&mut m).expect("该 sanitize 得了");
omgkit_io::stereo::perceive_bond_stereo(&mut m);
let r = omgkit_io::canon::classed_ranks(&m);
omgkit_chem::add_explicit_hs(&mut m, &r);
m
}
#[test]
fn every_pair_has_a_non_empty_interval() {
let mut checked = 0;
let mut n_mol = 0;
let smis = [
"CCO", "c1ccccc1", "C1CCCCC1", "C1CC2CCC1CC2", "C1CCC12CCC2", "c1ccc2ccccc2c1", "CC(=O)N(C)C", "C#N", "CC(Cl)=C=C(C)Cl", "FS(F)(F)(F)(F)F", "N[Co](N)(N)(N)(N)N", "C1CCCCCCCCCCC1", "C1CC1", "C1CS1", "C1CO1", "C1NN1", "C1COC1", "C1C(CS1)O", "C1CC2C1C2", "C1CCC2(C1)CC2", ];
for smi in smis {
let m = prep(smi);
let (b, _) = build(&m);
n_mol += 1;
for i in 0..b.len() {
for j in (i + 1)..b.len() {
assert!(
b.lower(i, j) <= b.upper(i, j),
"{smi}:第 {i}/{j} 对的下限 {} 超过了上限 {}",
b.lower(i, j),
b.upper(i, j)
);
checked += 1;
}
}
}
assert_eq!(n_mol, smis.len(), "有分子没验到");
assert!(checked > 300, "只验了 {checked} 对");
}
#[test]
fn a_bond_gets_the_measured_interval() {
let m = prep("CC");
let (b, st) = build(&m);
assert!(st.n12 >= 7, "乙烷该有 7 根键,写了 {}", st.n12);
let bd = &m.bonds()[0];
let want = params::bond_length(
m.atoms()[bd.begin as usize].atomic_num,
m.atoms()[bd.end as usize].atomic_num,
bd.order,
0,
);
let (i, j) = (bd.begin as usize, bd.end as usize);
assert!(
b.lower(i, j) <= want.value && want.value <= b.upper(i, j),
"区间 [{:.4}, {:.4}] 没包住表里的中位 {:.4}",
b.lower(i, j),
b.upper(i, j),
want.value
);
let w = b.upper(i, j) - b.lower(i, j);
assert!(w < 0.06, "键长区间宽 {w:.4} —— 太松,精修就没东西可依");
}
#[test]
fn a_one_three_distance_comes_from_the_law_of_cosines() {
let m = prep("O");
let (b, _) = build(&m);
let hs: Vec<usize> = (0..m.num_atoms())
.filter(|i| m.atoms()[*i].atomic_num == 1)
.collect();
assert_eq!(hs.len(), 2, "水该有两个氢");
let (lo, hi) = (b.lower(hs[0], hs[1]), b.upper(hs[0], hs[1]));
let bl = params::bond_length(8, 1, omgkit_core::BondOrder::Single, 0).value;
let ang = params::angle(8, 2, false, 0, 0).value.to_degrees();
let d = third_side(bl, bl, ang);
assert!(
lo <= d && d <= hi,
"H···H 区间 [{lo:.3}, {hi:.3}] 没包住表算的 {d:.3}(键长 {bl:.3}、键角 {ang:.1}°)"
);
assert!(hi - lo < 0.2, "区间 [{lo:.3}, {hi:.3}] 太宽了");
}
#[test]
fn a_hypervalent_centre_gets_a_polyhedral_envelope_not_an_arbitrary_floor() {
let m = prep("FS(F)(F)(F)(F)F");
let (b, st) = build(&m);
assert!(st.angle_envelope > 0, "六配位的硫该走包络");
let fs: Vec<usize> = (0..m.num_atoms())
.filter(|i| m.atoms()[*i].atomic_num == 9)
.collect();
assert_eq!(fs.len(), 6);
let (lo, hi) = (b.lower(fs[0], fs[1]), b.upper(fs[0], fs[1]));
assert!(
lo > 1.5,
"F···F 下限 {lo:.3} —— RDKit 那个无依据的 1.0 不该出现在这里"
);
assert!(hi < 3.6, "F···F 上限 {hi:.3} 该被 180° 压住");
assert!(lo < hi);
}
#[test]
fn a_ring_needs_no_special_case_the_smoothing_finds_the_short_way_round() {
let m = prep("c1ccccc1");
let (mut b, _) = build(&m);
crate::smooth::triangle_smooth(&mut b).expect("苯的界该可行");
let cs: Vec<usize> = (0..m.num_atoms())
.filter(|i| m.atoms()[*i].atomic_num == 6)
.collect();
assert_eq!(cs.len(), 6);
let mut worst = 0.0f64;
for x in 0..6 {
for y in (x + 1)..6 {
worst = crate::linalg::max_nan_wins(worst, b.upper(cs[x], cs[y]));
}
}
assert!(
worst < 6.0,
"苯里最远的一对碳上限 {worst:.3} —— 光滑化没把 MAX_UPPER 压下来"
);
}
#[test]
fn 三元环里成键那一对不许被_1_3_挤宽() {
for smi in ["C1CS1", "C1CO1", "C1CC1", "C1NN1"] {
let m = prep(smi);
let (b, _) = build(&m);
let z: Vec<u8> = m.atoms().iter().map(|a| a.atomic_num).collect();
for bd in m.bonds() {
let (i, j) = (bd.begin as usize, bd.end as usize);
if z[i] == 1 || z[j] == 1 {
continue;
}
let w = b.upper(i, j) - b.lower(i, j);
assert!(
w <= 2.0 * DIST12_TOL + 1e-9,
"{smi} 的键 {i}-{j} 宽 {w:.4},超过 2×DIST12_TOL —— 1-3 挤到键上了"
);
assert!(
b.lower(i, j) <= b.upper(i, j),
"{smi} 的键 {i}-{j} 区间空了"
);
}
}
}
#[test]
fn a_one_four_across_an_aromatic_ring_has_a_pinned_torsion() {
let m = prep("c1ccccc1");
let (b, _) = build(&m);
let n = m.num_atoms();
let mut widths = Vec::new();
for i in 0..n {
for j in (i + 1)..n {
if topo_dist(&m, i, j) != 3 {
continue;
}
widths.push(b.upper(i, j) - b.lower(i, j));
}
}
assert!(
widths.len() >= 6,
"苯该有不少 1-4 对,只找到 {}",
widths.len()
);
let worst = widths
.iter()
.fold(0.0f64, |a, x| crate::linalg::max_nan_wins(a, *x));
assert!(
worst <= 2.0 * DIST14_TOL + 1e-9,
"苯的 1-4 最宽 {worst:.4} Å,超过钉死宽度 {:.4} —— 芳环扭转没被钉住",
2.0 * DIST14_TOL
);
}
fn topo_dist(m: &MolBuilder, a: usize, b: usize) -> u8 {
let n = m.num_atoms();
let mut d = vec![u8::MAX; n];
d[a] = 0;
let mut q = std::collections::VecDeque::from([a]);
while let Some(x) = q.pop_front() {
if d[x] >= 4 {
continue;
}
for (y, _) in m.neighbors(u32::try_from(x).expect("下标")) {
let y = y as usize;
if d[y] == u8::MAX {
d[y] = d[x] + 1;
q.push_back(y);
}
}
}
d[b].min(4)
}
#[test]
fn cis_and_trans_get_different_bounds() {
let cis = prep(r"C/C=C\C");
let trans = prep(r"C/C=C/C");
let mut got = Vec::new();
for m in [&cis, &trans] {
let (b, _) = build(m);
let cs: Vec<usize> = (0..m.num_atoms())
.filter(|i| {
m.atoms()[*i].atomic_num == 6
&& m.neighbors(u32::try_from(*i).expect("下标"))
.filter(|(y, _)| m.atoms()[*y as usize].atomic_num > 1)
.count()
== 1
})
.collect();
assert_eq!(cs.len(), 2, "2-丁烯该有两个端甲基");
let (lo, hi) = (b.lower(cs[0], cs[1]), b.upper(cs[0], cs[1]));
assert!(hi - lo < 0.30, "区间 [{lo:.3}, {hi:.3}] 没被立体标记钉住");
got.push(((lo) + (hi)) / 2.0);
}
assert!(
got[1] - got[0] > 0.5,
"顺式 {:.3}、反式 {:.3} —— 立体那一支没生效,两个写法拿到了同一个区间",
got[0],
got[1]
);
}
#[test]
fn tiny_molecules_are_fine() {
for smi in ["C", "[H][H]", "O"] {
let m = prep(smi);
let (b, st) = build(&m);
assert_eq!(b.len(), m.num_atoms());
assert_eq!(st.atoms, m.num_atoms());
}
}
}