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

//! MMFF94 out-of-plane bending: E = 0.5*143.9325*koop*chi^2 (Wilson angle)

use crate::ff::forcefield::Params;
use crate::ff::potential::Potential;
use crate::ff::potential::geometry::{cross3, dot3, mag3, sub3, validate_coords};
use molrs::store::frame::Frame;
use molrs::types::F;

use crate::ff::constants::MDYNE_A_TO_KCAL;

pub struct MMFFOutOfPlane {
    atom_i: Vec<usize>,
    atom_j: Vec<usize>,
    atom_k: Vec<usize>,
    atom_l: Vec<usize>,
    koop: Vec<F>,
}

impl Potential for MMFFOutOfPlane {
    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()];
        let conv = MDYNE_A_TO_KCAL as F;

        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 a = sub3(coords, i, coords, j); // ji
            let b = sub3(coords, k, coords, j); // jk
            let c = sub3(coords, l, coords, j); // jl
            let n = cross3(a, b);
            let n_mag = mag3(n);
            let c_mag = mag3(c);
            if n_mag < 1e-12 as F || c_mag < 1e-12 as F {
                continue;
            }

            let sin_chi = (dot3(c, n) / (c_mag * n_mag)).clamp(-1.0, 1.0);
            let chi = sin_chi.asin();
            let cos_chi = chi.cos();
            energy += 0.5 * conv * self.koop[idx] * chi * chi;

            if cos_chi.abs() < 1e-12 as F {
                continue;
            }
            let prefactor = conv * self.koop[idx] * chi / cos_chi;

            // d(sin chi)/d(r_l)
            let inv_cn = 1.0 / (c_mag * n_mag);
            for dim in 0..3 {
                let dsdc = n[dim] * inv_cn - sin_chi * c[dim] / (c_mag * c_mag);
                forces[l * 3 + dim] += -prefactor * dsdc;
                forces[j * 3 + dim] -= -prefactor * dsdc;
            }
            // d(sin chi)/d(r_i)  (a = r_i - r_j, affects n = a x b)
            let bxc = cross3(b, c);
            let bxn = cross3(b, n);
            for dim in 0..3 {
                let dsda = bxc[dim] * inv_cn - sin_chi * bxn[dim] / (n_mag * n_mag);
                forces[i * 3 + dim] += -prefactor * dsda;
                forces[j * 3 + dim] -= -prefactor * dsda;
            }
            // d(sin chi)/d(r_k)  (b = r_k - r_j, affects n = a x b)
            let axc = cross3(a, c);
            let axn = cross3(a, n);
            for dim in 0..3 {
                let dsdb = -axc[dim] * inv_cn + sin_chi * axn[dim] / (n_mag * n_mag);
                forces[k * 3 + dim] += -prefactor * dsdb;
                forces[j * 3 + dim] -= -prefactor * dsdb;
            }
        }
        (energy, forces)
    }
}

pub fn mmff_oop_ctor(
    _sp: &Params,
    _tp: &[(&str, &Params)],
    frame: &Frame,
) -> Result<Box<dyn Potential>, String> {
    // Per-instance parameters: the MMFF typifier baked koop onto each improper.
    // This kernel only reads the column and evaluates.
    let block = frame
        .get("impropers")
        .ok_or("mmff_oop: missing \"impropers\"")?;
    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 koopc = block
        .get_float("koop")
        .ok_or("mmff_oop: missing \"koop\" column (typifier did not bake oop params)")?;

    let n = ic.len();
    let (mut ai, mut aj, mut ak, mut al, mut koop) = (
        Vec::with_capacity(n),
        Vec::with_capacity(n),
        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);
        koop.push(koopc[idx] as F);
    }
    Ok(Box::new(MMFFOutOfPlane {
        atom_i: ai,
        atom_j: aj,
        atom_k: ak,
        atom_l: al,
        koop,
    }))
}