Computational Chemistry · Solvent Molecular Dynamics

Solvent Molecular Dynamics Simulation Workflow

Pure solvent box preparation · equilibration · production MD · complete solvent structure and dynamics analysis using GROMACS, VMD, MDAnalysis, and Python

GROMACS SPC/E TIP3P TIP4P/2005 RDF MSD H-bond Lifetime Diffusion · Residence Time · FES

System
Prep

🧪

Select Solvent Model

Step 1 SPC/E · TIP4P

📦

Build Solvent Box

Step 2 solvate · Packmol

📄

Topology & Force Field

Step 3 topol.top · .itp

⚡

Add Ions / Conc. Control

Step 4 genion · conc

Sim. &
Analysis

📉

Energy Minimization

Step 5 mdrun EM

🌡️

NVT & NPT Equil.

Step 6 300 K · 1 bar

🚀

Production MD

Step 7 10–100 ns

📊

Solvent Analysis A–Z

Step 8 RDF·MSD·HB·FES

System Preparation

🧪

Select Solvent Model SPC/E · TIP3P · TIP4P · OPLS · GAFF

Choose water, organic, mixed, or ionic solvent model based on target property, force-field compatibility, and experimental reference

Water model Organic solvent Force-field compat. Density Diffusion Dielectric constant

▼

Why Solvent Model Selection Matters

The solvent model controls hydrogen bonding, diffusion, density, dielectric behavior, viscosity, surface tension, and solvation-shell structure. The correct model depends on the physical property most important for the simulation — not on which model is most popular.

Common Water Models

Model	Sites	Best For	Key Note
SPC	3-site	Simple GROMACS systems	Fast, widely used, less accurate for bulk physics
SPC/E	3-site + polarization	Bulk water, ion hydration	Better density and diffusion than SPC
TIP3P	3-site	CHARMM/AMBER biomolecular	May diffuse too fast in some settings
TIP4P	4-site massless	Improved electrostatics	Better water structure
TIP4P/2005	4-site	Physical chemistry, water properties	Best broad reproduction of water properties

Organic Solvent Checklist

Confirm partial charges, bonded parameters, Lennard-Jones parameters
Verify combination rules match force field
Confirm compatibility with solute force field
Check correct molecule, residue, and atom naming
Use OPLS-AA, GAFF, CGenFF, or ATB for parameterization

Water models

SPC/E · TIP3P · TIP4P/2005

Organic solvent

OPLS-AA · GAFF · CGenFF · ATB

Key validation

density · diffusion · RDF · dielectric

Main risk

incompatible FF parameters

Solvent Model Property Comparison

Normalized agreement with experimental density, diffusion, and dielectric constant

Higher bars indicate better agreement with experimental bulk water properties for that model

📦

Build Solvent Box gmx solvate · Packmol

Generate cubic, dodecahedral, or rectangular solvent box using GROMACS, Packmol, or CHARMM-GUI

Cubic box Dodecahedral Solvent packing PBC Packmol

▼

Water Box with GROMACS

gmx solvate -cs spc216.gro -o solvent_box.gro -p topol.top
# -box 4.0 4.0 4.0 for explicit cubic box dimensions

Organic / Mixed Solvent with Packmol

packmol < packmol_input.inp

# Example Packmol input:
tolerance 2.0
filetype pdb
output solvent_box.pdb

structure methanol.pdb
number 1000
inside box 0. 0. 0. 60. 60. 60.
end structure

Box Types and Use Cases

Cubic box — simple and universal; most common for pure solvent
Dodecahedral box — ~71% cubic volume for spherical solutes; reduces solvent count
Rectangular slab — interfaces, membranes, confined layers
Mixed solvent box — multiple components at desired mole/volume fraction

Box type

cubic / dodecahedral / rectangular

Builder

gmx solvate / Packmol / CHARMM-GUI

PBC

periodic in x, y, z

Quality check

no vacuum gap, no severe overlap

After building, always verify starting density visually (VMD) and numerically. A vacuum bubble in the initial box will cause explosive forces at the start of dynamics.

📄

Define Topology & Force Field topol.top · .itp

Prepare solvent topology, include molecule parameters, define correct solvent count, verify FF compatibility

topol.top Molecule count Solvent .itp Combination rules

▼

Pure Water Topology Example

#include "charmm36.ff/forcefield.itp"
#include "spce.itp"

[ system ]
Pure SPC/E Water Box

[ molecules ]
SOL 4096

Mixed Solvent Topology Example

#include "oplsaa.ff/forcefield.itp"
#include "water.itp"
#include "dmso.itp"

[ system ]
Water-DMSO Mixed Solvent

