Setting Up & Running MD Simulations

Adapted using tutorial from "http://www.mdtutorials.com/gmx/lysozyme/index.html"

All simulation parameterization is done using the HCC cluster.

Download required PDB file from "https://www.rcsb.org/"
Upload PDB file to desired folder on HCC cluster.
In a terminal or command line, navigate to the working folder you added the PDB file to.
Use "grep -v HOH prot.pdb > prot_clean.pdb to remove water from crystal structures.
1. If PDB structure has other atoms, such as SO4, ligands or PO4, run this command again, replacing HOH with other atom.
2. This can also be accomplished by removing the HOH lines in the PDB file using vi
Run the command "gmx pdb2gmx -f prot_clean.pdb -o prot_processed.gro -water spce"
1. This step will prompt you to chose a force field. A commonly used field is AMBER99SB
2. This will convert the pdb file into a file that GROMACS can use
Solvate the system by running "gmx editconf -f prot_processed.gro -o prot_newbox.gro -c -d 1.0 -bt cubic" and "gmx solvate -cp prot_newbox.gro -cs spc216.gro -o prot_solv.gro -p topol.top"
1. These commands define a 1nm box around the protein, and then fill that box with water molecules
At this point, it is beneficial to create the mdp files that GROMACS will use to run the simulations. These can be found at "http://www.mdtutorials.com/gmx/lysozyme/Files/"
Next, you want to add counter-ions to the solvent. Run the commands "gmx grompp -f ions.mdp -c prot_solv.gro -p topol.top -o ions.tpr" and "gmx genion -s ions.tpr -o prot_solv_ions.gro -p topol.top -pname NA -nname CL -neutral"
1. To run at physiological salt conditions, add the "-conc 0.1" flag before -pname
Once solvent and ions are added, the system is ready for minimization. Edit "minim.mdp" to the desired number of steps (50,000 is fairly standard for this step). Run "gmx grompp -f minim.mdp -c prot_solv_ions.gro -p topol.top -o em.tpr" and gmx mdrun -v -deffnm em".
1. Once minimization is complete, run command "gmx energy -f em.edr -o potential.xvg", inputting "10 0" when prompted and open .xvg file to check plot. Plot should show steady convergence of Epot.
The next step is to equilibrate the temperature of the system. Edit the nvt.mdp file to run for 100ps by changing nsteps to 50,000. The reference temperature you want to run the simulation at can also be edited here by changing ref_t.
1. Run commands "gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr" and "gmx mdrun -deffnm nvt". THE EQUILIBRATOIN MAY TAKE A FEW HOURS
2. Plot these results using "gmx energy -f nvt.edr -o temperature.xvg" inputting "16 0" when prompted
3. If temperature doesn't remain stable at the target temperature, a longer equilibration is necessary.
Next, equilibrate the density of the system. Edit the npt.mdp as you did the nvt.mdp.
1. Run commands "gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr" and "gmx mdrun -deffnm npt" THE EQUILIBRATOIN MAY TAKE A FEW HOURS
2. Plot results using "gmx energy -f npt.edr -o pressure.xvg" inputting "18 0" when prompted.
The system is now ready for the production simulation. Edit the md.mdp to run for the desired length (500,000 steps is a 1ns simulation). Don't forget to edit the reference temperature.
1. Run commands "gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md_0_1.tpr" and "gmx mdrun -deffnm md_0_1"
2. However, longer runs will need to run using GPU cores. To do this, you must submit a slurm job to HCC. An example script is shown below.

   #!/bin/sh
 #SBATCH --time=168:00:00   # walltime limit (HH:MM:SS)
 #SBATCH --nodes=3   # number of nodes
 #SBATCH --ntasks-per-node=32   # 16 processor core(s) per node 
 #SBATCH --mem=32G   # maximm memory per node
 #SBATCH --gres=gpu:1
 #SBATCH --partition=gpu    # gpu node(s
 #SBATCH --error=/work/powers/tessajoy17/SO237_2Pro/md.err
 cd /work/powers/tessajoy17/SO237_2Pro
 gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md_0_1.tpr
 gmx mdrun -deffnm md_0_1 -g md.log -x traj_comp.xtc -e ener.edr -nb gpu -cpi md_0_1.cpt -v
 hostname
 sleep 60