(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 2015

AMBER Advanced Tutorial 25
SECTION 1

GIST Analysis Example
Analysis of water thermodynamics using Grid Inhomogeneous Solvation Theory of Factor Xa active site.

By Romelia Salomon, Crystal Nguyen, Steven Ramsey, Jonathan D. Gough, Mike Gilson, Tom Kurtzman, & Ross Walker

1) Generating and Relaxing the Initial Structure

The first thing we need to do is generate an initial structure file that can be used as the input for the molecular dynamics simulation. Typically one obtains such a structure from the pdb databank. For simplicity we are including it here 1FJS.pdb. We then use the program Leap, included in Amber Tools, to generate the input files:

For this system several steps need to be followed to prepare the files.

  1. All of the water molecules, keeping the ions, need to be removed.
  2. The appropriate disulfide bonds need to be made.
  3. The Hydrogen atoms need to be added.
  4. The entire protein needs to be solvated with explicit water (in this case we have chosen to use the TIP3P water model).

The following leap commands are an example of how the pdb and ion parameters are loaded into leap, the disulfide bonds are made, the protein is solvated, and the parameters are saved:

$AMBERHOME/bin/tleap
> source leaprc.ff99SB
> loadamberparams frcmod.ff99SB
> loadoff Ca2.lib
> prot = loadpdb ChainA_IONS_H.pdb
> bond prot.7.SG prot.12.SG
> bond prot.27.SG prot.43.SG
> bond prot.156.SG prot.170.SG
> bond prot.181.SG prot.209.SG
> solvateBox prot TIP3PBOX 10
> saveamberparm prot ChA_IONS_disul_TIP3P.top ChA_IONS_disul_TIP3P.inpcrd
> quit

Leap outputs the following files which will be used as the paramater-topology and the initial coordinates for our MD simulation: ChA_IONS_disul_TIP3P.prmtop, ChA_IONS_disul_TIP3P.inpcrd.

Using these input files, we equilibrate the system before performing the production MD simulation. This process entails minimizing, heating, and relaxing the system before performing the production run. The equilibration was done in 6 sequential stages:

  1. Minimize only the water while restraining the protein (20000 cycles)
  2. Minimize the water and the protein hydrogen atoms (20000 cycles)
  3. Heat the system to 50K while restraining the protein heavy atoms (NVT 0 to 50K)
  4. Heat the system from 50K to 300K in increments of 50K for 200ps each
  5. Relax the system, restraining the protein heavy atoms (NPT, 300K, 10ns)
  6. Relax the system, restraining all atoms in the protein (NVT, 300K, 5ns)
min_sys-1.in min_sys-2.in md-heat-50K.in
minimize water
 &cntrl
   imin=1, ntmin=1, nmropt=0, drms=0.1,
   maxcyc=20000, ncyc=1500,
   ntx=1,	// not a restart
   irest=0,
   ntpr=100, ntwr=100, iwrap=1,
   ntf=1,	// SHAKE is off
   ntb=1,	 // constant volume pbc 
   cut=9.0,
   ntr=1, restraint_wt=100, // Cartesian restraint
   restraintmask=':1-234',
 /
minimize all
 &cntrl
   imin=1, ntmin=1, nmropt=0, drms=0.1,
   maxcyc=20000, ncyc=1500,
   ntx=1, irest=0,
   ntpr=100, ntwr=100, iwrap=1,
   ntf=1, ntb=1, cut=9.0,
   ntr=1,
   restraint_wt=100,
   restraintmask=':1-234 & !@H',
 /
heat NVT 200ps
 &cntrl
  imin = 0,
  irest = 0, ntx = 1, 
  ntb = 1, ntc = 2, ntf = 2,
  cut = 9.0,
  tempi = 0.0, temp0 = 50.0,
  ntt = 3, gamma_ln = 2.0,
  nstlim = 10000, dt = 0.002,
  iwrap = 1,
  ntpr = 1000, ntwx = 1000, ntwr = 10000,
  nscm = 1000, 
  ntr=1, restraint_wt=100,
  restraintmask=':1-234 & !@H',
  ig=-1,
  ioutfm=1,
 /
md-heat-300K.in md-eq-npt-10ns.in md-eq-nvt-5ns.in
heat NPT 200ps
 &cntrl
  imin = 0,
  irest = 1, ntx = 5,
  ntb = 2, ntc = 2, ntf = 2,
  cut = 9.0,
  temp0 = 300.0,
  ntt = 3, gamma_ln = 2.0,
  nstlim = 10000, dt = 0.002,
  ntp = 1, taup = 0.5,
  iwrap = 1,
  ntpr = 1000, ntwx = 1000, ntwr = 10000,
  nscm = 1000,
  ntr=1, restraint_wt=100,
  restraintmask=':1-234 & !@H',
  ig=-1,
  ioutfm=1,
 /
