Struct MdState 

pub struct MdState {

    pub cfg: MdConfig,
    pub atoms: Vec<AtomDynamics>,
    pub water: Vec<WaterMol>,
    pub adjacency_list: Vec<Vec<usize>>,
    pub time: f64,
    pub step_count: usize,
    pub snapshots: Vec<Snapshot>,
    pub cell: SimBox,
    pub neighbors_nb: NeighborsNb,
    pub water_pme_sites_forces: Vec<[Vec3; 3]>,
    pub kinetic_energy: f64,
    pub potential_energy: f64,
    pub potential_energy_nonbonded: f64,
    pub potential_energy_bonded: f64,
    pub potential_energy_between_mols: Vec<f64>,
    pub neighbor_rebuild_count: usize,
    pub computation_time: ComputationTimeSums,
    pub mol_start_indices: Vec<usize>,
    /* private fields */
}

Fields

cfg: MdConfig

atoms: Vec<AtomDynamics>

water: Vec<WaterMol>

adjacency_list: Vec<Vec<usize>>

Note: We don’t use bond structs once the simulation is set up; the adjacency list is the source of this.

time: f64

Current simulation time, in picoseconds.

step_count: usize

snapshots: Vec<Snapshot>

These are the snapshots we keep in memory, accumulating.

cell: SimBox

neighbors_nb: NeighborsNb

water_pme_sites_forces: Vec<[Vec3; 3]>

kinetic_energy: f64

kcal/mol

potential_energy: f64

potential_energy_nonbonded: f64

A newer, simpler approach for energy between molecules, compared to potential_energy_between_mols. This is simply the potential energy from non-bonded interactions, and excludes that from bonded.

potential_energy_bonded: f64

E.g. energy in covalent bonds, as modelled as oscillators.

potential_energy_between_mols: Vec<f64>

Every so many snapshots, write these to file, then clear from memory. Used to track which molecule each atom is associated with in our flattened structures. This is the potential energy between every pair of molecules.

neighbor_rebuild_count: usize

computation_time: ComputationTimeSums

mol_start_indices: Vec<usize>

Used to track which molecule each atom is associated with in our flattened structures.

Implementations

impl MdState

pub fn apply_angle_bending_forces(&mut self)

This maintains bond angles between sets of three atoms as they should be from hybridization. It reflects this hybridization, steric clashes, and partial double-bond character. This identifies deviations from the ideal angle, calculates restoring torque, and applies forces based on this to push the atoms back into their ideal positions in the molecule.

Valence angles, which are the angle formed by two adjacent bonds ba et bc in a same molecule; a valence angle tends to maintain constant the anglê abc. A valence angle is thus concerned by the positions of three atoms.

impl MdState

pub fn step(&mut self, dev: &ComputationDevice, dt: f32)

Perform one integration step. This is the entry point for running the simulation. One step of length dt is in picoseconds (10^-12), with typical values of 0.001, or 0.002ps (1 or 2fs). This method orchestrates the dynamics at each time step. Uses a Verlet Velocity base, with different thermostat approaches depending on configuration.

Examples found in repository

examples/minimal.rs (line 44)

fn main() {
    let dev = ComputationDevice::Cpu;
    let param_set = FfParamSet::new_amber().unwrap();

    let mut protein = MmCif::load(Path::new("1c8k.cif")).unwrap();
    let mol = Mol2::load(Path::new("CPB.mol2")).unwrap();

    // Add Hydrogens, force field type, and partial charge to atoms in the protein; these usually aren't
    // included from RSCB PDB. You can also call `populate_hydrogens_dihedrals()`, and
    // `populate_peptide_ff_and_q() separately. Add bonds.
    let (_bonds, _dihedrals) = prepare_peptide_mmcif(
        &mut protein,
        &param_set.peptide_ff_q_map.as_ref().unwrap(),
        7.0,
    )
    .unwrap();

    let mols = vec![
        MolDynamics::from_mol2(&mol, None),
        MolDynamics {
            ff_mol_type: FfMolType::Peptide,
            atoms: protein.atoms.clone(),
            static_: true,
            ..Default::default()
        },
    ];

    let mut md = MdState::new(&dev, &MdConfig::default(), &mols, &param_set).unwrap();

    let n_steps = 100;
    let dt = 0.002; // picoseconds.

    for _ in 0..n_steps {
        md.step(&dev, dt);
    }

    let snap = &md.snapshots[md.snapshots.len() - 1]; // A/R.
    println!(
        "KE: {}, PE: {}, Atom posits:",
        snap.energy_kinetic, snap.energy_potential
    );
    for posit in &snap.atom_posits {
        println!("Posit: {posit}");
        // Also keeps track of velocities, and water molecule positions/velocity
    }

    // Do something with snapshot data, like displaying atom positions in your UI.
    // You can save to DCD file, and adjust the ratio they're saved at using the `MdConfig.snapshot_setup`
    // field: See the example below.
    for snap in &md.snapshots {}
}

