Solution-Phase Molecular Dynamics

Goal

Set up and run molecular dynamics (MD) simulations of molecules in explicit solvent. This skill covers three stages: (1) building a solvation box with Packmol, (2) running NPT/NVT MD using MLIPs, and (3) analyzing the trajectory for radial distribution functions (RDFs), coordination numbers, density convergence, and mean-square displacement (MSD).

[!IMPORTANT] This skill bridges gas-phase chem-* skills and condensed-phase mat-* skills by providing workflows for solvation dynamics, liquid structure characterization, and dissolution studies.

1. Prerequisites

Packmol binary must be installed and on PATH in the base-agent environment.
RDKit must be available in the base-agent environment (for SMILES → 3D geometry).
An MLIP backend must be available via MCP tools (MACE, MatGL, or FairChem).

2. MLIP Selection

Refer to the foundation-potentials skill for model selection.

[!NOTE]

Organic solvents: Use MACE-MH-1 with omol head, or UMA with omol task.

Aqueous inorganic systems: Use MACE-MH-1 with omat_pbe head, or MatGL/CHGNet.

Mixed organic-inorganic: Use UMA which handles both.

3. Workflow

Step 1: Build Solvation Box

Use the box-building script to create a solvated system with Packmol:

# Env: base-agent
# Pure solvent box (64 water molecules)
python .agents/skills/chem-solution-md/scripts/build_solvation_box.py \
    --solvent water \
    --num_solvent 64 \
    --output_dir research/my_folder/solvation_box

# Solute in solvent (NaCl in 64 water molecules)
python .agents/skills/chem-solution-md/scripts/build_solvation_box.py \
    --solute_smiles "[Na+].[Cl-]" \
    --solvent water \
    --num_solvent 64 \
    --output_dir research/my_folder/solvation_box

Key Parameters:

| Argument | Description | |:---|:---| | --solvent | Pre-defined solvent name (see available solvents below) | | --solvent_smiles | SMILES string for custom solvent | | --solvent_file | Path to solvent structure file | | --solute_smiles | SMILES string for solute (optional) | | --solute_file | Path to solute structure file (optional) | | --num_solvent | Number of solvent molecules (default: 64) | | --box_size | Cubic box side in Å (auto-calculated from density if omitted) | | --tolerance | Minimum inter-molecular distance in Å (default: 2.0) | | --output_dir | Output directory |

Available pre-defined solvents: water, methanol, ethanol, acetonitrile, dmso, dmf, thf, toluene, acetone, dichloromethane, chloroform, hexane

Output files:

solvated_box.cif — Periodic structure for MD
solvated_box.xyz — Non-periodic XYZ for visualization
box_metadata.json — Box size, atom counts, solute indices

Step 2: Run MD with MLIP

Use MCP run_md tools for NPT equilibration followed by NVT production.

NPT Equilibration (stabilize density):

mcp_mace_load_model(
    model_name="MACE-MH-1",
    task_name="omol"
)
mcp_mace_run_md(
    structure_data="research/my_folder/solvation_box/solvated_box.cif",
    temperature=300,
    ensemble="npt",
    pressure=1.01325,        # 1 atm in bar
    steps=5000,              # 2.5 ps at 0.5 fs timestep
    timestep=0.5,            # 0.5 fs for systems with water (fast O-H vibrations)
    log_interval=10,
    monitor=True,
    monitor_type=["explosion", "volume"],
    output_dir="research/my_folder/npt_equilibration"
)

NVT Production (use the equilibrated structure):

mcp_mace_run_md(
    structure_data="research/my_folder/npt_equilibration/final_structure.cif",
    temperature=300,
    ensemble="nvt",
    steps=20000,             # 10 ps at 0.5 fs timestep
    timestep=0.5,
    log_interval=10,
    monitor=True,
    monitor_type="explosion",
    output_dir="research/my_folder/nvt_production"
)

Step 3: Analyze Trajectory

Run the analysis script on the production trajectory:

# Env: base-agent
python .agents/skills/chem-solution-md/scripts/analyze_solution_md.py \
    --trajectory research/my_folder/nvt_production/trajectory.traj \
    --rdf_pairs "Na-O,Cl-O,O-O" \
    --msd_elements "Na,Cl" \
    --log_interval_fs 5.0 \
    --output_dir research/my_folder/analysis

Key Parameters:

| Argument | Description | |:---|:---| | --trajectory | Path to ASE .traj trajectory file | | --rdf_pairs | Comma-separated element pairs for RDF, e.g. "Na-O,Cl-O" | | --rmax | Maximum RDF distance in Å (default: 8.0) | | --start_frame | First frame to include in analysis (default: 0) | | --stride | Frame stride (default: 1) | | --log_interval_fs | Time between frames in fs (default: 10.0) | | --msd_elements | Comma-separated elements for MSD (optional) |

Output files:

solution_analysis.json — Full results (RDF data, coordination numbers, density, MSD)
rdf_plots.png — RDF plots for each element pair
density_convergence.png — Density vs. time
msd_plot.png — MSD for specified elements (if requested)

4. Examples

Pure Water Box

# Env: base-agent
python .agents/skills/chem-solution-md/scripts/build_solvation_box.py \
    --solvent water --num_solvent 64 \
    --output_dir .agents/skills/chem-solution-md/examples/pure_water

Expected: 192 atoms (64 × 3), box ~12.4 Å

NaCl in Water

# Env: base-agent
python .agents/skills/chem-solution-md/scripts/build_solvation_box.py \
    --solute_smiles "[Na+].[Cl-]" --solvent water --num_solvent 64 \
    --output_dir .agents/skills/chem-solution-md/examples/NaCl_in_water

After MD + analysis, expected RDF peak positions:

Na–O first peak: ~2.4 Å
Cl–O first peak: ~3.2 Å
Na coordination number: ~5–6

5. Constraints

Timestep: Use 0.5 fs for water and systems with O–H/N–H bonds (fast vibrations). Can use 1.0 fs for heavier solvents without H.
Equilibration: NPT equilibration is critical. Verify density stabilization before production run.
System size: A minimum of 64 solvent molecules is recommended for reliable RDFs. Larger boxes (128–256) reduce finite-size effects.
PBC interactions: Ensure the box is large enough that periodic images do not interact (box side > 2 × rmax for RDF).
Environments:
- base-agent for box building and analysis scripts
- MCP tools for MD (any MLIP backend)

References

Martínez et al., "PACKMOL: A package for building initial configurations for molecular dynamics simulations", J. Comput. Chem., 2009. DOI
pymatgen PackmolBoxGen: pymatgen.io.packmol

Author: Bowen Deng Contact: GitHub @learningmatter-mit