molcrafts-molrs 0.7.0

//! MMFF94 parameter resolution: equivalence-level table lookups, type codes,
//! and the empirical-rule fallbacks used when explicit parameters are absent.
//!
//! Ported from RDKit (BSD-3, Paolo Tosco / RDKit contributors):
//! - `Code/ForceField/MMFF/Params.h` (the `operator()` / `getMMFF*Params`
//!   equivalence-level lookups + vdW combining rules).
//! - `Code/GraphMol/ForceFieldHelpers/MMFF/AtomTyper.cpp`
//!   (`getMMFFAngleType` / `getMMFFBondType` / `getMMFFStretchBendType` /
//!   `getMMFFTorsionType`, and the `getMMFF*EmpiricalRuleParams` fallbacks).
//!
//! Reference: Halgren, T. A. *J. Comput. Chem.* 1996, 17, 490-641 (MMFF.I-V).

use crate::ff::mmff::MmffVariant;
use crate::ff::mmff::charges::mmff_bond_type;
use crate::ff::mmff::tables::{
    MmffProp, mmff_angle, mmff_bndk, mmff_bond, mmff_cov_rad_pau_ele, mmff_def, mmff_dfsb,
    mmff_herschbach_laurie, mmff_oop, mmff_oop_s, mmff_prop, mmff_stbn, mmff_tor, mmff_tor_s,
    mmff_vdw,
};
use crate::ff::mmff::topo::{BondOrder, Topo};

// MMFF numeric constants live in the crate-level `constants` module so the
// `potential` adapters share one definition; re-exported here under the names
// the energy kernels already import (`COULOMB` keeps its short local name).
pub(super) use crate::ff::constants::{
    COULOMB_MMFF as COULOMB, DEG2RAD, ELE_BUFFER, MDYNE_A_TO_KCAL, RAD2DEG,
};

// vdW combining-rule globals (RDKit `MMFFVdWCollection`, Params.cpp:
// "power B Beta DARAD DAEPS" = "0.25 0.2 12. 0.8 0.5").
const VDW_B: f64 = 0.2;
const VDW_BETA: f64 = 12.0;
const VDW_DARAD: f64 = 0.8;
const VDW_DAEPS: f64 = 0.5;

const DA_DONOR: u8 = b'D';
const DA_ACCEPTOR: u8 = b'A';

#[inline]
fn is_zero(x: f64) -> bool {
    x.abs() < 1.0e-8
}

/// MMFF resolved bond parameters.
#[derive(Clone, Copy, Debug)]
pub(crate) struct BondParams {
    pub r0: f64,
    pub kb: f64,
}

/// MMFF resolved angle parameters.
#[derive(Clone, Copy, Debug)]
pub(crate) struct AngleParams {
    pub theta0: f64,
    pub ka: f64,
}

/// MMFF resolved stretch-bend force constants (`kba_ijk`, `kba_kji`).
#[derive(Clone, Copy, Debug)]
pub(crate) struct StbnParams {
    pub kba_ijk: f64,
    pub kba_kji: f64,
}

/// MMFF resolved torsion Fourier coefficients.
#[derive(Clone, Copy, Debug)]
pub(crate) struct TorParams {
    pub v1: f64,
    pub v2: f64,
    pub v3: f64,
}

/// Combined vdW pair parameters (already donor/acceptor-scaled).
#[derive(Clone, Copy, Debug)]
pub(super) struct VdwParams {
    pub r_star: f64,
    pub epsilon: f64,
}

// --- periodic-table rows (AtomTyper.cpp) ---------------------------------

pub(super) fn periodic_row(atno: u8) -> u8 {
    match atno {
        3..=10 => 1,
        11..=18 => 2,
        19..=36 => 3,
        37..=54 => 4,
        _ => 0,
    }
}

fn periodic_row_hl(atno: u8) -> u8 {
    let mut row = match atno {
        2 => 1,
        3..=10 => 2,
        11..=18 => 3,
        19..=36 => 4,
        37..=54 => 5,
        _ => 0,
    };
    if (21..=30).contains(&atno) || (39..=48).contains(&atno) {
        row *= 10;
    }
    row
}

// --- bond type (single sbmb/arom pair) -----------------------------------

/// RDKit `getMMFFBondType`. Delegates to the already-ported charge helper.
pub(super) fn bond_type(topo: &Topo, types: &[u8], i: usize, j: usize) -> u8 {
    mmff_bond_type(topo, types, i, j)
}

