(Note: These tutorials are meant to provide illustrative examples of how to use the AMBER software suite to carry out simulations that can be run on a simple workstation in a reasonable period of time. They do not necessarily provide the optimal choice of parameters or methods for the particular application area.)
Copyright Ross Walker 2006

TUTORIAL 2 - SECTION 3

Using dynamics simulations to estimate binding energetics

By Ross Walker

Stage 3 - Run Production MD and Estimate the Binding Energy

The final stage now we have equilibrated our system is to run a production simulation. During this time we should monitor the energy of our 3 systems and calculate the average energy over the production run. Fortunately sander does this for us automatically and writes out the average energy at the end of the output file. So, we will run our 3 production calculations. These are identical to the equilibration phase above except that we will state that we are doing a restart so that we read the positions and velocities from the restrt file produced by the equilibration phase. We will also only run our MD for 1ps. Here is the input file:

$AMBERHOME/exe/sander -O -i md2.in -o md2_complex.out -p complex.prmtop -c complex_md1.restrt -r complex_md2.restrt -x complex_md2.mdcrd &

$AMBERHOME/exe/sander -O -i md2.in -o md2_protein.out -p protein.prmtop -c protein_md1.restrt -r protein_md2.restrt -x protein_md2.mdcrd &

$AMBERHOME/exe/sander -O -i md2.in -o md2_peptide.out -p peptide.prmtop -c peptide_md1.restrt -r peptide_md2.restrt -x peptide_md2.mdcrd &

Here are the output files: md2_complex.out, complex_md2.restrt, complex_md2.mdcrd, md2_protein.out, protein_md2.restrt, protein_md2.mdcrd, md2_peptide.out, peptide_md2.restrt, peptide_md2.mdcrd

Once again let's take a look at the energies in the output files:

mkdir analysis_md2_complex
mkdir analysis_md2_protein
mkdir analysis_md2_peptide
cd  analysis_md2_complex
process_mdout.perl ../md2_complex.out
cd ../analysis_md2_protein
process_mdout.perl ../md2_protein.out
cd ../analysis_md2_peptide
process_mdout.perl ../md2_peptide.out
cd ..

Total Energy

xmgr analysis_md2_complex/summary.ETOT analysis_md2_protein/summary.ETOT analysis_md2_peptide/summary.ETOT

Potential Energy

xmgr analysis_md2_complex/summary.EPTOT analysis_md2_protein/summary.EPTOT analysis_md2_peptide/summary.EPTOT

Total Energy	Potential Energy

Again, we can see that we have not equilibrated sufficiently since the mean value of our energy still drifts over time. Thus we cannot expect our calculated binding energy to be accurate. For this tutorial though it will suffice. Now, look at the three output files and find the average energies. The relevant lines are:

COMPLEX

      A V E R A G E S   O V E R    1000 S T E P S


 NSTEP =     1000   TIME(PS) =       3.000  TEMP(K) =   279.63  PRESS =     0.0
 Etot   =      -512.9891  EKtot   =       726.8297  EPtot      =     -1239.8187
 BOND   =       293.4816  ANGLE   =       437.7059  DIHED      =       530.1570
 1-4 NB =       180.9881  1-4 EEL =      2439.6489  VDWAALS    =      -322.5055
 EELEC  =     -3795.3103  EGB     =     -1003.9845  RESTRAINT  =         0.0000
 ------------------------------------------------------------------------------

PROTEIN

      A V E R A G E S   O V E R    1000 S T E P S


 NSTEP =     1000   TIME(PS) =       3.000  TEMP(K) =   274.17  PRESS =     0.0
 Etot   =      -501.3990  EKtot   =       599.8667  EPtot      =     -1101.2657
 BOND   =       242.2721  ANGLE   =       352.3939  DIHED      =       435.6360
 1-4 NB =       151.4390  1-4 EEL =      2007.4174  VDWAALS    =      -266.9382
 EELEC  =     -3251.1970  EGB     =      -772.2890  RESTRAINT  =         0.0000
 ------------------------------------------------------------------------------

PEPTIDE

      A V E R A G E S   O V E R    1000 S T E P S


 NSTEP =     1000   TIME(PS) =       3.000  TEMP(K) =   260.67  PRESS =     0.0
 Etot   =       -27.2990  EKtot   =       105.6709  EPtot      =      -132.9699
 BOND   =        45.3439  ANGLE   =        71.0219  DIHED      =        86.3453
 1-4 NB =        28.5914  1-4 EEL =       432.7471  VDWAALS    =       -25.7689
 EELEC  =      -548.5866  EGB     =      -222.6639  RESTRAINT  =         0.0000
 ------------------------------------------------------------------------------

From these three files we can get the average potential energies of our three systems:

Complex = -1239.8187 KCal/mol

Protein = -1101.2657 KCal/mol

Peptide = -132.9699 KCal/mol

Our binding energy is then simply Ecomplex - Eprotein - Epeptide = -5.5831 KCal/mol

So, our complex is stabilised by about 5.6 KCal/mol. Note, however, that we did not run a long enough equilibration or production run and so our answer is on an approximation at best. Nevertheless you should now be familiar with how to run minimisation and molecular dynamics with Amber.

(Note: These tutorials are meant to provide illustrative examples of how to use the AMBER software suite to carry out simulations that can be run on a simple workstation in a reasonable period of time. They do not necessarily provide the optimal choice of parameters or methods for the particular application area.)
Copyright Ross Walker 2006