Amber Workshop - tutorial 8

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

Case Study: All Atom Structure Prediction and Folding Simulations
of a Stable Protein (Folding Trp-Cage Peptide)

Before we start running molecular dynamics we need to perform a short minimisation of our starting structure. This won't cause our system to fold up since minimisation remains in the local minima of the starting structure. It should cleanup the structure, fix up hydrogen positions etc., such that our MD, when we start it, will be stable.

min1.in

Stage 1 - minimisation of TC5b
 &cntrl
  imin=1, maxcyc=1000, ncyc=500,
  cut=999., rgbmax=999.,igb=1, ntb=0,
  ntpr=100
 /

We will be running a total of 1000 steps of minimisation; 500 of steepest descent (ncyc=500) followed by 500 of conjugate gradient (maxcyc-ncyc). This should be sufficient to clean up our structure. Note the large cutoff I am using (cut=999. angstroms). The reason for this is that we are running a non-periodic simulation (ntb=0) and so we do not have the Particle Mesh Ewald (PME) method to give us infinite electrostatics. With PME the recommended cutoff is 8 angstroms as this only truncates the very short range VDW interactions. However, without PME both the VDW and electrostatic interactions are truncated at the cutoff so ideally we need to make this as large as we can afford. Unfortunately the time taken for a simulation is proportional to the square of the cutoff - see tutorial 1 section 5.1.2 for a plot of simulation time vs cutoff size. Fortunately this system is small enough that we can afford to run with no cutoff. A 999 angstroms cutoff is sufficiently large that we will never truncate any of our nonbonded interactions during this simulation. Similarly we set rgbmax to 999 angstroms as well. This controls the maximum distance between atom pairs that will be considered when calculating the effective Born radii. As with the cutoff the bigger this value is the better but the bigger the value the longer the simulation takes. Since this system is only 20 residues in size we can afford to include all atom pairs in the effective Born radii calculation, so we set rgbmax to be significantly more than the largest extent of the system.

$AMBERHOME/exe/sander -O -i min1.in -o min1.out -p TC5b.prmtop
-c TC5b.inpcrd -r min1.rst

Open the two files (min1.pdb and TC5b_linear.pdb) in the viewing package of your choice, e.g. VMD

The original structure is shown in blue while the minimised structure is shown in yellow. As you can see the minimisation has not changed the backbone structure much, but it has moved the tryptophan and tyrosine residues considerably to remove steric clashes. These high energy hotspots would have caused problems had we just started out with molecular dynamics. If you don't believe me try running the MD from the inpcrd file we created and watch what happens.