// --- ring-size helpers (AtomTyper.cpp) -----------------------------------

/// RDKit `isAngleInRingOfSize3or4`.
fn angle_ring_size(topo: &Topo, i: usize, j: usize, k: usize) -> u8 {
    if topo.bond_order(i, j).is_none() || topo.bond_order(j, k).is_none() {
        return 0;
    }
    if topo.bond_order(k, i).is_some() {
        return 3;
    }
    // 4-ring: a shared neighbour of i and k (other than j).
    for &a in &topo.nbrs[i] {
        if a == j {
            continue;
        }
        if topo.nbrs[k].iter().any(|&b| b != j && b == a) {
            return 4;
        }
    }
    0
}

/// RDKit `isTorsionInRingOfSize4or5`.
fn torsion_ring_size(topo: &Topo, i: usize, j: usize, k: usize, l: usize) -> u8 {
    if topo.bond_order(i, j).is_none()
        || topo.bond_order(j, k).is_none()
        || topo.bond_order(k, l).is_none()
    {
        return 0;
    }
    if topo.bond_order(l, i).is_some() {
        return 4;
    }
    for &a in &topo.nbrs[i] {
        if a == j {
            continue;
        }
        if topo.nbrs[l].iter().any(|&b| b != k && b == a) {
            return 5;
        }
    }
    0
}

// --- angle / stretch-bend / torsion type codes (AtomTyper.cpp) -----------

/// RDKit `getMMFFAngleType`.
pub(super) fn angle_type(topo: &Topo, types: &[u8], i: usize, j: usize, k: usize) -> u8 {
    let bts = bond_type(topo, types, i, j) + bond_type(topo, types, j, k);
    let mut at = bts;
    let size = angle_ring_size(topo, i, j, k);
    if size != 0 {
        at = size;
        if bts != 0 {
            at += bts + size - 2;
        }
    }
    at
}

/// RDKit `getMMFFStretchBendType`.
fn stretch_bend_type(angle_type: u8, bt1: u8, bt2: u8) -> u8 {
    match angle_type {
        1 => {
            if bt1 != 0 || bt1 == bt2 {
                1
            } else {
                2
            }
        }
        2 => 3,
        4 => 4,
        3 => 5,
        5 => {
            if bt1 != 0 || bt1 == bt2 {
                6
            } else {
                7
            }
        }
        6 => 8,
        7 => {
            if bt1 != 0 || bt1 == bt2 {
                9
            } else {
                10
            }
        }
        8 => 11,
        _ => 0,
    }
}

/// RDKit `getMMFFTorsionType`. Returns `(principal, secondary)`.
pub(super) fn torsion_type(
    topo: &Topo,
    types: &[u8],
    i: usize,
    j: usize,
    k: usize,
    l: usize,
) -> (u8, u8) {
    let bt_ij = bond_type(topo, types, i, j);
    let bt_jk = bond_type(topo, types, j, k);
    let bt_kl = bond_type(topo, types, k, l);
    let jk_single = topo.bond_order(j, k) == Some(BondOrder::Single);
    let mut tt = bt_jk;
    let mut second = 0u8;
    if bt_jk == 0 && jk_single && (bt_ij == 1 || bt_kl == 1) {
        tt = 2;
    }
    let size = torsion_ring_size(topo, i, j, k, l);
    if size == 4 && topo.bond_order(i, k).is_none() && topo.bond_order(j, l).is_none() {
        second = tt;
        tt = 4;
    } else if size == 5 && (types[i] == 1 || types[j] == 1 || types[k] == 1 || types[l] == 1) {
        second = tt;
        tt = 5;
    }
    (tt, second)
}

// --- equivalence-level lookups (Params.h operator()) ---------------------

/// RDKit `MMFFDefCollection::operator()` equivalence-level lookup.
///
/// `MMFF_DEF` (regenerated from RDKit `defaultMMFFDef`) now covers every MMFF
/// atom type 1-82 / 87-99, so the 5-stage equivalence search for aromatic /
/// charged-N / metal types reaches its wild-card defaults uniformly. Falls
/// back to the type itself for any key absent from the table (RDKit returns a
/// null params pointer there, which the callers below treat as "no match").
fn eq_level(atom_type: u8, level: usize) -> u8 {
    mmff_def(atom_type)
        .map(|d| d.eq_level[level])
        .unwrap_or(atom_type)
}

