AMBER TUTORIAL-1: How to simulate a protein-basic steps for md simulation in amber (AMBER 16)
BLOG #1
This blog is for beginners in MD simulation. After installation of AMBER, you may start this tutorial by choosing any protein with no ligand or water i.e. APO protein.
Key points to consider:
- Prepare your protein in terms of bond order, terminal capping (add ACE and NME on the N and C- terminals to prevent terminal residues interaction with the solvent), add disulfide bonds if exists in a protein structure like in the case of an immunoglobulin protein family (eg. CD28, PD-1, PD-L1, CD80, etc.), fill the missing regions or breaks in your crystal, etc.
#if you have prepared the protein in Maestro (Schrodinger), change the NMA to NME, increment the residue number, convert CA to CH3 in the pdb file with help of a text editor .. see example below:)
file before editing:
ATOM 1011 CE1 HIS A 143 26.670 14.396 200.638 1.00 0.00 C
ATOM 1012 NE2 HIS A 143 25.957 14.123 201.724 1.00 0.00 N
ATOM 1013 N NMA A 143A 24.758 8.646 198.431 1.00 0.00 N
ATOM 1014 CA NMA A 143A 24.368 7.303 198.023 1.00 0.00 C
file after editing:
ATOM 1011 CE1 HIS A 143 26.670 14.396 200.638 1.00 0.00 C
ATOM 1012 NE2 HIS A 143 25.957 14.123 201.724 1.00 0.00 N
ATOM 1013 N NME A 144 24.758 8.646 198.431 1.00 0.00 N
ATOM 1014 CH3 NME A 144 24.368 7.303 198.023 1.00 0.00 C
PS:-save this file and cross-check/view this file in any visualization software. the pdb file architecture should be maintained (like spaces) otherwise your structure will appear distorted.
- Remove hydrogens from the protein file as amber itself adds hydrogen
PS: '#' is used here after commands to explain them.
export AMBERHOME=set path of the directory
(example: export AMBERHOME=/home/user/Desktop/amber16)
PS: one can add this above command in their .bashrc (example for centos7 terminal: type: vi ~/.bashrc... insert the above command... save the file and exit (wq!)... and type: source ~/.bashrc.. ur path is set for amberhome)
to cross-check the path after export type: echo $AMBERHOME
Methodology for a Protein (APO: not without any ligand or any binding partner) simulation
$AMBERHOME/bin/tleap -f oldff/leaprc.ff14SB
#here amberff14SB forcefield is used
rec=loadpdb pro_noh.pdb #load ur prepared and without hydrogens protein file into a variable 'rec'
set default PBRadii mbondi2
saveamberparm rec pro.prmtop pro.inpcrd #saving prmtop and inpcrd files
source leaprc.water.tip3p #using water model TIP3P
charge rec #checking charge in protein structure
#suppose your output looks like this:
#Total unperturbed charge: -2.999000
#Total perturbed charge: -2.999000
# then u need to add counterions to it.. here we will add 3 Na+ ions...
addions2 rec Na+ 3
#otherwise if charge was negative then add Cl- ions like this: addions2 rec Cl- 3
charge rec #to check if charges are neutralised
solvatebox rec TIP3PBOX 10.0 ##OR #solvateoct dna2 TIP3PBOX 8.0####
#adding solvent or water to the system around 10 Å from each atom of the solute (protein)
#the solvation box can be of any type either cube, octahedral, etc. ..depends upon you.. usually a octahedral box than cube is considered so a to add only sufficient solvent around the protein that could be useful for increasing computational efficiency.. i.e. to have a long md simulation faster.. for more details please check amber16 manual attached to our blog
set default PBRadii mbondi2
saveamberparm rec pro_solv.prmtop pro_solv.inpcrd #saving the solvated protein prmtop and inpcrd
savepdb rec pro_solv.pdb #saving solvated protein system as pdb file to visualise the system.
quit #exit from the tleap
The .in files for the upcoming commands are attached in the blog. These files are editable as per ur requirement.
Here the commands will be running sander and pmemd on multiple processors for faster output using mpirun feature (in centos) for minimisation, heating, density and equilibration steps. If you want to run on single CPU then remove mpi from the command: example:- $AMBERHOME/bin/sander -O -i min.in -o min.out -p comp_solv.prmtop -c comp_solv.inpcrd -r min.rst -ref comp_solv.inpcrd
#perform 1st Minimization ( lines in min.in can be changed as per ur requirement: maxcyc=2000, ncyc=1000] Here, the systems will be minimized for 2000 steps, out of which 1000 steps will be used by the steepest descent method followed by 1000 steps using the conjugate gradient to remove bad contacts in the systems.
mpirun -np 6 $AMBERHOME/bin/sander.MPI -O -i min.in -o min.out -p pro_solv.prmtop -c pro_solv.inpcrd -r min.rst -ref comp_solv.inpcrd
#output files: min.out and min.rst
#2nd Minimization [restraintmask=':1-254 & !@H='] please change the restrain mask as per the number solute amino acids.. example if ur protein contains 50 amino acids, then restraintmask=':1-50& !@H=' all non hydrogen atoms. The position restraints of 10 and 2 kcal mol-1 Å-2, respectively, are placed on the entire protein system to relax the solvent molecules in the first and second minimisation steps.
mpirun -np 6 $AMBERHOME/bin/sander.MPI -O -i min1.in -o min1.out -p pro_solv.prmtop -c min.rst -r min1.rst -ref min.rst
#output files: min1.out and min1.rst
#3rd Minimization [The third step of minimization will be carried out unrestrained. ]
mpirun -np 6 $AMBERHOME/bin/sander.MPI -O -i min2.in -o min2.out -p pro_solv.prmtop -c min1.rst -r min2.rst -ref min1.rst
#output files: min2.out and min2.rst
#Heating [The systems will be gradually heated for 20ps from 0K to 300K]..[nstlim=10000,
dt=0.002] the values in these can be changed if you want to run this step for more time. calculation is time of heating: nstlim*dt i.e. 10000*0.002 is 20ps.
mpirun -np 6 $AMBERHOME/bin/sander.MPI -O -i heat.in -o heat.out -p pro_solv.prmtop -c min2.rst -r heat.rst -ref min2.rst -x heat.nc
#output files: heat.out, heat.rst, heat.nc and/or mdinfo
#Density ( change the restraintmast values here in density.in file also like in min1.in file: restraintmask=':1-254',)
mpirun -np 6 $AMBERHOME/bin/sander.MPI -O -i density.in -o density.out -p comp_solv.prmtop -c heat.rst -r density.rst -x density.nc -ref heat.rst
#output files: density.out, density.rst, density.nc and/or mdinfo
#Equilibration [Equilibration of the system for 500ps at 300K and 1atm pressure.]
mpirun -np 6 $AMBERHOME/bin/pmemd.MPI -O -i equil.in -o equil.out -p pro_solv.prmtop -c density.rst -r equil.rst -x equil.nc
#output files: equil.out, equil.rst, equil.nc and/or mdinfo
#Production or MD simulation (the md.in file contains MD length of 1ns).
$AMBERHOME/bin/pmemd.cuda -O -i prod.in -o md.out -p comp_solv.prmtop -c equil.rst -r md1ns.rst -x md1ns.nc
#output files: md.out, md1ns.rst, md1ns.nc and/or mdinfo
.nc file is a NETCDF trajectory file
.out and mdinfo provides status of the running commands which can be seen in either text editors or by typing commands on terminal such as " tail -f mdinfo" [in another terminal to simultaneously view the status of the jobs.]
Some important calculations to remember:
- length of simulation: nstlim*dt [ i.e. Number of MD steps in run *Time step in picoseconds (ps)(the time length of each MD step) ]
- 1000ps =1ns
- to run for 100ns, change nstlim like : 50000000, then 0.002 (i.e. dt) is multiplied to get 100,000ps=100ns
- ntpr=5000, ntwx=5000: can be changed to write the trajectory
we have also provided a short automated script for running these minimisation, heating, density, equilibration and production commands step by step].. the script is job.sh... u can run it by command " sh job.sh. please don't forget to change the name of ur input and output files in this job.sh script before running it.
PS: to extend the simulations: please keep irest=1 in prod.in file, if u want to change the time, change the nstlim number, save]
use this command to extend the simulation from the last step of the md file:
$AMBERHOME/bin/pmemd.cuda -O -i prod.in -o md.out -p pro_solv.prmtop -c md1ns.rst -r md9ns.rst -x md9ns.nc
Min.in file:
&cntrl
imin=1,
ntx=1,
irest=0,
maxcyc=2000,
ncyc=1000,
ntpr=100,
ntwx=0,
cut=16.0,
/
References:
1. Amber 16 user manual
2. https://ambermd.org/tutorials/basic/tutorial2/section6.htm
3. http://ambermd.org/tutorials/basic/tutorial0/
Please provide your suggestions, queries, and feedback related to our blogs in the comments section below. Also do share and follow our blog for more updates related to Computational biology techniques.
"Spend less time on techniques, more time on analysis."
Thank you for the information .
ReplyDeletei hope it is helpful
ReplyDelete