Difference between revisions of "2024 AMBER tutorial 1 with PDBID 2ITO"
Stonybrook (talk | contribs) (→Production) |
Stonybrook (talk | contribs) (→RMSD) |
||
| Line 520: | Line 520: | ||
Now we want to look at the RMSD of the receptor. Include the residues that include only the receptor: | Now we want to look at the RMSD of the receptor. Include the residues that include only the receptor: | ||
| − | vi cpp_rmsd_receptor.in | + | vi cpp_rmsd_receptor.in |
with the following: | with the following: | ||
| Line 538: | Line 538: | ||
If you want to look at the histogram, you can graph it using python. | If you want to look at the histogram, you can graph it using python. | ||
| − | |||
==H-Bonds== | ==H-Bonds== | ||
Latest revision as of 02:23, 6 May 2024
Contents
Introduction
AMBER is a molecular dynamics program that simulate the behavior of a protein/ligand complex of interest. AMBER considers all atoms and models their interactions over time, continuously updating positions and forces in small time increments. For this tutorial, we will focus on 2ITO.
Getting Started
Update the directory structure to the following:
001_structure 002_forceFieldParameters 003_tleap 004_equilibration 005_productionRun 006a_RMSD 006b_hbonds 006c_mmgbsa
Structure
Similarly to the DOCK tutorial, we need to obtain the structure of the protein and ligand. For the protein, we want to fill in all loops even if it is not near the binding pocket, since AMBER considers all atoms in its calculations. Follow the DOCK tutorial for 2ITO to do this.
For the protein, do not add H or add charge. Save as 2ITO_protein_only_for_AMBER.pdb and 2ITO_protein_only_for_AMBER.mol2. For the ligand, add H and add charge. Save as 2ITO_ligand_only_for_AMBER.pdb and 2ITO_ligand_only_for_AMBER.mol2.
Force Field Parameters
Force field parameters are required to calculation the interactions between all atoms. This needs to be done for the ligand. Create the following slurm script as follows:
vi 2ITO_ffp.slurm
and have the following within the file:
Make sure to change LIG -nc 0 based off of your ligand. It is the charge on the ligand. For 2ITO, this charge is +1.
Once the slurm file has run, the following will show up in your directory:
2ITO_ligand_antechamber.mol2 ANTECHAMBER_AC.AC0 ANTECHAMBER_AM1BCC_PRE.AC ANTECHAMBER_BOND_TYPE.AC0 sqm.in sqm.pdb ANTECHAMBER_AC.AC ANTECHAMBER_AM1BCC.AC ANTECHAMBER_BOND_TYPE.AC ATOMTYPE.INF sqm.out
2ITO_ligand_antechamber.mol2 is the file with the parameters we generated.
In the command line, load amber 16:
module load amber/16
Run the following command to modify the parameters slightly:
parmchk2 -i 2ITO_ligand_antechamber.mol2 -f mol2 -o 2ITO_ligand.am1bcc.frcmod
You will see 2ITO_ligand.am1bcc.frcmod file in your directory.
TLeap
TLeap is a program with AMBER that will prepare the system for MD simulations. It will generate the gas-phase and solvated systems.
Create the file:
vi 2ITO_tleap.in
and add the following to the file:
Create the following slurm file:
vi 2ITO_tleap.slurm
and add the following to the file:
Run the slurm and you should get the following new files in your directory:
2ITO.complex.parm7 2ITO.gas.complex.rst7 2ITO.gas.ligand.rst7 2ITO.gas.receptor.rst7 leap.log 2ITO.gas.complex.pdb 2ITO.gas.ligand.parm7 2ITO.gas.receptor.parm7 2ITO.wet.complex.pdb 2ITO.wet.complex.prmtop 2ITO.wet.complex.rst7 tleap_output.txt
You can look at the files you created by moving the files to your local computer. Move 2ITO.wet.complex.rst7 and 2ITO.wet.complex.prmtop. Open chimera and follow the following commands:
Tools → MD/Ensemble Analysis → MD Movie
Click browse and put in the prmtop file. Add the rst7 trajectory file. You should be able to view your system.
Equilibration
Look at the 2ITO.gas.complex.pdb file that was just generated. Scroll all the way down to see the residues in the complex. In our example, the protein is residues 1-324 and the ligand is 325.
We need to equilibrate the system before running the MD simulations. Note, for files 1-7, set the restraint mask to the entire system (1-325). For files 8-9, set the restraint mask to cover only the receptor (1-324).
Create:
vi 01.min.mdin
With the following:
Minmize all the hydrogens &cntrl imin=1, ! Minimize the initial structure maxcyc=5000, ! Maximum number of cycles for minimization ntb=1, ! Constant volume ntp=0, ! No pressure scaling ntf=1, ! Complete force evaluation ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":1-325 & !@H=", ! atoms to be restrained restraint_wt=5.0, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format /
Create:
vi 02.equil.mdin
With the following:
MD simulation &cntrl imin=0, ! Perform MD nstlim=50000 ! Number of MD steps ntb=2, ! Constant Pressure ntc=1, ! No SHAKE on bonds between hydrogens dt=0.001, ! Timestep (ps) ntp=1, ! Isotropic pressure scaling barostat=1 ! Berendsen taup=0.5 ! Pressure relaxtion time (ps) ntf=1, ! Complete force evaluation ntt=3, ! Langevin thermostat gamma_ln=2.0 ! Collision Frequency for thermostat ig=-1, ! Random seed for thermostat temp0=298.15 ! Simulation temperature (K) ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":1-325 & !@H=", ! atoms to be restrained restraint_wt=5.0, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format iwrap=1, ! iwrap is turned on /
Create:
vi 03.min.mdin
With the following:
Minmize all the hydrogens &cntrl imin=1, ! Minimize the initial structure maxcyc=1000, ! Maximum number of cycles for minimization ntb=1, ! Constant volume ntp=0, ! No pressure scaling ntf=1, ! Complete force evaluation ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":1-325 & !@H=", ! atoms to be restrained restraint_wt=2.0, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format /
Create:
vi 04.min.mdin
With the following:
Minmize all the hydrogens &cntrl imin=1, ! Minimize the initial structure maxcyc=1000, ! Maximum number of cycles for minimization ntb=1, ! Constant volume ntp=0, ! No pressure scaling ntf=1, ! Complete force evaluation ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":1-325 & !@H=", ! atoms to be restrained restraint_wt=0.1, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format /
Create:
vi 05.min.mdin
With the following:
Minmize all the hydrogens &cntrl imin=1, ! Minimize the initial structure maxcyc=1000, ! Maximum number of cycles for minimization ntb=1, ! Constant volume ntp=0, ! No pressure scaling ntf=1, ! Complete force evaluation ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":1-325 & !@H=", ! atoms to be restrained restraint_wt=0.05, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format /
Create:
vi 06.equil.mdin
With the following:
MD simulation &cntrl imin=0, ! Perform MD nstlim=50000 ! Number of MD steps ntb=2, ! Constant Pressure ntc=1, ! No SHAKE on bonds between hydrogens dt=0.001, ! Timestep (ps) ntp=1, ! Isotropic pressure scaling barostat=1 ! Berendsen taup=0.5 ! Pressure relaxtion time (ps) ntf=1, ! Complete force evaluation ntt=3, ! Langevin thermostat gamma_ln=2.0 ! Collision Frequency for thermostat ig=-1, ! Random seed for thermostat temp0=298.15 ! Simulation temperature (K) ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":1-325 & !@H=", ! atoms to be restrained restraint_wt=1.0, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format iwrap=1, ! iwrap is turned on /
Create:
vi 07.equil.mdin
With the following:
MD simulation &cntrl imin=0, ! Perform MD nstlim=50000 ! Number of MD steps ntx=5, ! Positions and velocities read formatted irest=1, ! Restart calculation ntc=1, ! No SHAKE on for bonds with hydrogen dt=0.001, ! Timestep (ps) ntb=2, ! Constant Pressure ntp=1, ! Isotropic pressure scaling barostat=1 ! Berendsen taup=0.5 ! Pressure relaxtion time (ps) ntf=1, ! Complete force evaluation ntt=3, ! Langevin thermostat gamma_ln=2.0 ! Collision Frequency for thermostat ig=-1, ! Random seed for thermostat temp0=298.15 ! Simulation temperature (K) ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":1-325 & !@H=", ! atoms to be restrained restraint_wt=0.5, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format iwrap=1, ! iwrap is turned on /
Create:
vi 08.equil.mdin
With the following:
MD simulations &cntrl imin=0, ! Perform MD nstlim=50000 ! Number of MD steps ntx=5, ! Positions and velocities read formatted irest=1, ! Restart calculation ntc=1, ! No SHAKE on for bonds with hydrogen dt=0.001, ! Timestep (ps) ntb=2, ! Constant Pressure ntp=1, ! Isotropic pressure scaling barostat=1 ! Berendsen taup=0.5 ! Pressure relaxtion time (ps) ntf=1, ! Complete force evaluation ntt=3, ! Langevin thermostat gamma_ln=2.0 ! Collision Frequency for thermostat ig=-1, ! Random seed for thermostat temp0=298.15 ! Simulation temperature (K) ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":1-324@CA.C.N", ! atoms to be restrained, only the backbone restraint_wt=0.1, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format iwrap=1, ! iwrap is turned on /
Create:
vi 09.equil.mdin
With the following:
MD simulations &cntrl imin=0, ! Perform MD nstlim=50000 ! Number of MD steps ntx=5, ! Positions and velocities read formatted irest=1, ! Restart calculation ntc=1, ! No SHAKE on for bonds with hydrogen dt=0.001, ! Timestep (ps) ntb=2, ! Constant Pressure ntp=1, ! Isotropic pressure scaling barostat=1 ! Berendsen taup=0.5 ! Pressure relaxtion time (ps) ntf=1, ! Complete force evaluation ntt=3, ! Langevin thermostat gamma_ln=2.0 ! Collision Frequency for thermostat ig=-1, ! Random seed for thermostat temp0=298.15 ! Simulation temperature (K) ntwx= 1000, ! Write to trajectory file every ntwx steps ntpr= 1000, ! Print to mdout every ntpr steps ntwr= 1000, ! Write a restart file every ntwr steps cut= 8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":1-324@CA.C.N", ! atoms to be restrained restraint_wt=0.1, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format iwrap=1, ! iwrap is turned on /
Create the following slurm file with the following and run:
2ITO.equil.slurm
#!/bin/bash
#
#SBATCH --job-name=2ITO_equilibration
#SBATCH --output=equilibration_output.txt
#SBATCH --ntasks-per-node=24
#SBATCH --nodes=1
#SBATCH --time=48:00:00
#SBATCH -p long-24core
module load amber/16
do_parallel="mpirun -np 80 pmemd.MPI"
prmtop="../003_tleap/2ITO.wet.complex.prmtop"
coords="../003_tleap/2ITO.wet.complex"
MDINPUTS=(01.min 02.equil 03.min 04.min 05.min 06.equil 07.equil 08.equil 09.equil)
for input in ${MDINPUTS[@]}; do
$do_parallel -O -i ${input}.mdin -o ${input}.mdout -p $prmtop -c ${coords}.rst7 -ref ${coords}.rst7 -x ${input}.trj -inf ${input}.info -r ${input}.rst7
coords=$input
done
Production
The system is all equilibrated and now we can run production run.
vi 10.prod.mdin
and add the following:
MD simulations &cntrl imin=0, ! Perform MD nstlim=5000000, ! Number of MD steps ntx=5, ! Positions and velocities read formatted irest=1, ! Restart calculation ntc=2, ! SHAKE on for bonds with hydrogen dt=0.002, ! Timestep (ps) ntb=2, ! Constant Pressure ntp=1, ! Isotropic pressure scaling barostat=1 ! Berendsen taup=0.5 ! Pressure relaxtion time (ps) ntf=2, ! No force evaluation for bonds with hydrogen ntt=3, ! Langevin thermostat gamma_ln=2.0 ! Collision Frequency for thermostat ig=-1, ! Random seed for thermostat temp0=298.15 ! Simulation temperature (K) ntwx= 2500, ! Write to trajectory file every ntwx steps ntpr= 2500, ! Print to mdout every ntpr steps ntwr= 5000000, ! Write a restart file every ntwr steps cut=8.0, ! Nonbonded cutoff in Angstroms ntr=1, ! Turn on restraints restraintmask=":1-324@CA,C,N", ! atoms to be restrained restraint_wt=0.1, ! force constant for restraint ntxo=1, ! Write coordinate file in ASCII format ioutfm=0, ! Write trajectory file in ASCII format iwrap=1, ! iwrap is turned on $end
Create the slurm file:
production.slurm
and add the following:
#!/bin/bash
#
#SBATCH --job-name=2ITO_production
#SBATCH --output=2ITO_production.txt
#SBATCH --ntasks-per-node=28
#SBATCH --nodes=4
#SBATCH --time=48:00:00
#SBATCH -p long-28core
module load amber/16
do_parallel="mpirun pmemd.MPI"
prmtop="../003_tleap/2ITO.wet.complex.prmtop"
coords="../004_equil/09.equil"
MDINPUTS=(10.prod)
for input in ${MDINPUTS[@]}; do
$do_parallel -O -i ${input}.mdin -o ${input}.mdout -p $prmtop -c ${coords}.rst7 -ref ${coords}.rst7 -x ${input}.trj -inf ${input}.info -r ${input}.rst7
coords=$input
done
Run the slurm script
Analysis
Now we can analyze the simulation by looking at the RMSD, H-Bonds, and MM/GBSA results.
RMSD
First we want to strip all the water:
vi cpp_strip_water.in
with the following:
#!/usr/bin/sh parm ../003_tleap/2ITO.wet.complex.prmtop #read in trajectory trajin ../005_production/10.prod.trj #read in reference reference ../003_tleap/2ITO.wet.complex.rst7 #compute RMSD + align CA to crystal structure rmsd rms1 reference :1-325@CA #strip solvent strip :WAT:Na+:Cl- #create gas-phase trajectory trajout 2ITO_stripfit_restrained_gas.trj nobox cpp_strip_water.in (END)
Run it with the following command:
cpptraj -i cpp_strip_water.in
To look at the RMSD of the ligand, we only want to look at the residue that includes the ligand:
vi cpp_rmsd_lig.in
with the following:
#!/usr/bin/sh #trajin the trajectory trajin 2ITO_stripfit_restrained_gas.trj #read in reference reference ../003_tleap/2ITO.gas.complex.rst7 #compute RMSD (do not fit internal geometris first, included rigid body motions, convert frames to ns (framenum*.005) rmsd rms1 ":325&!(@H=)" nofit mass out 2ITO_lig_restrained_rmsd_nofit.dat time .005 #histogram the nofit rmsd histogram rms1,*,*,.1,* norm out 2ITO_lig_restrained_rmsd_nofit_histogram.dat
Run it with the following command:
cpptraj -p ../003_tleap/2ITO.complex.parm7 -i cpp_rmsd_lig.in
Now we want to look at the RMSD of the receptor. Include the residues that include only the receptor:
vi cpp_rmsd_receptor.in
with the following:
#!/usr/bin/sh #trajin the trajectory trajin 2ITO_stripfit_restrained_gas.trj #read in reference reference ../003_tleap/2ITO.gas.complex.rst7 #compute RMSD (do not fit internal geometries first, included rigid body motions, convert frames to ns (framenum*.005) rmsd rms1 ":1-324&!(@H=)" nofit mass out 2ITO_receptor_restrained_rmsd_nofit.dat time .005 #histogram the nofit rmsd histogram rms1,*,*,.1,* norm out 2ITO_receptor_restrained_rmsd_nofit_histogram.dat
Run it with the following command:
cpptraj -p ../003_tleap/2ITO.complex.parm7 -i cpp_rmsd_receptor.in
If you want to look at the histogram, you can graph it using python.
H-Bonds
To look at the hydrogen bonds that occur in the system during the simulation, make sure to include the receptor and ligand in the mask. Create the program:
vi cpp_hbonds.in
with the following:
!/usr/bin/sh #read in trajectory trajin ../004_production/10_prod.trj #wrap everything in one periodic cell - could cause problems, may comment out #autoimage if problems later autoimage #compute intra + water-mediated H-bonds hbond hb1 :1-325 out 2ITO_hbond.out avgout 2ITO_hbond_avg.dat solventdonor :WAT solventacceptor :WAT@-O nointramol bridgeout 2ITO_bridge-water.dat dist 3.0 angle 140
Run it with the following command:
cpptraj -p ../003_tleap/2ITO.wet.complex.prmtop -i cpp_hbonds.in
MM/GBSA
Now we will calcualte the free energy of the system. Create the program:
vi mmgbsa_2ITO.in
with the following:
mmgbsa 2ITO analysis &general interval=1, netcdf=1, keep_files=0, / &gb igb=8 saltcon=0.0, surften=0.0072, surfoff=0.0, molsurf=0, / &nmode drms=0.001, maxcyc=10000, nminterval=250, nmendframe=2000, nmode_igb=1, /
Create the slurm file:
vi 2ITO_mmgbsa.slurm
with the following:
#bin/bash
#SBATCH --time=2-00:00:00
#SBATCH --nodes=2
#SBATCH --ntasks=28
#SBATCH --job-name=2ITO_MMGBSA
#SBATCH --output=2ITO_MMGBSA.log
#SBATCH -p extended-28core
#Define topology files
solv_parm="../003_tleap/2ITO.wet.complex.prmtop"
complex_parm="../003_tleap/2ITO.complex.parm7"
receptor_parm="../003_tleap/2ITO.gas.receptor.parm7"
lig_parm="../003_tleap/2ITO.gas.ligand.parm7"
trajectory="../005_production/10_prod.trj"
MMPBSA.py -O -i mmgbsa.in \
-o 2ITO_mmgbsa_results.dat \
-eo 2ITO_mmgbsa_perframe.dat \
-sp ${solv_parm} \
-cp ${complex_parm} \
-rp ${receptor_parm} \
-lp ${lig_parm} \
-y ${trajectory}
and run the slurm file.
sbatch 2ITO_mmgbsa.slurm
