Crate dynamics 

See the Readme for a general overview, or Molchanica docs, MD section for more information about assumptions. Or see the examples folder on Github for how to use this in your application.

The textual information here is informal, and aimed at code maintenance; not library use.

This module contains high-level tools for running Newtonian molecular dynamics simulations.

Good article A summary on molecular dynamics

Amber Force Fields reference Small molucules using GAFF2 Amber RM 2025

To download .dat files (GAFF2), download Amber source (Option 2) here. Files are in dat -> leap -> parm

Base units: Å, ps (10^-12), Dalton (AMU), native charge units (derive from other base units; not a traditional named unit).

Amber: ff19SB for proteins, gaff2 for ligands. (Based on recommendations from https://ambermd.org/AmberModels.php).

We use the term “Non-bonded” interactions to refer to Coulomb, and Lennard Interactions, the latter of which is an approximation for both Van der Waals force and exclusion.

A broad list of components of this simulation:

• Water: Rigid OPC solvent molecules that have mutual non-bonded interactions with dynamic atoms and solvent
• Thermostat/barostat, with a way to specify temp, pressure, solvent density
• OPC solvent model
• Cell wrapping
• Velocity Verlet integration (Water and non-solvent)
• Amber parameters for mass, partial charge, VdW (via LJ), dihedral/improper, angle, bond len
• Optimizations for Coulomb: Ewald/SPME.
• Optimizations for LJ: Dist cutoff for now.
• Amber 1-2, 1-3 exclusions, and 1-4 scaling of covalently-bonded atoms.
• Rayon parallelization of non-bonded forces
• WIP SIMD and CUDA parallelization of non-bonded forces, depending on hardware availability. todo
• A thermostat and barostat
• An energy-measuring system.
• An integrated tool for inferring atom types, bonded-parameter overrides, and partial charges for arbitrary small organic molecules. (Similar to Amber’s Antechamber)

---

A timing test, using bond-stretching forces between two atoms only. Measure the period of oscillation for these atom combinations, e.g. using custom Mol2 files. c6-c6: 35fs (correct). os-os: 47fs nc-nc: 34fs hw-hw: 9fs Our measurements, 2025-08-04 c6-c6: 35fs os-os: 31fs nc-nc: 34fs (Correct) hw-hw: 6fs

---

We use traditional MD non-bonded terms to maintain geometry: Bond length, valence angle between 3 bonded atoms, dihedral angle between 4 bonded atoms (linear), and improper dihedral angle between each hub and 3 spokes. (E.g. at ring intersections). We also apply Coulomb force between atom-centered partial charges, and Lennard Jones potentials to simulate Van der Waals forces. These use spring-like forces to retain most geometry, while allowing for flexibility.

We use the OPC solvent model. (See water_opc.rs). For both maintaining the geometry of each solvent molecule, and for maintaining Hydrogen atom positions, we do not apply typical non-bonded interactions: We use SHAKE + RATTLE algorithms for these. In the case of solvent, it’s required for OPC compliance. For H, it allows us to maintain integrator stability with a greater timestep, e.g. 2fs instead of 1fs.

On f32 vs f64 floating point precision: f32 may be good enough for most things, and typical MD packages use mixed precision. Long-range electrostatics are a good candidate for using f64. Or, very long runs.

Note on performance: It appears that non-bonded forces dominate computation time. This is my observation, and it’s confirmed by an LLM. Both LJ and Coulomb take up most of the time; bonded forces are comparatively insignificant. Building neighbor lists are also significant. These are the areas we focus on for parallel computation (Thread pools, SIMD, CUDA)

Re-exports

pub use integrate::Integrator;

Modules

alchemical

Alchemical free-energy calculations for high-level solvation and LogP workflows. This computes a result after varying simulation interactions between solute and solvent, and measuring potential energy at each variation.

integrate

Contains integration code, including the primary time step.

minimize_energy

param_inference

Uses reasoning similar to AnteChamber’s to estimate force field parameters using DEF files. Specifically, we use DEF_GFF2 for Gaff2 force field names, and DEF_ABCG2 for FRCMOD bonded paraemter (generally dihedral) overrides. Reference source code

params

Data functionality for Forcefield params. Includes Amber parameters built in to binaries which use this library, and can load params for other sets as required.

partial_charge_inference

For inferring partial charge of small organic molecules using Amber’s GeoStd library as training data. Uses a neural net. Force field types must be assigned prior to running this, e.g. from the param_inference module.

snapshot

Related to storing snapshots (also known as trajectories) of MD runs.

Structs

AtomDynamics

A trimmed-down atom for use with molecular dynamics. Contains parameters for single-atom, but we use ParametersIndex for multi-atom parameters.

BarostatCfg

Dihedral

An amino acid in a protein structure, including all dihedral angles required to determine the conformation. Includes backbone and side chain dihedral angles. Doesn’t store coordinates, but coordinates can be generated using forward kinematics from the angles.

ForcesOnWaterMol

Per-solvent, per-site force accumulator. Used transiently when applying nonbonded forces. This is the force on each atom in the molecule.

MdConfig

This is the primary way of configurating an MD run. It’s passed at init, along with the molecule list and FF params.

MdOverrides

These are primarily used for debugging and testing, but may be used for specific scenarios as well, e.g. if wishing to speed up computations for real-time use by removing long range forces. These are not standard MD config parameters.

MdState

MolDynamics

Packages information required to perform dynamics on a Molecule. This is used to initialize the simulation with atoms and related; one or more of these is passed at init.

ParamError

Represents problems loading parameters. For example, if an atom is missing a force field type or partial charge, or has a force field type that hasn’t been loaded.

ShrinkingBoxCfg

Settings shared by template preparation and driven shrinking-box simulations.

ShrinkingBoxPackingCfg

Perhaps formally called gradual isotropic compression.

SimBox

This bounds the area where atoms are wrapped. For now at least, it is only used for solvent atoms. Its size and position should be such as to keep the system solvated. We may move it around during the sim.

WaterInitTemplate

For 3 and 4 point water models.

WaterMolOpc

Contains 4 atoms for each solvent molecules, at a given time step. Note that these are not independent, but are useful in our general MD APIs, for compatibility with non-solvent atoms.

Enums

ComMotionRemoval

GROMACS-style center-of-mass motion removal mode.

ComputationDevice

CustomSolventCount

FfMolType

This is used to assign the correct force field parameters to a molecule.

HydrogenConstraint

We use this variant in the configuration API. Deferrs to HydrogenConstraintInner for holding constraints.

SimBoxInit

Note: The shortest edge should be > 2(r_cutoff + r_skin), to prevent atoms from interacting with their own image in the real-space component.

Solvent

Used when configuring a MD Sim. We use OPC (rigid) water as a default, but can use custom solvents as well, from arbitrary molecules using standard MD forcefields.

SolventTemplateType

Contains variants of templates we have built into this library. These are included in the binary of applications which use this.

Constants

LANGEVIN_GAMMA_DEFAULT

LINCS_ITER_DEFAULT

LINCS_ORDER_DEFAULT

OCTANOL_WATER_TEMPLATE

PRESSURE_DEFAULT

SHAKE_TOL_DEFAULT

TAU_PRESSURE_DEFAULT

TAU_TEMP_DEFAULT

WATER_TEMPLATE_60A

Functions

compute_energy_snapshot

Set up with no solvent molecules or relaxation. Run one step to compute energies, then return the snapshot taken.

find_planar_posit

Given two positions, find the third planar position. Uses the central_posit_0 distance for the computed distances.

find_tetra_posit_final

Given 3 satellite atoms, find the 4th, in tetrahedral config.

find_tetra_posits

Given two satelite positions, find the third and fourth tetrahedral position. Uses the central, posit_0 distance for the computed distances.

make_octanol

Using PubChem data as a reference. Partial charges are computed using ORCA. We use this input: ! HF 6-31G* Opt TightSCF TightOpt RESP

make_water_mols_grid

Creates a regular lattice of water molecules. We use this as the first part of creating a solvent template. Use this, run a sim with thermostat and barostat, then store the result in a WaterInitTemplate. We can save and load this to disk as binary, or in .gro format.

merge_params

Add items from one parameter set to the other. If there are duplicates, the second set’s overrides the baseline.

pack_solvent_with_shrinking_box

A think wrapper over pack_solvent_with_shrinking_box_cfg, but with the default config.

pack_solvent_with_shrinking_box_cfg

populate_hydrogens_dihedrals

Adds hydrogens to a molecule, and populdates residue dihedral angles. This is useful in particular for mmCIF files from RCSB PDB, as they don’t have these. Uses Amber (or similar)-provided parameters as a guide.

water_mols_from_template

Create solvent molecules from a template, tiling it as many times as needed to fill the cell. Works for any cell size: smaller than, equal to, or larger than the template. Deconflcits with solute molecules, and adds the proper amount based on the free volume (Volume of the cell not taken up by solute).

water_mols_from_template_in_region

Create water molecules from a template inside a rectangular sub-region of a larger simulation cell.