Struct MdState 

pub struct MdState {

    pub cfg: MdConfig,
    pub atoms: Vec<AtomDynamics>,
    pub water: Vec<WaterMolOpc>,
    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>,
    pub alchemical: StateAlchemical,
    pub run_index: Option<usize>,
    /* private fields */
}

Fields

cfg: MdConfig

atoms: Vec<AtomDynamics>

water: Vec<WaterMolOpc>

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.

alchemical: StateAlchemical

run_index: Option<usize>

Index assigned at the start of each MD run. Trajectory files are named traj_N.dcd, traj_N.trr, etc. so that successive runs never overwrite each other. Chosen as the lowest N for which no such files exist yet. None until the first call to handle_ss_file_writes.

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, external_force: Option<Vec<Vec3>>, )

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.

External force allows injection of a specific force into the system. It’s indexed by atom.

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, None);
    }

    let snap = &md.snapshots[md.snapshots.len() - 1]; // A/R.
    let energy = snap.energy_data.as_ref().unwrap();
    println!(
        "KE: {}, PE: {}, Atom posits:",
        energy.energy_kinetic, energy.energy_potential
    );
    for posit in &snap.atom_posits {
        println!("Posit: {posit}");
        // Also keeps track of velocities, and solvent 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 {}
}

examples/lig_and_protein.rs (line 92)

pub fn build_dynamics(
    dev: &ComputationDevice,
    ligs: Vec<&mut Sdf>,
    peptide: &MmCif,
    param_set: &FfParamSet,
    cfg: &MdConfig,
    n_steps: u32,
    dt: f32,
) -> Result<MdState, ParamError> {
    println!("Setting up dynamics...");

    let mut mols = Vec::new();

    for lig in &ligs {
        mols.push(MolDynamics {
            ff_mol_type: FfMolType::SmallOrganic,
            atoms: lig.atoms.clone(),
            atom_posits: None,
            atom_init_velocities: None,
            bonds: lig.bonds.clone(),
            // Adj list is created atomatically, but can be passed from a cache.
            adjacency_list: None,
            static_: false,
            mol_specific_params: None,
            bonded_only: false,
            // Or: ..Default::default()
        });
    }

    // A demonstration of only selecting peptide atoms near a ligand.
    // We assume hetero atoms are ligands, solvent etc, and are not part of the protein.
    let atoms: Vec<_> = peptide
        .atoms
        .iter()
        .filter(|a| {
            let mut closest_dist = f64::MAX;
            for lig in &ligs {
                // You should probably use a special position set vs the native atom positions.
                for a in &lig.atoms {
                    let posit = a.posit;
                    let dist = (posit - a.posit).magnitude();
                    if dist < closest_dist {
                        closest_dist = dist;
                    }
                }
            }

            !a.hetero && closest_dist < STATIC_ATOM_DIST_THRESH
        })
        .map(|a| a.clone())
        .collect();

    let bonds = create_bonds(&atoms);

    mols.push(MolDynamics {
        ff_mol_type: FfMolType::Peptide,
        atoms,
        bonds,
        static_: true,
        ..Default::default()
    });

    // See also: `MolDynamics::from_sdf()`, `::from_mol2()`, and `::from_amber_geostd("CPB")`

    println!("Initializing MD state...");
    let (mut md_state, _) = MdState::new(dev, cfg, &mols, param_set)?;
    println!("Done.");

    let start = Instant::now();

    for _ in 0..n_steps {
        md_state.step(dev, dt, None);
    }

    let elapsed = start.elapsed();
    println!("MD complete in {:.2} s", elapsed.as_secs());

    change_snapshot(ligs, &md_state.snapshots[0]);

    Ok(md_state)
}

impl MdState

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

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

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

impl MdState

pub fn pack_atoms(&mut self)

impl MdState

pub fn flush_snapshot_queues(&mut self)

Flush any remaining snapshots in the DCD/TRR/XTC queues to disk. Call this at the end of a simulation run to ensure the last frames (those accumulated since the most recent TRAJ_FILE_SAVE_INTERVAL write) are not lost.

impl MdState

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

Use this to help initialize solvent molecules to realistic geometry of hydrogen bond networks, prior to the first proper simulation step. Runs MD on solvent 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 solvent 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 solvent is stable by the time the main sim starts.

impl MdState

pub fn shrink_cell_towards( &mut self, dev: &ComputationDevice, target_cell: SimBox, cfg: ShrinkingBoxCfg, ) -> bool

Scale each molecule’s centroid into the next smaller box while preserving its internal geometry. Water is treated as a rigid molecule as well.

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

Rebuild cell-dependent caches after a caller deliberately moves molecules.

pub fn redistribute_interleaved_opc_waters( &mut self, dev: &ComputationDevice, solvent_centers: &[Vec3F64], solvent_spacing: Vec3F64, ) -> bool

Redistribute the state’s existing OPC waters throughout a grid of co-packed solvent molecules. This avoids sparse explicit water counts occupying only the first region traversed by the template initializer.

impl MdState

