Molecular Dynamics

A Python and Rust library for molecular dynamics. Compatible with Linux, Windows, and Mac. Uses CPU with threadpools and SIMD, or an Nvidia GPU.

It uses traditional forcefield-based molecular dynamics, and is inspired by Amber. It does not use quantum-mechanics, nor ab-initio methods.

It uses the Bio-Files dependency to load molecule and force-field files.

Please reference the API documentation for details on the functionality of each data structure and function. You may with to reference the Bio Files API docs as well.

Note: We currently only support saving and loading snapshots/trajectories in a custom format.

We recommend running this on GPU; it's much faster. This requires an Nvidia GPU Rtx3 series or newer, with nvidia drivers 580 or newer.

Note: The Python version does not yet use GPU for long-range forces. We would like to fix this, but are having trouble linking the cuFFT dependency.

Use of this library

This is intended for integration into a Rust or Python program, another Rust or Python library, or in small scripts that describe a workflow. When scripting, you will likely load molecule files directly, use integrated force fields or load them from file, and save results to a reporter format like DCD. If incorporating into an application, you might do more in memory using the data structures we provide.

Our API goal is to both provide a default terse syntax with safe defaults, and allow customization and flexibility that facilitates integration into bigger systems.

Goals

• Runs Newtonian MD algorithms accurately
• Easy to install, learn, and use
• Fast
• Easy integration into workflows, scripting, and applications

Installation

Python: pip install mol-dynamics

Rust: Add dynamics to Cargo.toml. Likely bio_files as well.

For a GUI application that uses this library, download the Daedalus molecule viewer This provides an easy-to-use way to set up the simulation, and play back trajectories.

Input topology

The simulation accepts sets of AtomicGeneric and BondGeneric. You can get these by loading molecular file formats (mmCIF, Mol2, SDF, etc) using the Bio Files library (biology-files in Python), or by creating them directly. See examples below and in the examples folder, and the docs links above; those are structs of plain data that can be built from from arbitrary input sources. For example, if you're building an application, you might use a more complicated Atom format; you can create a function that converts between yours, and AtomGeneric.

Small organic molecule note

There is some trouble with small organic molecule file formats. Some tools like openff-toolkit discourage or disallow use of Mol2 due to its lack of well-defined specification. They recommend SDF files instead. This is surprising, as SDF files also don't have a standard way to define partial charge or force field type. For example, OpenMM, and PubChem use different formats for partial charges. We view this with context: the alternatives also have problems or equal or greater severity to this. For example, PDBQT isn't intended for general use; XYZ is minimal; CML uses XML which is not convenient.

We suggest using whichever files are convenient, and converting them to in-memory data structure using tools you like, or ones we include. Mol2's ability to store force field type and partial charge in column data, having a published spec variant of Mol, and being adopted by Amber's Geostd library makes it a good choice.

Requirements by molecule type:

• Proteins/amino acids chains: No special requirements. Uses mmCif. (also known as PdbX). Force field type and partial charges are inferred automatically.
• Small organic molecules: Must have Force Field name, and partial charge populated on all atoms. Mol2 files from
• Amber's GeoStd library have these. SDF files from Drugbank are generally missing these fields. SDF files from PubChem usually includes partial charges, but not forcefield names. Molecule-specific parameters (e.g. from .frcmod files) may be required; these can also be [automatically] loaded from Amber Geostd.
• Lipids: TBD
• Nucleic acids: TBD
• Carbohydrates: TBD
• Lipids: TBD

Parameters

Integrates the following Amber parameters:

• Small organic molecules, e.g. ligands: General Amber Force Fields: GAFF2
• Protein and amino acids: FF19SB
• Nucleic acids: Amber OL3 and RNA libraries
• Lipids: lipids21
• Carbohydrates: GLYCAM_06j
• Water: Explicit, with OPC

The algorithm

This library uses a traditional MD workflow. We use the following components:

Integrators, thermostats, barostats

