2013 AMBER Tutorial with UMP and OMP

From Rizzo_Lab
Revision as of 15:05, 27 March 2013 by Stonybrook (talk | contribs) (MM-GBSA Energy Calculation)
Jump to: navigation, search

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 - Assisted Model Building with Energy Refinement - 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 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, Amber Tools User's Manual serves as another reference while using Amber tools.

Here are some programs in Amber

  1. 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.
  2. ANTECHAMBER: in additional to LEap, this main Antechamber suite program is for preparing input files other than standard nucleic acids and proteins.
  3. 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.
  4. PMEMD: version of SANDER that has improved parallel scaling property and optimized speed.
  5. 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 mailing list you could sign-up for, as an additional resource.


For information of the UMP-OMP system, see 2013 DOCK tutorial with Orotodine Monophosphate Decarboxylase.

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:


II. Structural Preparation

Preparation in Chimera

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.

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

In this section, we will create three new files and save them in the 001.CHIMERA.MOL.PREP/ folder:

 1LOQ.dockprep.mol2      (complete system with hydrogens and charges)
 1LOQ.receptor.noH.mol2  (the receptor alone, without hydrogens)
 1LOQ.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:

 %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 (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, 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).

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

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 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.


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 H5 cannot be distinguished from each other. )

    22 H5'        40.0697   36.0506   37.6716 H         1 LIG210      0.0761
    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:


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

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:

antechamber -i ../001.CHIMERA.MOL.PREP/1LOQ.lig.mol2 -fi mol2 -o 1LOQ.lig.ante.pdb -fo 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 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:



Parmchk is another program included in the Antechamber package. 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:

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


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.

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

cp  ~lingling/AMS536/AMBER_Tutorial/002.ANTE.TLEAP/tleap.lig.in
cp  ~lingling/AMS536/AMBER_Tutorial/002.ANTE.TLEAP/tleap.rec.in
cp  ~lingling/AMS536/AMBER_Tutorial/002.ANTE.TLEAP/tleap.com.in

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

You will see the following output files: 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

III. Simulation using sander




Running jobs on the queue

To run equilibration and production, Write the tcsh file equil.produc.qsub.csh:

 #! /bin/tcsh
 #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

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 que and monitor progress.

IV. Simulation Analysis


RMSD Plots

Measuring h-bond distances

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
 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/csh
 #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/yechen1/AMBER-Tutorial/005.MMGBSA"
 set sander  = "/nfs/user03/wjallen/local/amber12/bin/sander"
 cd ${workdir}

 $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 \
 -r restrt.com -ref ../002.ANTE.TLEAP/1LOQ.com.gas.leap.rst7 -x mdcrd.com -inf mdinfo.com
 $sander -O -i gb.rescore.in -o gb.rescore.out.lig -p ../002.ANTE.TLEAP/1LOQ.lig.gas.leap.prm7 \
 -c ../002.ANTE.TLEAP/1LOQ.lig.gas.leap.rst7 -y ../004.PTRAJ/1LOQ.lig.trj.stripfit \
 -r restrt.lig -ref ../002.ANTE.TLEAP/1LOQ.lig.gas.leap.rst7 -x mdcrd.lig -inf mdinfo.lig
 $sander -O -i gb.rescore.in -o gb.rescore.out.rec -p ../002.ANTE.TLEAP/1LOQ.rec.gas.leap.prm7 \
 -c ../002.ANTE.TLEAP/1LOQ.rec.gas.leap.rst7 -y ../004.PTRAJ/1LOQ.rec.trj.stripfit \
 -r restrt.rec -ref ../002.ANTE.TLEAP/1LOQ.rec.gas.leap.rst7 -x mdcrd.rec -inf mdinfo.rec

~ ~ ~ ~ ~ ~ ~ ~ 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:


  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
1-4 VDW =      756.7767  1-4 EEL =     7260.2823  RESTRAINT  =        0.0000
ESURF   =     9752.0291

minimization completed, ENE= 0.36269092E+04 RMS= 0.187371E+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.              

Equations for analysis

Plotting Energy

V. Frequently Encountered Problems