Skip to main content

minimal/
minimal.rs

1//! A minimal example, demonstrating the most important syntax.
2
3use std::path::Path;
4
5use bio_files::{MmCif, Mol2};
6use dynamics::{
7    ComputationDevice, FfMolType, MdConfig, MdState, MolDynamics,
8    params::{FfParamSet, prepare_peptide_mmcif},
9};
10
11fn main() {
12    let dev = ComputationDevice::Cpu;
13    let param_set = FfParamSet::new_amber().unwrap();
14
15    let mut protein = MmCif::load(Path::new("1c8k.cif")).unwrap();
16    let mol = Mol2::load(Path::new("CPB.mol2")).unwrap();
17
18    // Add Hydrogens, force field type, and partial charge to atoms in the protein; these usually aren't
19    // included from RSCB PDB. You can also call `populate_hydrogens_dihedrals()`, and
20    // `populate_peptide_ff_and_q() separately. Add bonds.
21    let (_bonds, _dihedrals) = prepare_peptide_mmcif(
22        &mut protein,
23        &param_set.peptide_ff_q_map.as_ref().unwrap(),
24        7.0,
25    )
26    .unwrap();
27
28    let mols = vec![
29        MolDynamics::from_mol2(&mol, None),
30        MolDynamics {
31            ff_mol_type: FfMolType::Peptide,
32            atoms: protein.atoms.clone(),
33            static_: true,
34            ..Default::default()
35        },
36    ];
37
38    let (mut md, _) = MdState::new(&dev, &MdConfig::default(), &mols, &param_set).unwrap();
39
40    let n_steps = 100;
41    let dt = 0.002; // picoseconds.
42
43    for _ in 0..n_steps {
44        md.step(&dev, dt, None);
45    }
46
47    let snap = &md.snapshots[md.snapshots.len() - 1]; // A/R.
48    let energy = snap.energy_data.as_ref().unwrap();
49    println!(
50        "KE: {}, PE: {}, Atom posits:",
51        energy.energy_kinetic, energy.energy_potential
52    );
53    for posit in &snap.atom_posits {
54        println!("Posit: {posit}");
55        // Also keeps track of velocities, and solvent molecule positions/velocity
56    }
57
58    // Do something with snapshot data, like displaying atom positions in your UI.
59    // You can save to DCD file, and adjust the ratio they're saved at using the `MdConfig.snapshot_setup`
60    // field: See the example below.
61    for _snap in &md.snapshots {}
62}