molframe-pdb 0.6.0

Legacy PDB reading and writing for molframe
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

//! Splits MMTF inter-group bonds and maps bond endpoints onto each model's atoms.
//!
//! Complexity: O(bonds) time and space; O(1) per endpoint lookup.

use super::group_bonds::decoded_order;
use super::validation::{schema_error, usize_of};
use super::{Decoded, ModelRange};
use molframe_core::BondOrder;
use molframe_core::diagnostic::Diagnostic;

/// One global MMTF bond translated to model-local atom indices.
#[derive(Clone, Copy)]
pub(super) struct GlobalBond {
    pub(super) atom_a: u32,
    pub(super) atom_b: u32,
    pub(super) order: BondOrder,
}

/// Partitions global bonds once instead of scanning the complete bond table for
/// every model. Cross-model bonds retain the previous behavior and are ignored.
pub(super) fn partition_global_bonds(
    decoded: &mut Decoded,
    ranges: &[ModelRange],
) -> Result<Vec<Vec<GlobalBond>>, Diagnostic> {
    if !decoded.global_bonds.len().is_multiple_of(2) {
        return Err(schema_error("global bond endpoint list has odd length"));
    }
    let mut counts = vec![0usize; ranges.len()];
    for pair in decoded.global_bonds.as_chunks::<2>().0 {
        if let Some((model, _, _)) = local_bond_endpoints(pair, ranges)? {
            counts[model] = counts[model]
                .checked_add(1)
                .ok_or_else(|| schema_error("MMTF model bond count overflows"))?;
        }
    }
    let mut partitioned = counts
        .into_iter()
        .map(Vec::with_capacity)
        .collect::<Vec<_>>();
    for (bond_index, pair) in decoded.global_bonds.as_chunks::<2>().0.iter().enumerate() {
        let Some((model, atom_a, atom_b)) = local_bond_endpoints(pair, ranges)? else {
            continue;
        };
        partitioned[model].push(GlobalBond {
            atom_a,
            atom_b,
            order: decoded_order(
                &decoded.global_orders,
                &decoded.global_resonance,
                bond_index,
            ),
        });
    }
    decoded.global_bonds = Vec::new();
    decoded.global_orders = Vec::new();
    decoded.global_resonance = Vec::new();
    Ok(partitioned)
}

/// Resolves both global endpoints and returns their common model and offsets.
pub(super) fn local_bond_endpoints(
    pair: &[i32],
    ranges: &[ModelRange],
) -> Result<Option<(usize, u32, u32)>, Diagnostic> {
    let first = usize_of(pair[0])?;
    let second = usize_of(pair[1])?;
    let Some((first_model, first_local)) = model_atom(first, ranges)? else {
        return Ok(None);
    };
    let Some((second_model, second_local)) = model_atom(second, ranges)? else {
        return Ok(None);
    };
    if first_model == second_model {
        Ok(Some((first_model, first_local, second_local)))
    } else {
        Ok(None)
    }
}

/// Binary-searches the disjoint, ascending model atom ranges.
pub(super) fn model_atom(
    atom: usize,
    ranges: &[ModelRange],
) -> Result<Option<(usize, u32)>, Diagnostic> {
    let Some(model) = ranges
        .partition_point(|range| range.atom_start <= atom)
        .checked_sub(1)
    else {
        return Ok(None);
    };
    let range = ranges[model];
    let end = range
        .atom_start
        .checked_add(range.atom_count)
        .ok_or_else(|| schema_error("MMTF model atom range overflows"))?;
    if atom >= end {
        return Ok(None);
    }
    let local = atom - range.atom_start;
    let local =
        u32::try_from(local).map_err(|_| schema_error("MMTF local atom index exceeds u32"))?;
    Ok(Some((model, local)))
}