impl MdState

pub fn apply_nonbonded_forces(&mut self, dev: &ComputationDevice)

Run the appropriate force-computation function to get force on non-water atoms, force on water atoms, and virial sum for the barostat. Uses GPU if available.

Applies Coulomb and Van der Waals (Lennard-Jones) forces on non-water atoms, in place. We use the MD-standard [S]PME approach to handle approximated Coulomb forces. This function applies forces from non-water, and water sources.

impl MdState

pub fn pack_atoms(&mut self)

impl MdState

pub fn md_on_water_only(&mut self, dev: &ComputationDevice)

Use this to help initialize water molecules to realistic geometry of hydrogen bond networks, prior to the first proper simulation step. Runs MD on water only. Make sure to only run this after state is properly initialized, e.g. towards the end of init; not immediately after populating waters.

This will result in an immediate energy bump as water positions settle from their grid into position. As they settle, the thermostat will bring the velocities down to set the target temp. This sim should run long enough to the water is stable by the time the main sim starts.

impl MdState

pub fn zero_linear_momentum(&mut self)

Remove center-of-mass drift. This can help stabilize system energy. We perform the sums here as f64.

pub fn zero_angular_momentum(&mut self)

Remove rigid-body rotation. Computes ω from I ω = L about the atoms’ COM, then sets v’ = v - ω × (r - r_cm).

impl MdState

pub fn minimize_energy(&mut self, dev: &ComputationDevice, max_iters: usize)

Relaxes the molecules. Use this at the start of the simulation to control kinetic energy that arrises from differences between atom positions, and bonded parameters. It can also be called externally. It also stabilizes the water molecules, so that their hydrogen bond structure is correct at initialization.

Uses flexible bonds to hydrogen. (Not Shake/Rattle constraints)

We don’t apply this to water molecules, as we have a pre-sim set up for them that runs prior to this.

impl MdState

pub fn new( dev: &ComputationDevice, cfg: &MdConfig, mols: &[MolDynamics], param_set: &FfParamSet, ) -> Result<Self, ParamError>

Examples found in repository

examples/minimal.rs (line 38)

fn main() {
    let dev = ComputationDevice::Cpu;
    let param_set = FfParamSet::new_amber().unwrap();

    let mut protein = MmCif::load(Path::new("1c8k.cif")).unwrap();
    let mol = Mol2::load(Path::new("CPB.mol2")).unwrap();

    // Add Hydrogens, force field type, and partial charge to atoms in the protein; these usually aren't
    // included from RSCB PDB. You can also call `populate_hydrogens_dihedrals()`, and
    // `populate_peptide_ff_and_q() separately. Add bonds.
    let (_bonds, _dihedrals) = prepare_peptide_mmcif(
        &mut protein,
        &param_set.peptide_ff_q_map.as_ref().unwrap(),
        7.0,
    )
    .unwrap();

    let mols = vec![
        MolDynamics::from_mol2(&mol, None),
        MolDynamics {
            ff_mol_type: FfMolType::Peptide,
            atoms: protein.atoms.clone(),
            static_: true,
            ..Default::default()
        },
    ];

    let mut md = MdState::new(&dev, &MdConfig::default(), &mols, &param_set).unwrap();

    let n_steps = 100;
    let dt = 0.002; // picoseconds.

    for _ in 0..n_steps {
        md.step(&dev, dt);
    }

    let snap = &md.snapshots[md.snapshots.len() - 1]; // A/R.
    println!(
        "KE: {}, PE: {}, Atom posits:",
        snap.energy_kinetic, snap.energy_potential
    );
    for posit in &snap.atom_posits {
        println!("Posit: {posit}");
        // Also keeps track of velocities, and water molecule positions/velocity
    }

    // Do something with snapshot data, like displaying atom positions in your UI.
    // You can save to DCD file, and adjust the ratio they're saved at using the `MdConfig.snapshot_setup`
    // field: See the example below.
    for snap in &md.snapshots {}
}

pub fn computation_time(&self) -> Result<ComputationTime>

pub fn save_snapshots_to_file( &self, path: &Path, ratio: usize, ) -> Result<(), Error>

For calling by the applicastion. Saves in-memory snapshots to file, e.g. DCD, XTC, or MDT.

Trait Implementations

impl Default for MdState

fn default() -> MdState

Returns the “default value” for a type.

impl Display for MdState

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter.