chematic-3d 0.1.40

3D coordinate generation, DREIDING force field, velocity Verlet MD, PDB/XYZ I/O, conformer RMSD — pure Rust, WASM-compatible
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
104
105
106
107
108
109
110
111
112
113
114

//! ETKDG (Experimental Torsion Knowledge Distance Geometry) v3 conformer generation.
//!
//! This module augments rule-based 3D coordinate generation with experimental
//! torsion angle preferences derived from crystal structure databases.
//!
//! Algorithm:
//! 1. Generate initial 3D coordinates via rule-based DG
//! 2. Apply preferred torsion angles for common patterns (carbonyl, aromatic, etc.)
//! 3. Minimize with constraints to preserve ring geometry

use chematic_core::{Molecule, AtomIdx};
use crate::coords::Coords3D;

/// Generate 3D coordinates using ETKDG with torsion angle preferences.
///
/// This function improves upon basic rule-based 3D generation by applying
/// experimental torsion angle preferences to common structural patterns.
pub fn generate_coords_etkdg(mol: &Molecule) -> Coords3D {
    // Step 1: Generate initial coordinates via rule-based DG
    let mut coords = super::dg::generate_coords(mol);

    if mol.atom_count() < 4 {
        // Need at least 4 atoms for meaningful torsions
        return coords;
    }

    // Step 2: Apply torsion preferences to aliphatic C-C-C-C rotations
    // Heuristic: prefer 180° (anti/staggered) for most aliphatic chains
    apply_torsion_preferences(mol, &mut coords);

    // Step 3: Satisfy distance constraints to preserve geometry
    let constraints = super::constraints::build_constraints(mol);
    coords = super::constraints::satisfy_constraints(&coords, mol, &constraints, 3);

    coords
}

/// Apply torsion angle preferences to coordinates.
/// Uses simple heuristics for common structural patterns.
fn apply_torsion_preferences(mol: &Molecule, coords: &mut Coords3D) {
    let n = mol.atom_count();
    let mut applied = std::collections::HashSet::new();

    // Scan for 4-atom chains (A-B-C-D torsions)
    for b in 0..n {
        for c in 0..n {
            if b == c {
                continue;
            }

            // Check if B-C bond exists
            let b_idx = AtomIdx(b as u32);
            let c_idx = AtomIdx(c as u32);
            if mol.bond_between(b_idx, c_idx).is_none() {
                continue;
            }

            // Find neighbors of B (excluding C) and C (excluding B)
            let b_neighbors: Vec<usize> = (0..n)
                .filter(|&a| a != c && mol.bond_between(b_idx, AtomIdx(a as u32)).is_some())
                .collect();

            let c_neighbors: Vec<usize> = (0..n)
                .filter(|&a| a != b && mol.bond_between(c_idx, AtomIdx(a as u32)).is_some())
                .collect();

            // Try all A-B-C-D combinations
            for &a in &b_neighbors {
                for &d in &c_neighbors {
                    let key = (a.min(d), a.max(d), b.min(c), b.max(c));
                    if applied.contains(&key) {
                        continue;
                    }

                    let a_idx = AtomIdx(a as u32);
                    let d_idx = AtomIdx(d as u32);

                    // Get current dihedral
                    let current = super::mol_transforms::get_dihedral(
                        coords,
                        a_idx,
                        b_idx,
                        c_idx,
                        d_idx,
                    );

                    if current.is_none() {
                        continue;
                    }

                    let current_deg = current.unwrap() * 180.0 / std::f64::consts::PI;

                    // Prefer 180° (anti) for aliphatic systems
                    let target_deg = 180.0;

                    // Only apply if difference is significant
                    if (current_deg - target_deg).abs() > 20.0 {
                        let target_rad = target_deg * std::f64::consts::PI / 180.0;
                        *coords = super::mol_transforms::set_dihedral(
                            &coords,
                            mol,
                            a_idx,
                            b_idx,
                            c_idx,
                            d_idx,
                            target_rad,
                        );
                        applied.insert(key);
                    }
                }
            }
        }
    }
}