/// RDKit `MMFFAngleCollection::operator()` (5-stage equivalence search).
fn angle_lookup(angle_type: u8, i: u8, j: u8, k: u8) -> Option<AngleParams> {
    for iter in 0..4 {
        let (mut ci, mut ck) = (eq_level(i, iter), eq_level(k, iter));
        if ci > ck {
            std::mem::swap(&mut ci, &mut ck);
        }
        if let Some(a) = mmff_angle(angle_type, ci, j, ck) {
            return Some(AngleParams {
                theta0: a.theta0,
                ka: a.ka,
            });
        }
    }
    None
}

/// RDKit `MMFFOopCollection::operator()` (5-stage equivalence search).
///
/// RDKit selects the whole table at construction
/// (`defaultMMFFsOop` vs `defaultMMFFOop`) based on `isMMFFs`. The `_S` table
/// shares every key with the base table, so we look up `_S` first under
/// [`MmffVariant::Mmff94s`] and fall back to the base table for safety.
fn oop_lookup(variant: MmffVariant, i: u8, j: u8, k: u8, l: u8) -> Option<f64> {
    for iter in 0..4 {
        let mut ikl = [eq_level(i, iter), eq_level(k, iter), eq_level(l, iter)];
        ikl.sort_unstable();
        let hit = match variant {
            MmffVariant::Mmff94s => mmff_oop_s(ikl[0], j, ikl[1], ikl[2])
                .or_else(|| mmff_oop(ikl[0], j, ikl[1], ikl[2])),
            MmffVariant::Mmff94 => mmff_oop(ikl[0], j, ikl[1], ikl[2]),
        };
        if let Some(o) = hit {
            return Some(o.koop);
        }
    }
    None
}

/// RDKit `MMFFTorCollection::getMMFFTorParams`. Returns `(matched_type, params)`.
///
/// Like the Oop collection, RDKit picks `defaultMMFFsTor` vs `defaultMMFFTor`
/// wholesale on `isMMFFs`; the `_S` table shares every key with the base
/// table, so under [`MmffVariant::Mmff94s`] we query `_S` first and fall back
/// to the base table.
fn torsion_lookup(
    variant: MmffVariant,
    tor_type: (u8, u8),
    i: u8,
    j: u8,
    k: u8,
    l: u8,
) -> Option<(u8, TorParams)> {
    let mut iter: i32 = 0;
    let mut max_iter = 5i32;
    let mut can_tor = tor_type.0;
    // Mirrors the RDKit `while` guard: keep iterating while we have not yet
    // found a hit (we return on the first), with the special last-resort
    // restart when torType is (5, secondary).
    loop {
        let cont = (iter < max_iter) || (iter == 4 && tor_type.0 == 5 && tor_type.1 != 0);
        if !cont {
            break;
        }
        if max_iter == 5 && iter == 4 {
            max_iter = 4;
            iter = 0;
            can_tor = tor_type.1;
        }
        let (mut i_wild, mut l_wild) = (iter as usize, iter as usize);
        if iter == 1 {
            i_wild = 1;
            l_wild = 3;
        } else if iter == 2 {
            i_wild = 3;
            l_wild = 1;
        }
        let mut ci = eq_level(i, i_wild);
        let mut cj = j;
        let mut ck = k;
        let mut cl = eq_level(l, l_wild);
        if cj > ck {
            std::mem::swap(&mut cj, &mut ck);
            std::mem::swap(&mut ci, &mut cl);
        } else if cj == ck && ci > cl {
            std::mem::swap(&mut ci, &mut cl);
        }
        let hit = match variant {
            MmffVariant::Mmff94s => {
                mmff_tor_s(can_tor, ci, cj, ck, cl).or_else(|| mmff_tor(can_tor, ci, cj, ck, cl))
            }
            MmffVariant::Mmff94 => mmff_tor(can_tor, ci, cj, ck, cl),
        };
        if let Some(t) = hit {
            return Some((
                can_tor,
                TorParams {
                    v1: t.v1,
                    v2: t.v2,
                    v3: t.v3,
                },
            ));
        }
        iter += 1;
    }
    None
}

// --- bond stretch (explicit + empirical) ---------------------------------

