ChemServerHow to do a molecular dynamics simulation
The standard workflow for simulating a protein ligand complex in water, the settings most papers use, and how to tell whether your simulation says anything useful.
What MD tells you
Molecular dynamics solves Newton's equations of motion for every atom in the system, femtosecond by femtosecond. For drug design the main questions are: does the docked ligand stay in the pocket, which interactions survive, and how flexible are the protein and the ligand? Common programs are GROMACS, AMBER, OpenMM and NAMD; CHARMM-GUI is a popular web tool for building the input files.
Step 1. Start from a good complex
- Use a crystal structure of the complex or your best docked pose.
- Repair missing atoms and short loops, and set protonation states for pH 7.4. Look closely at histidines and catalytic residues in the pocket.
- Make sure the ligand has correct bond orders and charge; a SMILES string is the easiest way to define them.
Step 2. Pick a force field and ligand parameters
- Protein: AMBER ff14SB or ff19SB, or CHARMM36m.
- Ligand: GAFF2 with AM1-BCC charges, OpenFF (Sage) or CGenFF for CHARMM.
- Water: use the model that belongs with the protein force field, for example TIP3P with ff14SB or CHARMM36m and OPC with ff19SB.
Step 3. Solvate and add ions
- Put the complex in a box with at least 1.0 nm of water between the protein and the box edge.
- A truncated octahedron or rhombic dodecahedron needs fewer waters than a cube.
- Neutralize the charge and add about 0.15 M NaCl to mimic physiological salt.
Step 4. Minimize the energy
Remove bad contacts left by building and solvating, for example with steepest descent until the maximum force falls below about 1000 kJ/mol/nm.
Step 5. Equilibrate
- NVT, about 100 ps: heat to 300 K at constant volume with the protein and ligand heavy atoms restrained.
- NPT, about 100 ps to 1 ns: bring the pressure to 1 bar and the density to a stable value, still restrained.
Check that temperature, pressure and density have levelled off before you go on.
Step 6. Run the production simulation
- Remove the restraints.
- Time step of 2 fs with bonds to hydrogen constrained, or 4 fs with hydrogen mass repartitioning.
- PME for long range electrostatics and a 1.0 nm cut off for the rest.
- Length: 50 to 100 ns per run is common for checking complex stability; three shorter replicas with different starting velocities are often more convincing than one long run.
- Save a frame every 10 to 100 ps.
Speed depends on the system size and the GPU; a modern GPU typically runs a solvated protein ligand system at tens to hundreds of nanoseconds per day.
Step 7. Analyse
| Analysis | What to look for |
|---|---|
| Backbone RMSD | Levels off after a few nanoseconds: the protein is stable |
| Ligand RMSD (fitted on the protein) | Stays within about 2 to 3 Å: the pose holds |
| RMSF per residue | Which loops move; pocket residues should usually be calm |
| Hydrogen bonds and contacts | Which interactions last for most of the run |
| Radius of gyration | The protein stays compact |
| MM/GBSA | Ranks ligands against each other; not an absolute binding free energy |
Common mistakes
- Wrong ligand charge or bond orders, which no amount of simulation fixes.
- Skipping equilibration and starting production from a minimized structure.
- Calling a complex stable from a single short run.
- Measuring ligand RMSD without first fitting the trajectory on the protein.
Do it on ChemServer
Molecular dynamics on ChemServer covers every step above, from a docked pose to the final plots, on GPU servers. It is available on request: email chris@chemserver.online with your system.