Difference between revisions of "2014 AMBER tutorial with HIV Protease"

From Rizzo_Lab
Jump to: navigation, search
(Created page with "For additional Rizzo Lab tutorials see AMBER Tutorials. In this tutorial, we will learn how to run a molecular dynamics simulation of a protein-ligand complex. We will t...")
 
(→‎tleap)
 
(102 intermediate revisions by 3 users not shown)
Line 1: Line 1:
 +
 
For additional Rizzo Lab tutorials see [[AMBER Tutorials]].
 
For additional Rizzo Lab tutorials see [[AMBER Tutorials]].
  
Line 6: Line 7:
  
 
===AMBER===
 
===AMBER===
−
[http://ambermd.org/ Amber] - '''A'''ssisted '''M'''odel '''B'''uilding with '''E'''nergy '''R'''efinement - is a suite of about multiple programs  for perform macromolecular simulations. Amber11, the current version of Amber, includes newly released functionality such as PMEMD, particle mesh Ewald MD and soft-core Thermodynamics Integration MD. For the tutorial, we are using the newest version AMBER12.
 
  
−
The [http://ambermd.org/doc12/Amber12.pdf Amber 12 Manual] is the primary resource to get started with Amber12. (Tip: Using Adobe Acrobat to view the file, you can simply search the document for keywords such as the name of a simulation parameter, which saves much time.) In addition, [http://ambermd.org/doc12/AmberTools12.pdf Amber Tools User's Manual] serves as another reference while using Amber tools.  
+
[http://ambermd.org/ Amber] - '''A'''ssisted '''M'''odel '''B'''uilding with '''E'''nergy '''R'''efinement - is a multi-program suite for macromolecular simulations. Amber12 is the most recent version of the software and it includes new force fields such as ff12SB and Lipid11, expanded options for Poisson-Boltzmann solvation calculations, accelerated molecular dynamics, additional features in ''sander'' pmemd code, and expanded methods for free energy calculations. Our lab is set up with Amber12 and the latest update of AmberTools13 which contains the programs such as antechamber and tleap to set up your simulation.
 +
 
 +
The [http://ambermd.org/doc12/Amber12.pdf Amber 12 Manual] is available to get started with using Amber12 as well as a troubleshooting reference. (Tip: You can search the document for keywords if you use Adobe Acrobat to view the file which can save time.) Additionally, [http://ambermd.org/doc12/AmberTools13.pdf Amber Tools User's Manual] is another reference for the probrgrams available under Amber tools.  
 +
 
 +
Here some of programs availabe in Amber and AmberTools
  
−
Here are some programs in Amber
+
#'''LEaP''': is a preparatory program for constructing new or modified systems in Amber. It contains the functions: prep, link, edit, and parm from earlier version of Amber.
 +
#'''ANTECHAMBER''': is the primary program used to prepare input files excluding the original nucleic acid and protein information taken from the PDB.
 +
#'''SANDER''': is 'a basic energy minimizer and molecular dynamics program' that is used to minimize, equilibrate and sample molecular conformations. This is the program used to generate trajectory files of the system.
 +
#'''PMEMD''': is an improved version of SANDER optimized for periodic, PME simulations, and for GB simulations. It is faster and scales better on parallel machines than SANDER.
 +
#'''PTRAJ''': is an analysis program for processing trajectory files. ptraj can be used to rotate and translate the structures, evaluate geometrical features, calculate RMSDs among other analyses.
  
−
#'''LEaP''': an preparing program for constructing new or modified systems in Amber. It consists of the functions of prep, link, edit, and parm for earlier version of Amber.
+
There is a [http://lists.ambermd.org/mailman/listinfo/amber mailing list] available as an additional resource.
−
#'''ANTECHAMBER''': in additional to LEap, this main Antechamber suite program is for preparing input files other than standard nucleic acids and proteins.
 
−
#'''SANDER''': according to the Amber 12 manual, it is 'a basic energy minimizer and molecular dynamics program' that can be used to minimize, equilibrate and sample molecular conformations. And this is the program we mainly use in this tutorial to generate trajectory files of the molecular system.
 
−
#'''PMEMD''': version of SANDER that has improved parallel scaling property and optimized speed.
 
−
#'''PTRAJ''': an analysis program for processing trajectory files. One can use ptraj to rotate, translate the structures, evaluate geometrical features and so on.  
 
  
−
There is a [http://lists.ambermd.org/mailman/listinfo/amber mailing list] you could sign-up for, as an additional resource.
+
===HIV Protease===
  
−
===UMP and OMP===
+
Human Immunodeficiency Virus (HIV) is a retrovirus which causes Acquired Immune Deficiency Syndrome (AIDS).  HIV protease is a aspartyl protease which is a dimer, each monomer consisting of 99 amino acids.  This protein is essential for the lifecycle of the virus as it cleaves immature polyproteins into mature proteins which can then be used for proper virion assembly.  HIV protease has been validated as a drug target for the treatment for HIV as inhibition of this protein can slow HIV progression. However, retroviruses have high mutation rates and drug-resistance is a constant problem which needs to be addressed. There are a number of structure of HIV protease in the protein databank, we will use 1HVR for this tutorial.
−
For information of the UMP-OMP system, see [[2013 DOCK tutorial with Orotodine Monophosphate Decarboxylase]].
 
  
 
===Organizing Directories===
 
===Organizing Directories===
  
−
While performing MD simulations, it is convenient to adopt a standard directory structure / naming scheme, so that files are easy to find / identify. For this tutorial, we will use something similar to the following:
+
The following directory structure and naming scheme is a convienient way to organize your files for this tutorial. Create these directories first before doing anything else.
  
−
  ~username/AMS536/AMBER-Tutorial/001.CHIMERA.MOL.PREP/   
+
  ~username/AMS536/AMBER-Tutorial/001.MOL.PREP/   
 
                                 002.ANTE.TLEAP/  
 
                                 002.ANTE.TLEAP/  
−
                                 003.SANDER/       
+
                                 003.PMEMD/       
 
                                 004.PTRAJ/
 
                                 004.PTRAJ/
 
                                 005.MMGBSA/
 
                                 005.MMGBSA/
Line 35: Line 38:
 
==II. Structural Preparation==
 
==II. Structural Preparation==
  
−
===Preparation in Chimera===
+
====Downloading the PDB file====
  
−
In this AMBER tutorial, we will use the same system with previous DOCK part.  To begin with, we need three files under directory 001.CHIMERA.MOL.PREP.
+
Go to PDB homepage (http://www.rcsb.org/pdb/home/home.do ) enter the protein ID (1HVR) in the search bar, click Download Files in the top-right of the webpage, then select PDB File (text). In the new window, save the file.
−
 
 
−
====Downloading the PDB file (1LOQ)====
 
−
 
 
−
Since we need to edit the PDB before we use it in Chimera we should do this manually. Go to PDB homepage (http://www.rcsb.org/pdb/home/home.do ) enter the protein ID (1LOQ) in the search bar, click Download Files in the top-right of the webpage, then select PDB File (text). In the new window, save the file in Downloads.  
 
  
 
====Preparing the ligand and receptor in Chimera====
 
====Preparing the ligand and receptor in Chimera====
  
−
In this section, we will create three new files and save them in the 001.CHIMERA.MOL.PREP/ folder:
+
In this section, three new files would be created in the 001.CHIMERA.MOL.PREP/ folder:
  
−
   1LOQ.dockprep.mol2      (complete system with hydrogens and charges)
+
   1HVR.dockprep.mol2      (complete system with hydrogens and charges)
−
   1LOQ.receptor.noH.mol2 (the receptor alone, without hydrogens)
+
   1HVR.receptor.noH.pdb (the receptor alone, without hydrogens)
−
   1LOQ.ligand.mol2        (the ligand alone)
+
   1HVR.ligand.mol2        (the ligand alone)
  
−
To prepare these files, first copy the original PDB file into the 001.CHIMERA.MOL.PREP/ folder and open it with VIM ($ vim 1LOQ.pdb). Because the residue name of the ligand (U) will give us some problems when assigning charges, change the residue name "U" to "LIG" starting at line 2082. Here is an example command that will change all instances of "  U" into "LIG", while preserving the correct spacing:
+
To prepare these files, first copy the original PDB file into '''001.MOL.PREP''' folder and open it with VIM ($ vim 1HVR.pdb). Because the residue name of the ligand (U) will give us some problems when assigning charges, change the residue name "U" to "LIG". Here is an example command that will change all instances of "  U" into "LIG", while preserving the correct spacing (Note: there are two spaces before U):
  
 
   :%s/  U/LIG/gc
 
   :%s/  U/LIG/gc
Line 57: Line 56:
 
For this command, '''g'''  is short for global and '''c''' is short for check with the user before making the change.
 
For this command, '''g'''  is short for global and '''c''' is short for check with the user before making the change.
  
−
Next, open up the PDB file (1LOQ.pdb) in Chimera. To delete water molecules and other ligands, click ''Tools'' -> ''Structure Editing'' -> ''Dock Prep''. Check all boxes and click ''Okay'' to the end. Alternatively, the waters can be deleted manually by choosing ''Select'' -> ''Residue'' -> ''HOH'', then go to ''Actions'' -> ''Atoms/Bonds'' -> ''Delete''. Hydrogen atoms can be added manually by choosing ''Tools'' -> ''Structure Editing'' -> ''Add H''.
+
Next, open up the PDB file (1HVR.pdb) in Chimera. To generate '''1HVR.dockprep.mol2''' files, delete the water molecules; delete the original hydrogen atoms; add the charge and add the hydrogen atoms.  
  
−
Next, to add charges to the ligand and receptor,  go to ''Select'' -> ''Residue'' -> ''LIG'', then go to ''Tools'' -> ''Structure Editing'' -> ''Add Charge''. Choose ''AMBER ff99SB'' as the charge model, click ''Okay'', and when prompted chose AM1-BCC charges for the ligand, and make sure the ''Net Charge'' is set to ''-1''. (You must consider the chemistry of the ligand when assigning a charge state).
+
To delete water molecules and other ligands, click '''Select -> Residue -> HOH''', then go to '''Actions -> Atoms/Bonds -> Delete'''. Hydrogen atoms can be added manually by choosing '''Tools -> Structure Editing -> Add H'''. To add charges to the ligand and receptor,  go to '''Select -> Residue -> LIG''', then go to '''Tools -> Structure Editing -> Add Charge'''. Choose '''AMBER ff99SB''' as the charge model, click '''Okay''', and when prompted chose AM1-BCC charges for the ligand, and make sure the '''Net Charge''' is set to '''0'''. (The chemistry of the ligand should be considered when assigning a charge state).
  
−
Finally, save this file as '''1LOQ.dockprep.mol2'''.
+
Alternatively, you can do all of the above by clicking '''Tools -> Structure Editing -> Dock Prep'''. Note when adding the charge to the ligand, you can choose AMBER ff99SB as the charge model and chose gasteiger as the charge method. In this 1HVR case, we set Net Charge to 0.
 +
 
 +
Finally, save this file as '''1HVR.dockprep.mol2'''.
  
 
====Generating the final files====
 
====Generating the final files====
  
−
To create the receptor file with no hydrogen atoms: Open 1LOQ.dockprep.mol2, click ''Select'' -> ''Chemistry'' -> ''Element'' -> ''H'', then chose ''Actions'' -> ''Atoms/Bonds'' -> ''Delete''. Save the file as '''1LOQ.receptor.noH.mol2'''.
+
To create the receptor file with no hydrogen atoms: Open 1HVR.dockprep.mol2, click '''Select -> Chemistry -> Element -> H'', then chose '''Actions -> Atoms/Bonds -> Delete'''. Save the file as '''1HVR.receptor.noH.pdb'''.
  
−
To create the ligand file: Open 1LOQ.dockprep.mol2, click ''Select'' -> ''Residue'' -> ''LIG'', then click ''Select'' -> ''Invert'', then chose ''Actions'' -> ''Atoms/Bonds'' -> ''Delete''. Save the file as '''1LOQ.ligand.mol2'''.
+
To create the ligand file: Open 1HVR.dockprep.mol2, click '''Select -> Residue -> LIG''', then click '''Select -> Invert''', then chose '''Actions -> Atoms/Bonds -> Delete'''. Save the file as '''1HVR.ligand.mol2'''.
  
 
===antechamber===
 
===antechamber===
−
+
An antichamber input file need all the atom names to be unique and only the first three characters will be used. For example, H301 and H302 can not be distinguished with each other. So first we need to manually rename the atoms which have the same first three characters.
−
The antechamber program itself is the main program of Antechamber package. In most of cases, one should use this program instead of a series of separated programs to do molecular format convertion, atom type assignment and charge generation etc. An antechamber input file requires all the atom names to be unique. So if we use 1LOQ.ligand.mol2 as the input file, it will cause errors. The program can only recognize atom names of 3 characters ( In this case, H5' and <nowiki>H5'' </nowiki> cannot be distinguished from each other. )
+
         [[File:007.png]]
−
 
+
The file after renaming:
−
    22 H5'        40.0697  36.0506  37.6716 H        1 LIG210      0.0761
+
         [[File:009.png]]
−
    23 H5''        40.6060  37.0416  36.2883 H        1 LIG210      0.0349
 
−
    24 H4'        38.2510  37.5673  38.1082 H        1 LIG210      0.0967
 
−
    25 H3'        38.2723  38.5564  35.5613 H        1 LIG210      0.1052
 
−
    26 H2'        40.3056  39.7322  35.6920 H        1 LIG210      0.1188
 
−
    27 H1'        39.7587  40.5927  38.4795 H        1 LIG210      0.0753
 
−
    28 H6        41.3847  41.7815  39.3543 H        1 LIG210      0.0322
 
−
    29 H5        43.6843  42.3632  39.6188 H        1 LIG210      0.1693
 
−
    30 H3        44.6515  39.3570  36.8437 H        1 LIG210      0.3481
 
−
    31 HO3'      36.5607  38.9743  36.9843 H        1 LIG210      0.4373
 
−
    32 HP3        39.3989  32.7036  35.9199 H        1 LIG210      0.4245
 
−
    33 HO2'      39.5501  41.7087  35.4449 H        1 LIG210      0.4228
 
−
 
 
−
We need to manually rename the atoms.  
 
−
 
 
−
    22 H1         40.0697  36.0506  37.6716 H        1 LIG210      0.0761
 
−
    23 H2        40.6060  37.0416  36.2883 H        1 LIG210      0.0349
 
−
    24 H3         38.2510  37.5673  38.1082 H        1 LIG210      0.0967
 
−
    25 H4        38.2723  38.5564  35.5613 H        1 LIG210      0.1052
 
−
    26 H5        40.3056  39.7322  35.6920 H        1 LIG210      0.1188
 
−
    27 H6        39.7587  40.5927  38.4795 H        1 LIG210      0.0753
 
−
    28 H7        41.3847  41.7815  39.3543 H        1 LIG210      0.0322
 
−
    29 H8        43.6843  42.3632  39.6188 H        1 LIG210      0.1693
 
−
    30 H9        44.6515  39.3570  36.8437 H        1 LIG210      0.3481
 
−
    31 HO1        36.5607  38.9743  36.9843 H        1 LIG210      0.4373
 
−
    32 HP3        39.3989  32.7036  35.9199 H        1 LIG210      0.4245
 
−
    33 HO2        39.5501  41.7087  35.4449 H        1 LIG210      0.4228
 
−
 
 
−
To begin with, go to 002.ANTE.TLEAP directory.
 
−
To make sure we have access to the three programs that we want to run (antechamber, parmchk and tleap) and we are using the correct version of amber, we can use the which command, type:
 
−
 
 
−
which antechamber
 
−
which parmchk
 
−
which tleap
 
−
 
 
−
Your results should be similar to this:
 
−
 
−
/home/wjallen/AMS536/local/amber12/bin/antechamber
 
−
/home/wjallen/AMS536/local/amber12/bin/parmchk
 
−
/home/wjallen/AMS536/local/amber12/bin/tleap
 
  
 +
To run antichamber, we need make sure these three programs are available, (antechamber, parmchk and tleap) and we are using the correct version of amber, we can use the which command, type command:
 +
  which antechamber
 +
  which parmchk
 +
  which tleap
 +
The result should look like this:
 +
  /nfs/user03/wjallen/local/amber12/bin/antechamber
 +
  /nfs/user03/wjallen/local/amber12/bin/parmchk
 +
  /nfs/user03/wjallen/local/amber12/bin/tleap
 
Copy parameters of ions to your working directory from the following resource:
 
Copy parameters of ions to your working directory from the following resource:
−
 
+
  scp -r ~yuchzhou/AMS536/AMBER_Tutorial/002.ANTE.TLEAP/rizzo_amber7.ionparms ./
−
cp -r ~lingling/AMS536/AMBER_Tutorial/002.ANTE.TLEAP/rizzo_amber7.ionparms
 
−
 
 
 
Then we use antechamber to convert our input mol2 file into files ready for LEaP.Type command:
 
Then we use antechamber to convert our input mol2 file into files ready for LEaP.Type command:
  
−
antechamber -i ../001.CHIMERA.MOL.PREP/1LOQ.lig.mol2 -fi mol2 -o 1LOQ.lig.ante.pdb -fo pdb
+
  antechamber -i ../001.MOL.PREP/1HVR.lig.mol2 -fi mol2 -o 1HVR.lig.ante.pdb -fo pdb
−
 
+
Here, -i means input file name; -fi means input file format; -o means output file name; -fo output file format. You will have an output file:1HVR.lig.ante.pdb
−
Here, -i input file name; -fi input file format; -o output file name; -fo output file format. You will have an output file:1LOQ.lig.ante.pdb
+
Similarly, we can use antechamber to change the format of 1HVR.lig.mol2 file to prep file:
−
 
+
  antechamber -i ../001.MOL.PREP/1HVR.lig.mol2 -fi mol2 -o 1HVR.lig.ante.prep -fo prepi
−
Similarly, we can use antechamber to change the fomat of 1LOQ.lig.mol2 file to prep file:
 
−
 
 
−
antechamber -i ../001.CHIMERA.MOL.PREP/1LOQ.lig.mol2 -fi mol2 -o 1LOQ.lig.ante.prep -fo prepi
 
−
 
 
 
You will get a set of output files:
 
You will get a set of output files:
−
 
+
  ANTECHAMBER_AC.AC   
−
ANTECHAMBER_AC.AC   
+
  ANTECHAMBER_AC.AC0   
−
ANTECHAMBER_AC.AC0   
+
  ANTECHAMBER_BOND_TYPE.AC   
−
ANTECHAMBER_BOND_TYPE.AC   
+
  ANTECHAMBER_BOND_TYPE.AC0   
−
ANTECHAMBER_BOND_TYPE.AC0   
+
  ANTECHAMBER_PREP.AC   
−
ANTECHAMBER_PREP.AC   
+
  ANTECHAMBER_PREP.AC0
−
ANTECHAMBER_PREP.AC0
+
  ATOMTYPE.INF  
−
ATOMTYPE.INF  
+
  NEWPDB.PDB  
−
NEWPDB.PDB  
+
  PREP.INF  
−
PREP.INF  
+
  1HVR.lig.ante.prep
−
1LOQ.lig.ante.prep
 
  
 
===parmchk===
 
===parmchk===
  
−
Parmchk is another program in Antechamber. After running antechamber, we run parmchk to check the parameters. If there are missing parameters after antechamber is finished, the frcmod template generated by parmchk will help us in generating the needed parameters:
+
After running antechamber, we run parmchk to check the parameters. Parmchk is another program in Antechamber. If there are missing parameters or mistake after antechamber is finished, the frcmod template generated by parmchk will help us in generating the needed parameters:
 +
  parmchk -i 1HVR.lig.ante.prep -f prepi -o 1HVR.lig.ante.frcmod
 +
In the command,  -i input file name,-f input file format (prepi, prepc, ac, mol2),-o frcmod file name. When the programming is done, you should have a new file in the directory.
 +
  1HVR.lig.ante.frcmod
  
−
  parmchk -i 1LOQ.lig.ante.prep -f prepi -o 1LOQ.lig.ante.frcmod
+
===tleap===
 +
Then, three input files are needed to run TLEAP, they are “tleap.lig.in”, “tleap.rec.in” and “tleap.com.in”, for the ligand, receptor, and complex respectively. These input files will be used by tleaP to create parameter/topology files and initial coordinate files.
 +
Copy parameters to your working directory from the following resource:
 +
  scp ~yuchzhou/AMS536/AMBER_Tutorial/002.ANTE.TLEAP/tleap.lig.in
 +
  scp  ~yuchzhou/AMS536/AMBER_Tutorial/002.ANTE.TLEAP/tleap.rec.in
 +
  scp  ~yuchzhou/AMS536/AMBER_Tutorial/002.ANTE.TLEAP/tleap.com.in
  
−
===tleaP===
+
here are what the files should look like. 
 +
* tleap.lig.in
  
−
Next, we need 3 new input files – “tleap.lig.in”, “tleap.rec.in” and “tleap.com.in”, for the ligand, the receptor, and the complex, respectively. These input files will be used by LEaP to create parameter/topology files and initial coordinate files.
+
set default PBradii mbondi2                                      #set default PBradii                                       
 +
source leaprc.ff99SB                                            #load a force field
 +
loadoff rizzo_amber7.ionparms/ions.lib                          #Loading Rob Rizzo's ion parameters
 +
loadamberparams rizzo_amber7.ionparms/ions.frcmod                #load an AMBER format parameter file
 +
source leaprc.gaff
 +
loadamberparams rizzo_amber7.ionparms/gaff.dat.rizzo
 +
loadamberparams 1HVR.lig.ante.frcmod
 +
loadamberprep 1HVR.lig.ante.prep                                #load an AMBER PREP file
 +
lig = loadpdb 1HVR.lig.ante.pdb                                  #load the ligand pdb file
 +
saveamberparm lig 1HVR.lig.gas.leap.prm7 1HVR.lig.gas.leap.rst7  #save the ligand gas phase AMBER topology and coordinate files
 +
solvateBox lig TIP3PBOX 10.0                                    #solvate the receptor using TIP3P, solvent box radii 10 angstroms
 +
saveamberparm lig 1HVR.lig.wat.leap.prm7 1HVR.lig.wat.leap.rst7  #save the ligand water phase AMBER topology and coordinate files
 +
charge lig                                                      #this writes out the charge of the ligand
 +
quit
  
−
Copy parameters of ions to your working directory from the following resource:
+
Coments (text after the '#' symbol) are for edification and should not apper in the file.
  
−
  cp ~lingling/AMS536/AMBER_Tutorial/002.ANTE.TLEAP/tleap.lig.in
+
* tleap.rec.in
−
  cp ~lingling/AMS536/AMBER_Tutorial/002.ANTE.TLEAP/tleap.rec.in
+
  set default PBradii mbondi2
−
  cp ~lingling/AMS536/AMBER_Tutorial/002.ANTE.TLEAP/tleap.com.in
+
  source leaprc.ff99SB
 +
loadoff rizzo_amber7.ionparms/ions.lib
 +
loadamberparams rizzo_amber7.ionparms/ions.frcmod
 +
REC = loadpdb ../001.MOL.PREP/1HVR.rec.noh.pdb
 +
saveamberparm REC 1HVR.rec.gas.leap.prm7 1HVR.rec.gas.leap.rst7
 +
solvateBox REC TIP3PBOX 10.0                                   
 +
saveamberparm REC 1HVR.rec.wat.leap.prm7 1HVR.rec.wat.leap.rst7
 +
charge REC
 +
quit
 +
* tleap.com.in
 +
  set default PBradii mbondi2
 +
source leaprc.ff99SB
 +
  loadoff rizzo_amber7.ionparms/ions.lib
 +
loadamberparams rizzo_amber7.ionparms/ions.frcmod
 +
source leaprc.gaff
 +
loadamberparams rizzo_amber7.ionparms/gaff.dat.rizzo
 +
REC = loadpdb ../001.MOL.PREP/1HVR.rec.noh.pdb
 +
loadamberparams 1HVR.lig.ante.frcmod
 +
loadamberprep 1HVR.lig.ante.prep
 +
LIG = loadpdb 1HVR.lig.ante.pdb
 +
  COM = combine {REC LIG}
 +
  saveamberparm COM 1HVR.com.gas.leap.prm7 1HVR.com.gas.leap.rst7
 +
solvateBox COM TIP3PBOX 10.0
 +
saveamberparm COM 1HVR.com.wat.leap.prm7 1HVR.com.wat.leap.rst7
 +
charge COM
 +
quit
  
 
Use leap command:
 
Use leap command:
−
tleap -s -f tleap.lig.in > tleap.lig.out
+
  tleap -s -f tleap.lig.in > tleap.lig.out
−
tleap -s -f tleap.rec.in > tleap.rec.out
+
  tleap -s -f tleap.rec.in > tleap.rec.out
−
tleap -s -f tleap.com.in > tleap.com.out
+
  tleap -s -f tleap.com.in > tleap.com.out
 +
The following output files will be generated:
 +
  1HVR.com.gas.leap.prm7
 +
  1HVR.com.wat.leap.prm7
 +
  1HVR.lig.ante.frcmod
 +
  1HVR.lig.ante.prep
 +
  1HVR.lig.gas.leap.rst7
 +
  1HVR.lig.wat.leap.rst7
 +
  1HVR.rec.gas.leap.rst7
 +
  1HVR.rec.wat.leap.rst7
 +
  1HVR.com.gas.leap.rst7
 +
  1HVR.com.wat.leap.rst7
 +
  1HVR.lig.ante.pdb
 +
  1HVR.lig.gas.leap.prm7
 +
  1HVR.lig.wat.leap.prm7
 +
  1HVR.rec.gas.leap.prm7
 +
  1HVR.rec.wat.leap.prm7
 +
  tleap.lig.out
 +
  tleap.rec.out
 +
  tleap.com.out
 +
  leap.log
  
−
You will see the following output files:
+
===Visualization in VMD===
−
1LOQ.com.gas.leap.prm7
 
−
1LOQ.com.wat.leap.prm7
 
−
1LOQ.lig.ante.frcmod 
 
−
1LOQ.lig.ante.prep     
 
−
1LOQ.lig.gas.leap.rst7 
 
−
1LOQ.lig.wat.leap.rst7 
 
−
1LOQ.rec.gas.leap.rst7 
 
−
1LOQ.rec.wat.leap.rst7
 
−
1LOQ.com.gas.leap.rst7
 
−
1LOQ.com.wat.leap.rst7
 
−
1LOQ.lig.ante.pdb   
 
−
1LOQ.lig.gas.leap.prm7 
 
−
1LOQ.lig.wat.leap.prm7 
 
−
1LOQ.rec.gas.leap.prm7 
 
−
1LOQ.rec.wat.leap.prm7
 
−
tleap.lig.out
 
−
tleap.rec.out
 
−
tleap.com.out
 
−
leap.log
 
  
−
===Visualization in VMD===
+
Visualization is an important step in AMBER molecular dynamics simulation as it allows for the viewing of molecules and molecule movements within a specified field of view. Several files prepared in the previous two steps will be required for the visualization of the ligand and its movement in the protein binding pocket. The first step that must be completed is the copying of all necessary information from Seawulf to Herbie.
 +
  scp -r ~/AMS536/AMBER_Tutorial/002.ANTE.TLEAP username@herbie.mathlab.sunysb.edu:~/AMS536/AMBER_Tutorial/
 +
We can now visualize several files in herbie by VMD: the protein-ligand complex in the gas or water phase, and the ligand in the gas or water phase.In the "VMD Main" window, click File > New Molecule. Now in the new window titled "Molecule File Browser" a few files must be selected and loaded.Viewing the protein-ligand complex in the gas phase we select 1LOQ.com.gas.leap.prm7 and the file type AMBER7 Parm. This would allow you to load the parameter information of the milecule, but nothing will shoe up in the VMD visualization window because coordinates has not been loaded yet.Next you need to select and load the file 1LOQ.com.gas.leap.rst7 and the file type AMBER7 Restart. This can be done for the complex or ligand for either the gas or the water phase by selecting and loading the corresponding .parm7 files and .rst7 files.
  
−
Visualization is an important step in AMBER molecular dynamics simulation as it allows for the viewing of molecules and molecule movements within a specified field of view. Several files prepared in the previous two steps will be required for the visualization of the ligand and its movement in the protein binding pocket.  The first step that must be completed is the copying of all necessary information from Seawulf to Herbie.
+
Next, you can choose to edit the structure shown for better visualization. By clicking Graphics > Representations... a new window "Graphical Representations" will open, in which you can create new representations of the protein-ligand complex. It is possible to choose the entire complex or simply the protein or ligand, to color the structure in various ways, and to choose how it is best represented (e.g. lines, thicker bonds, full surface).
  
−
Before exiting Seawulf, you can type:
+
[[File:protein-ligand complex in the gas visualization.png]]
−
scp -r ~/AMS536/AMBER_Tutorial/002.ANTE.TLEAP username@herbie.mathlab.sunysb.edu:~/AMS536/AMBER_Tutorial/
 
  
−
Then, when in Herbie, type the command '''vmd''' in any directory and the program VMD will open on the computer. In the "VMD Main" window, click '''File > New Molecule'''.  Now in the new window titled "Molecule File Browser" a few files must be selected and loaded. To select and load a file follow these steps:
+
==III. Simulation using sander==
−
1. Click '''Browse...''' and select the file that you desire.
+
Junjie, Tianao and Kai
−
2. Click the down button in the "Determine file type:" field and select the proper file type.
+
=== Minimization ===
−
3. Click '''Load''' and view the molecule/system in the original "VMD 1.8.5 OpenGL Display" window.
+
Before running any equilibration and production, we have to minimize the energy of system. This process can remove any unfavorable interaction between atoms due to the low resolution of biomolecule. If an unfavorable interaction is present when doing equilibration and production, atoms may "fly apart" and "crash" into other part of the biomolecule.  
  
−
We can now visualize several files: the protein-ligand complex in the gas or water phase, and the ligand in the gas or water phase. Viewing the protein-ligand complex in the gas phase we select '''1LOQ.com.gas.leap.prm7''' and the file type '''AMBER7 Parm'''. Next you need to select and load the file '''1LOQ.com.gas.leap.rst7''' and the file type '''AMBER7 Restart'''.  This can be done for the complex or ligand for either the gas or the water phase by selecting and loading the corresponding '''.parm7''' files and '''.rst7''' files.
+
An example of a typical input file “01mi.in” for minimization is listed blow:
 +
  01mi.in: minimization
 +
  &cntrl
 +
  imin = 1, maxcyc = 1000, ntmin = 2,
 +
  ntx = 1, ntc = 1, ntf = 1,
 +
  ntb = 1, ntp = 0,
 +
  ntwx = 1000, ntwe = 0, ntpr = 1000,
 +
  cut = 8.0,
 +
  ntr = 1,
 +
  restraintmask = ':1-199 & !@H=',
 +
  restraint_wt = 5.0,
 +
  /
 +
You may have to minimize the system several times. Each time you release the constraint on heavy atoms slowly.
  
−
Next, you can choose to edit the structure shown for better visualization.  By clicking '''Graphics > Representations...''' a new window "Graphical Representations" will open, in which you can create new representations of the protein-ligand complex. It is possible to choose the entire complex or simply the protein or ligand, to color the structure in various ways, and to choose how it is best represented (e.g. lines, thicker bonds, full surface).
+
The meaning of each flag is:
  
−
==III. Simulation using sander==
+
imin = 1: do minimization only.
  
−
=== Minimization ===
+
maxcyc: cycles of minimization will be performed.
−
Energy minimization is first performed on the stucture before the equilibration and production runs may be performed. Model building often creates unwanted structural artifacts that must be removed before a molecular dynamics simulation is performed.  
 
  
−
To begin, create four files for minimization steps 1,3,4 and 5.  
+
ntmin: specify the number of cycles that use a deepest decent method, after that minimization will use the conjugate gradient method.
  
−
vim 01mi.in
+
ntx = 1: read coordinates only from rst file. There is no need to read velocity before minimization.
  
−
and follow the naming according to minimization run number, i.e. 03mi.in
+
ntc = 1: Minimization is performed on every bond.
  
−
All "#mi.in" file content will be identical except the last parameter, the restraint weight (restraint_wt). This restraint will decrease with increasing minimization runs. Run number 1,3,4,and 5 has restraint weight 5,2,0.1,and 0.05 respectively.
+
ntf = 1: Typically, ntf = ntc.
  
 +
ntb: periodic boundary is applied.
  
−
An example of minimization code for 01mi.in:
+
ntp = 0: NO constant pressure is applied.
  
−
01mi.in: equilibration
+
ntwx = 1000: saving snapshot to a trajectory file every ntwx steps.
−
&cntrl
 
−
  imin = 1, maxcyc = 1000, ntmin = 2,
 
−
  ntx = 1, ntc = 1, ntf = 1,
 
−
  ntb = 1, ntp = 0,
 
−
  ntwx = 1000, ntwe = 0, ntpr = 1000,
 
−
  ntr = 1,
 
−
  restraintmask = ':1-210 & !@H=',
 
−
  restraint_wt=5.0,
 
  
 +
ntwe = 0: human readable energy file is appended every ntwe steps.
  
−
The following parameters are defined as follows:
+
ntpr = 1000: newest calculated energy is stored in .info file every ntpr steps.
  
−
'''imin''':specifies minimization not molecular dynamics
+
cut: cutoff used in energy evaluation.
  
−
'''maxcyc''':total number of minimization cycles to be performed
+
ntr = 1: positional restraint method is applied.
  
−
'''ntmin''':how many cycles will use a deepest decent method, the remaining cycles use an approximation of this called the conjugate gradient method.
+
restraintmask = ':1-199 & !@H : atoms of residue 1-199 except H is restrained.
  
−
'''ntx''':only coordinates and not velocities are to be read from previous step
+
restraint_wt: restraint strength in unit of kcal/(mol*Angstrom^2).
  
−
'''ntc''':indicates level of constraint on bonds. if =1, SHAKE algorithm is off so no bonds are constrained. If =2, constrains any bonds with H atoms. If =3, constrains all bonds.
+
=== Equilibration ===
 +
System need to be heat up in a every careful way. An abrupt heating up may cause unwanted perturbation.  
  
−
'''ntf''':=1, all parts of the potential must be evaluated
+
An example of the content of "02md.in" is listed below:
 +
  02md.in: equilibration
 +
  &cntrl
 +
  imin = 0, ntx = 1, irest = 0, nstlim = 50000,
 +
  temp0 = 298.15, tempi = 298.15, ig = 71287,
 +
  ntc = 2, ntf = 1, ntt = 1, dt = 0.001,
 +
  ntb = 2, ntp = 1, tautp = 1.0, taup = 1.0,
 +
  ntwx = 1000, ntwe = 0, ntwr = 1000, ntpr = 1000,
 +
  cut = 8.0, iwrap = 1,
 +
  ntr = 1, nscm = 100,
 +
  restraintmask = ':1-199 & !@H=', restraint_wt = 5.0,
 +
  /
  
−
'''ntb''':periodic boundary to keep system at constant volume
+
The parameters are described as below:
  
−
'''ntp''':=0, NO constant pressure applied
+
imin = 0: now md simulation is performed instead of minimization.
  
−
The frequencies at which the program records data are in controlled by the paramenters ntwx, ntwe, and ntpr.
+
ntc = 2: SHAKE is applied to constraint bonds involve H.
  
−
'''ntwx''':=1000, atom coordinates saved into .trj file every 1000 cycles
+
ntb = 2 and ntp = 1: constant pressure is applied.
  
−
'''ntwe''':=0, no .en file is generated
+
irest = 0: we are not reading in any velocities so there is no need to use restart.
  
−
'''ntpr''':=1000, energy readins are written as .out and .info files every 1000 steps
+
ntt = 1: weak-coupling is applied to keep the temperature constant during equilibration.
  
−
'''ntr''':=1, positional restraint method applied
+
tautp = 1.0 and taup = 1.0: time constant in unit of ps and the pressure relaxation time in unit of ps.
  
−
'''restraintmask'''= ':1-119 & !@H : position of atom within residues 1-119 that is not a H atom is being restrained
+
iwrap = 1: all molecules are in the "primary box" when rst and traj files are written. Molecule may move to a neighboring box during simulation especially those molecules are close to the boundary, which is very annoying when we visualize the system.
  
−
'''restraint_wt''': restraint weight indicating how strong the restraint on the atoms is
+
nscm = 100: net translational and rotational velocity is cancelled every 100 steps. During a long simulation, the whole system may start moving and rotating which give an unwanted kinetic energy.
  
−
=== Equilibration ===
+
=== Production ===
−
To further prepare our complex for the molecular dynamic simulations, the subsequent step of energy minimization is equilibrate the system at some certain temperatures. We repeat the process of minimization and equilibration for twice in our case, of course with varied parameters and restraints put on our system.
+
Before the production runs, first we need to generate the input files 10md.in (as shown below).
  
−
Right after the first 1000 steps of minimization, we perform a 50000 steps('''nstlim = 50000''') X the step length of 1 fs('''dt = 0.001''') (that is 50 ps in total) equilibration at the temperature 298.15K('''temp0 = 298.15, tempi = 298.15'''), contining putting on the weight (in kcal/mol) of 5.0 for the positional restraints on all non hydrogen atoms('''restraintmask = ':1-210 & !@H=', restraint_wt = 5.0''').
+
  10md.in: production
 +
    &cntrl
 +
      imin = 0, ntx = 5, irest = 1, nstlim = 500000,
 +
      temp0 = 298.15, tempi = 298.15, ig = 71287,
 +
      ntc = 2, ntf = 1, ntt = 1, dt = 0.002,
 +
      ntb = 2, ntp = 1, tautp = 1.0, taup = 1.0,
 +
      ntwx = 500, ntwe = 0, ntwr = 500, ntpr = 500,
 +
      cut = 8.0, iwrap = 1,
 +
      ntr = 1, nscm = 100,
 +
      restraintmask = ':1-198@CA,C,N', restraint_wt = 0.1,
 +
 
 +
imin = 0
  
−
02md.in: equilibration (50ps)
+
ntx = 5: read coordinates and velocities
−
&cntrl
 
−
  imin = 0, ntx = 1, irest = 0, nstlim = 50000,
 
−
  temp0 = 298.15, tempi = 298.15,
 
−
  ntc = 2, ntf = 1, ntt = 1, dt = 0.001,
 
−
  ntb = 2, ntp = 1, tautp = 1.0, taup = 1.0,
 
−
  ntwx = 1000, ntwe = 0, ntwr = 1000, ntpr = 1000,
 
−
  ntr = 1, nscm = 100, iwrap = 1, cut = 8.0,
 
−
  restraintmask = ':1-210 & !@H=', restraint_wt = 5.0,
 
−
/
 
  
−
*Notice that some flags used in minimization have been changed:
+
irest = 1: now we are restarting. restarting requires reading in velocities.
  
−
'''imin = 0''' indicates that we are doing an MD simulation instead of doing energy minimization.
+
nstlim = 500000: 500000 steps.
  
−
'''ntc = 2 ''' indicates that all bonds involving H bonds are constrained by the SHAKE algorithm to eliminate high frequency oscillations in the system.
+
ig = 71287: random number seeds.
  
−
'''ntb = 2 '''and '''ntp = 1 '''indicate that constant pressure instead of constant volume is applied.
+
dt = 0.002: stepsize = 2fs.
  
−
*Some flags (including those explained above) are newly added;
+
restraintmask = ':1-198@CA,C,N': ligand (resid 199) is free to move around.
  
−
'''ntx = 1''' means coordinates, but no velocities, will be read.
+
=== Running jobs on the queue ===
 +
To perform minimization, equilibration and production, create the following script and submit through qsub:
  
−
'''irest = 0''' irest is a flag to restart a simulation, we always start from running a new simulation until the production part, so it is 0 from 02md.in, 06md.in to 09md.in.  
+
  #!/bin/tcsh
 +
  #PBS -l nodes=2:ppn=2
 +
  #PBS -l walltime=10:00:00
 +
  #PBS -N 1HVR.vs
 +
  #PBS -o 1HVR.output
 +
  #PBS -j oe
 +
  #PBS -V
 +
 
 +
  cd /nfs/user03/<<your username>>/CHE536/AMBER-Tutorial/003.PMEMD
 +
 
 +
  mpirun -n 4 pmemd.MPI -O -i 01mi.in -o 01mi.out -p ../002.ANTE.TLEAP/1HVR.com.wat.leap.prm7 \
 +
  -c ../002.ANTE.TLEAP/1HVR.com.wat.leap.rst7 -ref ../002.ANTE.TLEAP/1HVR.com.wat.leap.rst7 \
 +
  -x 01mi.trj -inf 01mi.info -r 01mi.rst7
 +
  mpirun -n 4 pmemd.MPI -O -i 02md.in -o 02md.out -p ../002.ANTE.TLEAP/1HVR.com.wat.leap.prm7 \
 +
  -c 01mi.rst7 -ref 01mi.rst7 -x 02md.trj -inf 02md.info -r 02md.rst7
 +
  mpirun -n 4 pmemd.MPI -O -i 03mi.in -o 03mi.out -p ../002.ANTE.TLEAP/1HVR.com.wat.leap.prm7 \
 +
  -c 02md.rst7 -ref 02md.rst7 -x 03mi.trj -inf 03mi.info -r 03mi.rst7
 +
  mpirun -n 4 pmemd.MPI -O -i 04mi.in -o 04mi.out -p ../002.ANTE.TLEAP/1HVR.com.wat.leap.prm7 \
 +
  -c 03mi.rst7 -ref 03mi.rst7 -x 04mi.trj -inf 04mi.info -r 04mi.rst7
 +
  mpirun -n 4 pmemd.MPI -O -i 05mi.in -o 05mi.out -p ../002.ANTE.TLEAP/1HVR.com.wat.leap.prm7 \
 +
  -c 04mi.rst7 -ref 04mi.rst7 -x 05mi.trj -inf 05mi.info -r 05mi.rst7
 +
  mpirun -n 4 pmemd.MPI -O -i 06md.in -o 06md.out -p ../002.ANTE.TLEAP/1HVR.com.wat.leap.prm7 \
 +
  -c 05mi.rst7 -ref 05mi.rst7 -x 06md.trj -inf 06md.info -r 06md.rst7
 +
  mpirun -n 4 pmemd.MPI -O -i 07md.in -o 07md.out -p ../002.ANTE.TLEAP/1HVR.com.wat.leap.prm7 \
 +
  -c 06md.rst7 -ref 05mi.rst7 -x 07md.trj -inf 07md.info -r 07md.rst7
 +
  mpirun -n 4 pmemd.MPI -O -i 08md.in -o 08md.out -p ../002.ANTE.TLEAP/1HVR.com.wat.leap.prm7 \
 +
  -c 07md.rst7 -ref 05mi.rst7 -x 08md.trj -inf 08md.info -r 08md.rst7
 +
  mpirun -n 4 pmemd.MPI -O -i 09md.in -o 09md.out -p ../002.ANTE.TLEAP/1HVR.com.wat.leap.prm7 \
 +
  -c 08md.rst7 -ref 05mi.rst7 -x 09md.trj -inf 09md.info -r 09md.rst7
 +
  mpirun -n 4 pmemd.MPI -O -i 10md.in -o 10md.out -p ../002.ANTE.TLEAP/1HVR.com.wat.leap.prm7 \
 +
  -c 09md.rst7 -ref 05mi.rst7 -x 10md.trj -inf 10md.info -r 10md.rst7
 +
  mpirun -n 4 pmemd.MPI -O -i 11md.in -o 11md.out -p ../002.ANTE.TLEAP/1HVR.com.wat.leap.prm7 \
 +
  -c 10md.rst7 -ref 05mi.rst7 -x 11md.trj -inf 11md.info -r 11md.rst7
  
−
'''ntt = 1''' switch the temperature scaling, 1 means constant temperature, using the weak-coupling algorithm.
+
You can also separate each run into several scripts.
  
−
'''tautp = 1.0''' and '''taup = 1.0''' defines the time constant (in ps) and the pressure relaxation time (in ps), respectively.
+
mpirun -n 4: run in parallel with 4 processors.
  
−
'''iwrap = 1''' indicate that the coordinates written to the .rst7 and .trj files will be wrapped into a primary box. This reduces the amount of coordinate output, thus preventing overflowing and making visualization more convenient.
+
pmemd.MPI: parallel version of pmemd.
  
−
'''cut = 8.0''' set the the non-bonded cutoff distance at 8.0 Angstroms.
+
-i: input file.
  
−
Continuing 4 simulations are carried out after the fourth minimization, there are only little differents in these input files:
+
-o: output file.
  
−
06md.in: equilibration (50ps)
+
-p: topology and parameter of the system.
−
  &cntrl
 
−
    imin = 0, ntx = 1, irest = 0, nstlim = 50000,
 
−
    temp0 = 298.15, tempi = 298.15, '''ig = 71287''',
 
−
    ntc = 2, ntf = 1, ntt = 1, dt = 0.001,
 
−
    ntb = 2, ntp = 1, tautp = 1.0, taup = 1.0,
 
−
    ntwx = 1000, ntwe = 0, ntwr = 1000, ntpr = 1000,
 
−
    ntr = 1, nscm = 100, iwrap = 1, cut = 8.0,
 
−
    restraintmask = ':1-210 & !@H=', '''restraint_wt = 1.0''',
 
−
/
 
  
−
*As we can see in 06md.in file, the only changes are setting '''ig''' as 71287 and reducing '''restraint_wt''' value from 5.0 to 1.0.  
+
-c: read coordinates and velocities from the restart file generated by last run.
  
−
'''ig''' is the seed for the pseudo-random number generator which the MD starting velocity depend on.
+
-x: trajectories.
  
−
To further reduce the restraint put on the whole system, we vary '''restraint_wt''' value to be '''0.05, 0.1, 0.1''', respectively in 07md, 08md and 09md input files.
+
-inf: info file (with most recent energy).
  
−
*Note that starting from 08md to 09md simulations, the '''restraintmask = ':1-210 & !@H='''' has been changed to '''restraintmask = ':1-209 & !@H='''' which means not at all restrains are put on the ligand atoms.
+
-r: restart file.
  
−
=== Production ===
+
==IV. Simulation Analysis==
−
Before the production runs (10 and 11), first we need to generate the input files 10md.in and 11md.in (as shown below).
+
Arkin, Jess and Mosavverul
  
−
10md.in: production (500000 = 1ns)
+
===Checking Output Files===
−
&cntrl
 
−
  imin = 0, ntx = 5, irest = 1, nstlim = 500000,
 
−
  temp0 = 298.15, tempi = 298.15, ig = 71287,
 
−
  ntc = 2, ntf = 1, ntt = 1, dt = 0.002,
 
−
  ntb = 2, ntp = 1, tautp = 1.0, taup = 1.0,
 
−
  ntwx = 500, ntwe = 0, ntwr = 500, ntpr = 500,
 
−
  ntr = 1, nscm = 100, iwrap = 1, cut = 8.0,
 
−
  restraintmask = ':1-209@CA,C,N', restraint_wt = 0.1,
 
−
/
 
  
−
11md.in: production (500000 = 1ns)
+
===Ptraj===
−
&cntrl
 
−
  imin = 0, ntx = 5, irest = 1, nstlim = 500000,
 
−
  temp0 = 298.15, tempi = 298.15, ig = 71287,
 
−
  ntc = 2, ntf = 1, ntt = 1, dt = 0.002,
 
−
  ntb = 2, ntp = 1, tautp = 1.0, taup = 1.0,
 
−
  ntwx = 500, ntwe = 0, ntwr = 500, ntpr = 500,
 
−
  ntr = 1, nscm = 100, iwrap = 1, cut = 8.0,
 
−
  restraintmask = ':1-209@CA,C,N', restraint_wt = 0.1,
 
−
/
 
  
−
Explanation for some of the parameter values:
+
"Ptraj" stands for Process Trajectory. The frames in a trajectory can be thought of as snapshots in the MD simulation process. Each snapshot is taken at a frequency designated in the input.  
−
ntx = 5, irest = 1 #using velocity information from the previous restart file#
 
−
nstlim = 500000    #running for 500,000 steps#
 
−
temp0 = 298.15    #the target temperature is 298.15 K (the default is 300 K)#
 
−
tempi = 298.15    #the initial temperature is 298.15 K#
 
−
dt = 0.002        #the time step dt is 0.002 ps (2 fs)#
 
−
ntwr = 500        #to print energy output every 500 steps# 
 
−
ntb = 2            #the constant pressure is switched on#
 
  
−
To run the production, we need to create a c-shell file named "run.production.csh" of which the contents are shown below:
+
For this step we created another work directory ('''004.PTRAJ''', for example).  
−
#! /bin/tcsh
 
−
#-l nodes=4:ppn=2
 
−
#PBS -l walltime=720:00:00
 
−
#PBS -o zzz.qsub.prod.out
 
−
#PBS -e zzz.qsub.prod.err
 
−
#PBS -V
 
−
 
−
set workdir = "~<username>/AMS536/amber-tutorial/03.sander"
 
−
cd ${workdir}
 
−
cat $PBS_NODEFILE | sort | uniq
 
−
 
−
mpirun -n 8 /nfs/user03/wjallen/local/amber12/bin/sander.MPI -O -i 10md.in -o 10md.out -p ../02.ante.tleap/1LOQ.com.wat.leap.prm7 \
 
−
-c 09md.rst7 -ref 09md.rst7 -x 10md.trj -inf 10md.info -r 10md.rst7
 
−
mpirun -n 8 /nfs/user03/wjallen/local/amber12/bin/sander.MPI -O -i 11md.in -o 11md.out -p ../02.ante.tleap/1LOQ.com.wat.leap.prm7 \
 
−
-c 10md.rst7 -ref 10md.rst7 -x 11md.trj -inf 11md.info -r 11md.rst7
 
−
exit
 
  
−
The last step of a production run is to submit the c-shell file to the queue:
+
Create the file '''ptraj.strip.wat.in''' containing the following input:  
−
qsub zzz.production.csh
 
−
It will take almost 24 hour to finish the production run. To check the status of your job, type:
 
−
qstat -u <username>
 
−
After the production run, 10 new files should be created:
 
−
10md.info
 
−
10md.out
 
−
10md.rst7
 
−
10md.trj
 
−
11md.info
 
−
11md.out
 
−
11md.rst7
 
−
11md.trj
 
−
zzz.qsub.prod.err (should be empty if nothing is wrong)
 
−
zzz.qsub.prod.out
 
  
−
=== Running jobs on the queue ===
+
    trajin /nfs/user03/yuchzhou/AMS536/zzz.test/003.PMEMD/10md.trj 1 1000 1
−
To run equilibration and production,
+
    trajin /nfs/user03/yuchzhou/AMS536/zzz.test/003.PMEMD/11md.trj 1 1000 1
−
Write the tcsh file ''equil.produc.qsub.csh'':
+
    trajout 1HVR.trj.strip nobox
 +
    strip :WAT
  
−
  #! /bin/tcsh
+
The numbers 1 1000 1 give the input information about which frames are being used for Ptraj. The first two numbers 1 and 1000 specify the starting and ending snapshots from your trajectory file. The last number 1 specifies the frequency of the snapshot saved, in this case, every frame of the trajectory file.
−
  #PBS -l nodes=4:ppn=2
 
−
  #PBS -l walltime=720:00:00
 
−
  #PBS -o zzz.qsub.out
 
−
  #PBS -e zzz.qsub.err
 
−
  #PBS -V
 
−
 
 
−
  set workdir = "~lingling/AMS536/AMBER_Tutorial/003.SANDER"
 
−
  cd ${workdir}
 
−
  cat $PBS_NODEFILE | sort | uniq
 
−
 
 
−
  mpirun -n 8 /nfs/user03/wjallen/local/amber12/bin/sander.MPI -O -i 01mi.in -o 01mi.out -p ../002.ANTE.TLEAP/1LOQ.com.wat.leap.prm7 \
 
−
  -c ../002.ANTE.TLEAP/1LOQ.com.wat.leap.rst7 -ref ../002.ANTE.TLEAP/1LOQ.com.wat.leap.rst7 \
 
−
  -x 01mi.trj -inf 01mi.info -r 01mi.rst7
 
−
  mpirun -n 8 /nfs/user03/wjallen/local/amber12/bin/sander.MPI -O -i 02md.in -o 02md.out -p ../002.ANTE.TLEAP/1LOQ.com.wat.leap.prm7 \
 
−
  -c 01mi.rst7 -ref 01mi.rst7 -x 02md.trj -inf 02md.info -r 02md.rst7
 
−
  mpirun -n 8 /nfs/user03/wjallen/local/amber12/bin/sander.MPI -O -i 03mi.in -o 03mi.out -p ../002.ANTE.TLEAP/1LOQ.com.wat.leap.prm7 \
 
−
  -c 02md.rst7 -ref 02md.rst7 -x 03mi.trj -inf 03mi.info -r 03mi.rst7
 
−
  mpirun -n 8 /nfs/user03/wjallen/local/amber12/bin/sander.MPI -O -i 04mi.in -o 04mi.out -p ../002.ANTE.TLEAP/1LOQ.com.wat.leap.prm7 \
 
−
  -c 03mi.rst7 -ref 03mi.rst7 -x 04mi.trj -inf 04mi.info -r 04mi.rst7
 
−
  mpirun -n 8 /nfs/user03/wjallen/local/amber12/bin/sander.MPI -O -i 05mi.in -o 05mi.out -p ../002.ANTE.TLEAP/1LOQ.com.wat.leap.prm7 \
 
−
  -c 04mi.rst7 -ref 04mi.rst7 -x 05mi.trj -inf 05mi.info -r 05mi.rst7
 
−
  mpirun -n 8 /nfs/user03/wjallen/local/amber12/bin/sander.MPI -O -i 06md.in -o 06md.out -p ../002.ANTE.TLEAP/1LOQ.com.wat.leap.prm7 \
 
−
  -c 05mi.rst7 -ref 05mi.rst7 -x 06md.trj -inf 06md.info -r 06md.rst7
 
−
  mpirun -n 8 /nfs/user03/wjallen/local/amber12/bin/sander.MPI -O -i 07md.in -o 07md.out -p ../002.ANTE.TLEAP/1LOQ.com.wat.leap.prm7 \
 
−
  -c 06md.rst7 -ref 06md.rst7 -x 07md.trj -inf 07md.info -r 07md.rst7
 
−
  mpirun -n 8 /nfs/user03/wjallen/local/amber12/bin/sander.MPI -O -i 08md.in -o 08md.out -p ../002.ANTE.TLEAP/1LOQ.com.wat.leap.prm7 \
 
−
  -c 07md.rst7 -ref 07md.rst7 -x 08md.trj -inf 08md.info -r 08md.rst7
 
−
  mpirun -n 8 /nfs/user03/wjallen/local/amber12/bin/sander.MPI -O -i 09md.in -o 09md.out -p ../002.ANTE.TLEAP/1LOQ.com.wat.leap.prm7 \
 
−
  -c 08md.rst7 -ref 08md.rst7 -x 09md.trj -inf 09md.info -r 09md.rst7
 
−
  mpirun -n 8 /nfs/user03/wjallen/local/amber12/bin/sander.MPI -O -i 10md.in -o 10md.out -p ../002.ANTE.TLEAP/1LOQ.com.wat.leap.prm7 \
 
−
  -c 09md.rst7 -ref 09md.rst7 -x 10md.trj -inf 10md.info -r 10md.rst7
 
−
  mpirun -n 8 /nfs/user03/wjallen/local/amber12/bin/sander.MPI -O -i 11md.in -o 11md.out -p ../002.ANTE.TLEAP/1LOQ.com.wat.leap.prm7 \
 
−
  -c 10md.rst7 -ref 09md.rst7 -x 11md.trj -inf 11md.info -r 11md.rst7
 
−
  exit
 
−
 
 
−
Note in the above script that for each run, the .rst7, .trj and .info files generated from the previous run provides the initial state to start from. The .prm7 file generated for the hydrated complex by TLEAP provides the force parameters.
 
  
−
Submit the file to the queue and monitor progress.
+
The "strip :WAT" will strip off the waters, and the output will be a new file called 1HVR.trj.strip.  And the last line of the input file will remove the water molecules.
  
−
==IV. Simulation Analysis==
+
To execute the first Ptraj analysis, run the following command specifying the input file '''ptraj.strip.wat.in''' just created and the designated output file '''ptraj.strip.wat.log''':
  
−
===Ptraj===
+
    ptraj  ../002.ANTE.TLEAP/1HVR.com.wat.leap.prm7  ptraj.strip.wat.in  > ptraj.strip.wat.log
  
−
You should create another work directory for this step (004.PTRAJ, for example), if you don't have one.
+
The output file '''1HVR.trj.strip''' will contain the 2000 snapshots of the trajectory.  
−
There will be five runs of analysis, each of which will require different input files.
 
  
−
'''1.''' At first we want to concatenate the two 1ns trajectories together, stripping off the waters, and creating a .strip-file as an output. Below is the input file which will allow us to do so.
+
The second input file for ptraj analysis was called '''ptraj.rec.rmsd.in''' and contained the following input:
  
−
'''ptraj.1.in'''
+
    trajin 1HVR.trj.strip 1 2000 1
−
trajin ../003.SANDER/10md.trj 1 1000 1
+
    trajout 1HVR.com.trj.stripfit
−
trajin ../003.SANDER/11md.trj 1 1000 1
+
    reference ../002.ANTE.TLEAP/1HVR.com.gas.leap.rst7
−
trajout 1LOQ.trj.strip nobox
+
    rms reference out 1HVR.rmsd.CA.dat :1-198@CA
−
strip :WAT
 
  
−
The two sets of numbers ''1 1000 1'' give the input information about which frames are used for the Ptraj. The first two numbers ''1'' and ''1000'' specify the starting and ending snapshots from the trajectory file. The ending number of the snapshot doesn't need to be accurate because if you actually don't have enough snapshots in your trajectory file, Ptraj will read up to the last one you have. The last number ''1'' specifies the frequency of the snapshot saved, in this case, we are saving every frame of the trajectory file. And the last line of the input file will take away all the water molecules.  
+
After the first step the 1HVR.com.trj.strip file is created which contains 2000 snapshots of the trajectory. This file will be compared to the reference file 1HVR.com.gas.leap.rst7 specified in the third line. All the water molecules have been removed in the first step, and it is ready for comparison to the gas phase complex trajectory. The last line says calculates the rmsd for alpha carbon number 1 to 198.
  
−
As the input file is prepared we can launch the first analysis as follows:
+
To execute the the analysis, run the following line.
−
ptraj  ../002.ANTE.TLEAP/1LOQ.com.wat.leap.parm  ptraj.1.in  > ptraj.1.log
 
  
−
As the output you will obtain '''1LOQ.trj.strip''' file which contains ''2000'' snapshots of the trajectory.
+
    ptraj  ../002.ANTE.TLEAP/1HVR.com.gas.leap.prm7  ptraj.rec.rmsd.in  > ptraj.rec.rmsd.log
  
−
'''2.''' Later on we want to compare the output file just obtained to our reference file '''1LOQ.com.gas.leap.rst7''', using the following input file.  
+
This will give a txt file with two columns in it. The first column contains the number of frames being compared and the second column is the rmsd.  
  
−
'''ptraj.2.in'''
+
Then a similar rmsd file for the ligand is generated using the following input file '''ptraj.lig.rmsd.in''':
−
trajin 1LOQ.trj.strip  1 2000 1
 
−
trajout 1LOQ.com.trj.stripfit
 
−
reference ../002.ANTE.TLEAP/1LOQ.com.gas.leap.rst7
 
−
rms reference out 1LOQ.rmsd.CA.txt :1-209@CA
 
  
−
Since we have just concatenated the two trajectories, we will have ''2000'' snapshots in '''1LOQ.trj.strip'''. The third line in the input specifies the reference file, we have taken away all the water molecules during the first step, hence we are working here with the gas phase complex. The last line says we are calculating the rmsd for alpha carbon number ''1'' to ''209''.  
+
    trajin 1HVR.com.trj.stripfit 1 2000 1
 +
    reference ../002.ANTE.TLEAP/1HVR.com.gas.leap.rst7
 +
    rms reference out 1HVR.lig.rmsd.dat :199@C*,N*,O*,S* nofit
  
−
And as we have this file filled out, we can run this step:
+
Similarly, this time the trajectory file will be compared to the reference file 1HVR.com.gas.leap.rst7, only this time the ligand is being compared instead of the receptor. The last line states that the rmsd values for carbon, nitrogen, oxygen and sulfur in residue 199 are being calculated. The <nofit> in the last line specifies that the rmsd of the ligand to the crystal structure is being calculated without any fitting.
−
ptraj  ../002.ANTE.TLEAP/1LOQ.com.gas.leap.parm ptraj.2.in > ptraj.2.log
 
  
−
The output file '''1LOQ.rmsd.CA.txt''' wil contain two columns, the first one is the number of frame and the second one stands for rmsd value.  
+
To execute the the analysis, run the following line.
  
−
'''3.''' Afterwards we will generate a similar file for ligand, using the following input file.   
+
    ptraj  ../002.ANTE.TLEAP/1HVR.com.gas.leap.prm7  ptraj.lig.rmsd.in > ptraj.lig.rmsd.log
  
−
'''ptraj.3.in'''
 
−
trajin 1LOQ.com.trj.stripfit 1 2000 1
 
−
reference ../002.ANTE.TLEAP/1LOQ.com.gas.leap.rst7
 
−
rms reference out 1LOQ.lig.rmsd.txt :210@C*,N*,O*,S* nofit
 
  
−
And then run this step:
+
====RMSD Plots====
−
ptraj  ../002.ANTE.TLEAP/1LOQ.com.gas.leap.parm  ptraj.3.in  > ptraj.3.log
 
  
−
By doing this we will compare the trajectory file to '''1LOQ.com.gas.leap.rst7''' as well, but working with the ligand instead of the receptor (we specified that by pointing that we want to calculate the rmsd for carbon, nitrogen, oxygen and sulfur in residue ''210''.
 
  
−
The last two steps are to obtain energetic information about the system. To do this we take a trajectory file of the gas phase complex '''1LOQ.com.trj.stripfit''', and want to create two more trajectory files containing the information on only receptor and only ligand correspondingly.
+
====Measuring h-bond distances====
  
−
'''4.''' At this step we consider receptor only. The input file is provided below:
+
To measure H-bond distances, first create a script to identify the bonds of interest:
  
−
'''ptraj.4.in'''
+
  trajin 1HVR.com.trj.stripfit                                       #The name of the complex file
−
  trajin 1LOQ.com.trj.stripfit 1 2000 1
+
  distance 124OD1_199O4 :124@OD1 :199@O4 out 124OD1_199O4.dat        #Bonds of interest (Multiple bonds can be listed here)
−
  trajout 1LOQ.rec.trj.stripfit
+
                                                                          -Residue@atom number : Residue@atom number
−
strip :210
 
  
−
And then run this step:
 
−
ptraj  ../002.ANTE.TLEAP/1LOQ.com.gas.leap.parm  ptraj.4.in  > ptraj.4.log
 
  
−
'''5.''' And the same procedure for the ligand, with the following input file:
+
It is important to know which atoms to put in the script.
  
−
'''ptraj.5.in'''
+
Next, enter the following command in your current directory to get an output(.dat) file.
−
trajin 1LOQ.com.trj.stripfit 1 2000 1
 
−
trajout 1LOQ.lig.trj.stripfit
 
−
strip :1-209
 
  
−
And then run this step:
 
−
ptraj  ../002.ANTE.TLEAP/1LOQ.com.gas.leap.parm  ptraj.5.in  > ptraj.5.log
 
  
−
'''6.''' ''(optional)'' Visualization
+
ptraj ../Parameter File Input File > Name.log
  
−
As we've gone through all these steps, the analysis is done. If you want to visualize the trajectories, you first need to copy the trajectory files to Herbie like this, for example (being a level above '''004.PTRAJ''' directory):
 
−
scp 004.PTRAJ/ your_username@herbie.mathlab.sunysb.edu:~AMS536/AMBER_tutorial/004.PTRAJ
 
  
−
Now, launch VMD, then open one of the .prm7 files in 002.ANTE.TLEAP. If you want to visualize the whole complex in gas state, you can open 1LOQ.com.gas.leap.prm7 with AMBER7 Parm from 002.ANTE.TLEAP and then 1LOQ.com.trj.stripfit with AMBER coordinates from 004.PTRAJ. With these files, you can look at the real-time movement of the complex in the gas state. You can repeat this procedure to observe the real-time movement of the complex in the water state. Just open 1LOQ.com.wat.leap.prm7 instead of 1LOQ.com.gas.leap.prm7.  
+
The .dat file you end up with can be opened using the program "xmgrace".
  
−
====RMSD Plots====
+
To do this, enter the command:
  
−
Once having finished running Ptraj, you should find two RMSD files in your directory: '''1LOQ.lig.rmsd.txt''' and '''1LOQ.rmsd.CA.txt'''
+
xmgrace Name.dat
−
The data in these files can be plotted to visualize the structural drift (ie change in RMSD) over the course of your simulation. Large fluctuations in RMSD may be an indication that the simulation was unstable.
 
  
−
Here is a plot for ligand RMSD:
+
Xmgrace allows you to open the .dat file and customize labels and axes.
  
−
[[Image:AMBER_RMSD_plot.jpg|500px|center]]
 
  
−
====Measuring h-bond distances====
+
[[File:1HVR.rmsd.CA.jpg]]
  
 
===MM-GBSA Energy Calculation===
 
===MM-GBSA Energy Calculation===
Line 536: Line 477:
  
 
Then create a csh script, run.sander.rescore.csh, that contains the following lines of command:  
 
Then create a csh script, run.sander.rescore.csh, that contains the following lines of command:  
−
   #! /bin/csh
+
   #! /bin/tcsh
 
   #PBS -l nodes=1:ppn=2
 
   #PBS -l nodes=1:ppn=2
 
   #PBS -l walltime=24:00:00
 
   #PBS -l walltime=24:00:00
Line 542: Line 483:
 
   #PBS -e zzz.qsub.err
 
   #PBS -e zzz.qsub.err
 
   #PBS -V
 
   #PBS -V
−
   set workdir = "/nfs/user03/yechen1/AMBER-Tutorial/005.MMGBSA"
+
 
−
   set sander  = "/nfs/user03/wjallen/local/amber12/bin/sander"
+
   set workdir = nfs/user03/username/AMS536/AMBER_Tutorial/005.MMGBSA
−
   cd ${workdir}
+
   cd /nfs/user03/username/AMS536/AMBER_Tutorial/005.MMGBSA
−
   $sander -O -i gb.rescore.in -o gb.rescore.out.com -p ../002.ANTE.TLEAP /1LOQ.com.gas.leap.prm7 \
+
    
−
   -c ../002.ANTE.TLEAP/1LOQ.com.gas.leap.rst7 -y ../004.PTRAJ/1LOQ.com.trj.stripfit \
+
   mpirun -n 2 sander.MPI -O -i gb.rescore.in \
−
   -r restrt.com -ref ../002.ANTE.TLEAP/1LOQ.com.gas.leap.rst7 -x mdcrd.com -inf mdinfo.com
+
  -o gb.rescore.out.com \
−
   $sander -O -i gb.rescore.in -o gb.rescore.out.lig -p ../002.ANTE.TLEAP/1LOQ.lig.gas.leap.prm7 \
+
  -p ../002.TLEAP/1HVR.com.gas.leap.prm7 \
−
   -c ../002.ANTE.TLEAP/1LOQ.lig.gas.leap.rst7 -y ../004.PTRAJ/1LOQ.lig.trj.stripfit \
+
   -c ../002.TLEAP/1HVR.com.gas.leap.rst7 \
−
   -r restrt.lig -ref ../002.ANTE.TLEAP/1LOQ.lig.gas.leap.rst7 -x mdcrd.lig -inf mdinfo.lig
+
  -y ../004.PTRAJ/1HVR.com.trj.stripfit \
−
   $sander -O -i gb.rescore.in -o gb.rescore.out.rec -p ../002.ANTE.TLEAP/1LOQ.rec.gas.leap.prm7 \
+
   -r restrt.com \
−
   -c ../002.ANTE.TLEAP/1LOQ.rec.gas.leap.rst7 -y ../004.PTRAJ/1LOQ.rec.trj.stripfit \
+
  -ref ../002.TLEAP/1HVR.com.gas.leap.rst7 \
−
   -r restrt.rec -ref ../002.ANTE.TLEAP/1LOQ.rec.gas.leap.rst7 -x mdcrd.rec -inf mdinfo.rec
+
  -x mdcrd.com \
 +
  -inf mdinfo.com
 +
    
 +
  mpirun -n 2 sander.MPI -O -i gb.rescore.in \
 +
  -o gb.rescore.out.lig \
 +
  -p ../002.TLEAP/1HVR.lig.gas.leap.prm7 \
 +
   -c ../002.TLEAP/1HVR.lig.gas.leap.rst7 \
 +
  -y ../004.PTRAJ/1HVR.lig.trj.stripfit \
 +
   -r restrt.lig \
 +
  -ref ../002.TLEAP/1HVR.lig.gas.leap.rst7 \
 +
  -x mdcrd.lig \
 +
  -inf mdinfo.lig
 +
    
 +
  mpirun -n 2 sander.MPI -O -i gb.rescore.in \
 +
  -o gb.rescore.out.test.rec \
 +
  -p ../002.TLEAP/1HVR.rec.gas.leap.prm7 \
 +
   -c ../002.TLEAP/1HVR.rec.gas.leap.rst7 \
 +
  -y ../004.PTRAJ/1HVR.rec.trj.stripfit \
 +
   -r restrt.rec \
 +
  -ref ../002.TLEAP/1HVR.rec.gas.leap.rst7 \
 +
  -x mdcrd.rec \
 +
  -inf mdinfo.rec
 +
 
 
   exit
 
   exit
 
~                                                                                                                                                   
 
~                                                                                                                                                   
Line 566: Line 529:
 
When the job is complete, you will obtain the following output files: <font color="#666666">'''gb.rescore.out.com''', '''gb.rescore.out.lig''', and '''gb.rescore.out.rec'''</font>   
 
When the job is complete, you will obtain the following output files: <font color="#666666">'''gb.rescore.out.com''', '''gb.rescore.out.lig''', and '''gb.rescore.out.rec'''</font>   
 
In these files, the single point energy calculation results will be written for each individual frame. It will be found in the results section and the output file will have an infrastrucutre that is similar to the following:
 
In these files, the single point energy calculation results will be written for each individual frame. It will be found in the results section and the output file will have an infrastrucutre that is similar to the following:
−
      FINAL RESULTS
+
                      FINAL RESULTS
−
    NSTEP      ENERGY          RMS            GMAX        NAME    NUMBER
+
 
−
      1      3.6269E+03    1.8737E+01    1.0472E+02    CB       585
+
 
−
   BOND    =      580.2786 ANGLE  =    1563.7704 DIHED      =    2161.5659
+
 
−
   VDWAALS =    -1684.2762 EEL    =  -13809.8494 EGB        =    -2953.6681
+
    NSTEP      ENERGY          RMS            GMAX        NAME    NUMBER
−
   1-4 VDW =      756.7767 1-4 EEL =    7260.2823 RESTRAINT  =        0.0000
+
        1      6.0197E+03    1.9231E+01    1.0478E+02    CA       950
−
   ESURF  =     9752.0291
+
 
−
   minimization completed, ENE= 0.36269092E+04 RMS= 0.187371E+02
+
   BOND    =      613.6427 ANGLE  =    1575.2750 DIHED      =    2052.4780
 +
   VDWAALS =    -1544.4789 EEL    =  -13857.3791 EGB        =    -2026.3402
 +
   1-4 VDW =      697.7785 1-4 EEL =    8382.6220 RESTRAINT  =        0.0000
 +
   ESURF  =   10126.1041
 +
   minimization completed, ENE= 0.60197022E+04 RMS= 0.192312E+02
 
   TRAJENE: Trajectory file ended
 
   TRAJENE: Trajectory file ended
 
   TRAJENE: Trajene complete.
 
   TRAJENE: Trajene complete.
Line 637: Line 604:
  
 
2013 AMBER MMGBSA plot.jpg
 
2013 AMBER MMGBSA plot.jpg
−
[[File:Example.png]]
+
[[File:MMGBSA.png]]
  
 
==V. Frequently Encountered Problems==
 
==V. Frequently Encountered Problems==
 +
===Mike===
 +
It is very important that you know which filetype you need.  For example, for minimizations be sure that you have the correct input .mol2 file.  The wrond mol2 file will result in file output with no information.
  
−
===Artem===
+
===Ivan===
−
 
 
−
 
 
−
===Brian===
 
−
 
 
−
 
 
−
===He===
 
−
 
 
−
 
 
−
===Jiahui===
 
−
 
 
−
 
 
−
===Jiaye===
 
−
 
 
−
 
 
−
===Koushik===
 
−
 
 
−
 
 
−
===Natalie===
 
  
 +
''Double'' and '''''triple''''' check your syntax before you submit a job to seawulf! It is incredibly disappointing to come back the next morning and realize your job failed simply because of a misplaced period or dash.
  
−
===Nikolay===
+
===Lu===
  
 +
===Yan===
 +
Pay attention to the names of your directories and files. Mostly, the reason why your commands fail to run is just a minor mirror in typing.
  
−
===Weiliang===
+
===Yao===
  
 +
===Junjie===
 +
Before you start minimization, equilibration and production, if the queue is not busy, you can type "qsub -I" in seawulf and "cd" to where you want to run sander. Then copy your sander command in your .csh file that you are going to submit and run it by simply paste to the command line and hit "return". By doing this, you are running your job interactively. If you can tell the job is running, your script is likely to be correct. It is easy to debug your .csh script in this way, since error is printed out right after you run commands instead of stored in some .er or .out files that have no context.
  
−
===Ye===
+
===Tianao===
  
 +
===Kai===
  
−
===Yuan===
+
===Arkin===
  
 +
===Jess===
  
−
===Zach===
+
===Mosavverul===

Latest revision as of 20:14, 23 April 2014

For additional Rizzo Lab tutorials see AMBER Tutorials.

In this tutorial, we will learn how to run a molecular dynamics simulation of a protein-ligand complex. We will then post-process that simulation by calculating structural fluctuations (with RMSD) and free energies of binding (MM-GBSA).

I. Introduction

AMBER

Amber - Assisted Model Building with Energy Refinement - is a multi-program suite for macromolecular simulations. Amber12 is the most recent version of the software and it includes new force fields such as ff12SB and Lipid11, expanded options for Poisson-Boltzmann solvation calculations, accelerated molecular dynamics, additional features in sander pmemd code, and expanded methods for free energy calculations. Our lab is set up with Amber12 and the latest update of AmberTools13 which contains the programs such as antechamber and tleap to set up your simulation.

The Amber 12 Manual is available to get started with using Amber12 as well as a troubleshooting reference. (Tip: You can search the document for keywords if you use Adobe Acrobat to view the file which can save time.) Additionally, Amber Tools User's Manual is another reference for the probrgrams available under Amber tools.

Here some of programs availabe in Amber and AmberTools

  1. LEaP: is a preparatory program for constructing new or modified systems in Amber. It contains the functions: prep, link, edit, and parm from earlier version of Amber.
  2. ANTECHAMBER: is the primary program used to prepare input files excluding the original nucleic acid and protein information taken from the PDB.
  3. SANDER: is 'a basic energy minimizer and molecular dynamics program' that is used to minimize, equilibrate and sample molecular conformations. This is the program used to generate trajectory files of the system.
  4. PMEMD: is an improved version of SANDER optimized for periodic, PME simulations, and for GB simulations. It is faster and scales better on parallel machines than SANDER.
  5. PTRAJ: is an analysis program for processing trajectory files. ptraj can be used to rotate and translate the structures, evaluate geometrical features, calculate RMSDs among other analyses.

There is a mailing list available as an additional resource.

HIV Protease

Human Immunodeficiency Virus (HIV) is a retrovirus which causes Acquired Immune Deficiency Syndrome (AIDS). HIV protease is a aspartyl protease which is a dimer, each monomer consisting of 99 amino acids. This protein is essential for the lifecycle of the virus as it cleaves immature polyproteins into mature proteins which can then be used for proper virion assembly. HIV protease has been validated as a drug target for the treatment for HIV as inhibition of this protein can slow HIV progression. However, retroviruses have high mutation rates and drug-resistance is a constant problem which needs to be addressed. There are a number of structure of HIV protease in the protein databank, we will use 1HVR for this tutorial.

Organizing Directories

The following directory structure and naming scheme is a convienient way to organize your files for this tutorial. Create these directories first before doing anything else. 

~username/AMS536/AMBER-Tutorial/001.MOL.PREP/  
                                002.ANTE.TLEAP/ 
                                003.PMEMD/       
                                004.PTRAJ/
                                005.MMGBSA/

II. Structural Preparation

Downloading the PDB file

Go to PDB homepage (http://www.rcsb.org/pdb/home/home.do ) enter the protein ID (1HVR) in the search bar, click Download Files in the top-right of the webpage, then select PDB File (text). In the new window, save the file.

Preparing the ligand and receptor in Chimera

In this section, three new files would be created in the 001.CHIMERA.MOL.PREP/ folder: 

 1HVR.dockprep.mol2      (complete system with hydrogens and charges)
 1HVR.receptor.noH.pdb  (the receptor alone, without hydrogens)
 1HVR.ligand.mol2        (the ligand alone)

To prepare these files, first copy the original PDB file into 001.MOL.PREP folder and open it with VIM ($ vim 1HVR.pdb). Because the residue name of the ligand (U) will give us some problems when assigning charges, change the residue name "U" to "LIG". Here is an example command that will change all instances of " U" into "LIG", while preserving the correct spacing (Note: there are two spaces before U): 

 :%s/  U/LIG/gc

For this command, g is short for global and c is short for check with the user before making the change.

Next, open up the PDB file (1HVR.pdb) in Chimera. To generate 1HVR.dockprep.mol2 files, delete the water molecules; delete the original hydrogen atoms; add the charge and add the hydrogen atoms.

To delete water molecules and other ligands, click Select -> Residue -> HOH, then go to Actions -> Atoms/Bonds -> Delete. Hydrogen atoms can be added manually by choosing Tools -> Structure Editing -> Add H. To add charges to the ligand and receptor, go to Select -> Residue -> LIG, then go to Tools -> Structure Editing -> Add Charge. Choose AMBER ff99SB as the charge model, click Okay, and when prompted chose AM1-BCC charges for the ligand, and make sure the Net Charge is set to 0. (The chemistry of the ligand should be considered when assigning a charge state).

Alternatively, you can do all of the above by clicking Tools -> Structure Editing -> Dock Prep. Note when adding the charge to the ligand, you can choose AMBER ff99SB as the charge model and chose gasteiger as the charge method. In this 1HVR case, we set Net Charge to 0.

Finally, save this file as 1HVR.dockprep.mol2.

Generating the final files

To create the receptor file with no hydrogen atoms: Open 1HVR.dockprep.mol2, click Select -> Chemistry -> Element -> H, then chose Actions -> Atoms/Bonds -> Delete'. Save the file as 1HVR.receptor.noH.pdb.

To create the ligand file: Open 1HVR.dockprep.mol2, click Select -> Residue -> LIG, then click Select -> Invert, then chose Actions -> Atoms/Bonds -> Delete. Save the file as 1HVR.ligand.mol2.

antechamber

An antichamber input file need all the atom names to be unique and only the first three characters will be used. For example, H301 and H302 can not be distinguished with each other. So first we need to manually rename the atoms which have the same first three characters. 

       007.png

The file after renaming: 

       009.png

To run antichamber, we need make sure these three programs are available, (antechamber, parmchk and tleap) and we are using the correct version of amber, we can use the which command, type command: 

  which antechamber
  which parmchk
  which tleap

The result should look like this: 

  /nfs/user03/wjallen/local/amber12/bin/antechamber
  /nfs/user03/wjallen/local/amber12/bin/parmchk
  /nfs/user03/wjallen/local/amber12/bin/tleap

Copy parameters of ions to your working directory from the following resource: 

  scp -r ~yuchzhou/AMS536/AMBER_Tutorial/002.ANTE.TLEAP/rizzo_amber7.ionparms ./

Then we use antechamber to convert our input mol2 file into files ready for LEaP.Type command: 

  antechamber -i ../001.MOL.PREP/1HVR.lig.mol2 -fi mol2 -o 1HVR.lig.ante.pdb -fo pdb

Here, -i means input file name; -fi means input file format; -o means output file name; -fo output file format. You will have an output file:1HVR.lig.ante.pdb Similarly, we can use antechamber to change the format of 1HVR.lig.mol2 file to prep file: 

  antechamber -i ../001.MOL.PREP/1HVR.lig.mol2 -fi mol2 -o 1HVR.lig.ante.prep -fo prepi

You will get a set of output files: 

  ANTECHAMBER_AC.AC  
  ANTECHAMBER_AC.AC0  
  ANTECHAMBER_BOND_TYPE.AC  
  ANTECHAMBER_BOND_TYPE.AC0  
  ANTECHAMBER_PREP.AC  
  ANTECHAMBER_PREP.AC0
  ATOMTYPE.INF 
  NEWPDB.PDB 
  PREP.INF 
  1HVR.lig.ante.prep

parmchk

After running antechamber, we run parmchk to check the parameters. Parmchk is another program in Antechamber. If there are missing parameters or mistake after antechamber is finished, the frcmod template generated by parmchk will help us in generating the needed parameters: 

  parmchk -i 1HVR.lig.ante.prep -f prepi -o 1HVR.lig.ante.frcmod

In the command, -i input file name,-f input file format (prepi, prepc, ac, mol2),-o frcmod file name. When the programming is done, you should have a new file in the directory. 

  1HVR.lig.ante.frcmod

tleap

Then, three input files are needed to run TLEAP, they are “tleap.lig.in”, “tleap.rec.in” and “tleap.com.in”, for the ligand, receptor, and complex respectively. These input files will be used by tleaP to create parameter/topology files and initial coordinate files. Copy parameters to your working directory from the following resource: 

  scp  ~yuchzhou/AMS536/AMBER_Tutorial/002.ANTE.TLEAP/tleap.lig.in
  scp  ~yuchzhou/AMS536/AMBER_Tutorial/002.ANTE.TLEAP/tleap.rec.in
  scp  ~yuchzhou/AMS536/AMBER_Tutorial/002.ANTE.TLEAP/tleap.com.in

here are what the files should look like.

  • tleap.lig.in
set default PBradii mbondi2                                      #set default PBradii                                         
source leaprc.ff99SB                                             #load a force field
loadoff rizzo_amber7.ionparms/ions.lib                           #Loading Rob Rizzo's ion parameters
loadamberparams rizzo_amber7.ionparms/ions.frcmod                #load an AMBER format parameter file
source leaprc.gaff
loadamberparams rizzo_amber7.ionparms/gaff.dat.rizzo
loadamberparams 1HVR.lig.ante.frcmod
loadamberprep 1HVR.lig.ante.prep                                 #load an AMBER PREP file
lig = loadpdb 1HVR.lig.ante.pdb                                  #load the ligand pdb file
saveamberparm lig 1HVR.lig.gas.leap.prm7 1HVR.lig.gas.leap.rst7  #save the ligand gas phase AMBER topology and coordinate files
solvateBox lig TIP3PBOX 10.0                                     #solvate the receptor using TIP3P, solvent box radii 10 angstroms
saveamberparm lig 1HVR.lig.wat.leap.prm7 1HVR.lig.wat.leap.rst7  #save the ligand water phase AMBER topology and coordinate files
charge lig                                                       #this writes out the charge of the ligand
quit

Coments (text after the '#' symbol) are for edification and should not apper in the file.

  • tleap.rec.in
set default PBradii mbondi2
source leaprc.ff99SB
loadoff rizzo_amber7.ionparms/ions.lib
loadamberparams rizzo_amber7.ionparms/ions.frcmod
REC = loadpdb ../001.MOL.PREP/1HVR.rec.noh.pdb
saveamberparm REC 1HVR.rec.gas.leap.prm7 1HVR.rec.gas.leap.rst7
solvateBox REC TIP3PBOX 10.0                                     
saveamberparm REC 1HVR.rec.wat.leap.prm7 1HVR.rec.wat.leap.rst7
charge REC
quit
  • tleap.com.in
set default PBradii mbondi2
source leaprc.ff99SB
loadoff rizzo_amber7.ionparms/ions.lib
loadamberparams rizzo_amber7.ionparms/ions.frcmod
source leaprc.gaff
loadamberparams rizzo_amber7.ionparms/gaff.dat.rizzo
REC = loadpdb ../001.MOL.PREP/1HVR.rec.noh.pdb
loadamberparams 1HVR.lig.ante.frcmod
loadamberprep 1HVR.lig.ante.prep
LIG = loadpdb 1HVR.lig.ante.pdb
COM = combine {REC LIG}
saveamberparm COM 1HVR.com.gas.leap.prm7 1HVR.com.gas.leap.rst7
solvateBox COM TIP3PBOX 10.0
saveamberparm COM 1HVR.com.wat.leap.prm7 1HVR.com.wat.leap.rst7
charge COM
quit

Use leap command: 

  tleap -s -f tleap.lig.in > tleap.lig.out
  tleap -s -f tleap.rec.in > tleap.rec.out
  tleap -s -f tleap.com.in > tleap.com.out

The following output files will be generated: 

  1HVR.com.gas.leap.prm7 
  1HVR.com.wat.leap.prm7	 
  1HVR.lig.ante.frcmod 
  1HVR.lig.ante.prep 
  1HVR.lig.gas.leap.rst7 
  1HVR.lig.wat.leap.rst7 
  1HVR.rec.gas.leap.rst7 
  1HVR.rec.wat.leap.rst7 
  1HVR.com.gas.leap.rst7	 
  1HVR.com.wat.leap.rst7 
  1HVR.lig.ante.pdb 
  1HVR.lig.gas.leap.prm7 
  1HVR.lig.wat.leap.prm7 
  1HVR.rec.gas.leap.prm7 
  1HVR.rec.wat.leap.prm7 
  tleap.lig.out 
  tleap.rec.out 
  tleap.com.out 
  leap.log

Visualization in VMD

Visualization is an important step in AMBER molecular dynamics simulation as it allows for the viewing of molecules and molecule movements within a specified field of view. Several files prepared in the previous two steps will be required for the visualization of the ligand and its movement in the protein binding pocket. The first step that must be completed is the copying of all necessary information from Seawulf to Herbie. 

  scp -r ~/AMS536/AMBER_Tutorial/002.ANTE.TLEAP username@herbie.mathlab.sunysb.edu:~/AMS536/AMBER_Tutorial/

We can now visualize several files in herbie by VMD: the protein-ligand complex in the gas or water phase, and the ligand in the gas or water phase.In the "VMD Main" window, click File > New Molecule. Now in the new window titled "Molecule File Browser" a few files must be selected and loaded.Viewing the protein-ligand complex in the gas phase we select 1LOQ.com.gas.leap.prm7 and the file type AMBER7 Parm. This would allow you to load the parameter information of the milecule, but nothing will shoe up in the VMD visualization window because coordinates has not been loaded yet.Next you need to select and load the file 1LOQ.com.gas.leap.rst7 and the file type AMBER7 Restart. This can be done for the complex or ligand for either the gas or the water phase by selecting and loading the corresponding .parm7 files and .rst7 files.

Next, you can choose to edit the structure shown for better visualization. By clicking Graphics > Representations... a new window "Graphical Representations" will open, in which you can create new representations of the protein-ligand complex. It is possible to choose the entire complex or simply the protein or ligand, to color the structure in various ways, and to choose how it is best represented (e.g. lines, thicker bonds, full surface).

Protein-ligand complex in the gas visualization.png

III. Simulation using sander

Junjie, Tianao and Kai

Minimization

Before running any equilibration and production, we have to minimize the energy of system. This process can remove any unfavorable interaction between atoms due to the low resolution of biomolecule. If an unfavorable interaction is present when doing equilibration and production, atoms may "fly apart" and "crash" into other part of the biomolecule.

An example of a typical input file “01mi.in” for minimization is listed blow: 

  01mi.in: minimization
  &cntrl
  imin = 1, maxcyc = 1000, ntmin = 2,
  ntx = 1, ntc = 1, ntf = 1,
  ntb = 1, ntp = 0,
  ntwx = 1000, ntwe = 0, ntpr = 1000,
  cut = 8.0,
  ntr = 1,
  restraintmask = ':1-199 & !@H=',
  restraint_wt = 5.0,
  /

You may have to minimize the system several times. Each time you release the constraint on heavy atoms slowly.

The meaning of each flag is:

imin = 1: do minimization only.

maxcyc: cycles of minimization will be performed.

ntmin: specify the number of cycles that use a deepest decent method, after that minimization will use the conjugate gradient method.

ntx = 1: read coordinates only from rst file. There is no need to read velocity before minimization.

ntc = 1: Minimization is performed on every bond.

ntf = 1: Typically, ntf = ntc.

ntb: periodic boundary is applied.

ntp = 0: NO constant pressure is applied.

ntwx = 1000: saving snapshot to a trajectory file every ntwx steps.

ntwe = 0: human readable energy file is appended every ntwe steps.

ntpr = 1000: newest calculated energy is stored in .info file every ntpr steps.

cut: cutoff used in energy evaluation.

ntr = 1: positional restraint method is applied.

restraintmask = ':1-199 & !@H : atoms of residue 1-199 except H is restrained.

restraint_wt: restraint strength in unit of kcal/(mol*Angstrom^2).

Equilibration

System need to be heat up in a every careful way. An abrupt heating up may cause unwanted perturbation.

An example of the content of "02md.in" is listed below: 

  02md.in: equilibration
  &cntrl
  imin = 0, ntx = 1, irest = 0, nstlim = 50000,
  temp0 = 298.15, tempi = 298.15, ig = 71287,
  ntc = 2, ntf = 1, ntt = 1, dt = 0.001,
  ntb = 2, ntp = 1, tautp = 1.0, taup = 1.0,
  ntwx = 1000, ntwe = 0, ntwr = 1000, ntpr = 1000,
  cut = 8.0, iwrap = 1,
  ntr = 1, nscm = 100,
  restraintmask = ':1-199 & !@H=', restraint_wt = 5.0,
  /

The parameters are described as below:

imin = 0: now md simulation is performed instead of minimization.

ntc = 2: SHAKE is applied to constraint bonds involve H.

ntb = 2 and ntp = 1: constant pressure is applied.

irest = 0: we are not reading in any velocities so there is no need to use restart.

ntt = 1: weak-coupling is applied to keep the temperature constant during equilibration.

tautp = 1.0 and taup = 1.0: time constant in unit of ps and the pressure relaxation time in unit of ps.

iwrap = 1: all molecules are in the "primary box" when rst and traj files are written. Molecule may move to a neighboring box during simulation especially those molecules are close to the boundary, which is very annoying when we visualize the system.

nscm = 100: net translational and rotational velocity is cancelled every 100 steps. During a long simulation, the whole system may start moving and rotating which give an unwanted kinetic energy.

Production

Before the production runs, first we need to generate the input files 10md.in (as shown below). 

  10md.in: production
   &cntrl
     imin = 0, ntx = 5, irest = 1, nstlim = 500000,
     temp0 = 298.15, tempi = 298.15, ig = 71287,
     ntc = 2, ntf = 1, ntt = 1, dt = 0.002,
     ntb = 2, ntp = 1, tautp = 1.0, taup = 1.0,
     ntwx = 500, ntwe = 0, ntwr = 500, ntpr = 500,
     cut = 8.0, iwrap = 1,
     ntr = 1, nscm = 100,
     restraintmask = ':1-198@CA,C,N', restraint_wt = 0.1,

imin = 0

ntx = 5: read coordinates and velocities

irest = 1: now we are restarting. restarting requires reading in velocities.

nstlim = 500000: 500000 steps.

ig = 71287: random number seeds.

dt = 0.002: stepsize = 2fs.

restraintmask = ':1-198@CA,C,N': ligand (resid 199) is free to move around.

Running jobs on the queue

To perform minimization, equilibration and production, create the following script and submit through qsub: 

  #!/bin/tcsh
  #PBS -l nodes=2:ppn=2
  #PBS -l walltime=10:00:00
  #PBS -N 1HVR.vs
  #PBS -o 1HVR.output
  #PBS -j oe
  #PBS -V
  
  cd /nfs/user03/<<your username>>/CHE536/AMBER-Tutorial/003.PMEMD
  
  mpirun -n 4 pmemd.MPI -O -i 01mi.in -o 01mi.out -p ../002.ANTE.TLEAP/1HVR.com.wat.leap.prm7 \
  -c ../002.ANTE.TLEAP/1HVR.com.wat.leap.rst7 -ref ../002.ANTE.TLEAP/1HVR.com.wat.leap.rst7 \
  -x 01mi.trj -inf 01mi.info -r 01mi.rst7
  mpirun -n 4 pmemd.MPI -O -i 02md.in -o 02md.out -p ../002.ANTE.TLEAP/1HVR.com.wat.leap.prm7 \
  -c 01mi.rst7 -ref 01mi.rst7 -x 02md.trj -inf 02md.info -r 02md.rst7
  mpirun -n 4 pmemd.MPI -O -i 03mi.in -o 03mi.out -p ../002.ANTE.TLEAP/1HVR.com.wat.leap.prm7 \
  -c 02md.rst7 -ref 02md.rst7 -x 03mi.trj -inf 03mi.info -r 03mi.rst7
  mpirun -n 4 pmemd.MPI -O -i 04mi.in -o 04mi.out -p ../002.ANTE.TLEAP/1HVR.com.wat.leap.prm7 \
  -c 03mi.rst7 -ref 03mi.rst7 -x 04mi.trj -inf 04mi.info -r 04mi.rst7
  mpirun -n 4 pmemd.MPI -O -i 05mi.in -o 05mi.out -p ../002.ANTE.TLEAP/1HVR.com.wat.leap.prm7 \
  -c 04mi.rst7 -ref 04mi.rst7 -x 05mi.trj -inf 05mi.info -r 05mi.rst7
  mpirun -n 4 pmemd.MPI -O -i 06md.in -o 06md.out -p ../002.ANTE.TLEAP/1HVR.com.wat.leap.prm7 \
  -c 05mi.rst7 -ref 05mi.rst7 -x 06md.trj -inf 06md.info -r 06md.rst7
  mpirun -n 4 pmemd.MPI -O -i 07md.in -o 07md.out -p ../002.ANTE.TLEAP/1HVR.com.wat.leap.prm7 \
  -c 06md.rst7 -ref 05mi.rst7 -x 07md.trj -inf 07md.info -r 07md.rst7
  mpirun -n 4 pmemd.MPI -O -i 08md.in -o 08md.out -p ../002.ANTE.TLEAP/1HVR.com.wat.leap.prm7 \
  -c 07md.rst7 -ref 05mi.rst7 -x 08md.trj -inf 08md.info -r 08md.rst7
  mpirun -n 4 pmemd.MPI -O -i 09md.in -o 09md.out -p ../002.ANTE.TLEAP/1HVR.com.wat.leap.prm7 \
  -c 08md.rst7 -ref 05mi.rst7 -x 09md.trj -inf 09md.info -r 09md.rst7
  mpirun -n 4 pmemd.MPI -O -i 10md.in -o 10md.out -p ../002.ANTE.TLEAP/1HVR.com.wat.leap.prm7 \
  -c 09md.rst7 -ref 05mi.rst7 -x 10md.trj -inf 10md.info -r 10md.rst7
  mpirun -n 4 pmemd.MPI -O -i 11md.in -o 11md.out -p ../002.ANTE.TLEAP/1HVR.com.wat.leap.prm7 \
  -c 10md.rst7 -ref 05mi.rst7 -x 11md.trj -inf 11md.info -r 11md.rst7

You can also separate each run into several scripts.

mpirun -n 4: run in parallel with 4 processors.

pmemd.MPI: parallel version of pmemd.

-i: input file.

-o: output file.

-p: topology and parameter of the system.

-c: read coordinates and velocities from the restart file generated by last run.

-x: trajectories.

-inf: info file (with most recent energy).

-r: restart file.

IV. Simulation Analysis

Arkin, Jess and Mosavverul

Checking Output Files

Ptraj

"Ptraj" stands for Process Trajectory. The frames in a trajectory can be thought of as snapshots in the MD simulation process. Each snapshot is taken at a frequency designated in the input.

For this step we created another work directory (004.PTRAJ, for example).

Create the file ptraj.strip.wat.in containing the following input: 

    trajin /nfs/user03/yuchzhou/AMS536/zzz.test/003.PMEMD/10md.trj 1 1000 1
    trajin /nfs/user03/yuchzhou/AMS536/zzz.test/003.PMEMD/11md.trj 1 1000 1
    trajout 1HVR.trj.strip nobox
    strip :WAT

The numbers 1 1000 1 give the input information about which frames are being used for Ptraj. The first two numbers 1 and 1000 specify the starting and ending snapshots from your trajectory file. The last number 1 specifies the frequency of the snapshot saved, in this case, every frame of the trajectory file.

The "strip :WAT" will strip off the waters, and the output will be a new file called 1HVR.trj.strip. And the last line of the input file will remove the water molecules.

To execute the first Ptraj analysis, run the following command specifying the input file ptraj.strip.wat.in just created and the designated output file ptraj.strip.wat.log: 

    ptraj  ../002.ANTE.TLEAP/1HVR.com.wat.leap.prm7  ptraj.strip.wat.in  > ptraj.strip.wat.log

The output file 1HVR.trj.strip will contain the 2000 snapshots of the trajectory.

The second input file for ptraj analysis was called ptraj.rec.rmsd.in and contained the following input: 

    trajin 1HVR.trj.strip 1 2000 1
    trajout 1HVR.com.trj.stripfit
    reference ../002.ANTE.TLEAP/1HVR.com.gas.leap.rst7
    rms reference out 1HVR.rmsd.CA.dat :1-198@CA

After the first step the 1HVR.com.trj.strip file is created which contains 2000 snapshots of the trajectory. This file will be compared to the reference file 1HVR.com.gas.leap.rst7 specified in the third line. All the water molecules have been removed in the first step, and it is ready for comparison to the gas phase complex trajectory. The last line says calculates the rmsd for alpha carbon number 1 to 198.

To execute the the analysis, run the following line. 

    ptraj  ../002.ANTE.TLEAP/1HVR.com.gas.leap.prm7  ptraj.rec.rmsd.in  > ptraj.rec.rmsd.log

This will give a txt file with two columns in it. The first column contains the number of frames being compared and the second column is the rmsd.

Then a similar rmsd file for the ligand is generated using the following input file ptraj.lig.rmsd.in: 

   trajin 1HVR.com.trj.stripfit 1 2000 1
   reference ../002.ANTE.TLEAP/1HVR.com.gas.leap.rst7
   rms reference out 1HVR.lig.rmsd.dat :199@C*,N*,O*,S* nofit

Similarly, this time the trajectory file will be compared to the reference file 1HVR.com.gas.leap.rst7, only this time the ligand is being compared instead of the receptor. The last line states that the rmsd values for carbon, nitrogen, oxygen and sulfur in residue 199 are being calculated. The <nofit> in the last line specifies that the rmsd of the ligand to the crystal structure is being calculated without any fitting.

To execute the the analysis, run the following line. 

    ptraj  ../002.ANTE.TLEAP/1HVR.com.gas.leap.prm7  ptraj.lig.rmsd.in  > ptraj.lig.rmsd.log


RMSD Plots

Measuring h-bond distances

To measure H-bond distances, first create a script to identify the bonds of interest: 

trajin 1HVR.com.trj.stripfit                                       #The name of the complex file
distance 124OD1_199O4 :124@OD1 :199@O4 out 124OD1_199O4.dat        #Bonds of interest (Multiple bonds can be listed here)
                                                                          -Residue@atom number : Residue@atom number


It is important to know which atoms to put in the script.

Next, enter the following command in your current directory to get an output(.dat) file. 


ptraj ../Parameter File Input File > Name.log


The .dat file you end up with can be opened using the program "xmgrace".

To do this, enter the command: 

xmgrace Name.dat

Xmgrace allows you to open the .dat file and customize labels and axes.


1HVR.rmsd.CA.jpg

MM-GBSA Energy Calculation

MM/GBSA is the acronym for Molecular Mechanics/Generalized Born Surface Area. This part of AMBER combines molecular mechanics (MM) with both the electrostatic (PB) and nonpolar (SA) contribution to solvation . Topology files are needed for the receptor, ligand, and receptor-ligand complex. The trajectory files generate the coordinates. Therefore, molecular dynamics is used to generate a set of snapshots taken at fixed intervals from the trajectories. These snapshots are processed to remove solvent and generate the free energy of binding.

In the AMBER_Tutorial directory, create a new directory: 

 mkdir 005.MMGBSA

In this new directory, create the file gb.rescore.in containing: 

 Single point GB energy calc
&cntrl
 ntf    = 1,        ntb    = 0,        ntc    = 2,
 idecomp= 0,
 igb    = 5,        saltcon= 0.00,
 gbsa   = 2,        surften= 1.0,
 offset = 0.09,     extdiel= 78.5,
 cut    = 99999.0,  nsnb   = 99999,
 imin   = 5,        maxcyc = 1,        ncyc   = 0,
 /

Then create a csh script, run.sander.rescore.csh, that contains the following lines of command: 

 #! /bin/tcsh
 #PBS -l nodes=1:ppn=2
 #PBS -l walltime=24:00:00
 #PBS -o zzz.qsub.out
 #PBS -e zzz.qsub.err
 #PBS -V
 
 set workdir = nfs/user03/username/AMS536/AMBER_Tutorial/005.MMGBSA
 cd /nfs/user03/username/AMS536/AMBER_Tutorial/005.MMGBSA
 
 mpirun -n 2 sander.MPI -O -i gb.rescore.in \
 -o gb.rescore.out.com \
 -p ../002.TLEAP/1HVR.com.gas.leap.prm7 \
 -c ../002.TLEAP/1HVR.com.gas.leap.rst7 \
 -y ../004.PTRAJ/1HVR.com.trj.stripfit \
 -r restrt.com \
 -ref ../002.TLEAP/1HVR.com.gas.leap.rst7 \
 -x mdcrd.com \
 -inf mdinfo.com
 
 mpirun -n 2 sander.MPI -O -i gb.rescore.in \
 -o gb.rescore.out.lig \
 -p ../002.TLEAP/1HVR.lig.gas.leap.prm7 \
 -c ../002.TLEAP/1HVR.lig.gas.leap.rst7 \
 -y ../004.PTRAJ/1HVR.lig.trj.stripfit \
 -r restrt.lig \
 -ref ../002.TLEAP/1HVR.lig.gas.leap.rst7 \
 -x mdcrd.lig \
 -inf mdinfo.lig
 
 mpirun -n 2 sander.MPI -O -i gb.rescore.in \
 -o gb.rescore.out.test.rec \
 -p ../002.TLEAP/1HVR.rec.gas.leap.prm7 \
 -c ../002.TLEAP/1HVR.rec.gas.leap.rst7 \
 -y ../004.PTRAJ/1HVR.rec.trj.stripfit \
 -r restrt.rec \
 -ref ../002.TLEAP/1HVR.rec.gas.leap.rst7 \
 -x mdcrd.rec \
 -inf mdinfo.rec
 
 exit

~ ~ ~ ~

Then this script should be sent to the queue, i.e., qsub the script using the commands: 

 qsub run.sander.rescore.csh

You can monitor your progress by typing 

 qstat -u username

When the job is complete, you will obtain the following output files: gb.rescore.out.com, gb.rescore.out.lig, and gb.rescore.out.rec In these files, the single point energy calculation results will be written for each individual frame. It will be found in the results section and the output file will have an infrastrucutre that is similar to the following: 

                     FINAL RESULTS
 
 
 
    NSTEP       ENERGY          RMS            GMAX         NAME    NUMBER  
       1       6.0197E+03     1.9231E+01     1.0478E+02     CA        950
 
  BOND    =      613.6427  ANGLE   =     1575.2750  DIHED      =     2052.4780
  VDWAALS =    -1544.4789  EEL     =   -13857.3791  EGB        =    -2026.3402
  1-4 VDW =      697.7785  1-4 EEL =     8382.6220  RESTRAINT  =        0.0000
  ESURF   =    10126.1041
 minimization completed, ENE= 0.60197022E+04 RMS= 0.192312E+02
 TRAJENE: Trajectory file ended
 TRAJENE: Trajene complete.

In the command line, type: 

 grep VDWAALS gb.rescore.out.com > vdw.com.txt.
 grep ESURF gb.rescore.out.com > surf.com.txt.

You can take these text files, import them into Excel, and do the rest of your modifications there.

Equations for analysis

Remember that to obtain the Gvdw term, you need to take the SASA (which is ESURF) and input into equation 1:

ΔGnonpolar = SASA*0.00542 + 0.92

Also, the mmgbsa of a given system can be determined by equation 2:

ΔGmmgbsa = ΔGvdw + ΔGcoul + ΔGpolar + ΔGnonpolar

From the output file:

VDWAALS = ΔGvdw

EELS = ΔGcoul

EGB = ΔGpolar

You can then easily calculate the ΔΔGbind by using equation 3:

ΔΔGbind = ΔGmmgbsa,complex – (ΔGmmgbsa,lig + ΔGmmgbsa,rec) You will want to careful when doing your analysis that the results from frame 1 for the receptor and ligand are subtracted from the results from frame 1 for your complex. By doing this in excel, you should have 2000 frames for each, and the values should cleanly line up. Finally, you will want to plot your ΔΔGbind and examine if you see corresponding changes in the ligand position and the ΔΔGbind. Also, you should determine the mean and standard deviation for your ΔΔGbind.

Plotting Energy

When your rescoring calculation finishes, you have the following three output files: "gb.rescore.out.com", "gb.rescore.out.lig", and "gb.rescore.out.rec".

Use the following script, entitled get.mmgbsa.bash, to extract data and calculate MMGBSA energy for each snap shot. 

 #! /bin/bash
 # by Haoquan
 echo com lig rec > namelist
 LIST=`cat namelist`
 for i in $LIST ; do
 grep VDWAALS gb.rescore.out.$i | awk '{print $3}' > $i.vdw
 grep EGB     gb.rescore.out.$i | awk '{print $9}' > $i.polar
 grep EELS    gb.rescore.out.$i | awk '{print $6}' > $i.coul
 grep ESURF   gb.rescore.out.$i | awk '{print $3 * 0.00542 + 0.92}' > $i.surf
 paste -d " " $i.vdw $i.polar $i.surf $i.coul | awk '{print $1 + $2 + $3 + $4}' > data.$i
 rm $i.*
 done
 paste -d " " data.com data.lig data.rec | awk '{print $1 - $2 - $3}' > data.all 
 for ((j=1; j<=`wc -l data.all | awk '{print $1}'`; j+=1)) do
 echo $j , >> time
 done
 paste -d " " time data.all > MMGBSA_vs_time.dat  
 rm namelist time data.*

To run this script do: 

 bash get.mmgbsa.sh

This will create a text file called MMGBSA_vs_time.dat with x and y values separated by a space and comma. These values can be imported to Excel or Origin or to XMGRACE if you are using Linux: 

 xmgrace MMGBSA_vs_time.dat

2013 AMBER MMGBSA plot.jpg MMGBSA.png

V. Frequently Encountered Problems

Mike

It is very important that you know which filetype you need. For example, for minimizations be sure that you have the correct input .mol2 file. The wrond mol2 file will result in file output with no information.

Ivan

Double and triple check your syntax before you submit a job to seawulf! It is incredibly disappointing to come back the next morning and realize your job failed simply because of a misplaced period or dash.

Lu

Yan

Pay attention to the names of your directories and files. Mostly, the reason why your commands fail to run is just a minor mirror in typing.

Yao

Junjie

Before you start minimization, equilibration and production, if the queue is not busy, you can type "qsub -I" in seawulf and "cd" to where you want to run sander. Then copy your sander command in your .csh file that you are going to submit and run it by simply paste to the command line and hit "return". By doing this, you are running your job interactively. If you can tell the job is running, your script is likely to be correct. It is easy to debug your .csh script in this way, since error is printed out right after you run commands instead of stored in some .er or .out files that have no context.

Tianao

Kai

Arkin

Jess

Mosavverul

Retrieved from "https://ringo.ams.stonybrook.edu/index.php?title=2014_AMBER_tutorial_with_HIV_Protease&oldid=8661"