pub fn initialize_velocities(&mut self, target_k: f32, zero_com_drift: bool)

Assign Maxwell-Boltzmann velocities at the requested temperature.

Solute atoms are sampled independently. OPC water is sampled as a rigid body, with a translational COM velocity and angular momentum, so the initial velocities respect the same rigid-water model used by integration.

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, external_force: Option<Vec<Vec3>>, )

Relaxes the molecules using a steepest-descent energy minimizer. 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 solvent 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 solvent molecules, as we have a pre-sim set up for them that runs prior to this.

pub fn minimize_energy_setup( &mut self, dev: &ComputationDevice, external_force: &Option<Vec<Vec3>>, ) -> (Vec<Vec3>, f32, f64, Vec<Vec3>, bool)

Separate, so can be called `separately by an application, e.g. if it needs to apply a new external force each step.

pub fn minimize_energy_cleanup( &mut self, dev: &ComputationDevice, prev_long_range: bool, initial_velocities: &[Vec3], )

See the note on minimize_energy_setup; this is broken out so it can be called separately by an application.

pub fn step_energy_min( &mut self, dev: &ComputationDevice, last_step: &mut [Vec3], alpha: &mut f32, e_prev: &mut f64, external_force: &Option<Vec<Vec3>>, ) -> bool

One iteration of energy minimization. Returns true if the energy is converged, indicating to abort further steps.

impl MdState

pub fn configure_alchemical_window( &mut self, dev: &ComputationDevice, mol_idx: usize, lambda: f64, ) -> Result<(), AlchemicalError>

Enable alchemical decoupling for one molecule at a fixed λ value. This is the entry point to enable an alchemical computation for an MD run. The application sequences the λ values by setting up MD runs, and calling this with the appropraite λ.

It validates the molecule index, stores the λ value, clears cached reciprocal data, and rebuilds non-bonded pairs so cross interactions with the selected molecule use alchemical LJ/Coulomb force handling.

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

Disable alchemical scaling and rebuild non-bonded pairs.

impl MdState

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

Also returns any explicit solvent molecules added. This may be needed by applications in order to create molecule sets for rendering the trajectories. These are placed, in the trajectory, after all solute atoms.

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, None);
    }

    let snap = &md.snapshots[md.snapshots.len() - 1]; // A/R.
    let energy = snap.energy_data.as_ref().unwrap();
    println!(
        "KE: {}, PE: {}, Atom posits:",
        energy.energy_kinetic, energy.energy_potential
    );
    for posit in &snap.atom_posits {
        println!("Posit: {posit}");
        // Also keeps track of velocities, and solvent 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 {}
}

examples/lig_and_protein.rs (line 86)

pub fn build_dynamics(
    dev: &ComputationDevice,
    ligs: Vec<&mut Sdf>,
    peptide: &MmCif,
    param_set: &FfParamSet,
    cfg: &MdConfig,
    n_steps: u32,
    dt: f32,
) -> Result<MdState, ParamError> {
    println!("Setting up dynamics...");

    let mut mols = Vec::new();

    for lig in &ligs {
        mols.push(MolDynamics {
            ff_mol_type: FfMolType::SmallOrganic,
            atoms: lig.atoms.clone(),
            atom_posits: None,
            atom_init_velocities: None,
            bonds: lig.bonds.clone(),
            // Adj list is created atomatically, but can be passed from a cache.
            adjacency_list: None,
            static_: false,
            mol_specific_params: None,
            bonded_only: false,
            // Or: ..Default::default()
        });
    }

    // A demonstration of only selecting peptide atoms near a ligand.
    // We assume hetero atoms are ligands, solvent etc, and are not part of the protein.
    let atoms: Vec<_> = peptide
        .atoms
        .iter()
        .filter(|a| {
            let mut closest_dist = f64::MAX;
            for lig in &ligs {
                // You should probably use a special position set vs the native atom positions.
                for a in &lig.atoms {
                    let posit = a.posit;
                    let dist = (posit - a.posit).magnitude();
                    if dist < closest_dist {
                        closest_dist = dist;
                    }
                }
            }

            !a.hetero && closest_dist < STATIC_ATOM_DIST_THRESH
        })
        .map(|a| a.clone())
        .collect();

    let bonds = create_bonds(&atoms);

    mols.push(MolDynamics {
        ff_mol_type: FfMolType::Peptide,
        atoms,
        bonds,
        static_: true,
        ..Default::default()
    });

    // See also: `MolDynamics::from_sdf()`, `::from_mol2()`, and `::from_amber_geostd("CPB")`

    println!("Initializing MD state...");
    let (mut md_state, _) = MdState::new(dev, cfg, &mols, param_set)?;
    println!("Done.");

    let start = Instant::now();

    for _ in 0..n_steps {
        md_state.step(dev, dt, None);
    }

    let elapsed = start.elapsed();
    println!("MD complete in {:.2} s", elapsed.as_secs());

    change_snapshot(ligs, &md_state.snapshots[0]);

    Ok(md_state)
}

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

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.