Amber Workshop - tutorial 2 - Section 2: Minimize and Equilibrate our systems

(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 next stage is to minimise our complex and then use molecular dynamics to sample conformational space. Note, the procedure given here is very short in order to make the simulations compatible with the timescale of this tutorial. In a real "production" simulation you would typically run a much longer simulation (ns) in order to obtain good statistical convergence.

Okay, step 1 let's minimise it to remove any bad contacts created by the hydrogenation.

Here's our input file, we will do a total of 500 steps of minimisation (MAXCYC) with the first 200 being steepest descent (NCYC), the remainder will be conjugate gradient (MAXCYC-NCYC). We will use a reasonably large cut off of 16 Angstroms since this is not going to be a periodic simulation and we want to deal with our electrostatics accurately (NTB=0,CUT=16). For our implicit solvent model we will use the GB model of Hawkins, Cramer and Truhlar, see the AMBER manual for a full description (IGB=1):

Initial minimisation of our complex
&cntrl
imin=1, maxcyc=500, ncyc=200,
cut=16, ntb=0, igb=1,
/

min.in

$AMBERHOME/exe/sander -O -i min.in -o min_complex.out -p complex.prmtop -c complex.inpcrd -r complex_min.crd &

The ampersand will background the job so that we get the command line back. We can monitor the progress by tailing the output file:

Now we wait for it to run. On a 3GHz P4 this will take around 80 seconds to run. While you wait you could take a look at the starting structures in vmd (or go get a coffee). If you can't wait the output files are here:

We now repeat the process for the protein and the peptide (if our workstation was a multiprocessor machine we could run these simulations all at the same time and save us time. If we had a cluster with a queuing system we could also submit all three jobs to the queue and they would run in parallel. It is also possible to use the parallel version of sander to use multiple processors for a single job in order to decrease the time required. This however is beyond the scope of this tutorial):

$AMBERHOME/exe/sander -O -i min.in -o min_protein.out -p protein.prmtop -c protein.inpcrd -r protein_min.crd &

$AMBERHOME/exe/sander -O -i min.in -o min_peptide.out -p peptide.prmtop -c peptide.inpcrd -r peptide_min.crd &

You can use ambpdb to generate a pdb of the final minimised structures if you want:

It is always a good idea to visualise your results since the human eye is very good at spotting anomolies.

The next step if to heat and equilibrate our three systems. For speed we will do this very rapidly over 2ps. Ideally you should do this for much longer.

Here's our input file, we will run MD (imin=0) and this is not a restart (irest=0). In this example we will not use shake since it is possible that the hydrogen motion may effect the binding energy (probably not, but it serves as an example here). As we are not using shake we will need a time step smaller than normal. Here I will use a time step of 1 fs and run for 2000 steps [2 ps] (dt = 0.001, nstlim=2000, ntc=1). We will also write to our output file every 20 steps and to our trajectory [mdcrd] file every 20 steps (ntpr=20,ntwx=20). For temperature control we will use a Langevin dynamics approach with a collision frequency of 1 ps^-1. We will start our system at 0K and we want a target temperature of 300K (ntt=3, gamma_ln=1.0, tempi=0.0, temp0=300.0):

Initial MD equilibration
&cntrl
imin=0, irest=0,
nstlim=2000,dt=0.001, ntc=1,
ntpr=20, ntwx=20,
cut=16, ntb=0, igb=1,
ntt=3, gamma_ln=1.0,
tempi=0.0, temp0=300.0,
/

md1.in

We run this 2ps of equilibration for all 3 of our systems (using the final structure from the minimisation as our input structure). On a P4 3GHz machine the complex md sim takes about 322 seconds.:

$AMBERHOME/exe/sander -O -i md1.in -o md1_complex.out -p complex.prmtop -c complex_min.crd -r complex_md1.restrt -x complex_md1.mdcrd &

$AMBERHOME/exe/sander -O -i md1.in -o md1_protein.out -p protein.prmtop -c protein_min.crd -r protein_md1.restrt -x protein_md1.mdcrd &

$AMBERHOME/exe/sander -O -i md1.in -o md1_peptide.out -p peptide.prmtop -c peptide_min.crd -r peptide_md1.restrt -x peptide_md1.mdcrd &

At this point we should check that our simulations are acceptable and haven't "blown up" due to instabilities, poor hydrogenation etc. First off we can visualise the trajectory in VMD as we did in tutorial 1. Try this now. We are looking for a nice stable, equilibrated trajectory. We can also look at the energies and temperatures to see if these have stabilised. For this we will use the process_mdout.perl script we used in tutorial 1. We should also compute the RMSD using ptraj and see this is stable as we did in tutorial one but in the interests of speed we will skip this here and just look at the energies.

mkdir analysis_md1_complex
mkdir analysis_md1_protein
mkdir analysis_md1_peptide
cd analysis_md1_complex
process_mdout.perl ../md1_complex.out
cd ../analysis_md1_protein
process_mdout.perl ../md1_protein.out
cd ../analysis_md1_peptide
process_mdout.perl ../md1_peptide.out
cd ..

You can now plot these data files using your favourite plotting program. I will use xmgr as we did in tutorial one:

xmgr analysis_md1_complex/summary.ETOT analysis_md1_protein/summary.ETOT analysis_md1_peptide/summary.ETOT

xmgr analysis_md1_complex/summary.EPTOT analysis_md1_protein/summary.EPTOT analysis_md1_peptide/summary.EPTOT

The two graphs above show that our simulation seems reasonable, there are no sudden jumps in our energy. However, as we can see we have not yet reached an equilibrium. This is because our simulation was way too short. In reality this simulation should be extended much further. For our purposes we will ignore this for the moment and just bare in mind that the final binding energy we get will not be accurate.


Total Energy	Potential Energy