Step-by-Step MD Simulation of a Protein–Ligand Complex with GROMACS
This is a recorded CHARMM36 recipe, not a universal GROMACS map. The protein uses CHARMM36, the ligand uses CGenFF, and the write-up calls the solvent TIP3P. The commands below are the ones to copy. Menu numbers and index groups belong to the system they were answered on.
em.mdp, nvt.mdp, npt.mdp, and md.mdp are not in this note. The production run is the 100 ns job written as md_0_100. Ions use the short ions.mdp that is included.
1 | 13 in make_ndx are the answers printed for this system. Solvation uses spc216.gro even though the force field named here is CHARMM36 with TIP3P. That coordinate file is left as written.0. Repair the PDB
Missing atoms get filled before the force field sees the file.
conda activate pdbfixer
pdbfixer protein.pdb --output=protein_fixed.pdb
1. Separate the ligand
Unknown residues and cofactors should not go through pdb2gmx as part of the protein. Two variants are recorded. One pulls UNK only. The other also drops FAD.
grep UNK protein.pdb > unk.pdb
grep -v "UNK" protein.pdb > clean.pdb
grep -v -e "UNK" -e "FAD" protein_fixed.pdb > clean.pdb
conda activate pdbfixer
pdbfixer clean.pdb --output=clean.pdb
2. Protein topology
pdb2gmx is run with -ter and -ignh. The six lines are the menu answers recorded for this CHARMM36 setup. On another GROMACS build the same digit can be a different force field or a different terminus. Read the menu.
gmx pdb2gmx -f clean.pdb -o processed.gro -ter -ignh << EOF
1
1
0
0
0
0
EOF
3. Ligand MOL2
Hydrogens are added and a 3D guess is generated. The sed line sets the residue name to UNK.
obabel unk.pdb -O unk.mol2 --addh --gen3d
sed -i '2s/.*/UNK/; s/UNK1/UNK/g' unk.mol2
4. Bond order in the MOL2
perl sort_mol2_bonds.pl unk.mol2 unk_fix.mol2
5. CGenFF
The binary path is the one used on this machine. The conversion script targets charmm36-jul2022.ff. Penalty scores in the CGenFF log are a reason to stop. This recipe does not inspect them for you.
~/mdrun/silcsbio.2024.1/cgenff/cgenff unk_fix.mol2 -f unk.str
chmod +x cgenff_charmm2gmx.py
conda activate cgenff
python cgenff_charmm2gmx.py UNK unk_fix.mol2 unk.str charmm36-jul2022.ff
rm cgenff_charmm2gmx.py
6. Ligand coordinates as GRO
The next step expects unk_ini.pdb and, when a cofactor is present, fad_ini.pdb. The command that writes those two names is not in this note. The MOL2 above is the file that went to CGenFF.
gmx editconf -f unk_ini.pdb -o unk.gro
gmx editconf -f fad_ini.pdb -o fad.gro
7. Merge into one GRO
Protein plus ligand. The second line of a GRO file is the atom count, so it is rewritten after the append. The box line is dropped from the protein file and a blank box line is added back at the end.
cp processed.gro complex.gro && sed -i '$d' complex.gro && tail -n +3 unk.gro >> complex.gro && echo " " >> complex.gro
num1=$(sed -n '2p' complex.gro)
num2=$(sed -n '2p' unk.gro)
sum=$((num1 + num2))
sed -i "2s/.*/ $sum/" complex.gro
Protein plus FAD plus UNK:
cp processed.gro complex.gro && \
sed -i '$d' complex.gro && \
tail -n +3 fad.gro >> complex.gro && \
tail -n +3 unk.gro >> complex.gro && \
echo " " >> complex.gro
num1=$(sed -n '2p' processed.gro)
num2=$(sed -n '2p' fad.gro)
num3=$(sed -n '2p' unk.gro)
sum=$((num1 + num2 + num3))
sed -i "2s/.*/ $sum/" complex.gro
8. Topology includes
These sed edits are brittle. They assume the CHARMM36 comments and the molecule name Protein_chain_A are still in topol.top. Open the topology after they run.
Ligand only:
sed -i '/; Include Position restraint file/{:a;N;/#endif/!ba;s/\n#endif\n/\n#endif/;s/#endif/#endif\n\n; Include ligand topology\n#include "unk.itp"/}' topol.top
sed -i '/; Include forcefield parameters/{:a;N;/#include ".\/charmm36-jul2022.ff\/forcefield.itp"/!ba;s/#include ".\/charmm36-jul2022.ff\/forcefield.itp"/&\n\n; Include ligand parameters\n#include "unk.prm"/}' topol.top
sed -i '/\[ molecules \]/,/Protein_chain_A/ {/Protein_chain_A/ a\\nUNK 1\n}' topol.top
FAD and UNK:
sed -i '/; Include Position restraint file/{:a;N;/#endif/!ba;s/#endif/#endif\n\n; Include ligand topologies\n#include "fad.itp"\n#include "unk.itp"/}' topol.top
sed -i '/; Include forcefield parameters/{:a;N;/#include ".\/charmm36-jul2022.ff\/forcefield.itp"/!ba;s|#include ".\/charmm36-jul2022.ff\/forcefield.itp"|&\n\n; Include ligand parameters\n#include "fad.prm"\n#include "unk.prm"|}' topol.top
sed -i '/\[ molecules \]/,/Protein_chain_A/ {/Protein_chain_A/ a\\nFAD 1\\nUNK 1\n}' topol.top
9. Box and solvent
The buffer is 1.4 nm and the box is cubic. Solvation is spc216.gro. That is the SPC216 coordinate file, used here with the CHARMM36 topology. It is not rewritten as tip3p.gro.
gmx editconf -f complex.gro -o newbox.gro -bt cubic -c -d 1.4
gmx solvate -cp newbox.gro -cs spc216.gro -p topol.top -o solv.gro
10. Ions
The ion .mdp is only a steep minimizer so genion has a .tpr. Neutralize, then add 0.15 M NaCl. The heredoc answers group 15. On this system that was solvent. genion prints the group list. Choose solvent, whatever number it has.
cat <<EOF > ions.mdp
integrator = steep
emtol = 1000.0
nsteps = 50000
cutoff-scheme = Verlet
coulombtype = PME
rcoulomb = 1.0
rvdw = 1.0
constraints = none
EOF
gmx grompp -f ions.mdp -c solv.gro -p topol.top -o ions.tpr -maxwarn 3
gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -neutral -conc 0.15 <<EOF
15
EOF
11. Energy minimization
em.mdp itself is not listed.
gmx grompp -f em.mdp -c solv_ions.gro -p topol.top -o em.tpr
gmx mdrun -v -deffnm em
12. Ligand position restraints
Non-hydrogen atoms of the ligand are restrained with -fc 1000 1000 1000. Group 3 is the selection recorded after 0 & ! a H*. Check the index make_ndx prints.
gmx make_ndx -f unk.gro -o index_unk.ndx << EOF
0 & ! a H*
q
EOF
gmx genrestr -f unk.gro -n index_unk.ndx -o posre_unk.itp -fc 1000 1000 1000 <<EOF
3
EOF
sed -i '/; Include water topology/{x;s/.*/\n; Ligand position restraints\n#ifdef POSRES\n#include "posre_unk.itp"\n#endif\n/;G}' topol.top
13. Complex index
1 | 13 unions two groups that were protein and ligand on this system. Use the numbers your own index file prints.
gmx make_ndx -f solv_ions.gro -o index.ndx << EOF
1 | 13
q
EOF
14. NVT
Temperature equilibration at 300 K, with position restraints. The thermostat settings live in nvt.mdp, which is not pasted here.
gmx grompp -f nvt.mdp -c solv_ions.gro -r solv_ions.gro -p topol.top -n index.ndx -o nvt.tpr
gmx mdrun -deffnm nvt
15. NPT
Pressure equilibration at 1 bar with the Parrinello–Rahman barostat.
gmx grompp -f npt.mdp -c nvt.gro -t nvt.cpt -r nvt.gro -p topol.top -n index.ndx -o npt.tpr
gmx mdrun -deffnm npt
16. Production
The file names are the 100 ns run.
gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -n index.ndx -o md_0_100.tpr
gmx mdrun -deffnm md_0_100
Before you plot anything
Fix periodic boundary crossings before RMSD, radius of gyration, or SASA. Drop the first 5–10 ns. A flat RMSD, radius of gyration, and SASA is the stability check this note uses. PCA and a free-energy landscape are the next look at slow motions. Those analysis commands are a separate post.