/// RDKit `getMMFFBondStretchParams` + `getMMFFBondStretchEmpiricalRuleParams`.
pub(crate) fn bond_params(topo: &Topo, types: &[u8], i: usize, j: usize) -> Option<BondParams> {
    let bt = bond_type(topo, types, i, j);
    let (ti, tj) = (types[i].min(types[j]), types[i].max(types[j]));
    if let Some(b) = mmff_bond(bt, ti, tj) {
        return Some(BondParams { r0: b.r0, kb: b.kb });
    }
    bond_empirical(topo, i, j)
}

/// RDKit `getMMFFBondStretchEmpiricalRuleParams` (MMFF.V eq. 18/19 + HL rule).
fn bond_empirical(topo: &Topo, i: usize, j: usize) -> Option<BondParams> {
    let (an1, an2) = (topo.atno[i], topo.atno[j]);
    let cr1 = mmff_cov_rad_pau_ele(an1)?;
    let cr2 = mmff_cov_rad_pau_ele(an2)?;
    let c = if an1 == 1 || an2 == 1 { 0.050 } else { 0.085 };
    let n = 1.4;
    let r0 = cr1.r0 + cr2.r0 - c * (cr1.chi - cr2.chi).abs().powf(n);
    let kb = if let Some(bndk) = mmff_bndk(an1, an2) {
        let coeff = bndk.r0 / r0;
        let coeff2 = coeff * coeff;
        bndk.kb * coeff2 * coeff2 * coeff2
    } else {
        let hl = mmff_herschbach_laurie(periodic_row_hl(an1), periodic_row_hl(an2))?;
        10f64.powf(-(r0 - hl.a_ij) / hl.d_ij)
    };
    Some(BondParams { r0, kb })
}

// --- angle bend (explicit + empirical) -----------------------------------

/// RDKit `getMMFFAngleBendParams`.
pub(crate) fn angle_params(
    topo: &Topo,
    types: &[u8],
    i: usize,
    j: usize,
    k: usize,
) -> Option<AngleParams> {
    let at = angle_type(topo, types, i, j, k);
    let explicit = angle_lookup(at, types[i], types[j], types[k]);
    match explicit {
        Some(p) if !is_zero(p.ka) => Some(p),
        _ => {
            let b1 = bond_params(topo, types, i, j)?;
            let b2 = bond_params(topo, types, j, k)?;
            angle_empirical(topo, types, explicit, (&b1, &b2), [i, j, k])
        }
    }
}

/// RDKit `getMMFFAngleBendEmpiricalRuleParams` (MMFF.V eq. 20 + Table VI).
fn angle_empirical(
    topo: &Topo,
    types: &[u8],
    old: Option<AngleParams>,
    bonds: (&BondParams, &BondParams),
    idx: [usize; 3],
) -> Option<AngleParams> {
    let [i, j, k] = idx;
    let (b1, b2) = bonds;
    let atno = [topo.atno[i], topo.atno[j], topo.atno[k]];
    let ring = angle_ring_size(topo, i, j, k);
    let prop_j = mmff_prop(types[j])?;

    let theta0 = if let Some(o) = old {
        o.theta0
    } else {
        let mut t = 120.0;
        match prop_j.crd {
            4 => t = 109.45,
            2 => {
                if atno[1] == 8 {
                    t = 105.0;
                } else if prop_j.linh == 1 {
                    t = 180.0;
                }
            }
            3 if prop_j.val == 3 && prop_j.mltb == 0 => {
                t = if atno[1] == 7 { 107.0 } else { 92.0 };
            }
            _ => {}
        }
        if ring == 3 {
            t = 60.0;
        } else if ring == 4 {
            t = 90.0;
        }
        t
    };

    // Table VI (MMFF.V p.628) Z and C constants.
    let zc = |an: u8| -> (f64, f64) {
        match an {
            1 => (1.395, 0.0),
            6 => (2.494, 1.016),
            7 => (2.711, 1.113),
            8 => (3.045, 1.337),
            9 => (2.847, 0.0),
            14 => (2.350, 0.811),
            15 => (2.350, 1.068),
            16 => (2.980, 1.249),
            17 => (2.909, 1.078),
            35 => (3.017, 0.0),
            53 => (3.086, 0.0),
            _ => (0.0, 0.0),
        }
    };
    let (z0, _) = zc(atno[0]);
    let (_, c1) = zc(atno[1]);
    let (z2, _) = zc(atno[2]);

    let mut beta = 1.75;
    let r0ij = b1.r0;
    let r0jk = b2.r0;
    let d = (r0ij - r0jk) * (r0ij - r0jk) / ((r0ij + r0jk) * (r0ij + r0jk));
    let theta0_rad = DEG2RAD * theta0;
    if ring == 4 {
        beta *= 0.85;
    } else if ring == 3 {
        beta *= 0.05;
    }
    let ka = beta * z0 * c1 * z2 / ((r0ij + r0jk) * theta0_rad * theta0_rad * (2.0 * d).exp());
    Some(AngleParams { theta0, ka })
}

