sci-form 0.15.2

High-performance 3D molecular conformer generation using ETKDG distance geometry
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

//! Debug: compare analytical and numerical UFF gradients on a test molecule
//! to verify gradient correctness.
//!
//! Category: debug

use nalgebra::DMatrix;
use sci_form::forcefield::{
    energy::calculate_total_energy, gradients::compute_analytical_gradient, FFParams,
};
use sci_form::graph::{Atom, Molecule};

fn main() {
    let mut mol = Molecule::new("test");
    // Central atom 0, neighbors 1, 2, 3
    let a0 = mol.add_atom(Atom::new(6, 0.0, 0.0, 0.0));
    let a1 = mol.add_atom(Atom::new(6, 1.0, 0.0, 0.0));
    let a2 = mol.add_atom(Atom::new(6, 0.0, 1.0, 0.0));
    let a3 = mol.add_atom(Atom::new(6, 0.02, 0.02, 0.05));

    // Add bonds so they are neighbors
    mol.add_bond(
        a0,
        a1,
        sci_form::graph::Bond {
            order: sci_form::graph::BondOrder::Single,
            stereo: sci_form::graph::BondStereo::None,
        },
    );
    mol.add_bond(
        a0,
        a2,
        sci_form::graph::Bond {
            order: sci_form::graph::BondOrder::Single,
            stereo: sci_form::graph::BondStereo::None,
        },
    );
    mol.add_bond(
        a0,
        a3,
        sci_form::graph::Bond {
            order: sci_form::graph::BondOrder::Single,
            stereo: sci_form::graph::BondStereo::None,
        },
    );

    // Planar coordinates with very small volume (~0.05)
    let mut coords = DMatrix::from_element(4, 3, 0.0);
    coords[(0, 0)] = 0.0;
    coords[(0, 1)] = 0.0;
    coords[(0, 2)] = 0.0;
    coords[(1, 0)] = 1.0;
    coords[(1, 1)] = 0.0;
    coords[(1, 2)] = 0.0;
    coords[(2, 0)] = 0.0;
    coords[(2, 1)] = 1.0;
    coords[(2, 2)] = 0.0;
    coords[(3, 0)] = 0.02;
    coords[(3, 1)] = 0.02;
    coords[(3, 2)] = 0.05;

    let params = FFParams {
        kb: 0.0,
        k_theta: 0.0,
        k_omega: 0.0,
        k_oop: 0.0,
        k_bounds: 0.0,
        k_chiral: 100.0,
        k_vdw: 0.0,
    };

    let n = 4;
    let bounds_matrix = DMatrix::from_element(n, n, 0.0);

    let e0 = calculate_total_energy(&coords, &mol, &params, &bounds_matrix);
    println!("Initial Total Energy: {}", e0);

    let analytical_grad = compute_analytical_gradient(&coords, &mol, &params, &bounds_matrix);

    let h = 1e-4;
    let mut numerical_grad = DMatrix::from_element(n, 3, 0.0);

    for i in 0..n {
        for dim in 0..3 {
            let mut c_plus = coords.clone();
            c_plus[(i, dim)] += h;
            let mut c_minus = coords.clone();
            c_minus[(i, dim)] -= h;

            let e_plus = calculate_total_energy(&c_plus, &mol, &params, &bounds_matrix);
            let e_minus = calculate_total_energy(&c_minus, &mol, &params, &bounds_matrix);

            numerical_grad[(i, dim)] = (e_plus - e_minus) / (2.0 * h);
        }
    }

    println!("Analytical vs Numerical Chiral Gradients:");
    let mut max_diff: f32 = 0.0;
    for i in 0..n {
        for dim in 0..3 {
            let a = analytical_grad[(i, dim)];
            let num = numerical_grad[(i, dim)];
            let diff = (a - num).abs();
            if diff > max_diff {
                max_diff = diff;
            }
            println!(
                "Atom {}, dim {}: An: {:.5}, Num: {:.5}, Diff: {:.5}",
                i, dim, a, num, diff
            );
        }
    }
    println!("Max diff: {}", max_diff);
}