[ molecules ]
SOL 3000
DMSO 500

Topology Checklist (9 points)

Force field included correctly
Solvent .itp included correctly
Ion .itp included if ions present
Molecule names match coordinate file exactly
Molecule count matches coordinate file
Water model constraints match .mdp settings
Nonbonded settings match force-field recommendation
Long-range electrostatics configured (PME)
Temperature and pressure coupling groups defined sensibly

Most common fatal error: residue or molecule name mismatch between coordinate file and [ molecules ] directive. GROMACS will fail with "number of coordinates does not match topology."

⚡

Add Ions / Concentration Control genion · concentration

Add ions or control solvent composition for electrolyte, physiological salt, buffer, ionic liquid, or mixed-solvent systems

NaCl 0.15 M Ionic strength Counterions Ion hydration

▼

Add Ions at Physiological Concentration

gmx grompp -f ions.mdp -c solvent_box.gro -p topol.top -o ions.tpr
gmx genion -s ions.tpr -o solvent_ions.gro -p topol.top -pname NA -nname CL -conc 0.15

# For neutralization only:
gmx genion -s ions.tpr -o solvent_neutral.gro -p topol.top -pname NA -nname CL -neutral

Important Ion Rules

Use -neutral for charge neutralization only; use -conc for specific molar concentration
Ion names must match the force-field .itp exactly
Never mix ion parameters from a different force field without validation
Check ion-water RDF after simulation to confirm correct hydration shell geometry
Check ion-pair RDF to detect unrealistic salt aggregation

Salt conc.

0.15 M NaCl

Ion tool

gmx genion

Main validation

ion-water RDF

Main risk

ion clustering or wrong parameters

Equilibration

📉

Energy Minimization mdrun EM

Remove steric clashes, high-energy contacts, and bad ion placements before dynamics

Steepest descent Fmax < 1000 Epot negative em.edr

▼

gmx grompp -f minim.mdp -c solvent_ions.gro -p topol.top -o em.tpr
gmx mdrun -v -deffnm em
gmx energy -f em.edr -o potential.xvg # "11 0" to select Potential

Output Files

File	Contents
em.gro	Minimized coordinates → input for NVT
em.edr	Energy trajectory — extract with gmx energy
em.log	Convergence log — check for warnings

Troubleshooting Failures

Positive or non-converging Epot → inspect structure visually in VMD for overlaps
Wrong topology molecule count or residue naming
Bad ion placement near charged headgroups → explosive Coulomb force
If using Packmol, reduce initial packing density or increase tolerance

Potential Energy Minimization

Steepest descent convergence

Sharp decrease → stable minimum confirms system is ready for dynamics

Maximum Force (Fmax)

Convergence below threshold

Fmax must fall below the emtol threshold (typically 1000 kJ mol⁻¹ nm⁻¹)

🌡️

NVT & NPT Equilibration thermostat · barostat

Stabilize solvent temperature, pressure, density, and box volume before production MD

300 K 1 bar V-rescale Parrinello-Rahman density

▼

NVT — Temperature Stabilization

gmx grompp -f nvt.mdp -c em.gro -p topol.top -o nvt.tpr
gmx mdrun -deffnm nvt
gmx energy -f nvt.edr -o temperature.xvg # "16 0"

NPT — Pressure and Density Stabilization

gmx grompp -f npt.mdp -c nvt.gro -t nvt.cpt -p topol.top -o npt.tpr
gmx mdrun -deffnm npt
gmx energy -f npt.edr -o pressure.xvg # "17 0"
gmx energy -f npt.edr -o density.xvg # "23 0"
gmx energy -f npt.edr -o volume.xvg # "24 0"

Key Parameter Choices

Thermostat: V-rescale for NVT equilibration; Nosé-Hoover for production MD
Barostat: Berendsen for early NPT equilibration; Parrinello-Rahman for production
Density convergence: the most reliable indicator that NPT equilibration is complete
Pressure: use running average — instantaneous pressure fluctuates by ±hundreds of bar normally

NVT target

300 K

NPT target

1 bar

SPC/E density

~997 kg/m³

Pressure avg

use running avg

Temperature — NVT

Stabilization at 300 K

Temperature plateau confirms NVT equilibration is complete

Pressure Running Average — NPT

1 bar reference

Running average converges to 1 bar — instantaneous fluctuations are normal

Solvent Density Convergence — NPT

kg m⁻³ vs time

Stable density confirms correct pressure-volume equilibration

Box Volume Stabilization — NPT

nm³ vs time

Volume converges as barostat adjusts box to target pressure

Production MD

🚀

