Quick Start

MDCraft was developed to be an all-encompassing tool that helps simplify the entire research workflow, from setting up and running molecular dynamics simulations to analyzing, modeling, and plotting data from the trajectories for publication in scientific journals.

Of special note, MDCraft offers unprecedented flexibility in how data analysis can be carried out. It not only has high-level analysis classes in the analysis submodule with carefully thought-out selections of keyword-only arguments to control exactly what is calculated, but also provides direct low-level access to the optimized algorithms used in the analysis classes in the algorithm submodule for more seasoned programmers.

The following section will largely focus on the analysis, openmm, and plot submodules since they are the most likely to be used by end users.

openmm submodule

The openmm submodule contains helper functions for setting up OpenMM simulations. Notably, it enables coarse-grained simulations with reduced units in a simulation toolkit largely meant for atomistic simulations with real units.

The following code snippet sets up and simulates a simple charged Lennard-Jones fluid without the need for any external topology or force field files.

import sys

from mdcraft.openmm.pair import lj_coul, ljts
from mdcraft.openmm.reporter import NetCDFReporter
from mdcraft.openmm.system import register_particles
from mdcraft.openmm.topology import create_atoms
from mdcraft.openmm.unit import get_lj_scale_factors
import numpy as np
import openmm
from openmm import app, unit

# Define constants and parameters
m = 39.948 * unit.amu
sigma = 3.405 * unit.angstrom
epsilon = 119.8 * unit.kelvin * unit.BOLTZMANN_CONSTANT_kB
N = 10_000
rho_reduced = 0.8
T = 300 * unit.kelvin
varepsilon_r = 78
dt_reduced = 0.0025
every = 100
timesteps = 2_000

# Get scale factors for reduced Lennard-Jones units
scales = get_lj_scale_factors({
    "energy": epsilon,
    "length": sigma,
    "mass": m
})

# Determine system dimensions
dimensions = (N * scales["length"] ** 3 / rho_reduced) ** (1 / 3) * np.ones(3)

# Initialize simulation system and topology
system = openmm.System()
system.setDefaultPeriodicBoxVectors(*(dimensions * np.diag(np.ones(3))))
topology = app.Topology()
topology.setUnitCellDimensions(dimensions)

# Set up excluded volume (Lennard-Jones) and electrostatic (Coulomb) pair potentials
cutoff = 2.5 * scales["length"]
pair_lj_cut = ljts(cutoff, shift=False)
pair_coul = lj_coul(cutoff)
system.addForce(pair_lj_cut)
system.addForce(pair_coul)

# Register particles to pair potentials
for q, name, element in zip(
        (-1, 1),
        ("ANI", "CAT"), 
        (app.Element.getBySymbol(e) for e in ("Cl", "Na"))  # arbitrary
    ):
    register_particles(
        system, 
        topology, 
        N // 2, 
        scales["mass"],
        element=element, 
        name=name, 
        nbforce=pair_coul, 
        charge=q / np.sqrt(varepsilon_r), 
        cnbforces={pair_lj_cut: (scales["length"], scales["molar_energy"])}
    )

# Generate initial particle positions
positions = create_atoms(dimensions, N)

# Set up simulation
dt = dt_reduced * scales["time"]
friction = 1 / scales["time"]
platform = openmm.Platform.getPlatformByName("CPU")
integrator = openmm.LangevinIntegrator(T, friction, dt)
simulation = app.Simulation(topology, system, integrator, platform)
context = simulation.context
context.setPositions(positions)

# Minimize energy
simulation.minimizeEnergy()

# Initialize velocities using Maxwell-Boltzmann distribution
context.setVelocitiesToTemperature(T)

# Write topology file
with open("topology.cif", "w") as f:
    app.PDBxFile.writeFile(
        topology, 
        context.getState(getPositions=True).getPositions(asNumpy=True), 
        f, 
        keepIds=True
    )

# Register topology and state data reporters
simulation.reporters.append(NetCDFReporter("trajectory.nc", every))
simulation.reporters.append(
    app.StateDataReporter(
        sys.stdout, 
        reportInterval=every, 
        step=True, 
        temperature=True,
        volume=True, 
        potentialEnergy=True, 
        kineticEnergy=True,
        totalEnergy=True, 
        remainingTime=True, 
        speed=True,
        totalSteps=timesteps
    )
)