We provide a Velocity-Verlet integrator. It can be used with a Berendsen barostat, and either a CSVR/Bussi, or Langevin thermostat (Middle or traditional). These continuously update atom velocities (for molecules and solvents) to match target pressure and temperatures. The Langevin Middle thermostat is a good starting point.

Solvation

We use an explicit water solvation model: A 4-point rigid OPC model, with a point partial charge on each Hydrogen, and a M (or EP) point offset from the Oxygen. We use the SETTLE algorithm to maintain rigidity, while applying forces to each atom. Only the Oxygen atom carries a Lennard Jones(LJ) force.

Bonded forces

We use Amber-style spring-based forces to maintain covalent bonds. We maintain the following parameters:

• Bond length between each covalently-bonded atom pair
• Angle (sometimes called valence angle) between each 3-atom line of covalently-bonded atoms. (4atoms, 2 bonds)
• Dihedral (aka torsion) angles between each 4-atom line of covalently-bonded atoms. (4 atoms, 3 bonds). These usually have rotational symmetry of 2 or 3 values which are stable.
• Improper Dihedral angles between 4-atoms in a hub-and-spoke configuration. These, for example, maintain stability
• where rings meet other parts of the molecule, or other rings.

Non-bonded forces

These are Coulomb and Lennard Jones (LJ) interactions. These make up the large majority of computational effort. Coulomb forces represent electric forces occurring from dipoles and similar effects, or ions. We use atom-centered pre-computed partial charges for these. They occur within a molecule, between molecules, and between molecules and solvents.

We use a neighbors list (Sometimes called Verlet neighbors; not directly related to the Verlet integrator) to reduce computational effort. We use the SPME Ewald approximation to reduce computation time. This algorithm is suited for periodic boundary conditions, which we use for the solvent.

We use Amber's scaling and exclusion rules: LJ and Coulomb force is reduced between atoms separated by 1 and 2 covalent bonds, and skipped between atoms separated by 3 covalent bonds.

We have two modes of handling Hydrogen in bonded forces: The same as other atoms, and rigid, with position maintained using SHAKE and RATTLE algorithms. The latter allows for stability under higher timesteps. (e.g. 2ps)

Initial relaxation

We run a relaxation / energy-minimization function prior to starting each simulation. This adjusts atom positions to reduce the amount of energy that comes from initial conditions deviating from bonded parameters.

Floating point precision

Mixed precision: 32-bit floating points for most operations. We use 64-bit accumulators, and in thermostat and barostat computations.

How pH adjustment works

pH in proteins is represented by the protenation state of certain amino acids. In particular, His, Asp, Cys, Glu, and Lys are affected. These changes are affected in utility functions we provide that add Hydrogen atoms.

Saving results

Snapshots of results can be returned in memory, or saved to disk in DCD format.

More info

We plan to support carbohydrates and lipids later. If you're interested in these, please add a Github Issue.

These general parameters do not need to be loaded externally; they provide the information needed to perform MD with any amino acid sequence, and provide a baseline for dynamics of small organic molecules. You may wish to load frcmod data over these that have overrides for specific small molecules.

This program can automatically load ligands with Amber parameters, for the Amber Geostd set. This includes many common small organic molecules with force field parameters, and partial charges included. It can infer these from the protein loaded, or be queried by identifier.

You can load these molecules with parameters directly from the GUI by typing the identifier. If you load an SDF molecule, the program may be able to automatically update it using Amber parameters and partial charges.

For details on how dynamics using this parameterized approach works, see the Amber Reference Manual. Section 3 and 15 are of particular interest, regarding force field parameters.

Molecule-specific overrides to these general parameters can be loaded from .frcmod and .dat files. We delegate this to the bio files library.

We load partial charges for ligands from mol2, PDBQT etc files. Protein dynamics and water can be simulated using parameters built-in to the program (The Amber one above). Simulating ligands requires the loaded file (e.g. mol2) include partial charges. we recommend including ligand-specific override files as well, e.g. to load dihedral angles from .frcmod that aren't present in Gaff2.

