#![allow(
unused_imports,
unused_variables,
dead_code,
clippy::unnecessary_cast,
clippy::needless_range_loop,
clippy::manual_repeat_n,
clippy::manual_str_repeat,
clippy::manual_is_multiple_of,
clippy::redundant_field_names,
clippy::useless_vec,
clippy::single_range_in_vec_init
)]
use nalgebra::{DMatrix, Vector3};
use petgraph::visit::EdgeRef;
use rand::SeedableRng;
use serde::Deserialize;
#[derive(Deserialize)]
struct OracleAtom {
element: u8,
x: f32,
y: f32,
z: f32,
formal_charge: i8,
hybridization: String,
}
#[derive(Deserialize)]
struct OracleBond {
start: usize,
end: usize,
order: String,
}
#[derive(Deserialize)]
struct OracleMolecule {
smiles: String,
atoms: Vec<OracleAtom>,
bonds: Vec<OracleBond>,
}
fn build_mol(mol: &OracleMolecule) -> sci_form::graph::Molecule {
let mut our_mol = sci_form::graph::Molecule::new(&mol.smiles);
let mut node_indices = Vec::new();
for atom in &mol.atoms {
let hybridization = match atom.hybridization.as_str() {
"SP" => sci_form::graph::Hybridization::SP,
"SP2" => sci_form::graph::Hybridization::SP2,
"SP3" => sci_form::graph::Hybridization::SP3,
"SP3D" => sci_form::graph::Hybridization::SP3D,
"SP3D2" => sci_form::graph::Hybridization::SP3D2,
_ => sci_form::graph::Hybridization::Unknown,
};
let new_atom = sci_form::graph::Atom {
element: atom.element,
position: nalgebra::Vector3::zeros(),
charge: 0.0,
formal_charge: atom.formal_charge,
hybridization,
chiral_tag: sci_form::graph::ChiralType::Unspecified,
explicit_h: if atom.element == 1 || atom.element == 0 {
1
} else {
0
},
};
node_indices.push(our_mol.add_atom(new_atom));
}
for bond in &mol.bonds {
let order = match bond.order.as_str() {
"DOUBLE" => sci_form::graph::BondOrder::Double,
"TRIPLE" => sci_form::graph::BondOrder::Triple,
"AROMATIC" => sci_form::graph::BondOrder::Aromatic,
_ => sci_form::graph::BondOrder::Single,
};
our_mol.add_bond(
node_indices[bond.start],
node_indices[bond.end],
sci_form::graph::Bond {
order,
stereo: sci_form::graph::BondStereo::None,
},
);
}
our_mol
}
fn energy_breakdown(
coords: &DMatrix<f32>,
mol: &sci_form::graph::Molecule,
params: &sci_form::forcefield::FFParams,
bounds_matrix: &DMatrix<f64>,
) -> (f32, f32, f32, f32, f32, f32) {
let n = mol.graph.node_count();
let mut e_bonds = 0.0f32;
let mut e_angles = 0.0f32;
let mut e_bounds = 0.0f32;
let mut e_oop = 0.0f32;
let mut e_torsion_uff = 0.0f32;
let mut e_torsion_etkdg = 0.0f32;
for edge in mol.graph.edge_references() {
let idx1 = edge.source().index();
let idx2 = edge.target().index();
let p1 = Vector3::new(coords[(idx1, 0)], coords[(idx1, 1)], coords[(idx1, 2)]);
let p2 = Vector3::new(coords[(idx2, 0)], coords[(idx2, 1)], coords[(idx2, 2)]);
let r_eq = sci_form::distgeom::get_bond_length(mol, edge.source(), edge.target()) as f32;
e_bonds += sci_form::forcefield::energy::bond_stretch_energy(&p1, &p2, params.kb, r_eq);
}
for i in 0..n {
let ni = petgraph::graph::NodeIndex::new(i);
let nbs: Vec<_> = mol.graph.neighbors(ni).collect();
for j in 0..nbs.len() {
for k in (j + 1)..nbs.len() {
let p1 = Vector3::new(
coords[(nbs[j].index(), 0)],
coords[(nbs[j].index(), 1)],
coords[(nbs[j].index(), 2)],
);
let pc = Vector3::new(coords[(i, 0)], coords[(i, 1)], coords[(i, 2)]);
let p2 = Vector3::new(
coords[(nbs[k].index(), 0)],
coords[(nbs[k].index(), 1)],
coords[(nbs[k].index(), 2)],
);
let ideal =
sci_form::graph::get_corrected_ideal_angle(mol, ni, nbs[j], nbs[k]) as f32;
e_angles += sci_form::forcefield::energy::angle_bend_energy(
&p1,
&pc,
&p2,
params.k_theta,
ideal,
);
}
}
}
for i in 0..n {
for j in (i + 1)..n {
let p1 = Vector3::new(coords[(i, 0)], coords[(i, 1)], coords[(i, 2)]);
let p2 = Vector3::new(coords[(j, 0)], coords[(j, 1)], coords[(j, 2)]);
let upper = bounds_matrix[(i, j)] as f32;
let lower = bounds_matrix[(j, i)] as f32;
e_bounds += sci_form::forcefield::energy::bounds_energy(
&p1,
&p2,
lower,
upper,
params.k_bounds,
);
}
}
if params.k_oop.abs() > 1e-8 {
for i in 0..n {
let ni = petgraph::graph::NodeIndex::new(i);
if mol.graph[ni].hybridization != sci_form::graph::Hybridization::SP2 {
continue;
}
let nbs: Vec<_> = mol.graph.neighbors(ni).collect();
if nbs.len() != 3 {
continue;
}
let pc = Vector3::new(coords[(i, 0)], coords[(i, 1)], coords[(i, 2)]);
let p1 = Vector3::new(
coords[(nbs[0].index(), 0)],
coords[(nbs[0].index(), 1)],
coords[(nbs[0].index(), 2)],
);
let p2 = Vector3::new(
coords[(nbs[1].index(), 0)],
coords[(nbs[1].index(), 1)],
coords[(nbs[1].index(), 2)],
);
let p3 = Vector3::new(
coords[(nbs[2].index(), 0)],
coords[(nbs[2].index(), 1)],
coords[(nbs[2].index(), 2)],
);
let v1 = p1 - pc;
let v2 = p2 - pc;
let v3 = p3 - pc;
let vol = v1.dot(&v2.cross(&v3));
e_oop += params.k_oop * vol * vol;
}
}
if params.k_omega.abs() > 1e-8 && n >= 4 {
for edge in mol.graph.edge_references() {
let u = edge.source();
let v = edge.target();
let hyb_u = mol.graph[u].hybridization;
let hyb_v = mol.graph[v].hybridization;
if hyb_u == sci_form::graph::Hybridization::SP
|| hyb_v == sci_form::graph::Hybridization::SP
{
continue;
}
let (n_fold, gamma, weight) =
sci_form::forcefield::energy::torsion_params(hyb_u, hyb_v);
let neighbors_u: Vec<_> = mol.graph.neighbors(u).filter(|&x| x != v).collect();
let neighbors_v: Vec<_> = mol.graph.neighbors(v).filter(|&x| x != u).collect();
for &nu in &neighbors_u {
for &nv in &neighbors_v {
let p1 = Vector3::new(
coords[(nu.index(), 0)],
coords[(nu.index(), 1)],
coords[(nu.index(), 2)],
);
let p2 = Vector3::new(
coords[(u.index(), 0)],
coords[(u.index(), 1)],
coords[(u.index(), 2)],
);
let p3 = Vector3::new(
coords[(v.index(), 0)],
coords[(v.index(), 1)],
coords[(v.index(), 2)],
);
let p4 = Vector3::new(
coords[(nv.index(), 0)],
coords[(nv.index(), 1)],
coords[(nv.index(), 2)],
);
e_torsion_uff += sci_form::forcefield::energy::torsional_energy(
&p1,
&p2,
&p3,
&p4,
params.k_omega * weight,
n_fold,
gamma,
);
}
}
}
}
if n >= 4 {
for edge in mol.graph.edge_references() {
let u = edge.source();
let v = edge.target();
if sci_form::graph::min_path_excluding2(mol, u, v, u, v, 7).is_some() {
continue;
}
let m6 =
sci_form::forcefield::etkdg_lite::infer_etkdg_parameters(mol, u.index(), v.index());
if m6.v.iter().all(|&x| x.abs() < 1e-6) {
continue;
}
let neighbors_u: Vec<_> = mol.graph.neighbors(u).filter(|&x| x != v).collect();
let neighbors_v: Vec<_> = mol.graph.neighbors(v).filter(|&x| x != u).collect();
if neighbors_u.is_empty() || neighbors_v.is_empty() {
continue;
}
let nu = neighbors_u[0];
let nv = neighbors_v[0];
let p1 = Vector3::new(
coords[(nu.index(), 0)],
coords[(nu.index(), 1)],
coords[(nu.index(), 2)],
);
let p2 = Vector3::new(
coords[(u.index(), 0)],
coords[(u.index(), 1)],
coords[(u.index(), 2)],
);
let p3 = Vector3::new(
coords[(v.index(), 0)],
coords[(v.index(), 1)],
coords[(v.index(), 2)],
);
let p4 = Vector3::new(
coords[(nv.index(), 0)],
coords[(nv.index(), 1)],
coords[(nv.index(), 2)],
);
e_torsion_etkdg +=
sci_form::forcefield::etkdg_lite::calc_torsion_energy_m6(&p1, &p2, &p3, &p4, &m6);
}
}
(
e_bonds,
e_angles,
e_bounds,
e_oop,
e_torsion_uff,
e_torsion_etkdg,
)
}
fn gradient_norms(
coords: &DMatrix<f32>,
mol: &sci_form::graph::Molecule,
params: &sci_form::forcefield::FFParams,
bounds_matrix: &DMatrix<f64>,
) -> (f32, f32, f32, f32, f32, f32) {
let n = mol.graph.node_count();
let bond_params = sci_form::forcefield::FFParams {
kb: params.kb,
k_theta: 0.0,
k_omega: 0.0,
k_oop: 0.0,
k_bounds: 0.0,
k_chiral: 0.0,
k_vdw: 0.0,
};
let g_bond = sci_form::forcefield::gradients::compute_analytical_gradient(
coords,
mol,
&bond_params,
bounds_matrix,
);
let angle_params = sci_form::forcefield::FFParams {
kb: 0.0,
k_theta: params.k_theta,
k_omega: 0.0,
k_oop: 0.0,
k_bounds: 0.0,
k_chiral: 0.0,
k_vdw: 0.0,
};
let g_angle = sci_form::forcefield::gradients::compute_analytical_gradient(
coords,
mol,
&angle_params,
bounds_matrix,
);
let bounds_params = sci_form::forcefield::FFParams {
kb: 0.0,
k_theta: 0.0,
k_omega: 0.0,
k_oop: 0.0,
k_bounds: params.k_bounds,
k_chiral: 0.0,
k_vdw: 0.0,
};
let g_bounds = sci_form::forcefield::gradients::compute_analytical_gradient(
coords,
mol,
&bounds_params,
bounds_matrix,
);
let oop_params = sci_form::forcefield::FFParams {
kb: 0.0,
k_theta: 0.0,
k_omega: 0.0,
k_oop: params.k_oop,
k_bounds: 0.0,
k_chiral: 0.0,
k_vdw: 0.0,
};
let g_oop = sci_form::forcefield::gradients::compute_analytical_gradient(
coords,
mol,
&oop_params,
bounds_matrix,
);
let torsion_params = sci_form::forcefield::FFParams {
kb: 0.0,
k_theta: 0.0,
k_omega: params.k_omega,
k_oop: 0.0,
k_bounds: 0.0,
k_chiral: 0.0,
k_vdw: 0.0,
};
let g_torsion = sci_form::forcefield::gradients::compute_analytical_gradient(
coords,
mol,
&torsion_params,
bounds_matrix,
);
let g_full = sci_form::forcefield::gradients::compute_analytical_gradient(
coords,
mol,
params,
bounds_matrix,
);
let g_etkdg = &g_full - &g_bond - &g_angle - &g_bounds - &g_oop - &g_torsion;
(
g_bond.norm(),
g_angle.norm(),
g_bounds.norm(),
g_oop.norm(),
g_torsion.norm(),
g_etkdg.norm(),
)
}
#[test]
fn debug_ff_energy_breakdown() {
let data = std::fs::read_to_string("tests/fixtures/reference_coords_500.json").unwrap();
let mut molecules: Vec<OracleMolecule> = serde_json::from_str(&data).unwrap();
use rand::seq::SliceRandom;
let mut rng = rand::rngs::StdRng::seed_from_u64(123);
molecules.shuffle(&mut rng);
molecules.truncate(100);
let params = sci_form::forcefield::FFParams {
kb: 300.0,
k_theta: 200.0,
k_omega: 10.0,
k_oop: 20.0,
k_bounds: 10.0,
k_chiral: 0.0,
k_vdw: 0.0,
};
let test_indices = vec![62, 48, 28, 92, 0, 10];
for &mol_idx in &test_indices {
if mol_idx >= molecules.len() {
continue;
}
let mol_data = &molecules[mol_idx];
let our_mol = build_mol(mol_data);
let n = our_mol.graph.node_count();
let bounds = sci_form::distgeom::calculate_bounds_matrix(&our_mol);
let smoothed = sci_form::distgeom::smooth_bounds_matrix(bounds);
let mut rng = rand::rngs::StdRng::seed_from_u64(13);
let dists = sci_form::distgeom::pick_random_distances(&mut rng, &smoothed);
let metric = sci_form::distgeom::compute_metric_matrix(&dists);
let mut coords4d = sci_form::distgeom::generate_nd_coordinates(&mut rng, &metric, 4);
sci_form::forcefield::minimize_bounds_lbfgs_ex(
&mut coords4d,
&smoothed,
&[],
400,
1e-4,
5.0,
0.1,
);
sci_form::forcefield::minimize_bounds_lbfgs_ex(
&mut coords4d,
&smoothed,
&[],
200,
1e-4,
5.0,
1.0,
);
let coords_before = coords4d.columns(0, 3).into_owned();
println!("\n=== Mol {} ({}) n={} ===", mol_idx, &mol_data.smiles, n);
let (e_b, e_a, e_bnd, e_oop, e_t_uff, e_t_etkdg) =
energy_breakdown(&coords_before, &our_mol, ¶ms, &smoothed);
println!("BEFORE FF:");
println!(" Bonds: {:.4}", e_b);
println!(" Angles: {:.4}", e_a);
println!(" Bounds: {:.4}", e_bnd);
println!(" OOP: {:.4}", e_oop);
println!(" Tors UFF: {:.4}", e_t_uff);
println!(" Tors ETKDG: {:.4}", e_t_etkdg);
println!(
" TOTAL: {:.4}",
e_b + e_a + e_bnd + e_oop + e_t_uff + e_t_etkdg
);
let (gn_b, gn_a, gn_bnd, gn_oop, gn_t_uff, gn_t_etkdg) =
gradient_norms(&coords_before, &our_mol, ¶ms, &smoothed);
println!("Gradient norms:");
println!(" Bonds: {:.4}", gn_b);
println!(" Angles: {:.4}", gn_a);
println!(" Bounds: {:.4}", gn_bnd);
println!(" OOP: {:.4}", gn_oop);
println!(" Tors UFF: {:.4}", gn_t_uff);
println!(" Tors ETKDG: {:.4}", gn_t_etkdg);
let g_full = sci_form::forcefield::gradients::compute_analytical_gradient(
&coords_before,
&our_mol,
¶ms,
&smoothed,
);
println!(" Full norm: {:.4}", g_full.norm());
println!(" Full max: {:.4}", g_full.abs().max());
let coords_after =
sci_form::forcefield::minimize_etkdg_3d(&our_mol, &coords_before, &smoothed, 200, 1e-4);
if mol_idx == 62 {
let ff = sci_form::forcefield::build_etkdg_3d_ff(
&our_mol,
&coords_before.map(|v| v as f64),
&smoothed,
);
let anal_g = sci_form::forcefield::etkdg_3d_gradient(&coords_before, &our_mol, &ff);
let e0 = sci_form::forcefield::etkdg_3d_energy(&coords_before, &our_mol, &ff);
let eps = 1e-4f32;
println!(" COMPONENT GRADIENT CHECK (Mol 62):");
let mut g_dist = nalgebra::DMatrix::zeros(n, 3);
for c in ff
.dist_12
.iter()
.chain(ff.dist_13.iter())
.chain(ff.dist_long.iter())
{
let p1 = nalgebra::Vector3::new(
coords_before[(c.i, 0)],
coords_before[(c.i, 1)],
coords_before[(c.i, 2)],
);
let p2 = nalgebra::Vector3::new(
coords_before[(c.j, 0)],
coords_before[(c.j, 1)],
coords_before[(c.j, 2)],
);
sci_form::forcefield::gradients::analytical_grad_distance_constraint(
&p1,
&p2,
c.min_len as f32,
c.max_len as f32,
c.k as f32,
&mut g_dist,
c.i,
c.j,
);
}
for atom in 0..2 {
for d in 0..3 {
let mut cp = coords_before.clone();
let mut cm = coords_before.clone();
cp[(atom, d)] += eps;
cm[(atom, d)] -= eps;
let mut ep = 0.0f32;
let mut em = 0.0f32;
for c in ff
.dist_12
.iter()
.chain(ff.dist_13.iter())
.chain(ff.dist_long.iter())
{
let p1p = nalgebra::Vector3::new(cp[(c.i, 0)], cp[(c.i, 1)], cp[(c.i, 2)]);
let p2p = nalgebra::Vector3::new(cp[(c.j, 0)], cp[(c.j, 1)], cp[(c.j, 2)]);
ep += sci_form::forcefield::energy::distance_constraint_energy(
&p1p,
&p2p,
c.min_len as f32,
c.max_len as f32,
c.k as f32,
);
let p1m = nalgebra::Vector3::new(cm[(c.i, 0)], cm[(c.i, 1)], cm[(c.i, 2)]);
let p2m = nalgebra::Vector3::new(cm[(c.j, 0)], cm[(c.j, 1)], cm[(c.j, 2)]);
em += sci_form::forcefield::energy::distance_constraint_energy(
&p1m,
&p2m,
c.min_len as f32,
c.max_len as f32,
c.k as f32,
);
}
let num_g = (ep - em) / (2.0 * eps);
let anal = g_dist[(atom, d)];
let axis = ["x", "y", "z"][d];
println!(
" DIST atom={} {}: anal={:.4} num={:.4} diff={:.4}",
atom,
axis,
anal,
num_g,
(anal - num_g).abs()
);
}
}
println!(" --- SINGLE TORSION TEST ---");
'outer: for edge in our_mol.graph.edge_references() {
let u = edge.source();
let v = edge.target();
let hyb_u = our_mol.graph[u].hybridization;
let hyb_v = our_mol.graph[v].hybridization;
if hyb_u == sci_form::graph::Hybridization::SP
|| hyb_v == sci_form::graph::Hybridization::SP
{
continue;
}
let (n_fold, gamma, weight) =
sci_form::forcefield::energy::torsion_params(hyb_u, hyb_v);
let neighbors_u: Vec<_> = our_mol.graph.neighbors(u).filter(|&x| x != v).collect();
let neighbors_v: Vec<_> = our_mol.graph.neighbors(v).filter(|&x| x != u).collect();
if neighbors_u.is_empty() || neighbors_v.is_empty() {
continue;
}
let nu = neighbors_u[0];
let nv = neighbors_v[0];
let vv = 10.0 * weight;
println!(
" Torsion: {}-{}-{}-{} hyb={:?}-{:?} n={} gamma={:.2} V={:.2}",
nu.index(),
u.index(),
v.index(),
nv.index(),
hyb_u,
hyb_v,
n_fold,
gamma,
vv
);
let mk = |coords: &nalgebra::DMatrix<f32>, idx: usize| {
nalgebra::Vector3::new(coords[(idx, 0)], coords[(idx, 1)], coords[(idx, 2)])
};
let mut single_grad = nalgebra::DMatrix::zeros(n, 3);
sci_form::forcefield::gradients::analytical_grad_torsion(
&mk(&coords_before, nu.index()),
&mk(&coords_before, u.index()),
&mk(&coords_before, v.index()),
&mk(&coords_before, nv.index()),
vv,
n_fold,
gamma,
&mut single_grad,
nu.index(),
u.index(),
v.index(),
nv.index(),
);
let indices = [nu.index(), u.index(), v.index(), nv.index()];
for &idx in &indices {
for d in 0..3 {
let mut cp = coords_before.clone();
let mut cm = coords_before.clone();
cp[(idx, d)] += eps;
cm[(idx, d)] -= eps;
let ep = sci_form::forcefield::energy::torsional_energy(
&mk(&cp, nu.index()),
&mk(&cp, u.index()),
&mk(&cp, v.index()),
&mk(&cp, nv.index()),
vv,
n_fold,
gamma,
);
let em = sci_form::forcefield::energy::torsional_energy(
&mk(&cm, nu.index()),
&mk(&cm, u.index()),
&mk(&cm, v.index()),
&mk(&cm, nv.index()),
vv,
n_fold,
gamma,
);
let num_g = (ep - em) / (2.0 * eps);
let anal = single_grad[(idx, d)];
let axis = ["x", "y", "z"][d];
let ok = if num_g.abs() > 0.01 {
(anal / num_g - 1.0).abs() < 0.1
} else {
(anal - num_g).abs() < 0.1
};
if !ok {
println!(
" FAIL atom={} {} anal={:.6} num={:.6} ratio={:.4}",
idx,
axis,
anal,
num_g,
if num_g.abs() > 1e-6 {
anal / num_g
} else {
0.0
}
);
}
}
}
println!(" Single torsion check done");
break 'outer;
}
println!(" --- ETKDG M6 ---");
let no_params = sci_form::forcefield::FFParams {
kb: 0.0,
k_theta: 0.0,
k_omega: 0.0,
k_oop: 0.0,
k_bounds: 0.0,
k_chiral: 0.0,
k_vdw: 0.0,
};
let g_etkdg = sci_form::forcefield::gradients::compute_analytical_gradient(
&coords_before,
&our_mol,
&no_params,
&smoothed,
);
for atom in 0..2 {
for d in 0..3 {
let mut cp = coords_before.clone();
let mut cm = coords_before.clone();
cp[(atom, d)] += eps;
cm[(atom, d)] -= eps;
let mut ep = 0.0f32;
let mut em = 0.0f32;
for edge in our_mol.graph.edge_references() {
let u = edge.source();
let v = edge.target();
if sci_form::graph::min_path_excluding2(&our_mol, u, v, u, v, 7).is_some() {
continue;
}
let m6 = sci_form::forcefield::etkdg_lite::infer_etkdg_parameters(
&our_mol,
u.index(),
v.index(),
);
if m6.v.iter().all(|&x| x.abs() < 1e-6) {
continue;
}
let neighbors_u: Vec<_> =
our_mol.graph.neighbors(u).filter(|&x| x != v).collect();
let neighbors_v: Vec<_> =
our_mol.graph.neighbors(v).filter(|&x| x != u).collect();
if neighbors_u.is_empty() || neighbors_v.is_empty() {
continue;
}
let nu = neighbors_u[0];
let nv = neighbors_v[0];
let mk = |coords: &nalgebra::DMatrix<f32>, idx: usize| {
nalgebra::Vector3::new(
coords[(idx, 0)],
coords[(idx, 1)],
coords[(idx, 2)],
)
};
ep += sci_form::forcefield::etkdg_lite::calc_torsion_energy_m6(
&mk(&cp, nu.index()),
&mk(&cp, u.index()),
&mk(&cp, v.index()),
&mk(&cp, nv.index()),
&m6,
);
em += sci_form::forcefield::etkdg_lite::calc_torsion_energy_m6(
&mk(&cm, nu.index()),
&mk(&cm, u.index()),
&mk(&cm, v.index()),
&mk(&cm, nv.index()),
&m6,
);
}
let num_g = (ep - em) / (2.0 * eps);
let anal = g_etkdg[(atom, d)];
let axis = ["x", "y", "z"][d];
println!(
" M6 atom={} {}: anal={:.4} num={:.4} diff={:.4}",
atom,
axis,
anal,
num_g,
(anal - num_g).abs()
);
}
}
println!(" --- OOP ---");
let oop_params = sci_form::forcefield::FFParams {
kb: 0.0,
k_theta: 0.0,
k_omega: 0.0,
k_oop: 10.0,
k_bounds: 0.0,
k_chiral: 0.0,
k_vdw: 0.0,
};
let g_oop = sci_form::forcefield::gradients::compute_analytical_gradient(
&coords_before,
&our_mol,
&oop_params,
&smoothed,
);
for atom in 0..2 {
for d in 0..3 {
let mut cp = coords_before.clone();
let mut cm = coords_before.clone();
cp[(atom, d)] += eps;
cm[(atom, d)] -= eps;
let mut ep = 0.0f32;
let mut em = 0.0f32;
for i in 0..n {
let ni = petgraph::graph::NodeIndex::new(i);
if our_mol.graph[ni].hybridization != sci_form::graph::Hybridization::SP2 {
continue;
}
let nbs: Vec<_> = our_mol.graph.neighbors(ni).collect();
if nbs.len() != 3 {
continue;
}
let mk = |coords: &nalgebra::DMatrix<f32>, idx: usize| {
nalgebra::Vector3::new(
coords[(idx, 0)],
coords[(idx, 1)],
coords[(idx, 2)],
)
};
let (pcp, p1p, p2p, p3p) = (
mk(&cp, i),
mk(&cp, nbs[0].index()),
mk(&cp, nbs[1].index()),
mk(&cp, nbs[2].index()),
);
let v1p = p1p - pcp;
let v2p = p2p - pcp;
let v3p = p3p - pcp;
let volp = v1p.dot(&v2p.cross(&v3p));
ep += 10.0 * volp * volp;
let (pcm, p1m, p2m, p3m) = (
mk(&cm, i),
mk(&cm, nbs[0].index()),
mk(&cm, nbs[1].index()),
mk(&cm, nbs[2].index()),
);
let v1m = p1m - pcm;
let v2m = p2m - pcm;
let v3m = p3m - pcm;
let volm = v1m.dot(&v2m.cross(&v3m));
em += 10.0 * volm * volm;
}
let num_g = (ep - em) / (2.0 * eps);
let anal = g_oop[(atom, d)];
let axis = ["x", "y", "z"][d];
println!(
" OOP atom={} {}: anal={:.4} num={:.4} diff={:.4}",
atom,
axis,
anal,
num_g,
(anal - num_g).abs()
);
}
}
}
let (e_b2, e_a2, e_bnd2, e_oop2, e_t_uff2, e_t_etkdg2) =
energy_breakdown(&coords_after, &our_mol, ¶ms, &smoothed);
println!("AFTER FF:");
println!(" Bonds: {:.4} (delta {:.4})", e_b2, e_b2 - e_b);
println!(" Angles: {:.4} (delta {:.4})", e_a2, e_a2 - e_a);
println!(" Bounds: {:.4} (delta {:.4})", e_bnd2, e_bnd2 - e_bnd);
println!(" OOP: {:.4} (delta {:.4})", e_oop2, e_oop2 - e_oop);
println!(
" Tors UFF: {:.4} (delta {:.4})",
e_t_uff2,
e_t_uff2 - e_t_uff
);
println!(
" Tors ETKDG: {:.4} (delta {:.4})",
e_t_etkdg2,
e_t_etkdg2 - e_t_etkdg
);
println!(
" TOTAL: {:.4} (delta {:.4})",
e_b2 + e_a2 + e_bnd2 + e_oop2 + e_t_uff2 + e_t_etkdg2,
(e_b2 + e_a2 + e_bnd2 + e_oop2 + e_t_uff2 + e_t_etkdg2)
- (e_b + e_a + e_bnd + e_oop + e_t_uff + e_t_etkdg)
);
let diff = &coords_after - &coords_before;
let mut max_atom_move = 0.0f32;
let mut total_move = 0.0f32;
for i in 0..n {
let d = (diff[(i, 0)] * diff[(i, 0)]
+ diff[(i, 1)] * diff[(i, 1)]
+ diff[(i, 2)] * diff[(i, 2)])
.sqrt();
max_atom_move = max_atom_move.max(d);
total_move += d;
}
println!("Coord changes:");
println!(" Max atom move: {:.6} Å", max_atom_move);
println!(" Avg atom move: {:.6} Å", total_move / n as f32);
}
}