# Run simulation
simulation.step(timesteps)
Hide code cell output 
#"Step","Potential Energy (kJ/mole)","Kinetic Energy (kJ/mole)","Total Energy (kJ/mole)","Temperature (K)","Box Volume (nm^3)","Speed (ns/day)","Time Remaining"
100,-75601.6852557223,26927.301295265785,-48674.383960456515,215.90732944844424,493.47068906249956,0,--
200,-74193.81594809999,30262.741981328138,-43931.07396677185,242.6514165429806,493.47068906249956,7.76,1:48
300,-73117.93772773554,32243.894130384106,-40874.04359735143,258.53660552064184,493.47068906249956,7.74,1:42
400,-72233.83446452131,33913.88150719029,-38319.952957331014,271.9268265006488,493.47068906249956,7.77,1:35
500,-71442.52074629178,34642.327800210915,-36800.19294608086,277.76762324623104,493.47068906249956,7.76,1:30
600,-70737.0710529994,35143.88899084772,-35593.182062151674,281.7892196192926,493.47068906249956,7.76,1:24
700,-70758.194443823,36045.05277412124,-34713.141669701756,289.0148923188605,493.47068906249956,7.77,1:17
800,-70202.8354512464,36083.57705376152,-34119.25839748488,289.32378604143787,493.47068906249956,7.77,1:11
900,-70204.36327199747,36573.54617995029,-33630.81709204718,293.25243542176787,493.47068906249956,7.77,1:05
1000,-70046.13616483023,36400.44373608599,-33645.69242874424,291.86447284928715,493.47068906249956,7.77,0:59
1100,-69866.24208708217,36317.88150849092,-33548.360578591244,291.202475945925,493.47068906249956,7.76,0:53
1200,-70383.65860280504,37171.54663495196,-33212.111967853074,298.04729695774614,493.47068906249956,7.77,0:47
1300,-69957.16718398889,36410.26030133744,-33546.90688265145,291.943183610702,493.47068906249956,7.77,0:41
1400,-69956.69164840986,37119.81874988217,-32836.87289852769,297.6325346538312,493.47068906249956,7.77,0:35
1500,-69828.29763404389,37063.09316721193,-32765.204466831958,297.17770002590316,493.47068906249956,7.77,0:29
1600,-69574.79441599657,37210.0596371545,-32364.734778842067,298.3561002560571,493.47068906249956,7.76,0:24
1700,-69821.9706814598,37364.527530686944,-32457.44315077286,299.5946480782995,493.47068906249956,7.76,0:18
1800,-69833.66501844038,37497.877817126115,-32335.78720131427,300.66387161134355,493.47068906249956,7.76,0:12
1900,-69632.24554152896,37434.56784013212,-32197.67770139684,300.156242275959,493.47068906249956,7.76,0:06
2000,-69766.6486352097,37055.78557621977,-32710.863058989926,297.11910661401544,493.47068906249956,7.76,0:00

/home/docs/checkouts/readthedocs.org/user_builds/mdcraft/conda/v1.3.2/lib/python3.12/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

analysis submodule

Now that we have the topology.cif and trajectory.nc files from the simulation, we can perform data analysis to get structural and dynamic properties that we are interested in, like density profiles, structure factors, and self-diffusion coefficients.

As an illustrative example, the following code snippet first reads the topology and trajectory files using MDAnalysis and then calculates the radial distribution functions for the unique species pairs.

from itertools import combinations_with_replacement

import MDAnalysis as mda
from mdcraft.analysis.structure import RadialDistributionFunction

# Load topology and trajectory
universe = mda.Universe(app.PDBxFile("topology.cif"), "trajectory.nc")

# Select groups containing particles of the same species
groups = [universe.select_atoms(f"name {name}") for name in ("ANI", "CAT")]

# Determine unique species pairs
pairs = list(combinations_with_replacement(range(len(groups)), 2))