You can load (and save) combined atom and forcefield data from Amber PRMTOP files; these combine these two data types into one file.

Use the code below, the Examples folder on Github, and the API documentation to learn how to use it. General workflow:

-Create a MdState struct with MdState::new(). This accepts a configuration, the molecules to simulate, and force field parameters.

Run a simulation step by calling MdState::step(). This accepts an enum which defines the computation devices (CPU/GPU), and the time step in picoseconds. This step can be called as required for your application. For example you can call it repeatedly in a loop, or as required, e.g. to not block a GUI, or for interactive MD.

Example use (Python):

from mol_dynamics import *

def setup_dynamics(mol: Mol2, protein: MmCif, param_set: FfParamSet, lig_specific: ForceFieldParams) -> MdState:
    """
    Set up dynamics between a small molecule we treat with full dynamics, and a rigid one 
    which acts on the system, but doesn't move.
    """

    # Or, consider using these terse helpers instead for small organic molecules.
    # MolDynamics.from_amber_geostd("CPB")  # Can use with a PubChem CID as well.
    # MolDynamics.from_mol2(mol, lig_specific)
    # MolDynamics.from_sdf(mol, lig_specific)

    mols = [
        MolDynamics(
            ff_mol_type=FfMolType.SmallOrganic,
            atoms=mol.atoms,
            # Pass a [Vec3] of starting atom positions. If absent,
            # will use the positions stored in atoms.
            atom_posits=None,
            atom_init_velocities=None,
            bonds=mol.bonds,
            # Pass your own from cache if you want, or it will build.
            adjacency_list=None,
            static_=False,
            # This is usually mandatory for small organic molecules. Provided, for example,
            # in Amber FRCMOD files. Overrides general params.
            mol_specific_params=lig_specific,
            bonded_only=False,
        ),
        MolDynamics(
            ff_mol_type=FfMolType.Peptide,
            atoms=protein.atoms,
            atom_posits=None,
            atom_init_velocities=None,
            bonds=[],  # Not required if static.
            adjacency_list=None,
            static_=True,
            mol_specific_params=None,
            bonded_only=False,
        ),
    ]

    return MdState(
        MdConfig(),
        mols,
        param_set,
    )