Production Molecular Dynamics mdrun production

Run unrestrained solvent MD to generate trajectory for structure, diffusion, H-bonding, and solvation analysis

md.xtc md.edr 10–100 ns Nosé-Hoover Parrinello-Rahman

▼

gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md.tpr
gmx mdrun -deffnm md

Production Length Guide

Property	Minimum Length	Preferred
RDF, density	1–5 ns	10 ns+
Diffusion coefficient (MSD)	10 ns	50 ns+
H-bond lifetime	10 ns	50 ns+
Residence time near solute	50 ns	100 ns+
Dielectric constant	100 ns	very long
Viscosity	100 ns	very long + NEMD
Mixed solvent phase behavior	100 ns	200 ns+

Key Output Files

Trajectory

md.xtc / md.trr

Energy file

md.edr

Final structure

md.gro

Checkpoint

md.cpt

Convergence principle: Production length should be determined by property convergence, not a fixed number. Always verify whether the measured property reaches a stable average value before reporting results.

Solvent Structure & Dynamics Analysis

📊

Complete A–Z Solvent Analysis RDF · MSD · hbond · density · FES

Structure · solvation shell · H-bond network · diffusion · rotational dynamics · residence time · ion hydration · spatial density · clustering · free energy surfaces

RDF Coord. number H-bond lifetime MSD Residence time Dipole orient. Cluster size FES

▼

1. Radial Distribution Function (RDF)

gmx rdf -s md.tpr -f md.xtc -o rdf_OW_OW.xvg -ref 'name OW' -sel 'name OW'
# Also: O-H RDF, ion-water RDF, solute-water RDF, cosolvent-water RDF

Radial Distribution Function g(r)

SPC/E water O–O · first solvation shell

First peak ~0.28 nm = first solvation shell. First minimum defines coordination cutoff

Coordination Number

Running integral of RDF · cumulative N(r)

Plateau at first RDF minimum gives first-shell coordination number (~4.5 for SPC/E water)

2. Hydrogen Bond Network

gmx hbond -s md.tpr -f md.xtc -num hbnum.xvg # H-bond count vs time
gmx hbond -s md.tpr -f md.xtc -ac hbac.xvg # Autocorrelation C(t)

Hydrogen Bond Network Stability

H-bond count vs time

Stable average H-bond count indicates equilibrated solvent network

H-bond Lifetime Autocorrelation C(t)

MDAnalysis WaterDynamics · gmx hbond -ac

Slow decay = long-lived H-bonds and slower solvent rearrangement

3. Translational Diffusion (MSD)

gmx msd -s md.tpr -f md.xtc -o msd.xvg
# D = slope / 6 (3D Einstein relation) from linear MSD region

Mean Square Displacement

Translational motion · linear region → diffusion

Slope of linear MSD region = 6D. Early nonlinear region must not be used for fitting

Diffusion Coefficient Comparison

Across solvent models and conditions

TIP4P/2005 reproduces experimental D(water) ≈ 2.3 × 10⁻⁵ cm²/s best

4. Residence Time & Dipole Orientation

python -m MDAnalysis WaterDynamics survival_probability # residence time
python dipole_orientation.py # MDAnalysis dipole vector analysis

Residence Time — Survival Probability

First solvation shell retention

Slow decay = long-lived hydration shell. Fast decay = rapid solvent exchange

Dipole Orientation Distribution

cos(θ) probability near charged environment

Biased distribution indicates electrostatic ordering of solvent dipoles

5. Spatial Density & Cluster Analysis

gmx density -s md.tpr -f md.xtc -o density_z.xvg -d Z # axial density
python cluster_analysis.py # MDAnalysis distance-cutoff clustering

Solvent Occupancy Heatmap

2D spatial density projection · hydration hotspots

High-density regions reveal persistent hydration sites and preferred solvent positions

Cluster Size Distribution

Aggregation and phase separation

Large clusters indicate aggregation or phase separation in mixed solvent systems

6. Free Energy Surface (FES)

python fes_2d.py # 2D histogram of solvent coordinates → -kT ln P(x,y)

Solvent Free Energy Surface

2D FES: distance vs H-bond number · kcal mol⁻¹

FES minima identify preferred solvent configurations. Barriers show exchange difficulty between states

Additional Analysis Methods (A–Z)

