GROMACS #01 – simple protein in simple water

1. Process the protein to generate a topology
gmx pdb2gmx -f protein.pdb -o protein_solo_proc.gro -water spce -ter -p topology.top
2. Create a box centred on the protein
gmx editconf -f protein_solo_proc.gro -o protein_solo_newbox.gro -c -d 1.0 -bt cubic

NOTE: The most important question to ask now is why choose 1.2 nm? The short answer is that you don't want the protein to 'see' its periodic image across the boundary of the box. If you refer to the original paper for the AMBER 03 force field (Duan, et al., J. Comput. Chem. 2003, 24, 1999-2012), they use short-range VDW and electrostatic cut-offs of 0.8 nm - when you begin simulating, it is very important to do the same. In order for the protein to avoid seeing its image across the periodic boundary, it must be at least twice the cut-off distance from the next nearest image of itself. That being said, the space between the protein and the edge of the box only really needs to be slightly larger than 0.8 nm, but it is still a good idea to include extra space especially if there is a chance the protein might unfold a little bit or if is an oblong shape. It is better to have a slightly larger box size now than to find out later that your protein was interacting with its periodic image during the simulation.

3. Solvate the box with SPC216 water
gmx solvate -cp protein_solo_newbox.gro -cs spc216.gro -o protein_solo_solv.gro -p topology.top
4. Add ions to neutralize the system charges
gmx grompp -f ions_solo.mdp -c protein_solo_solv.gro -p topology.top -o ions_solo.tpr
gmx genion -s ions_solo.tpr -o protein_solo_solv_ions.gro -p topology.top -pname NA -nname CL -neutral -conc 0.15
$ Select SOL group
5. Prepare a minimization run
gmx grompp -f minim_solo.mdp -c protein_solo_solv_ions.gro -p topology.top -o em_solo.tpr
6. Minimize
gmx mdrun -v -deffnm em_solo

NOTE: We choose a standard cut-off of 1.0 nm, both for the neighborlist generation and the coulomb and Lennard-Jones interactions. nstlist=10 means it is updated at least every 10 steps, but for energy minimization it will usually be every step. Energies and other statistical data are stored every 10 steps (nstenergy), and we have chosen the more expensive Particle-Mesh-Ewald summation for electrostatic interactions. "e treatment of nonbonded interactions is frequently bordering to religion. One camp advocates standard cutoffs are fine, another swears by switched-off interactions, while the third wouldn’t even consider anything but PME. One argument in this context is that ‘true’ interactions should conserve energy, which is violated by sharp cut-offs since the force is no longer the exact derivative of the potential. On the other hand, just because an interaction conserves energy doesn’t mean it describes nature accurately. In practice, the difference is most pronounced for systems that are very small or with large charges, but the key lesson is really that it is a trade-off. PME is great, but also clearly slower than cut-offs. Longer cut-offs are always better than short ones (but slower), and while switched interactions improve energy conservation they introduce artificially large forces. Using PME is the safe option, but if that’s not fast enough it is worth investigating reaction-field or cut-off interactions. It is also a good idea to check and follow the recommended settings for the force field used.

7. Equilibrate temperature and pressure
gmx grompp -f nvt_solo.mdp -c em_solo.gro -p topology.top -o nvt_solo.tpr
gmx mdrun -v -deffnm nvt_solo

NOTE: The 'define = -DPOSRES' line initiates the backbone restraints on the protein, which should remain on until production MD. Temperature coupling is performed using the Berendsen method. In order to be sure velocities (and therefore temperature) are evenly distributed across all of the molecule types, it is a good idea to couple protein atoms and non-protein atoms to separate baths. Without separate coupling, you run the risk weird phenomena such as a system where the water molecules are at 350 K and the protein is at 200 K. Finally, all bonds are constrained using the linear constraint solver (LINCS) algorithm. Review the GROMACS manual if you are unfamiliar with the implementation or purpose of any of these parameters.

gmx grompp -f npt_solo.mdp -c nvt_solo.gro -t nvt_solo.cpt -p topol_solo.top -o npt_solo.tpr
gmx mdrun -v -deffnm npt_solo

  1. Preform the production run
    gmx grompp -f md_solo.mdp -c npt_solo.gro -t npt_solo.cpt -p topology.top -o md_solo.tpr
    gmx mdrun –v –deffnm md_solo