Amber Advanced Tutorials - Tutorial 3 - MM-PBSA - Section 2 : Run the production simulation and obtain an ensemble of snapshots.

(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

The production phase of the simulation should be run using the same conditions as the final phase of equilibration to prevent an abrupt jump in the potential energy due to a change in simulation conditions.

We will run a total of 2 ns or production recording the coordinates every 10 ps. This should be sufficiently far apart that the structures are uncorrelated. Depending on your system you might obtain good results with snapshots taken closer together. As long as all the structures you obtain are uncorrelated the more snapshots you have the lower the statistical error of your results should be. Note for system such as the RAS-RAF complex we have here a simulation time of 2 ns is most likely too short to obtain a set of uncorrelated snapshots that adequately sample the equilibrium ensemble. A value of 20 ns or so would probably be more appropriate. However, this will suffice for the purposes of this tutorial.

prod.in

prod ras-raf
 &cntrl
  imin=0,irest=1,ntx=5,
  nstlim=250000,dt=0.002,
  ntc=2,ntf=2,
  cut=8.0, ntb=2, ntp=1, taup=2.0,
  ntpr=5000, ntwx=5000,
  ntt=3, gamma_ln=2.0,
  temp0=300.0,
 /

This should then be run 4 times to obtain 2 ns of simulation time. Since this is a simple periodic boundary PME simulation one can used PMEMD to do the simulation if required. This will typical offer better performance and scaling in parallel. Below is an example script I used on San Diego Supercomputer Center's Teragrid Cluster to run this job on 96 processors. The calculation took a total of 10 hours.

run.x

#SDSC Teragrid PBS Script
#PBS -j oe
#PBS -l nodes=48:ppn=2
#PBS -l walltime=12:00:00
#PBS -q dque
#PBS -V
#PBS -M name@email.com
#PBS -A account_no
#PBS -N run_pmemd_96

cd /gpfs/projects/prod/
mpirun -v -machinefile $PBS_NODEFILE -np 96 /usr/local/apps/amber9/exe/pmemd -O -i prod.in -o prod1.out \
-p ras-raf_solvated.prmtop -c equil.rst -r prod1.rst -x prod1.mdcrd
mpirun -v -machinefile $PBS_NODEFILE -np 96 /usr/local/apps/amber9/exe/pmemd -O -i prod.in -o prod2.out \
-p ras-raf_solvated.prmtop -c prod1.rst -r prod2.rst -x prod2.mdcrd
mpirun -v -machinefile $PBS_NODEFILE -np 96 /usr/local/apps/amber9/exe/pmemd -O -i prod.in -o prod3.out \
-p ras-raf_solvated.prmtop -c prod2.rst -r prod3.rst -x prod3.mdcrd
mpirun -v -machinefile $PBS_NODEFILE -np 96 /usr/local/apps/amber9/exe/pmemd -O -i prod.in -o prod4.out \
-p ras-raf_solvated.prmtop -c prod3.rst -r prod4.rst -x prod4.mdcrd
gzip -9 prod*.mdcrd

It is essential, for good results, that our system still be exploring equilibrium phase space during the production phase. We will check this in the same fashion as we did for the last equilibration step by plotting the density, temperature, total energy and backbone RMSD.

Note the production RMSD does not look to be truly equilibrated while the other properties are essentially constant (note the small scales). Ideally we should probably run a much longer production run (ca. 20 ns). For the purposes of this tutorial, however, we will continue with what we have.

We now need to extract snapshots (without the water) from our production runs for use in the MM-PBSA calculation. The mm_pbsa.pl script (in the amber9 exe directory) automates this extraction process for us. Note if you don't have gzcat installed you will need to unzip the trajectory files before running mm_pbsa. We have to provide the following input file (see the linked file for an explanation of each of the various terms):

extract_coords.mmpbsa

@GENERAL
PREFIX                snapshot
PATH                  ./
COMPLEX               1
RECEPTOR              1
LIGAND                1
COMPT                 ./ras-raf.prmtop
RECPT                 ./ras.prmtop
LIGPT                 ./raf.prmtop
GC                    1
AS                    0
DC                    0
MM                    0
GB                    0
PB                    0
MS                    0
NM                    0
@MAKECRD
BOX                   YES
NTOTAL                42193
NSTART                1
NSTOP                 200
NFREQ                 1
NUMBER_LIG_GROUPS     1
LSTART                2622
LSTOP                 3862
NUMBER_REC_GROUPS     1
RSTART                1
RSTOP                 2621
@TRAJECTORY
TRAJECTORY            ./prod1.mdcrd
TRAJECTORY            ./prod2.mdcrd
TRAJECTORY            ./prod3.mdcrd
TRAJECTORY            ./prod4.mdcrd
@PROGRAMS

This just specifies which atoms are part of the receptor, ligand and complex as well as specifying the prmtop files corresponding to the unsolvated structures, the total number of snapshots in the trajectores, the stride length and the names of the trajectory files. We have specified that each of the output files should be written with the prefix 'snapshot'. In this setup we have defined RAS to be the receptor and RAF to be the ligand. This is purely a naming convention.

This will take a few minutes to run. Here are the relevant output files: extract_coords.tar.gz (14 MB)

You should check the log file for any errors. Also make sure the box sizes look reasonable. If there are any strange box sizes this is likely a problem with the atom numbers you selected or you a corruption in one of your trajectory files.

We can now proceed to section 3 where we will calculate the binding free energy.

DENSITY	TEMPERATURE
TOTAL ENERGY	BACKBONE RMSD