molcrafts-molrs 0.7.0

Molecular simulation toolkit: core data structures, IO, trajectory analysis, force fields, SMILES, and 3D conformer generation (feature-gated modules)
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
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137

//! MMFF94 torsional rotation: E = 0.5*(V1*(1+cos phi) + V2*(1-cos 2phi) + V3*(1+cos 3phi))

use crate::ff::forcefield::Params;
use crate::ff::potential::Potential;
use crate::ff::potential::geometry::{
    accumulate_dihedral_forces, compute_dihedral, validate_coords,
};
use molrs::store::frame::Frame;
use molrs::types::F;

pub struct MMFFTorsion {
    atom_i: Vec<usize>,
    atom_j: Vec<usize>,
    atom_k: Vec<usize>,
    atom_l: Vec<usize>,
    v1: Vec<F>,
    v2: Vec<F>,
    v3: Vec<F>,
}

impl Potential for MMFFTorsion {
    fn calc_energy_forces(&self, coords: &[F]) -> (F, Vec<F>) {
        let _n = validate_coords(coords);
        let mut energy: F = 0.0;
        let mut forces = vec![0.0 as F; coords.len()];

        for idx in 0..self.atom_i.len() {
            let (i, j, k, l) = (
                self.atom_i[idx],
                self.atom_j[idx],
                self.atom_k[idx],
                self.atom_l[idx],
            );
            let phi = compute_dihedral(coords, i, j, k, l);

            let (s1, c1) = phi.sin_cos();
            let (s2, c2) = (2.0 * phi).sin_cos();
            let (s3, c3) = (3.0 * phi).sin_cos();

            energy += 0.5
                * (self.v1[idx] * (1.0 + c1)
                    + self.v2[idx] * (1.0 - c2)
                    + self.v3[idx] * (1.0 + c3));

            let de_dphi =
                0.5 * (-self.v1[idx] * s1 + 2.0 * self.v2[idx] * s2 - 3.0 * self.v3[idx] * s3);
            accumulate_dihedral_forces(coords, i, j, k, l, de_dphi, &mut forces);
        }
        (energy, forces)
    }
}

pub fn mmff_torsion_ctor(
    _sp: &Params,
    _tp: &[(&str, &Params)],
    frame: &Frame,
) -> Result<Box<dyn Potential>, String> {
    // Per-instance parameters: the MMFF typifier baked v1/v2/v3 onto each
    // dihedral (table → empirical). This kernel only reads the columns and
    // evaluates — no force-field-specific resolution lives here.
    let block = frame
        .get("dihedrals")
        .ok_or("mmff_torsion: missing \"dihedrals\"")?;
    let ic = block.get_uint("atomi").ok_or("missing atomi")?;
    let jc = block.get_uint("atomj").ok_or("missing atomj")?;
    let kc = block.get_uint("atomk").ok_or("missing atomk")?;
    let lc = block.get_uint("atoml").ok_or("missing atoml")?;
    let v1c = block
        .get_float("v1")
        .ok_or("mmff_torsion: missing \"v1\" column (typifier did not bake torsion params)")?;
    let v2c = block
        .get_float("v2")
        .ok_or("mmff_torsion: missing \"v2\" column (typifier did not bake torsion params)")?;
    let v3c = block
        .get_float("v3")
        .ok_or("mmff_torsion: missing \"v3\" column (typifier did not bake torsion params)")?;

    let n = ic.len();
    let (mut ai, mut aj, mut ak, mut al) = (
        Vec::with_capacity(n),
        Vec::with_capacity(n),
        Vec::with_capacity(n),
        Vec::with_capacity(n),
    );
    let (mut v1, mut v2, mut v3) = (
        Vec::with_capacity(n),
        Vec::with_capacity(n),
        Vec::with_capacity(n),
    );

    for idx in 0..n {
        ai.push(ic[idx] as usize);
        aj.push(jc[idx] as usize);
        ak.push(kc[idx] as usize);
        al.push(lc[idx] as usize);
        v1.push(v1c[idx] as F);
        v2.push(v2c[idx] as F);
        v3.push(v3c[idx] as F);
    }
    Ok(Box::new(MMFFTorsion {
        atom_i: ai,
        atom_j: aj,
        atom_k: ak,
        atom_l: al,
        v1,
        v2,
        v3,
    }))
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_mmff_torsion() {
        let pot = MMFFTorsion {
            atom_i: vec![0],
            atom_j: vec![1],
            atom_k: vec![2],
            atom_l: vec![3],
            v1: vec![0.0],
            v2: vec![0.0],
            v3: vec![0.3],
        };
        let coords: Vec<F> = vec![0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, -1.0, 0.0];
        let (e, forces) = pot.calc_energy_forces(&coords);
        assert!(e.is_finite());
        let fx: F = forces.iter().step_by(3).sum();
        let fy: F = forces.iter().skip(1).step_by(3).sum();
        let fz: F = forces.iter().skip(2).step_by(3).sum();
        assert!(
            (fx.abs() + fy.abs() + fz.abs()) < 0.1,
            "force sum too large"
        );
    }
}