Analysis	Tool	Interpretation
SASA	gmx sasa	Solute surface exposure to solvent
Velocity ACF	Python/MDAnalysis	Short-time dynamics, cage effects
Rotational relaxation C₂(t)	MDAnalysis WaterDynamics	Reorientation of water dipole/OH vectors
Tetrahedral order q	Python custom	Local ice-like vs liquid structure of water
Dielectric constant	Dipole fluctuations	Requires long simulation; slow convergence
Density profile (z-axis)	gmx density	Interface layering and confinement effects
Ion pairing RDF	gmx rdf	Na-Cl, ion-water hydration shell geometry
Solvation shell exchange	RDF cutoff + Python	Rate of solvent entry/exit from first shell
Preferential solvation	Kirkwood-Buff / RDF	Water vs cosolvent enrichment near solute
Electrostatic potential	Charge density integration	Interface dipole and charge separation
H-bond network graph	NetworkX + MDAnalysis	Connectivity, fragmentation by salt/temp
Solvation free energy	TI / FEP / BAR	Thermodynamic cost of solvation

🧾

Common Solvent Models

Solvent / Model	Type	Best Use	Key Notes
SPC	3-site water	Simple GROMACS systems	Fast and common; less accurate for bulk physics
SPC/E	3-site + polarization correction	Bulk water and ion hydration	Better density and diffusion than SPC
TIP3P	3-site water	CHARMM/AMBER biomolecular	May diffuse too fast; highly compatible
TIP4P	4-site, massless charge	Improved electrostatics	Better water structure than 3-site models
TIP4P/2005	4-site water	Physical chemistry, water properties	Best broad reproduction of water properties
Methanol / Ethanol	Organic polar solvent	Cosolvent and solvation studies	Use OPLS, GAFF, CHARMM, or CGenFF
DMSO	Polar aprotic solvent	Protein cosolvent, drug solubility	Strongly affects hydrogen-bond network
Acetonitrile	Polar aprotic solvent	Electrochemistry, organic solvation	Check dielectric and diffusion properties
Ionic liquids	Charged solvent	Electrochemistry, catalysis	Requires carefully validated ion-pair parameters

⚙

Simulation Parameters

Force field

CHARMM36 / AMBER / OPLS-AA / GAFF

Water model

SPC/E · TIP3P · TIP4P/2005

Temperature

298–310 K

Pressure

1 bar

Thermostat

V-rescale / Nosé-Hoover

Barostat

Berendsen → Parrinello-Rahman

Electrostatics

PME

Constraints

LINCS (h-bonds)

Time step

2 fs

Production length

10–100 ns+

📈

Analysis Metrics

Density

Solvent packing / NPT quality

RDF

Solvation shell structure g(r)

Coord. N

Shell population from RDF integral

H-bonds

Network strength and stability

HB lifetime

C(t) autocorrelation decay

MSD

Translational motion → D

Diffusion D

Solvent mobility (cm²/s)

Residence τ

Survival probability decay

Dipole

Orientation P(cos θ)

Density map

3D hydration hotspots

Clusters

Aggregation / phase separation

FES

Preferred solvent states

Solvent Molecular Dynamics Simulation Workflow

Why Solvent Model Selection Matters

Common Water Models

Organic Solvent Checklist

Water Box with GROMACS

Organic / Mixed Solvent with Packmol

Box Types and Use Cases

Pure Water Topology Example

Mixed Solvent Topology Example

Topology Checklist (9 points)

Add Ions at Physiological Concentration

Important Ion Rules

Output Files

Troubleshooting Failures

NVT — Temperature Stabilization

NPT — Pressure and Density Stabilization

Key Parameter Choices

Production Length Guide

Key Output Files

1. Radial Distribution Function (RDF)

2. Hydrogen Bond Network

3. Translational Diffusion (MSD)

4. Residence Time & Dipole Orientation

5. Spatial Density & Cluster Analysis

6. Free Energy Surface (FES)

Additional Analysis Methods (A–Z)

Company

Links

Support

Recommend

Solvent Molecular Dynamics Simulation Workflow

Why Solvent Model Selection Matters

Common Water Models

Organic Solvent Checklist

Water Box with GROMACS

Organic / Mixed Solvent with Packmol

Box Types and Use Cases

Pure Water Topology Example

Mixed Solvent Topology Example

Topology Checklist (9 points)

Add Ions at Physiological Concentration

Important Ion Rules

Output Files

Troubleshooting Failures

NVT — Temperature Stabilization

NPT — Pressure and Density Stabilization

Key Parameter Choices

Production Length Guide

Key Output Files

1. Radial Distribution Function (RDF)

2. Hydrogen Bond Network

3. Translational Diffusion (MSD)

4. Residence Time & Dipole Orientation

5. Spatial Density & Cluster Analysis

6. Free Energy Surface (FES)

Additional Analysis Methods (A–Z)

Company

Links

Support

Recommend

Login with your site account