// --- stretch-bend (explicit + dfsb fallback) -----------------------------

/// RDKit `getMMFFStretchBendParams`. Returns the resolved params plus the
/// two bond rest lengths and the angle theta0 (needed by the energy term).
pub(crate) fn stretch_bend_params(
    topo: &Topo,
    types: &[u8],
    i: usize,
    j: usize,
    k: usize,
) -> Option<(StbnParams, f64, f64, f64)> {
    let prop_j = mmff_prop(types[j])?;
    if prop_j.linh != 0 {
        return None;
    }
    let b1 = bond_params(topo, types, i, j)?;
    let b2 = bond_params(topo, types, j, k)?;
    let angle = angle_params(topo, types, i, j, k)?;

    let bt1 = bond_type(topo, types, i, j);
    let bt2 = bond_type(topo, types, j, k);
    let at = angle_type(topo, types, i, j, k);
    let (ti, tk) = (types[i], types[k]);
    let sbt = stretch_bend_type(
        at,
        if ti <= tk { bt1 } else { bt2 },
        if ti < tk { bt2 } else { bt1 },
    );

    // explicit STBN lookup with the same swap convention as RDKit.
    let (swap, stbn) = stbn_lookup(sbt, bt1, bt2, ti, types[j], tk);
    let params = if let Some(s) = stbn {
        if swap {
            StbnParams {
                kba_ijk: s.1,
                kba_kji: s.0,
            }
        } else {
            StbnParams {
                kba_ijk: s.0,
                kba_kji: s.1,
            }
        }
    } else {
        // default-by-period-row (mmff_dfsb)
        let (swap_d, d) = dfsb_lookup(topo.atno[i], topo.atno[j], topo.atno[k]);
        let d = d?;
        if swap_d {
            StbnParams {
                kba_ijk: d.1,
                kba_kji: d.0,
            }
        } else {
            StbnParams {
                kba_ijk: d.0,
                kba_kji: d.1,
            }
        }
    };
    if is_zero(params.kba_ijk) && is_zero(params.kba_kji) {
        return None;
    }
    Some((params, b1.r0, b2.r0, angle.theta0))
}

/// RDKit `MMFFStbnCollection::getMMFFStbnParams`. Returns `(swap, (ijk,kji))`.
fn stbn_lookup(sbt: u8, bt1: u8, bt2: u8, i: u8, j: u8, k: u8) -> (bool, Option<(f64, f64)>) {
    let (mut ci, mut ck) = (i, k);
    let swap = if i > k {
        std::mem::swap(&mut ci, &mut ck);
        true
    } else if i == k {
        bt1 < bt2
    } else {
        false
    };
    let p = mmff_stbn(sbt, ci, j, ck).map(|s| (s.kba_ijk, s.kba_kji));
    (swap, p)
}

/// RDKit `MMFFDfsbCollection::getMMFFDfsbParams`. Returns `(swap, (ijk,kji))`.
fn dfsb_lookup(an_i: u8, an_j: u8, an_k: u8) -> (bool, Option<(f64, f64)>) {
    let (mut r1, mut r3) = (periodic_row(an_i), periodic_row(an_k));
    let swap = if r1 > r3 {
        std::mem::swap(&mut r1, &mut r3);
        true
    } else {
        false
    };
    let p = mmff_dfsb(r1, periodic_row(an_j), r3).map(|d| (d.kba_ijk, d.kba_kji));
    (swap, p)
}

// --- out-of-plane --------------------------------------------------------

/// RDKit `getMMFFOopBendParams` (no empirical fallback; term excluded if absent).
pub(crate) fn oop_koop(
    variant: MmffVariant,
    types: &[u8],
    i: usize,
    j: usize,
    k: usize,
    l: usize,
) -> Option<f64> {
    oop_lookup(variant, types[i], types[j], types[k], types[l])
}

