We first perform a short MD simulation with harmonic position restraints on the heavy protein atoms. This allows the solvent to equilibrate around the protein without disturbing the protein structure. In addition, we use “weak coupling” temperature and pressure coupling algorithms to obtain the desired temperatue and pressure.
This step requires the top/protein_posres.itp file with the default value for the harmonic force constants of 1000 kJ mol-1 nm-2. The position restraints are switched on by setting the -DPOSRES flag in the posres.mdp file (see mdp options).
Create the run input (TPR) file, using the energy minimized system as the starting structure:
grompp -f posres.mdp -o posres.tpr -p ../top/4ake.top -c ../emin/em.pdb
The mdp file contains cut-off settings that approximate the native CHARMM values (in the CHARMM program).
Weak (Berendsen) coupling is used for both temperature and pressure to quickly equilibrate. The protein and the solvent (water and ions) are coupled as separate groups. Gromacs provides a range of groups automatically (run make_ndx -f TPR to see them) and we use the groups Protein and non-Protein (these particularly groups work roughly since Gromacs 4.5.3). If the standard groups do not work then you will have to create the groups yourself using make_ndx -f TPR -o md.ndx (which would save them in a file md.ndx) and supply it to grompp -n md.ndx.
Run the simulation locally:
mdrun -v -stepout 10 -s posres.tpr -deffnm posres -c posres.pdb
(If this is too slow on your workstation, submit to saguaro using 8 cores.)
Note
Here the runtime of 10 ps is too short for real production use; typically 1 - 5 ns are used.
In order to visually check your system, create trajectory with all molecules in the primary unitcell (-ur compact; see also below the more extensive notes on Trajectory visualization):
echo "System" | trjconv -ur compact -s posres.tpr -f posres.xtc -pbc mol -o posres_ur.xtc
Check visually:
vmd emin/em.pdb posres_ur.xtc
(If you don’t have a vmd command available on the command line then launch VMD, load the emin/em.pdb file (File ‣ New Molecule...), highlight your molecule 1 (“em.pdb”) and load the posres/posres_ur.xtc trajectory into your molecule 1”, File ‣ Load Data Into Molecule. You should see that the first frame (from the energy minimization) looks as if the water is in a distorted box shape whereas all further frames show a roughly spherical unit cell (the rhombic dodecahedron).)