def main():
    mol = Mol2.load("CPB.mol2")
    protein = MmCif.load("1c8k.cif")

    param_set = FfParamSet.new_amber()
    lig_specific = ForceFieldParams.load_frcmod("CPB.frcmod")
    
    # Or, instead of loading atoms and mol-specific params separately:
    # mol, lig_specific = load_prmtop("my_mol.prmtop")
    
    # 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.
    _bonds, _dihedrals = prepare_peptide_mmcif(
        protein,
        param_set.peptide_ff_q_map,
        7.0,
    )
    # A variant of that function called `prepare_peptide` takes separate atom, residue, and chain
    # lists, for flexibility.
    
    md = setup_dynamics(mol, protein, param_set, lig_specific)
    
    n_steps = 100
    dt = 0.002  # picoseconds.
    
    for _ in range(n_steps):
        md.step(dt)
    
    snap = md.snapshots[len(md.snapshots) - 1]  # A/R.
    print(f"KE: {snap.energy_kinetic}, PE: {snap.energy_potential}, Atom posits:")
    for posit in snap.atom_posits:
        print(f"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:
        pass
        
        
main()

Example use (Rust):

use std::path::Path;

use bio_files::{MmCif, Mol2, md_params::ForceFieldParams};
use dynamics::{
    ComputationDevice, FfMolType, MdConfig, MdState, MolDynamics,
    params::{FfParamSet, prepare_peptide},
    populate_hydrogens_dihedrals,
};

/// Set up dynamics between a small molecule we treat with full dynamics, and a rigid one
/// which acts on the system, but doesn't move.
fn setup_dynamics(
    dev: &ComputationDevice,
    mol: &Mol2,
    protein: &MmCif,
    param_set: &FfParamSet,
    lig_specific: ForceFieldParams,
) -> MdState {
    
    // Or, consider using these terse helpers instead for small organic molecules.
    // MolDynamics::from_amber_geostd("CPB").unwrap();  // Can use with a PubChem CID as well.
    // MolDynamics::from_mol2(&mol, Some(lig_specific)).unwrap();
    // MolDynamics::from_sdf(&mol, Some(lig_specific)).unwrap();
    
    let mols = vec![
        MolDynamics {
            ff_mol_type: FfMolType::SmallOrganic,
            atoms: mol.atoms.clone(),
            // Pass a &[Vec3] of starting atom positions. If absent,
            // will use the positions stored in atoms.
            atom_posits: None,
            atom_init_velocities: None,
            bonds: mol.bonds.clone(),
            // Pass your own from cache if you want, or it will build.
            adjacency_list: None,
            static_: false,
            // This is usually mandatory for small organic molecules. Provided, for example,
            // in Amber FRCMOD files. Overrides general params.
            mol_specific_params: Some(lig_specific),
            bonded_only: false,
        },
        MolDynamics {
            ff_mol_type: FfMolType::Peptide,
            atoms: protein.atoms.clone(),
            atom_posits: None,
            atom_init_velocities: None,
            bonds: Vec::new(), // Not required if static.
            adjacency_list: None,
            static_: true,
            mol_specific_params: None,
            bonded_only: false,
        },
    ];

    MdState::new(dev, &MdConfig::default(), &mols, param_set).unwrap()
}

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();
    let mol_specific = ForceFieldParams::load_frcmod(Path::new("CPB.frcmod")).unwrap();

    // Or, instead of loading atoms and mol-specific params separately:
    // let (mol, lig_specific) = load_prmtop("my_mol.prmtop");
    
    // Or, if you have a small molecule available in Amber Geostd, load it remotely:
    // let data = bio_apis::amber_geostd::load_mol_files("CPB");
    // let mol = Mol2::new(&data.mol2);
    // let mol_specific = ForceFieldParams::from_frcmod(&data.frcmod);
    
    // 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();
    // A variant of that function called `prepare_peptide` takes separate atom, residue, and chain
    // lists, for flexibility.

    let mut md = setup_dynamics(&dev, &mol, &protein, &param_set, mol_specific);

    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 {}
}

Example of loading your own parameter files:

    param_paths = ParamGeneralPaths(
        peptide="parm19.dat",
        peptide_mod="frcmod.ff19SB",
        peptide_ff_q="amino19.lib",
        peptide_ff_q_c="aminoct12.lib",
        peptide_ff_q_n=None,
        small_organic="gaff2.dat",
        lipids="lipid21.dat",
    )
    
    param_set = FfParamSet(param_paths)

    let param_paths = ParamGeneralPaths {
        peptide: Some(&Path::new("parm19.dat")),
        peptide_mod: Some(&Path::new("frcmod.ff19SB")),
        peptide_ff_q: Some(&Path::new("amino19.lib")),
        peptide_ff_q_c: Some(&Path::new("aminoct12.lib")),
        peptide_ff_q_n: Some(&Path::new("aminont12.lib")),
        small_organic: Some(&Path::new("gaff2.dat")),
        lipids: Some(&Path::new("lipid21.dat")),
        ..default()
    };
    
    let param_set = FfParamSet::new(&param_paths);

An overview of configuration parameters. You may wish to (Rust) use a baseline of the Default implementation, then override specific fields you wish to change.