equil NTP 10ns
 &cntrl
  imin=0,
  ntpr   = 10000,    ntwx   = 10000,
  ntwr   = -1000000,
  iwrap=1,
  ntx=5, irest=1,
  ntf    = 2,       ntb    = 2,
  ntc    = 2,       cut    = 9.0,
  nstlim = 5000000,
  dt     = 0.002,
  nscm   = 1000,
  temp0  = 300.0, tempi = 300.0,
  ntt    = 3, gamma_ln = 2.0,
  ntp    = 1, taup   = 0.5,
  ntr=1, restraint_wt=100,
  restraintmask=':1-234 & !@H',
  ig=-1,
  ioutfm=1,
 /
equil NVT 5ns
 &cntrl
  imin=0,
  ntpr   = 10000,    ntwx   = 10000,
  ntwr   = -500000,
  iwrap=1,
  ntx=5, irest=1,
  ntf    = 2,       ntb    = 1,
  ntc    = 2,       cut    = 9.0,
  nstlim = 2500000,
  dt     = 0.002,
  nscm   = 1000,
  temp0  = 300.0,
  ntt    = 3, gamma_ln = 2.0,
  ntr=1, restraint_wt=100,
  restraintmask=':1-234',
  ig=-1,
  ioutfm=1,
 /

The following is an example of the commands to run the equilibration, using pmemd, on 2 processors of a dual cpu computer:

  1. $AMBERHOME/bin/pmemd -O -i min_sys-1.in -o min1.out -p ChA_IONS_disul_TIP3P.prmtop -c ChA_IONS_disul_TIP3P.inpcrd -r min_fxa-1.rst7


  2. $AMBERHOME/bin/pmemd -O -i min_sys-2.in -o min2.out -p ChA_IONS_disul_TIP3P.prmtop -c min_fxa-1.rst7 -r min_fxa-2.rst7 -ref min_fxa-1.rst7


  3. mpirun -np 2 $AMBERHOME/bin/pmemd.MPI -O -i md-heat-50K.in -o heat1.out -p ChA_IONS_disul_TIP3P.prmtop -c min_fxa-2.rst7 -r md_heat-50K.rst7 -x heat1.nc -ref min_fxa-2.rst7


  4. mpirun -np 2 $AMBERHOME/bin/pmemd.MPI -O -i md-heat-300K.in -o heat2.out -p ChA_IONS_disul_TIP3P.prmtop -c md_heat-50K.rst7 -r md_heat-300K.rst7 -x heat2.nc -ref md_heat-50K.rst7


  5. mpirun -np 2 $AMBERHOME/bin/pmemd.MPI -O -i md-eq-npt-10ns.in -o equil-npt.out -p ChA_IONS_disul_TIP3P.prmtop -c md_heat-300K.rst7 -r md_eq-npt.rst7 -x equil-npt.nc -ref md_heat-300K.rst7


  6. mpirun -np 2 $AMBERHOME/bin/pmemd.MPI -O -i md-eq-nvt-5ns.in -o equil-nvt.out -p ChA_IONS_disul_TIP3P.prmtop -c md_eq-npt.rst7 -r md_eq-nvt.rst7 -x equil-nvt.nc -ref md_eq-npt.rst7


After running these simulations, we get the following files for the different stages:

  1. Minimizing only the water gives the output file min1.out and the restart file min_fxa-1.rst7
  2. Minimizing the water and protein Hydrogen atoms give the output file min2.out and the restart file min_fxa-2.rst7
  3. Heating under NPT conditions from 0-50K gives the output file heat1.out the trajectory heat1.nc and the restart file md_heat-50K.rst7
  4. Heating under NPT conditions from 250-300K gives the output heat2.out the trajectory heat2.nc and the restart file md_heat-300K.rst7
  5. Relaxing the system under NPT conditions gives the output file equil-npt.out the trajectory equil-npt.nc and the restart file md_eq-npt.rst7
  6. Relaxing the system under NVT conditions gives the output file equil-nvt.out the trajectory equil-nvt.nc and the restart file md_eq-nvt.rst7

If you run the calculations as stated above, you won't get exactly the same files, due to the use of different random seeds, but an equivalent initial configuration

IMPORTANT NOTE: After the NPT run, the final structure is re-aligned to the reference structure to account for any potential translation of the protein during the equilibration. Slight movement of the host can cause the system to crash shortly after the NVT production starts. To accomplish this, a basic cpptraj script is used to as follows. image.ptraj

cpptraj -i image.ptraj

The final, re-aligned structure used as the input coordinates for the production run are as these prod_100ns.rst7


CLICK HERE TO GO TO SECTION 2


(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 2015