// --- torsion (explicit + empirical) --------------------------------------

/// RDKit `getMMFFTorsionParams`. Returns `None` when all coefficients vanish.
pub(crate) fn torsion_params(
    variant: MmffVariant,
    topo: &Topo,
    types: &[u8],
    i: usize,
    j: usize,
    k: usize,
    l: usize,
) -> Option<TorParams> {
    let tt = torsion_type(topo, types, i, j, k, l);
    let p = match torsion_lookup(variant, tt, types[i], types[j], types[k], types[l]) {
        Some((_, p)) => p,
        None => torsion_empirical(topo, types, j, k)?,
    };
    if is_zero(p.v1) && is_zero(p.v2) && is_zero(p.v3) {
        None
    } else {
        Some(p)
    }
}

/// RDKit `getMMFFTorsionEmpiricalRuleParams` (MMFF.V rules a-h, p.632).
fn torsion_empirical(topo: &Topo, types: &[u8], j: usize, k: usize) -> Option<TorParams> {
    let jp = mmff_prop(types[j])?;
    let kp = mmff_prop(types[k])?;
    let jt = types[j];
    let kt = types[k];
    let atno = [topo.atno[j], topo.atno[k]];
    let bond = topo.bond_order(j, k);
    let aromatic = crate::ff::mmff::tables::mmff_is_arom(jt)
        && crate::ff::mmff::tables::mmff_is_arom(kt)
        && topo.is_aromatic[j]
        && topo.is_aromatic[k]
        && bond == Some(BondOrder::Aromatic);

    let mut u = [0.0f64; 2];
    let mut v = [0.0f64; 2];
    let mut w = [0.0f64; 2];
    for (idx, &an) in atno.iter().enumerate() {
        match an {
            6 => {
                u[idx] = 2.0;
                v[idx] = 2.12;
            }
            7 => {
                u[idx] = 2.0;
                v[idx] = 1.5;
            }
            8 => {
                u[idx] = 2.0;
                v[idx] = 0.2;
                w[idx] = 2.0;
            }
            14 => {
                u[idx] = 1.25;
                v[idx] = 1.22;
            }
            15 => {
                u[idx] = 1.25;
                v[idx] = 2.40;
            }
            16 => {
                u[idx] = 1.25;
                v[idx] = 0.49;
                w[idx] = 8.0;
            }
            _ => {}
        }
    }

    let n_jk = ((jp.crd as i32 - 1) * (kp.crd as i32 - 1)) as f64;
    let mut tor = TorParams {
        v1: 0.0,
        v2: 0.0,
        v3: 0.0,
    };
    let (beta, pi_jk);
    let row = |a: u8| periodic_row(a);

    if jp.linh != 0 || kp.linh != 0 {
        // rule (a): all zero
    } else if aromatic {
        // rule (b)
        beta = if (jp.val == 3 && kp.val == 4) || (jp.val == 4 && kp.val == 3) {
            3.0
        } else {
            6.0
        };
        pi_jk = if jp.pilp == 0 && kp.pilp == 0 {
            0.5
        } else {
            0.3
        };
        tor.v2 = beta * pi_jk * (u[0] * u[1]).sqrt();
    } else if bond == Some(BondOrder::Double) {
        // rule (c)
        beta = 6.0;
        pi_jk = if jp.mltb == 2 && kp.mltb == 2 {
            1.0
        } else {
            0.4
        };
        tor.v2 = beta * pi_jk * (u[0] * u[1]).sqrt();
    } else if jp.crd == 4 && kp.crd == 4 {
        // rule (d)
        tor.v3 = (v[0] * v[1]).sqrt() / n_jk;
    } else if jp.crd == 4 && kp.crd != 4 {
        // rule (e)
        if (kp.crd == 3 && ((kp.val == 4 || kp.val == 34) || kp.mltb != 0))
            || (kp.crd == 2 && (kp.val == 3 || kp.mltb != 0))
        {
            // zero
        } else {
            tor.v3 = (v[0] * v[1]).sqrt() / n_jk;
        }
    } else if kp.crd == 4 && jp.crd != 4 {
        // rule (f)
        if (jp.crd == 3 && ((jp.val == 4 || jp.val == 34) || jp.mltb != 0))
            || (jp.crd == 2 && (jp.val == 3 || jp.mltb != 0))
        {
            // zero
        } else {
            tor.v3 = (v[0] * v[1]).sqrt() / n_jk;
        }
    } else if (bond == Some(BondOrder::Single) && jp.mltb != 0 && kp.mltb != 0)
        || (jp.mltb != 0 && kp.pilp != 0)
        || (jp.pilp != 0 && kp.mltb != 0)
    {
        // rule (g)
        if jp.pilp != 0 && kp.pilp != 0 {
            // case (1): zero
        } else if jp.pilp != 0 && kp.mltb != 0 {
            // case (2)
            let pi = if jp.mltb == 1 {
                0.5
            } else if row(atno[0]) == 2 && row(atno[1]) == 2 {
                0.3
            } else {
                0.15
            };
            tor.v2 = 6.0 * pi * (u[0] * u[1]).sqrt();
        } else if kp.pilp != 0 && jp.mltb != 0 {
            // case (3)
            let pi = if kp.mltb == 1 {
                0.5
            } else if row(atno[0]) == 2 && row(atno[1]) == 2 {
                0.3
            } else {
                0.15
            };
            tor.v2 = 6.0 * pi * (u[0] * u[1]).sqrt();
        } else if (jp.mltb == 1 || kp.mltb == 1) && (atno[0] != 6 || atno[1] != 6) {
            // case (4)
            tor.v2 = 6.0 * 0.4 * (u[0] * u[1]).sqrt();
        } else {
            // case (5)
            tor.v2 = 6.0 * 0.15 * (u[0] * u[1]).sqrt();
        }
    } else {
        // rule (h)
        if (atno[0] == 8 || atno[0] == 16) && (atno[1] == 8 || atno[1] == 16) {
            tor.v2 = -(w[0] * w[1]).sqrt();
        } else {
            tor.v3 = (v[0] * v[1]).sqrt() / n_jk;
        }
    }
    Some(tor)
}

