MD Simulations Tutorial

Learn how to run molecular dynamics simulations for ion exchange membranes using GROMACS through OHMind.

Table of Contents

Overview

Difficulty: 🟡 Intermediate
Time: 60 minutes
Requirements: GROMACS installed and on PATH

In this tutorial, you will:

  1. Build a polymer system from SMILES
  2. Parameterize molecules with force fields
  3. Set up and run MD simulations
  4. Analyze trajectories for IEM properties

What is MD in OHMind?

OHMind integrates with GROMACS to perform molecular dynamics simulations for Ion Exchange Membranes (IEMs). This enables:

  • Property prediction - Water uptake, conductivity, swelling
  • Structure analysis - Morphology, ion distribution
  • Dynamics - Diffusion coefficients, transport mechanisms
  • Temperature effects - Property changes with temperature

Prerequisites

GROMACS Installation

Verify GROMACS is available:

# Check GROMACS
gmx --version

# If using modules
module load gromacs
gmx --version

Workspace Setup 

# Check MD workspace
echo $MD_WORK_DIR
# Should show: /path/to/workspace/GROMACS

# Ensure directory exists
ls -la $MD_WORK_DIR

Additional Tools

The MD workflow uses:

  • Antechamber (AmberTools) - For parameterization
  • PACKMOL - For system packing
  • tleap - For topology building
# Check Antechamber
which antechamber

# Check PACKMOL
which packmol

Start the Interface 

cd OHMind
./start_OHMind_cli.sh

Part 1: System Preparation

Step 1.1: Analyze Monomer

Start by analyzing your monomer structure.

Prompt:

Analyze this AEM monomer SMILES for ion exchange groups:
c1ccc(cc1)C[N+](C)(C)C

Identify:
1. Number of ionizable sites
2. Suggested counterion
3. Recommended water model

Expected Output:

Monomer Analysis:

Structure: Benzyltrimethylammonium
- Ionizable sites: 1 (quaternary ammonium)
- Charge: +1
- Counterion needed: OH⁻ or Cl⁻

Recommendations:
- Water model: SPC/E (good for ion transport)
- Counterion: OH⁻ for AEM simulation

Step 1.2: Build Polymer

Create a polymer chain from the monomer.

Prompt:

Build a polymer from this monomer SMILES:
c1ccc(cc1)C[N+](C)(C)C

Settings:
- Degree of polymerization: 20
- Number of chains: 5

Generate the polymer PDB structure.

Expected Output:

Polymer Built:
- Monomer: Benzyltrimethylammonium
- DP: 20 units per chain
- Chains: 5
- Total atoms: ~2500
- Output: polymer.pdb

Structure saved to: $MD_WORK_DIR/polymer.pdb

Step 1.3: Parameterize Molecule

Generate force field parameters.

Prompt:

Parameterize the polymer using:
1. Antechamber for atom types and charges
2. GAFF2 force field
3. AM1-BCC charge method

Generate a GROMACS-compatible topology file.

What Happens:

  1. Antechamber assigns atom types
  2. AM1-BCC calculates partial charges
  3. tleap builds topology
  4. Conversion to GROMACS format

Expected Duration: 5-10 minutes

Expected Output:

Parameterization Complete:

Force Field: GAFF2
Charge Method: AM1-BCC
Total Charge: +5 (5 chains × +1)

Files Generated:
- polymer.itp (GROMACS topology)
- forcefield.itp (force field parameters)
- polymer_charged.pdb (structure with charges)

Step 1.4: Create System

Build the complete simulation system.

Prompt:

Create an AEM simulation system with:
- 5 polymer chains (from previous step)
- OH⁻ counterions (5 total for charge neutrality)
- Water content: 20 water molecules per ion
- Water model: SPC/E

Use PACKMOL to pack the system in a cubic box.

Expected Output:

System Created:

Components:
- Polymer chains: 5
- OH⁻ ions: 5
- Water molecules: 100
- Total atoms: ~2800

Box dimensions: 4.5 × 4.5 × 4.5 nm
Density: ~1.1 g/cm³

Files:
- system.pdb (packed structure)
- system.top (system topology)

Part 2: Running Simulations

Step 2.1: Energy Minimization

Minimize the system energy to remove bad contacts.

Prompt:

Run energy minimization on the AEM system:
- Algorithm: Steepest descent
- Maximum steps: 5000
- Force tolerance: 100 kJ/mol/nm

Report the initial and final potential energy.

Expected Output:

Energy Minimization Complete:

Algorithm: Steepest descent
Steps: 2341 (converged)

Energy:
- Initial: 1.23e+06 kJ/mol
- Final: -4.56e+04 kJ/mol

Maximum force: 89.3 kJ/mol/nm (< 100, converged)
Output: em.gro, em.edr

Step 2.2: NVT Equilibration

Equilibrate temperature at constant volume.

Prompt:

Run NVT equilibration:
- Temperature: 400 K
- Duration: 100 ps
- Thermostat: V-rescale
- Time step: 1 fs

Monitor temperature stability.

Expected Duration: 5-10 minutes

Expected Output:

NVT Equilibration Complete:

Settings:
- Temperature target: 400 K
- Duration: 100 ps (100,000 steps)
- Thermostat: V-rescale (τ = 0.1 ps)

Results:
- Average temperature: 399.8 ± 5.2 K
- Temperature stable after ~20 ps

Output: nvt.gro, nvt.edr, nvt.xtc

Step 2.3: NPT Equilibration

Equilibrate pressure at constant temperature.

Prompt:

Run NPT equilibration:
- Temperature: 400 K
- Pressure: 1 bar
- Duration: 200 ps
- Barostat: Parrinello-Rahman

Monitor density and box volume.

Expected Output:

NPT Equilibration Complete:

Settings:
- Temperature: 400 K
- Pressure: 1 bar
- Duration: 200 ps
- Barostat: Parrinello-Rahman (τ = 2.0 ps)

Results:
- Average density: 1.08 ± 0.02 g/cm³
- Average pressure: 1.0 ± 50 bar
- Box volume stable after ~100 ps

Output: npt.gro, npt.edr, npt.xtc

Step 2.4: Production MD

Run the production simulation.

Prompt:

Run production MD simulation:
- Temperature: 400 K
- Pressure: 1 bar
- Duration: 1 ns
- Save trajectory every 1 ps

This will be used for property analysis.

Expected Duration: 30-60 minutes (depends on system size)

Expected Output:

Production MD Complete:

Settings:
- Duration: 1 ns (1,000,000 steps)
- Trajectory frames: 1000
- Temperature: 400.1 ± 4.8 K
- Pressure: 1.2 ± 48 bar

Output:
- md.gro (final structure)
- md.xtc (trajectory, 1000 frames)
- md.edr (energy data)
- md.log (simulation log)

Part 3: Analysis

Step 3.1: Mean Square Displacement

Calculate diffusion coefficients.

Prompt:

Calculate the mean square displacement (MSD) for:
1. OH⁻ ions
2. Water molecules

From the MSD, estimate diffusion coefficients.

Expected Output:

MSD Analysis:

OH⁻ Ions:
- MSD at 1 ns: 45.2 Ų
- Diffusion coefficient: 7.5 × 10⁻⁶ cm²/s

Water:
- MSD at 1 ns: 120.5 Ų
- Diffusion coefficient: 2.0 × 10⁻⁵ cm²/s

Interpretation:
- OH⁻ diffusion is ~3× slower than water
- Values consistent with hydrated membrane

Step 3.2: Radial Distribution Function

Analyze ion-water structure.

Prompt:

Calculate radial distribution functions (RDF) for:
1. N (ammonium) - O (water)
2. OH⁻ - O (water)

Identify hydration shell distances.

Expected Output:

RDF Analysis:

N-O(water):
- First peak: 2.85 Å (first hydration shell)
- Coordination number: 4.2 waters

OH⁻-O(water):
- First peak: 2.65 Å
- Coordination number: 5.8 waters

Interpretation:
- Ammonium groups are well-hydrated
- OH⁻ has strong hydration shell

Step 3.3: Energy Analysis

Analyze system energetics.

Prompt:

Analyze the energy components from the MD simulation:
1. Total energy stability
2. Temperature fluctuations
3. Pressure fluctuations
4. Potential energy breakdown

Expected Output:

Energy Analysis:

Stability:
- Total energy drift: < 0.1% (stable)
- Temperature: 400.1 ± 4.8 K
- Pressure: 1.2 ± 48 bar

Potential Energy Components:
- Bond: 1234 ± 45 kJ/mol
- Angle: 2345 ± 67 kJ/mol
- Dihedral: 567 ± 23 kJ/mol
- Coulomb: -45678 ± 234 kJ/mol
- LJ: 3456 ± 123 kJ/mol

Step 3.4: Property Estimation

Estimate IEM-relevant properties.

Prompt:

From the MD simulation, estimate:
1. Water uptake (λ = waters per ion)
2. Ionic conductivity (from diffusion)
3. Density
4. Swelling ratio (if reference available)

Expected Output:

Property Estimates:

Water Uptake:
- λ = 20 H₂O/ion (as designed)
- Hydration level: Moderate

Ionic Conductivity (Nernst-Einstein):
- σ ≈ 45 mS/cm at 400 K
- Note: Overestimate due to correlation effects

Density:
- ρ = 1.08 g/cm³

Comparison with Experiment:
- Typical AEM conductivity: 20-80 mS/cm
- Estimate is within expected range

Part 4: Advanced Workflows

Step 4.1: Temperature Sweep

Study temperature dependence.

Prompt:

Run a temperature sweep for the AEM system:
- Temperatures: 300 K, 350 K, 400 K
- Duration: 500 ps each
- Analyze OH⁻ diffusion at each temperature

Calculate activation energy for ion transport.

Expected Output:

Temperature D(OH⁻) (cm²/s) ln(D)
300 K 2.1 × 10⁻⁶ -13.07
350 K 4.5 × 10⁻⁶ -12.31
400 K 7.5 × 10⁻⁶ -11.80

Activation Energy: ~15 kJ/mol (Arrhenius fit)

Step 4.2: Water Content Study

Study effect of hydration level.

Prompt:

Compare AEM properties at different water contents:
- λ = 10 (low hydration)
- λ = 20 (moderate)
- λ = 30 (high hydration)

For each, run 500 ps MD at 400 K and compare conductivity.

Expected Output:

λ (H₂O/ion) Density (g/cm³) σ (mS/cm)
10 1.25 25
20 1.08 45
30 0.98 58

Step 4.3: Complete IEM Workflow

Run the full automated workflow.

Prompt:

Run a complete IEM MD workflow from SMILES:

Monomer: c1ccc(cc1)C[N+](C)(C)C
Settings:
- DP: 15
- Chains: 3
- λ: 15
- Temperature: 373 K
- Production: 2 ns

Analyze and report key IEM properties.

Expected Duration: 1-2 hours

Expected Outputs

File Locations

MD simulations create files in:

$MD_WORK_DIR/
├── system_001/
│   ├── polymer.pdb           # Initial polymer
│   ├── polymer.itp           # Polymer topology
│   ├── system.pdb            # Packed system
│   ├── system.top            # System topology
│   ├── em.mdp, em.gro        # Energy minimization
│   ├── nvt.mdp, nvt.gro      # NVT equilibration
│   ├── npt.mdp, npt.gro      # NPT equilibration
│   ├── md.mdp                # Production parameters
│   ├── md.tpr                # Run input
│   ├── md.xtc                # Trajectory
│   ├── md.edr                # Energies
│   └── md.gro                # Final structure
└── analysis/
    ├── msd.xvg               # MSD data
    ├── rdf.xvg               # RDF data
    └── energy.xvg            # Energy data

Interpreting Results

Property Good for AEM Concerning
OH⁻ diffusion > 5×10⁻⁶ cm²/s < 1×10⁻⁶ cm²/s
Conductivity > 50 mS/cm < 20 mS/cm
Water uptake 15-25 λ < 10 or > 40 λ

Troubleshooting

Common Issues

Issue Cause Solution
“GROMACS not found” Not on PATH module load gromacs or check installation
“Topology error” Missing parameters Check parameterization step
“System exploded” Bad initial structure Re-minimize, use smaller timestep
“LINCS warning” Constraint issues Reduce timestep, check topology

Checking Simulation Health 

Show me the temperature and pressure from the last 100 ps 
of my NPT simulation. Are they stable?

Restarting Failed Simulations 

My MD simulation crashed at 500 ps. Restart it from the 
last checkpoint and continue to 1 ns.

Next Steps

After completing this tutorial:

  1. Analyze with Multiwfn - Extract snapshots for QM analysis
  2. Compare systems - Study different polymer architectures
  3. Validate predictions - Compare with experimental data

Suggested Follow-up Prompts 

Extract 5 representative snapshots from my MD trajectory and 
run QM calculations on the ion-water clusters to validate 
the force field.

Compare the conductivity predictions for three different 
cation types attached to the same backbone.

See Also

*Last updated: 2025-12-23 OHMind v1.0.0*

PolyAI Team
Copyright © 2009-2025 Changchun Institute of Applied Chemistry, Chinese Academy of Sciences
Address: No. 5625, Renmin Street, Changchun, Jilin, China. Postal Code: 130022