# Compute radial distribution functions
rdfs = []
for i, j in pairs:
    rdfs.append(
        RadialDistributionFunction(
            groups[i], 
            groups[j], 
            1_000, 
            (0, 5 * scales["length"].value_in_unit(unit.angstrom)),
            exclusion=(1, 1)
        ).run(start=timesteps // (2 * every))
    )
Hide code cell output 
0%|          | 0/10 [00:00<?, ?it/s]
10%|█         | 1/10 [00:02<00:18,  2.02s/it]
20%|██        | 2/10 [00:03<00:14,  1.86s/it]
30%|███       | 3/10 [00:05<00:12,  1.81s/it]
40%|████      | 4/10 [00:07<00:10,  1.79s/it]
50%|█████     | 5/10 [00:09<00:08,  1.77s/it]
60%|██████    | 6/10 [00:10<00:07,  1.76s/it]
70%|███████   | 7/10 [00:12<00:05,  1.76s/it]
80%|████████  | 8/10 [00:14<00:03,  1.76s/it]
90%|█████████ | 9/10 [00:16<00:01,  1.75s/it]
100%|██████████| 10/10 [00:17<00:00,  1.75s/it]
100%|██████████| 10/10 [00:17<00:00,  1.78s/it]
0%|          | 0/10 [00:00<?, ?it/s]
10%|█         | 1/10 [00:01<00:15,  1.75s/it]
20%|██        | 2/10 [00:03<00:13,  1.75s/it]
30%|███       | 3/10 [00:05<00:12,  1.75s/it]
40%|████      | 4/10 [00:06<00:10,  1.74s/it]
50%|█████     | 5/10 [00:08<00:08,  1.74s/it]
60%|██████    | 6/10 [00:10<00:06,  1.74s/it]
70%|███████   | 7/10 [00:12<00:05,  1.74s/it]
80%|████████  | 8/10 [00:13<00:03,  1.73s/it]
90%|█████████ | 9/10 [00:15<00:01,  1.74s/it]
100%|██████████| 10/10 [00:17<00:00,  1.73s/it]
100%|██████████| 10/10 [00:17<00:00,  1.74s/it]
0%|          | 0/10 [00:00<?, ?it/s]
10%|█         | 1/10 [00:01<00:15,  1.74s/it]
20%|██        | 2/10 [00:03<00:13,  1.74s/it]
30%|███       | 3/10 [00:05<00:12,  1.74s/it]
40%|████      | 4/10 [00:06<00:10,  1.74s/it]
50%|█████     | 5/10 [00:08<00:08,  1.74s/it]
60%|██████    | 6/10 [00:10<00:06,  1.74s/it]
70%|███████   | 7/10 [00:12<00:05,  1.73s/it]
80%|████████  | 8/10 [00:13<00:03,  1.74s/it]
90%|█████████ | 9/10 [00:15<00:01,  1.74s/it]
100%|██████████| 10/10 [00:17<00:00,  1.74s/it]
100%|██████████| 10/10 [00:17<00:00,  1.74s/it]

plot submodule

Finally, we can visualize the results from the data analysis in clean, aesthetic, and publication-ready figures using Matplotlib.

The following code snippet plots the radial distribution functions computed in the previous step according to the ACS journal guidelines.

import matplotlib.pyplot as plt
from mdcraft.plot.rcparam import update

# Update Matplotlib rcParams to adhere to the ACS journal guidelines
update("acs", font_scaling=5 / 4, size_scaling=3 / 2)

# Plot and display radial distribution functions
labels = [f"$g_{{{chr(43 + 2 * i)}{chr(43 + 2 * j)}}}(r)$" for i, j in pairs]
_, ax = plt.subplots()
for rdf, label in zip(rdfs, labels):
    ax.plot(
        rdf.results.bins * unit.angstrom / scales["length"], 
        rdf.results.rdf,
        label=label
    )
ax.legend()
ax.set_xlabel("$r/\sigma$")
ax.set_ylabel("$g_{ij}(r)$")
ax.text(-0.2, 0.959, " ", transform=ax.transAxes)
plt.show()

<>:17: SyntaxWarning: invalid escape sequence '\s'
<>:17: SyntaxWarning: invalid escape sequence '\s'
/tmp/ipykernel_3134/318846508.py:17: SyntaxWarning: invalid escape sequence '\s'
  ax.set_xlabel("$r/\sigma$")
Matplotlib is building the font cache; this may take a moment.
../../_images/2bb3cbb38c6d467a7c5f0fc8ebcf520a75d34517ff12599f041b75ecf792aa03.png