// --- van der Waals (combining rules) -------------------------------------

/// RDKit `getMMFFVdWParams` + `calcUnscaledVdWMinimum`/`WellDepth` + scaling.
pub(super) fn vdw_params(types: &[u8], i: usize, j: usize) -> Option<VdwParams> {
    let pi = mmff_vdw(types[i])?;
    let pj = mmff_vdw(types[j])?;

    let gamma = (pi.r_star - pj.r_star) / (pi.r_star + pj.r_star);
    let donor = pi.da == DA_DONOR || pj.da == DA_DONOR;
    let mut r_star = 0.5
        * (pi.r_star + pj.r_star)
        * (1.0
            + if donor {
                0.0
            } else {
                VDW_B * (1.0 - (-VDW_BETA * gamma * gamma).exp())
            });

    let r2 = r_star * r_star;
    let c4 = 181.16;
    let mut epsilon = c4 * pi.g_i * pj.g_i * pi.alpha_i * pj.alpha_i
        / (((pi.alpha_i / pi.n_i).sqrt() + (pj.alpha_i / pj.n_i).sqrt()) * r2 * r2 * r2);

    // donor/acceptor scaling (RDKit `scaleVdWParams`)
    let da_pair =
        (pi.da == DA_DONOR && pj.da == DA_ACCEPTOR) || (pi.da == DA_ACCEPTOR && pj.da == DA_DONOR);
    if da_pair {
        r_star *= VDW_DARAD;
        epsilon *= VDW_DAEPS;
    }
    Some(VdwParams { r_star, epsilon })
}

/// Periodic-distance class between atoms (RDKit `buildNeighborMatrix`):
/// number of bond hops, capped — 1 (1-2), 2 (1-3), 3 (1-4), or `u8::MAX`.
pub(super) fn relation(topo: &Topo, a: usize, b: usize) -> u8 {
    // BFS up to depth 3.
    use std::collections::VecDeque;
    let n = topo.n_atoms();
    let mut dist = vec![u8::MAX; n];
    dist[a] = 0;
    let mut q = VecDeque::new();
    q.push_back(a);
    while let Some(u) = q.pop_front() {
        if dist[u] >= 3 {
            continue;
        }
        for &v in &topo.nbrs[u] {
            if dist[v] == u8::MAX {
                dist[v] = dist[u] + 1;
                q.push_back(v);
            }
        }
    }
    dist[b]
}

/// Helper exposing the central-atom prop (linear flag) for angle bend.
pub(super) fn central_prop(types: &[u8], j: usize) -> Option<&'static MmffProp> {
    mmff_prop(types[j])
}