let cfg = MdConfig {
    // Defaults to Langevin middle.
    integrator: dynamics::Integrator::VelocityVerlet,
    // If enabled, zero the drift in center of mass of the system.
    zero_com_drift: true,
    // Kelvin. Defaults to 310 K.
    temp_target: 310.,
    // Bar (Pa/100). Defaults to 1 bar.
    pressure_target: 1.,
    // Allows constraining Hydrogens to be rigid with their bonded atom, using SHAKE and RATTLE
    // algorithms. This allows for higher time steps.
    hydrogen_constraint: dynamics::HydrogenConstraint::Fixed,
    // Deafults to in-memory, every step
    snapshot_handlers: vec![
        SnapshotHandler {
            save_type: SaveType::Memory,
            ratio: 1,
        },
        SnapshotHandler {
            save_type: SaveType::Dcd(PathBuf::from("output.dcd")),
            ratio: 10,
        },
    ],
    sim_box: SimBoxInit::Pad(10.),
    // Or sim_box: SimBoxInit::Fixed((Vec3::new(-10., -10., -10.), Vec3::new(10., 10., 10.)),
};

Python config syntax:

cfg = MdConfig() // Initializes with defaults.

cfg.integrator = dynamics.Integrator.VelocityVerlet
cfg.temp_target = 310.
# etc

You can run md_state.computation_time() after running to get a breakdown of how long each computation component took to run, averaged per step.

Using with GPU

We use the Cudarc library for GPU (CUDA) integration. In the python binding, it should be transparent. We've exposed a slightly lower level API in rust, where you use setup a Stream and modules with Cudarc in your application, and pass them to the library.

Rust setup example with Cudarc. Pass dev, defined below, to the step function.

let ctx = CudaContext::new(0).unwrap();
let stream = ctx.default_stream();

let dev = ComputationDevice::Gpu(GpuModules(stream);

To use with an Nvidia GPU, enable the cuda feature in Cargo.toml. The library will generate PTX instructions as a publicly exposed string. Set up your application to use it from dynamics::PTX. It requires CUDA 13 support, which requires Nvidia driver version 580 or higher.

On unflattening trajactory data

If you passed multiple molecules, these will be flattened during runtime, and in snapshots. You need to unflatten them if placing back into their original data structures.

Why this when OpenMM exists?

This library exists as part of a larger Rust biology infrastructure effort. It's not possible to use OpenMM there due to the language barrier. This library currently only has a limited subset of the functionality of OpenMM. It's unfortunate that, as a society, we've embraced a model of computing replete with obstacles. In this case, the major one is the one placed between programming languages.

While going around this obstacle, we attempt to jump over others, to make molecular dynamics more accessible. This includes operating systems, software distribution, and user experience. We hope that this is easier to install and use than OpenMM; it can be used on any Operating system, and any Python version >= 3.10, installable using pip or cargo.

This library is intended to just work. OpenMM does not natively work with molecules from online databases like RCSB PDB, PubChem, and Drugbank. It doesn't work with Amber GeoStd Mol2 files. OpenMM itself is easy to install with Pip, but the additional libraries it requires to load molecules and force fields are higher-friction. Getting a functional OpenMM configuration for a given system involves work which we hope to eschew.

Compiling from source

It requires these Amber parameter files to be present under the project's resources folder at compile time. These are available in Amber tools. Download, unpack, then copy these files from dat/leap/parm and dat/leap/lib:

• amino19.lib
• aminoct12.lib
• aminont12.lib
• parm19.dat
• frcmod.ff19SB
• gaff2.dat
• ff-nucleic-OL24.lib
• ff-nucleic-OL24.frcmod
• RNA.lib

We provide a copy of these files for convenience; this is a much smaller download than the entire Amber package, and prevents needing to locate the specific files. Unpack, and place these under resources prior to compiling.

To build the Python library wheel, from the python subdirectory, run maturin build. You can load the library locally for testing, once built, by running pip install .

Eratta

• GPU operations are slower than they should, as we're passing all data between CPU and GPU each time step.
• CPU SIMD unsupported

References

• Amber forcefields
• Amber reference manual
• Ewald Summation/SPME
• OPC water model
• Tripos Mol2 format